From bb2af8cf67cd8c760cab8c72c022013ab43e3df0 Mon Sep 17 00:00:00 2001 From: "anning.cheng" Date: Wed, 28 Jun 2023 23:28:58 -0400 Subject: [PATCH] optimized the code and decrease the stratosphere warm bias for mraerosol=T --- physics/GFS_rrtmgp_cloud_mp.F90 | 7 +------ physics/module_mp_thompson.F90 | 20 ++++++++++++-------- physics/mp_thompson.F90 | 20 +++++++++++--------- 3 files changed, 24 insertions(+), 23 deletions(-) diff --git a/physics/GFS_rrtmgp_cloud_mp.F90 b/physics/GFS_rrtmgp_cloud_mp.F90 index 32104b7f8..79ae1559a 100644 --- a/physics/GFS_rrtmgp_cloud_mp.F90 +++ b/physics/GFS_rrtmgp_cloud_mp.F90 @@ -875,17 +875,12 @@ subroutine cmp_reff_Thompson(nLev, nCol, i_cldliq, i_cldice, i_cldsnow, i_cldice qi_mp(iCol,iLay) = tracer(iCol,iLay,i_cldice) / (1.-q_lay(iCol,iLay)) qs_mp(iCol,iLay) = tracer(iCol,iLay,i_cldsnow) / (1.-q_lay(iCol,iLay)) ni_mp(iCol,iLay) = tracer(iCol,iLay,i_cldice_nc) / (1.-q_lay(iCol,iLay)) - if (ltaerosol) then + if (ltaerosol .or. mraerosol) then nc_mp(iCol,iLay) = tracer(iCol,iLay,i_cldliq_nc) / (1.-q_lay(iCol,iLay)) nwfa(iCol,iLay) = tracer(iCol,iLay,i_twa) if (qc_mp(iCol,iLay) > 1.e-12 .and. nc_mp(iCol,iLay) < 100.) then nc_mp(iCol,iLay) = make_DropletNumber(qc_mp(iCol,iLay)*rho, nwfa(iCol,iLay)*rho) * orho endif - elseif (mraerosol) then - nc_mp(iCol,iLay) = tracer(iCol,iLay,i_cldliq_nc) / (1.-q_lay(iCol,iLay)) - if (qc_mp(iCol,iLay) > 1.e-12 .and. nc_mp(iCol,iLay) < 100.) then - nc_mp(iCol,iLay) = make_DropletNumber(qc_mp(iCol,iLay)*rho, nwfa(iCol,iLay)*rho) * orho - endif else if (nint(lsmask(iCol)) == 1) then !land nc_mp(iCol,iLay) = nt_c_l*orho diff --git a/physics/module_mp_thompson.F90 b/physics/module_mp_thompson.F90 index 6a4ef5e02..0c708cb3d 100644 --- a/physics/module_mp_thompson.F90 +++ b/physics/module_mp_thompson.F90 @@ -3400,8 +3400,8 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & tcond(k) = (5.69 + 0.0168*tempc)*1.0E-5 * 418.936 ocp(k) = 1./(Cp*(1.+0.887*qv(k))) lvt2(k)=lvap(k)*lvap(k)*ocp(k)*oRv*otemp*otemp - - nwfa(k) = MAX(11.1E6*rho(k), (nwfa1d(k) + nwfaten(k)*DT)*rho(k)) + if (is_aerosol_aware) & + nwfa(k) = MAX(11.1E6*rho(k), (nwfa1d(k) + nwfaten(k)*DT)*rho(k)) enddo do k = kts, kte @@ -3654,7 +3654,8 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & qvten(k) = qvten(k) - prw_vcd(k) qcten(k) = qcten(k) + prw_vcd(k) ncten(k) = ncten(k) + pnc_wcd(k) - nwfaten(k) = nwfaten(k) - pnc_wcd(k) + if (is_aerosol_aware) & + nwfaten(k) = nwfaten(k) - pnc_wcd(k) tten(k) = tten(k) + lvap(k)*ocp(k)*prw_vcd(k)*(1-IFDRY) rc(k) = MAX(R1, (qc1d(k) + DT*qcten(k))*rho(k)) if (rc(k).eq.R1) L_qc(k) = .false. @@ -3743,7 +3744,8 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & qrten(k) = qrten(k) - prv_rev(k) qvten(k) = qvten(k) + prv_rev(k) nrten(k) = nrten(k) - pnr_rev(k) - nwfaten(k) = nwfaten(k) + pnr_rev(k) + if (is_aerosol_aware) & + nwfaten(k) = nwfaten(k) + pnr_rev(k) tten(k) = tten(k) - lvap(k)*ocp(k)*prv_rev(k)*(1-IFDRY) rr(k) = MAX(R1, (qr1d(k) + DT*qrten(k))*rho(k)) @@ -4232,10 +4234,12 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & qv1d(k) = MAX(1.E-10, qv1d(k) + qvten(k)*DT) qc1d(k) = qc1d(k) + qcten(k)*DT nc1d(k) = MAX(2./rho(k), MIN(nc1d(k) + ncten(k)*DT, Nt_c_max)) - nwfa1d(k) = MAX(11.1E6, MIN(9999.E6, & - (nwfa1d(k)+nwfaten(k)*DT))) - nifa1d(k) = MAX(naIN1*0.01, MIN(9999.E6, & - (nifa1d(k)+nifaten(k)*DT))) + if (is_aerosol_aware) then + nwfa1d(k) = MAX(11.1E6, MIN(9999.E6, & + (nwfa1d(k)+nwfaten(k)*DT))) + nifa1d(k) = MAX(naIN1*0.01, MIN(9999.E6, & + (nifa1d(k)+nifaten(k)*DT))) + end if if (qc1d(k) .le. R1) then qc1d(k) = 0.0 nc1d(k) = 0.0 diff --git a/physics/mp_thompson.F90 b/physics/mp_thompson.F90 index e62e8a596..6a95a706c 100644 --- a/physics/mp_thompson.F90 +++ b/physics/mp_thompson.F90 @@ -151,6 +151,10 @@ subroutine mp_thompson_init(ncol, nlev, con_g, con_rd, con_eps, & !> - Convert specific humidity to water vapor mixing ratio. !> - Also, hydrometeor variables are mass or number mixing ratio !> - either kg of species per kg of dry air, or per kg of (dry + vapor). + if (merra2_aerosol_aware) then + call get_niwfa(aerfld, nifa, nwfa, ncol, nlev) + end if + qv = spechum/(1.0_kind_phys-spechum) @@ -163,7 +167,7 @@ subroutine mp_thompson_init(ncol, nlev, con_g, con_rd, con_eps, & ni = ni/(1.0_kind_phys-spechum) nr = nr/(1.0_kind_phys-spechum) - if (is_aerosol_aware) then + if (is_aerosol_aware .or. merra2_aerosol_aware) then nc = nc/(1.0_kind_phys-spechum) nwfa = nwfa/(1.0_kind_phys-spechum) nifa = nifa/(1.0_kind_phys-spechum) @@ -208,8 +212,6 @@ subroutine mp_thompson_init(ncol, nlev, con_g, con_rd, con_eps, & nwfa(i,k) = naCCN1+naCCN0*exp(-((hgt(i,k)-hgt(i,1))/1000.)*niCCN3) enddo enddo - else if (merra2_aerosol_aware) then - call get_niwfa(aerfld, nifa, nwfa, ncol, nlev) else if (mpirank==mpiroot) write(*,*) ' Apparently initial CCN aerosols are present.' if (MAXVAL(nwfa2d) .lt. eps) then @@ -555,6 +557,9 @@ subroutine mp_thompson_run(ncol, nlev, con_g, con_rd, & else dtstep = dtp end if + if (merra2_aerosol_aware) then + call get_niwfa(aerfld, nifa, nwfa, ncol, nlev) + end if !> - Convert specific humidity to water vapor mixing ratio. !> - Also, hydrometeor variables are mass or number mixing ratio @@ -574,7 +579,7 @@ subroutine mp_thompson_run(ncol, nlev, con_g, con_rd, & ni = ni/(1.0_kind_phys-spechum) nr = nr/(1.0_kind_phys-spechum) - if (is_aerosol_aware) then + if (is_aerosol_aware .or. merra2_aerosol_aware) then nc = nc/(1.0_kind_phys-spechum) nwfa = nwfa/(1.0_kind_phys-spechum) nifa = nifa/(1.0_kind_phys-spechum) @@ -681,9 +686,6 @@ subroutine mp_thompson_run(ncol, nlev, con_g, con_rd, & ncten3 => diag3d(:,:,36:36) qcten3 => diag3d(:,:,37:37) end if set_extended_diagnostic_pointers - if (merra2_aerosol_aware) then - call get_niwfa(aerfld, nifa, nwfa, ncol, nlev) - end if !> - Call mp_gt_driver() with or without aerosols, with or without effective radii, ... if (is_aerosol_aware .or. merra2_aerosol_aware) then call mp_gt_driver(qv=qv, qc=qc, qr=qr, qi=qi, qs=qs, qg=qg, ni=ni, nr=nr, & @@ -921,8 +923,8 @@ subroutine get_niwfa(aerfld, nifa, nwfa, ncol, nlev) aerfld(:,:,4)/1011.5142+ aerfld(:,:,5)/5683.3501)*1.e15 nwfa=((aerfld(:,:,6)/0.0045435214+aerfld(:,:,7)/0.2907854+aerfld(:,:,8)/12.91224+ & - aerfld(:,:,9)/206.2216+ aerfld(:,:,10)/4326.23)*1.+aerfld(:,:,11)/0.3053104*5+ & - aerfld(:,:,15)/0.3232698*1)*1.e15 + aerfld(:,:,9)/206.2216+ aerfld(:,:,10)/4326.23)*9.+aerfld(:,:,11)/0.3053104*5+ & + aerfld(:,:,15)/0.3232698*8)*1.e15 end subroutine get_niwfa end module mp_thompson