diff --git a/physics/GFS_PBL_generic.F90 b/physics/GFS_PBL_generic.F90 index 401084479..aefcb9cf4 100644 --- a/physics/GFS_PBL_generic.F90 +++ b/physics/GFS_PBL_generic.F90 @@ -12,7 +12,7 @@ module GFS_PBL_generic_common contains subroutine set_aerosol_tracer_index(imp_physics, imp_physics_wsm6, & - imp_physics_thompson, ltaerosol,mraerosol, & + imp_physics_thompson, ltaerosol,mraerosol, & imp_physics_mg, ntgl, imp_physics_gfdl, & imp_physics_zhao_carr, kk, & errmsg, errflg) @@ -191,7 +191,7 @@ subroutine GFS_PBL_generic_pre_run (im, levs, nvdiff, ntrac, rtg_ozone_index, vdftra(i,k,10) = qgrs(i,k,ntoz) enddo enddo - rtg_ozone_index = 8 + rtg_ozone_index = 10 else do k=1,levs do i=1,im diff --git a/physics/GFS_rrtmg_pre.F90 b/physics/GFS_rrtmg_pre.F90 index bb5af6afe..516d1a300 100644 --- a/physics/GFS_rrtmg_pre.F90 +++ b/physics/GFS_rrtmg_pre.F90 @@ -197,7 +197,7 @@ subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, n_var_lndp, & real(kind=kind_phys), dimension(im,lm+LTP+1) :: tem2db, hz real(kind=kind_phys), dimension(im,lm+LTP,min(4,ncnd)) :: ccnd - real(kind=kind_phys), dimension(im,lm+LTP,2:ntrac+2) :: tracer1 + real(kind=kind_phys), dimension(im,lm+LTP,2:ntrac) :: tracer1 real(kind=kind_phys), dimension(im,lm+LTP,NF_CLDS) :: clouds real(kind=kind_phys), dimension(im,lm+LTP,NF_VGAS) :: gasvmr real(kind=kind_phys), dimension(im,lm+LTP,NBDSW,NF_AESW) :: faersw @@ -651,7 +651,6 @@ subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, n_var_lndp, & nwfa (i,k) = tracer1(i,k,ntwa) enddo enddo - elseif (imp_physics == imp_physics_thompson) then do k=1,LMK do i=1,IM diff --git a/physics/GFS_suite_interstitial.F90 b/physics/GFS_suite_interstitial.F90 index 7c9c33cda..0023cf686 100644 --- a/physics/GFS_suite_interstitial.F90 +++ b/physics/GFS_suite_interstitial.F90 @@ -513,8 +513,7 @@ end subroutine GFS_suite_interstitial_3_finalize !! \htmlinclude GFS_suite_interstitial_3_run.html !! subroutine GFS_suite_interstitial_3_run (im, levs, nn, cscnv, & - satmedmf, trans_trac, do_shoc, & - ltaerosol, mraerosol, ntrac, ntcw, & + satmedmf, trans_trac, do_shoc, ntrac, ntcw, & ntiw, ntclamt, ntrw, ntsw, ntrnc, ntsnc, ntgl, ntgnc, & xlon, xlat, gt0, gq0, imp_physics, imp_physics_mg, & imp_physics_zhao_carr, imp_physics_zhao_carr_pdf, & @@ -534,8 +533,7 @@ subroutine GFS_suite_interstitial_3_run (im, levs, nn, cscnv, & ntrnc, ntsnc, ntgl, ntgnc, imp_physics, imp_physics_mg, imp_physics_zhao_carr, imp_physics_zhao_carr_pdf, & imp_physics_gfdl, imp_physics_thompson, imp_physics_wsm6,imp_physics_fer_hires, me, index_of_process_conv_trans integer, intent(in ), dimension(:) :: islmsk, kpbl, kinver - logical, intent(in ) :: cscnv, satmedmf, trans_trac, do_shoc, ltaerosol, ras - logical, intent(in ) :: mraerosol + logical, intent(in ) :: cscnv, satmedmf, trans_trac, do_shoc, ras integer, intent(in) :: ntinc, ntlnc logical, intent(in) :: ldiag3d, qdiag3d @@ -658,11 +656,11 @@ subroutine GFS_suite_interstitial_3_run (im, levs, nn, cscnv, & save_tcp(i,k) = gt0(i,k) enddo enddo - if(ltaerosol .or. mraerosol) then + if (ntinc>0) then save_qi(:,:) = clw(:,:,1) + end if + if (ntlnc>0) then save_qc(:,:) = clw(:,:,2) - else - save_qi(:,:) = clw(:,:,1) endif elseif (imp_physics == imp_physics_wsm6 .or. imp_physics == imp_physics_mg .or. imp_physics == imp_physics_fer_hires) then do k=1,levs @@ -699,10 +697,10 @@ end subroutine GFS_suite_interstitial_4_finalize !> \section arg_table_GFS_suite_interstitial_4_run Argument Table !! \htmlinclude GFS_suite_interstitial_4_run.html !! - subroutine GFS_suite_interstitial_4_run (im, levs, ltaerosol, mraerosol, tracers_total, ntrac, ntcw, ntiw, ntclamt, & - ntrw, ntsw, ntrnc, ntsnc, ntgl, ntgnc, ntlnc, ntinc, nn, imp_physics, imp_physics_gfdl, imp_physics_thompson, & - imp_physics_zhao_carr, imp_physics_zhao_carr_pdf, convert_dry_rho, dtf, save_qc, save_qi, con_pi, dtidx, dtend,& - index_of_process_conv_trans, gq0, clw, prsl, save_tcp, con_rd, con_eps, nwfa, spechum, ldiag3d, & + subroutine GFS_suite_interstitial_4_run (im, levs, tracers_total, ntrac, ntcw, ntiw, ntclamt, & + ntrw, ntsw, ntrnc, ntsnc, ntgl, ntgnc, ntlnc, ntinc, nn, imp_physics, imp_physics_gfdl, imp_physics_thompson, & + imp_physics_zhao_carr, imp_physics_zhao_carr_pdf, convert_dry_rho, dtf, save_qc, save_qi, con_pi, dtidx, dtend, & + index_of_process_conv_trans, gq0, clw, prsl, save_tcp, con_rd, con_eps, nwfa, spechum, ldiag3d, & qdiag3d, save_lnc, save_inc, ntk, ntke, errmsg, errflg) use machine, only: kind_phys @@ -716,7 +714,7 @@ subroutine GFS_suite_interstitial_4_run (im, levs, ltaerosol, mraerosol, tracers ntsw, ntrnc, ntsnc, ntgl, ntgnc, ntlnc, ntinc, nn, imp_physics, imp_physics_gfdl, imp_physics_thompson, & imp_physics_zhao_carr, imp_physics_zhao_carr_pdf - logical, intent(in) :: ltaerosol, convert_dry_rho, mraerosol + logical, intent(in) :: convert_dry_rho real(kind=kind_phys), intent(in ) :: con_pi, dtf real(kind=kind_phys), intent(in ), dimension(:,:) :: save_qc diff --git a/physics/GFS_suite_interstitial.meta b/physics/GFS_suite_interstitial.meta index 3735b887c..2340c8770 100644 --- a/physics/GFS_suite_interstitial.meta +++ b/physics/GFS_suite_interstitial.meta @@ -1223,22 +1223,6 @@ type = logical intent = in optional = F -[ltaerosol] - standard_name = flag_for_aerosol_physics - long_name = flag for aerosol physics - units = flag - dimensions = () - type = logical - intent = in - optional = F -[mraerosol] - standard_name = flag_for_merra2_aerosol_aware_for_thompson - long_name = flag for merra2 aerosol-aware physics for thompson microphysics - units = flag - dimensions = () - type = logical - intent = in - optional = F [ntrac] standard_name = number_of_tracers long_name = number of tracers @@ -1695,22 +1679,6 @@ type = integer intent = in optional = F -[ltaerosol] - standard_name = flag_for_aerosol_physics - long_name = flag for aerosol physics - units = flag - dimensions = () - type = logical - intent = in - optional = F -[mraerosol] - standard_name = flag_for_merra2_aerosol_aware_for_thompson - long_name = flag for merra2 aerosol-aware physics for thompson microphysics - units = flag - dimensions = () - type = logical - intent = in - optional = F [tracers_total] standard_name = number_of_total_tracers long_name = total number of tracers diff --git a/physics/module_MYNNPBL_wrapper.F90 b/physics/module_MYNNPBL_wrapper.F90 index ae9c03a65..2ebb006bf 100644 --- a/physics/module_MYNNPBL_wrapper.F90 +++ b/physics/module_MYNNPBL_wrapper.F90 @@ -449,6 +449,8 @@ SUBROUTINE mynnedmf_wrapper_run( & qnc(i,k) = qgrs_cloud_droplet_num_conc(i,k) qni(i,k) = qgrs_cloud_ice_num_conc(i,k) ozone(i,k) = qgrs_ozone(i,k) + qnwfa(i,k) = 0. + qnifa(i,k) = 0. enddo enddo else diff --git a/physics/module_mp_thompson.F90 b/physics/module_mp_thompson.F90 index 4e376c6eb..08f11b61c 100644 --- a/physics/module_mp_thompson.F90 +++ b/physics/module_mp_thompson.F90 @@ -460,14 +460,20 @@ SUBROUTINE thompson_init(is_aerosol_aware_in, & ! Set module variable is_aerosol_aware/merra2_aerosol_aware is_aerosol_aware = is_aerosol_aware_in merra2_aerosol_aware = merra2_aerosol_aware_in + if (is_aerosol_aware .and. merra2_aerosol_aware) then + errmsg = 'Logic error in thompson_init: only one of the two options can be true, ' // & + 'not both: is_aerosol_aware or merra2_aerosol_aware' + errflg = 1 + return + end if if (mpirank==mpiroot) then - if (is_aerosol_aware) then - write (0,'(a)') 'Using aerosol-aware version of Thompson microphysics' - else if(merra2_aerosol_aware) then - write (0,'(a)') 'Using merra2 aerosol-aware version of Thompson microphysics' - else - write (0,'(a)') 'Using non-aerosol-aware version of Thompson microphysics' - end if + if (is_aerosol_aware) then + write (0,'(a)') 'Using aerosol-aware version of Thompson microphysics' + else if(merra2_aerosol_aware) then + write (0,'(a)') 'Using merra2 aerosol-aware version of Thompson microphysics' + else + write (0,'(a)') 'Using non-aerosol-aware version of Thompson microphysics' + end if end if micro_init = .FALSE. @@ -1026,9 +1032,7 @@ SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, nc, & REAL, DIMENSION(ims:ime, kms:kme, jms:jme), OPTIONAL, INTENT(IN):: & pii REAL, DIMENSION(ims:ime, kms:kme, jms:jme), OPTIONAL, INTENT(INOUT):: & - nc - REAL, DIMENSION(ims:ime, kms:kme, jms:jme), OPTIONAL, INTENT(INOUT):: & - nwfa, nifa + nc, nwfa, nifa REAL, DIMENSION(ims:ime, jms:jme), OPTIONAL, INTENT(IN):: nwfa2d, nifa2d REAL, DIMENSION(ims:ime, kms:kme, jms:jme), OPTIONAL, INTENT(INOUT):: & re_cloud, re_ice, re_snow @@ -1054,7 +1058,6 @@ SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, nc, & REAL, DIMENSION(ims:ime, kms:kme, jms:jme), OPTIONAL, INTENT(INOUT):: & vt_dbz_wt LOGICAL, INTENT(IN) :: first_time_step - REAL, INTENT(IN):: dt_in, dt_inner ! To support subcycling: current step and maximum number of steps INTEGER, INTENT (IN) :: istep, nsteps @@ -1179,9 +1182,9 @@ SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, nc, & else stop end if - else if (merra2_aerosol_aware .and. (.not.present(nc) .or. & - .not.present(nwfa) .or. & - .not.present(nifa) )) then + else if (merra2_aerosol_aware .and. (.not.present(nc) .or. & + .not.present(nwfa) .or. & + .not.present(nifa) )) then if (present(errmsg)) then write(errmsg, '(*(a))') 'Logic error in mp_gt_driver: provide nc, nwfa, and nifa', & ' for merra2 aerosol-aware version of Thompson microphysics' diff --git a/physics/mp_thompson.F90 b/physics/mp_thompson.F90 index 2a423fbdd..828791482 100644 --- a/physics/mp_thompson.F90 +++ b/physics/mp_thompson.F90 @@ -64,18 +64,17 @@ subroutine mp_thompson_init(ncol, nlev, con_g, con_rd, con_eps, & ! Aerosols logical, intent(in ) :: is_aerosol_aware logical, intent(in ) :: merra2_aerosol_aware - real(kind_phys), optional, intent(inout) :: nc(:,:) real(kind_phys), optional, intent(inout) :: nwfa(:,:) real(kind_phys), optional, intent(inout) :: nifa(:,:) real(kind_phys), optional, intent(inout) :: nwfa2d(:) real(kind_phys), optional, intent(inout) :: nifa2d(:) + real(kind_phys), intent(in) :: aerfld(:,:,:) ! State variables real(kind_phys), intent(in ) :: tgrs(:,:) real(kind_phys), intent(in ) :: prsl(:,:) real(kind_phys), intent(in ) :: phil(:,:) real(kind_phys), intent(in ) :: area(:) - real (kind=kind_phys), dimension(:,:,:),intent(in):: aerfld ! Cloud effective radii real(kind_phys), optional, intent( out) :: re_cloud(:,:) real(kind_phys), optional, intent( out) :: re_ice(:,:) @@ -169,7 +168,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 .or. merra2_aerosol_aware) then + if (is_aerosol_aware) then nc = nc/(1.0_kind_phys-spechum) nwfa = nwfa/(1.0_kind_phys-spechum) nifa = nifa/(1.0_kind_phys-spechum) @@ -384,18 +383,17 @@ subroutine mp_thompson_run(ncol, nlev, con_g, con_rd, & ! Aerosols logical, intent(in) :: is_aerosol_aware, reset_dBZ logical, intent(in) :: merra2_aerosol_aware - ! The following arrays are not allocated if is_aerosol_aware is false real(kind_phys), optional, intent(inout) :: nc(:,:) real(kind_phys), optional, intent(inout) :: nwfa(:,:) real(kind_phys), optional, intent(inout) :: nifa(:,:) real(kind_phys), optional, intent(in ) :: nwfa2d(:) real(kind_phys), optional, intent(in ) :: nifa2d(:) + real(kind_phys), intent(in) :: aerfld(:,:,:) ! State variables and timestep information real(kind_phys), intent(inout) :: tgrs(:,:) real(kind_phys), intent(in ) :: prsl(:,:) real(kind_phys), intent(in ) :: phii(:,:) real(kind_phys), intent(in ) :: omega(:,:) - real (kind=kind_phys), dimension(:,:,:),intent(in):: aerfld real(kind_phys), intent(in ) :: dtp logical, intent(in ) :: first_time_step integer, intent(in ) :: istep, nsteps @@ -518,6 +516,7 @@ subroutine mp_thompson_run(ncol, nlev, con_g, con_rd, & errflg = 1 return end if + ! Set reduced time step if subcycling is used if (nsteps>1) then dtstep = dtp/real(nsteps, kind=kind_phys) @@ -541,7 +540,7 @@ subroutine mp_thompson_run(ncol, nlev, con_g, con_rd, & ' nc, nwfa, nifa, nwfa2d, nifa2d' errflg = 1 return - else if (merra2_aerosol_aware .and. .not. (present(nc) .and. & + else if (merra2_aerosol_aware .and. .not. (present(nc) .and. & present(nwfa) .and. & present(nifa) )) then write(errmsg,fmt='(*(a))') 'Logic error in mp_thompson_run:', & @@ -552,7 +551,6 @@ subroutine mp_thompson_run(ncol, nlev, con_g, con_rd, & end if end if - !> - 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). @@ -699,8 +697,8 @@ subroutine mp_thompson_run(ncol, nlev, con_g, con_rd, & if (merra2_aerosol_aware) then call get_niwfa(aerfld, nifa, nwfa, ncol, nlev) end if - !> - Call mp_gt_driver() with or without aerosols - if (is_aerosol_aware .or. merra2_aerosol_aware) then + !> - Call mp_gt_driver() with or without aerosols, with or without effective radii, ... + if (is_aerosol_aware) then if (do_effective_radii) then call mp_gt_driver(qv=qv, qc=qc, qr=qr, qi=qi, qs=qs, qg=qg, ni=ni, nr=nr, & nc=nc, nwfa=nwfa, nifa=nifa, nwfa2d=nwfa2d, nifa2d=nifa2d, & @@ -781,7 +779,88 @@ subroutine mp_thompson_run(ncol, nlev, con_g, con_rd, & qiten3=qiten3, niten3=niten3, nrten3=nrten3, ncten3=ncten3, & qcten3=qcten3) end if - else + else if (merra2_aerosol_aware) then + if (do_effective_radii) then + call mp_gt_driver(qv=qv, qc=qc, qr=qr, qi=qi, qs=qs, qg=qg, ni=ni, nr=nr, & + nc=nc, nwfa=nwfa, nifa=nifa, & + tt=tgrs, p=prsl, w=w, dz=dz, dt_in=dtstep, dt_inner=dt_inner, & + rainnc=rain_mp, rainncv=delta_rain_mp, & + snownc=snow_mp, snowncv=delta_snow_mp, & + icenc=ice_mp, icencv=delta_ice_mp, & + graupelnc=graupel_mp, graupelncv=delta_graupel_mp, sr=sr, & + refl_10cm=refl_10cm, & + diagflag=diagflag, do_radar_ref=do_radar_ref_mp, & + re_cloud=re_cloud, re_ice=re_ice, re_snow=re_snow, & + has_reqc=has_reqc, has_reqi=has_reqi, has_reqs=has_reqs, & + rand_perturb_on=rand_perturb_on, kme_stoch=kme_stoch, & + ! DH* 2020-06-05 not passing this optional argument, see + ! comment in module_mp_thompson.F90 / mp_gt_driver + !rand_pert=rand_pert, & + ids=ids, ide=ide, jds=jds, jde=jde, kds=kds, kde=kde, & + ims=ims, ime=ime, jms=jms, jme=jme, kms=kms, kme=kme, & + its=its, ite=ite, jts=jts, jte=jte, kts=kts, kte=kte, & + reset_dBZ=reset_dBZ, istep=istep, nsteps=nsteps, & + first_time_step=first_time_step, errmsg=errmsg, errflg=errflg, & + ! Extended diagnostics + ext_diag=ext_diag, & + ! vts1=vts1, txri=txri, txrc=txrc, & + prw_vcdc=prw_vcdc, & + prw_vcde=prw_vcde, tpri_inu=tpri_inu, tpri_ide_d=tpri_ide_d, & + tpri_ide_s=tpri_ide_s, tprs_ide=tprs_ide, & + tprs_sde_d=tprs_sde_d, & + tprs_sde_s=tprs_sde_s, tprg_gde_d=tprg_gde_d, & + tprg_gde_s=tprg_gde_s, tpri_iha=tpri_iha, & + tpri_wfz=tpri_wfz, tpri_rfz=tpri_rfz, tprg_rfz=tprg_rfz, & + tprs_scw=tprs_scw, tprg_scw=tprg_scw, tprg_rcs=tprg_rcs, & + tprs_rcs=tprs_rcs, & + tprr_rci=tprr_rci, tprg_rcg=tprg_rcg, tprw_vcd_c=tprw_vcd_c, & + tprw_vcd_e=tprw_vcd_e, tprr_sml=tprr_sml, tprr_gml=tprr_gml, & + tprr_rcg=tprr_rcg, tprr_rcs=tprr_rcs, & + tprv_rev=tprv_rev, tten3=tten3, & + qvten3=qvten3, qrten3=qrten3, qsten3=qsten3, qgten3=qgten3, & + qiten3=qiten3, niten3=niten3, nrten3=nrten3, ncten3=ncten3, & + qcten3=qcten3) + else + call mp_gt_driver(qv=qv, qc=qc, qr=qr, qi=qi, qs=qs, qg=qg, ni=ni, nr=nr, & + nc=nc, nwfa=nwfa, nifa=nifa, & + tt=tgrs, p=prsl, w=w, dz=dz, dt_in=dtstep, dt_inner=dt_inner, & + rainnc=rain_mp, rainncv=delta_rain_mp, & + snownc=snow_mp, snowncv=delta_snow_mp, & + icenc=ice_mp, icencv=delta_ice_mp, & + graupelnc=graupel_mp, graupelncv=delta_graupel_mp, sr=sr, & + refl_10cm=refl_10cm, & + diagflag=diagflag, do_radar_ref=do_radar_ref_mp, & + has_reqc=has_reqc, has_reqi=has_reqi, has_reqs=has_reqs, & + rand_perturb_on=rand_perturb_on, kme_stoch=kme_stoch, & + ! DH* 2020-06-05 not passing this optional argument, see + ! comment in module_mp_thompson.F90 / mp_gt_driver + !rand_pert=rand_pert, & + ids=ids, ide=ide, jds=jds, jde=jde, kds=kds, kde=kde, & + ims=ims, ime=ime, jms=jms, jme=jme, kms=kms, kme=kme, & + its=its, ite=ite, jts=jts, jte=jte, kts=kts, kte=kte, & + reset_dBZ=reset_dBZ, istep=istep, nsteps=nsteps, & + first_time_step=first_time_step, errmsg=errmsg, errflg=errflg, & + ! Extended diagnostics + ext_diag=ext_diag, & + ! vts1=vts1, txri=txri, txrc=txrc, & + prw_vcdc=prw_vcdc, & + prw_vcde=prw_vcde, tpri_inu=tpri_inu, tpri_ide_d=tpri_ide_d, & + tpri_ide_s=tpri_ide_s, tprs_ide=tprs_ide, & + tprs_sde_d=tprs_sde_d, & + tprs_sde_s=tprs_sde_s, tprg_gde_d=tprg_gde_d, & + tprg_gde_s=tprg_gde_s, tpri_iha=tpri_iha, & + tpri_wfz=tpri_wfz, tpri_rfz=tpri_rfz, tprg_rfz=tprg_rfz, & + tprs_scw=tprs_scw, tprg_scw=tprg_scw, tprg_rcs=tprg_rcs, & + tprs_rcs=tprs_rcs, & + tprr_rci=tprr_rci, tprg_rcg=tprg_rcg, tprw_vcd_c=tprw_vcd_c, & + tprw_vcd_e=tprw_vcd_e, tprr_sml=tprr_sml, tprr_gml=tprr_gml, & + tprr_rcg=tprr_rcg, tprr_rcs=tprr_rcs, & + tprv_rev=tprv_rev, tten3=tten3, & + qvten3=qvten3, qrten3=qrten3, qsten3=qsten3, qgten3=qgten3, & + qiten3=qiten3, niten3=niten3, nrten3=nrten3, ncten3=ncten3, & + qcten3=qcten3) + end if + else ! neither is_aerosol_aware nor merra2_aerosol_aware if (do_effective_radii) then call mp_gt_driver(qv=qv, qc=qc, qr=qr, qi=qi, qs=qs, qg=qg, ni=ni, nr=nr, & tt=tgrs, p=prsl, w=w, dz=dz, dt_in=dtstep, dt_inner=dt_inner, & @@ -970,15 +1049,17 @@ subroutine mp_thompson_finalize(errmsg, errflg) end subroutine mp_thompson_finalize subroutine get_niwfa(aerfld, nifa, nwfa, ncol, nlev) - implicit none - integer, intent(in)::ncol, nlev - real (kind=kind_phys), dimension(:,:,:),intent(in):: aerfld - real (kind=kind_phys), dimension(:,:),intent(out)::nifa, nwfa - nifa=(aerfld(:,:,1)/4.0737762+aerfld(:,:,2)/30.459203+aerfld(:,:,3)/153.45048+ & - 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 + implicit none + integer, intent(in)::ncol, nlev + real (kind=kind_phys), dimension(:,:,:), intent(in) :: aerfld + real (kind=kind_phys), dimension(:,:), intent(out ):: nifa, nwfa + + nifa=(aerfld(:,:,1)/4.0737762+aerfld(:,:,2)/30.459203+aerfld(:,:,3)/153.45048+ & + 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 end subroutine get_niwfa end module mp_thompson