Merge pull request #88 from grantfirl/strat_warm_bias_fix_cheng

fixed stratosphere warm bias and code optimization for MERRA2
ufs-community · Jul 24, 2023 · 9b69974 · 9b69974
2 parents 5dc968e + a69c8a3
commit 9b69974
Show file tree

Hide file tree

Showing 3 changed files with 24 additions and 23 deletions.
diff --git 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
@@ -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
@@ -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