| bl_mynn_cloudpdf | mynnedmf_wrapper | flag to determine which cloud PDF to use \n
- 0: use Sommeria-Deardorff subgrid cloud PDF
@@ -384,6 +387,9 @@ show some variables in the namelist that must match the SDF.
| t_sub | gfdl_cloud_microphys | temperature threshold for sublimation. Cloud ice, snow or graupel stops(starts) sublimation when temperature is lower(higher) then \p t_sub | 184.
| | mp_print | gfdl_cloud_microphys | \a .true. to turn on GFDL cloud microphysics debugging print out. (not supported in GFS physics) | .false.
| | ltaerosol | mp_thompson | flag for using aerosol climotology in Thompson MP scheme | .false.
+ | | mraerosol | | flag for merra2 aerosol aware | .false.
+ | | lradar | | flag for radar reflectivity | .false.
+ | | nsfullradar_diag | | seconds between resetting radar reflectivity calculation | .-999.0
| | ttendlim | mp_thompson | temperature tendency limiter per time step in K/s, set to < 0 to deactivate | -999.0
| | ext_diag_thompson | mp_thompson | flag for extended diagnostic output from Thompson MP | .false.
| | thompson_ext_ndiag3d | mp_thompson | number of 3d arrays for extended diagnostic output from Thompson MP | 37
@@ -394,9 +400,11 @@ show some variables in the namelist that must match the SDF.
| | cnvcld | see \a GFS_typedefs.F90 | flag for convective cloud | .false.
| | lgfdlmprad | gfs_rrtmg_pre | flag for GFDL mp scheme and radiation consistency | .false.
| | nssl_cccn | mp_nssl | CCN concentration (m^-3) | 0.6e9
- | | nssl_alphar | mp_nssl | rain PSD shape parameter | 0.0
- | | nssl_alphah | mp_nssl | graupel PSD shape parameter | 0.0
- | | nssl_alphahl | mp_nssl | hail PSD shape parameter | 1.0
+ | | nssl_alphah | mp_nssl | graupel shape parameter | 0.0
+ | | nssl_alphahl | mp_nssl | hail shape parameter | 1.0
+ | | nssl_alphar | mp_nssl | shape paramter for rain (imurain=1 only) | 0.0
+ | | nssl_ehw0_in | mp_nssl | constant or max assumed graupel-droplet collection efficiency | 0.9
+ | | nssl_ehlw0_in | mp_nssl | constant or max assumed hail-droplet collection efficiency | 0.9
| | nssl_hail_on | mp_nssl | NSSL flag to activate the hail category | .false.
| | nssl_ccn_on | mp_nssl | NSSL flag to activate the CCN category | .true.
| | nssl_invertccn | mp_nssl | NSSL flag to treat CCN as activated or unactivated | .true.
@@ -508,15 +516,19 @@ show some variables in the namelist that must match the SDF.
- 0: USGS
- 1: IGBP(20 category): IGBP must be selected if NoahMP is used
- 2: UMD (13 category)
+
- 3: NLCD40 (40 category, NOAH WRFv4 only)
+
- 4: USGS-RUC (28 category, NOAH WRFv4 only)
+
- 5: MODIS-RUC (21 category, NOAH WRFv4 only)
| 2
| | isot | lsm_noah, lsm_ruc, \ref noahmpdrv | flag for soil type dataset choice:\n
- 0: Zobler soil type (9 category)
- 1: STATSGO soil type (19 category): STATSGO must be selected if NoahMP is used
+
- 2: STAS-RUC soil type (19 category, NOAH WRFv4 only
| 0
-
+ | | exticeden | lsm_noah,lsm_ruc | flag for calculating frozen precip ice density outside of the LSM | .false.
| | iopt_dveg | \ref noahmpdrv | options for dynamic vegetation \n
- 1: off (use table LAI; use FVEG = SHDFAC from input)
@@ -606,14 +618,42 @@ show some variables in the namelist that must match the SDF.
| 1
| | iopt_trs | \ref noahmpdrv | options for thermal roughness scheme:\n
-
- - 1: z0h=z0m
-
- 2: canopy based czil
-
- 3: ec
-
- 4: kb inversed
-
- | 2
+
+ - 1: z0h=z0m
+
- 2: canopy based czil
+
- 3: ec
+
- 4: kb inversed
+
+ | 2
+ | | mosaic_lu | \ref lsm_ruc | control for use of fractional landuse in RUC land surface model
+
+ - 1: use of fractional landuse in RUC LSM
+
+ | 0
+ | | mosaic_soil | \ref lsm_ruc | control for use of fractional soil in RUC land surface model
+
+ - 1: use of fractional soil in RUC LSM
+
+ | 0
+ | | isncond_opt | \ref lsm_ruc | control for soil thermal conductivity option in RUC land surface model
+
+ - 1: constant
+
- 2: Sturm et al.(1997) \cite sturm_1997
+
+ | 1
+ | | isncovr_opt | \ref lsm_ruc | control for snow cover fraction option in RUC land surface model
+
+ - 1: original formulation using threshold snow depth to compute snow fraction
+
- 2: Niu-Yang(2007) \cite Niu_2007
+
- 3: updated Niu-Yang similar to Noah MP
+
+ | 1
| | \b Parameters \b related \b to \b other \b surface \b scheme \b options
+ | | lkm | \ref clm_lake | 0: no lake; 1: lake; 2: lake&nsst | 0
+ | | iopt_lake | \ref clm_lake | 1: flake; 2: CLM lake | 2
+ | | lakedepth_threshold | \ref clm_lake | lakedepth must be greater than this value to enable a lake model | 1.0
+ | | lakefrac_threshold | \ref clm_lake | lakefrac must be greater than this value to enable a lake model | 0.0
+ | | use_lake2m | \ref clm_lake | use 2m T&Q from CLM lake model | .false.
| | nstf_name(5) | sfc_nst | NSST related paramters:\n
- nstf_name(1): 0=NSST off, 1= NSST on but uncoupled, 2= NSST on and coupled
diff --git a/physics/module_mp_thompson.F90 b/physics/module_mp_thompson.F90
index 0c708cb3d..271db11d0 100644
--- a/physics/module_mp_thompson.F90
+++ b/physics/module_mp_thompson.F90
@@ -708,9 +708,9 @@ SUBROUTINE thompson_init(is_aerosol_aware_in, &
dtc(n) = (Dc(n) - Dc(n-1))
enddo
-!> - Create bins of cloud ice (from min diameter up to 5x min snow size)
+!> - Create bins of cloud ice (from min diameter up to 2x min snow size)
xDx(1) = D0i*1.0d0
- xDx(nbi+1) = 5.0d0*D0s
+ xDx(nbi+1) = 2.0d0*D0s
do n = 2, nbi
xDx(n) = DEXP(DFLOAT(n-1)/DFLOAT(nbi) &
*DLOG(xDx(nbi+1)/xDx(1)) +DLOG(xDx(1)))
@@ -1046,7 +1046,7 @@ SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, nc, &
INTEGER, INTENT(IN) :: rand_perturb_on, kme_stoch, n_var_spp
REAL, DIMENSION(:,:), INTENT(IN) :: rand_pert
REAL, DIMENSION(:), INTENT(IN) :: spp_prt_list, spp_stddev_cutoff
- CHARACTER(len=3), DIMENSION(:), INTENT(IN) :: spp_var_list
+ CHARACTER(len=10), DIMENSION(:), INTENT(IN) :: spp_var_list
INTEGER, INTENT(IN):: has_reqc, has_reqi, has_reqs
#if ( WRF_CHEM == 1 )
REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT):: &
@@ -1509,6 +1509,14 @@ SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, nc, &
enddo
endif
+ if (merra2_aerosol_aware) then
+ do k = kts, kte
+ nc(i,k,j) = nc1d(k)
+ nwfa(i,k,j) = nwfa1d(k)
+ nifa(i,k,j) = nifa1d(k)
+ enddo
+ endif
+
do k = kts, kte
qv(i,k,j) = qv1d(k)
qc(i,k,j) = qc1d(k)
@@ -2822,7 +2830,7 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, &
prr_rcg(k) = MIN(DBLE(rg(k)*odts), prr_rcg(k))
prg_rcg(k) = -prr_rcg(k)
!> - Put in explicit drop break-up due to collisions.
- pnr_rcg(k) = -5.*tnr_gacr(idx_g1,idx_g,idx_r1,idx_r) ! RAIN2M
+ pnr_rcg(k) = -1.5*tnr_gacr(idx_g1,idx_g,idx_r1,idx_r) ! RAIN2M
endif
endif
endif
@@ -3053,16 +3061,15 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, &
if (prr_sml(k) .gt. 0.) then
prr_sml(k) = prr_sml(k) + 4218.*olfus*tempc &
* (prr_rcs(k)+prs_scw(k))
- endif
- prr_sml(k) = MIN(DBLE(rs(k)*odts), MAX(0.D0, prr_sml(k)))
- pnr_sml(k) = smo0(k)/rs(k)*prr_sml(k) * 10.0**(-0.25*tempc) ! RAIN2M
- pnr_sml(k) = MIN(DBLE(smo0(k)*odts), pnr_sml(k))
-
- if (ssati(k).lt. 0.) then
- prs_sde(k) = C_cube*t1_subl*diffu(k)*ssati(k)*rvs &
- * (t1_qs_sd*smo1(k) &
- + t2_qs_sd*rhof2(k)*vsc2(k)*smof(k))
- prs_sde(k) = MAX(DBLE(-rs(k)*odts), prs_sde(k))
+ prr_sml(k) = MIN(DBLE(rs(k)*odts), prr_sml(k))
+ pnr_sml(k) = smo0(k)/rs(k)*prr_sml(k) * 10.0**(-0.25*tempc) ! RAIN2M
+ pnr_sml(k) = MIN(DBLE(smo0(k)*odts), pnr_sml(k))
+ elseif (ssati(k).lt. 0.) then
+ prr_sml(k) = 0.0
+ prs_sde(k) = C_cube*t1_subl*diffu(k)*ssati(k)*rvs &
+ * (t1_qs_sd*smo1(k) &
+ + t2_qs_sd*rhof2(k)*vsc2(k)*smof(k))
+ prs_sde(k) = MAX(DBLE(-rs(k)*odts), prs_sde(k))
endif
endif
@@ -3070,17 +3077,16 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, &
prr_gml(k) = (tempc*tcond(k)-lvap0*diffu(k)*delQvs(k)) &
* N0_g(k)*(t1_qg_me*ilamg(k)**cge(10) &
+ t2_qg_me*rhof2(k)*vsc2(k)*ilamg(k)**cge(11))
-!-GT prr_gml(k) = prr_gml(k) + 4218.*olfus*tempc &
-!-GT * (prr_rcg(k)+prg_gcw(k))
- prr_gml(k) = MIN(DBLE(rg(k)*odts), MAX(0.D0, prr_gml(k)))
- pnr_gml(k) = N0_g(k)*cgg(2)*ilamg(k)**cge(2) / rg(k) & ! RAIN2M
- * prr_gml(k) * 10.0**(-0.5*tempc)
-
- if (ssati(k).lt. 0.) then
- prg_gde(k) = C_cube*t1_subl*diffu(k)*ssati(k)*rvs &
- * N0_g(k) * (t1_qg_sd*ilamg(k)**cge(10) &
- + t2_qg_sd*vsc2(k)*rhof2(k)*ilamg(k)**cge(11))
- prg_gde(k) = MAX(DBLE(-rg(k)*odts), prg_gde(k))
+ if (prr_gml(k) .gt. 0.) then
+ prr_gml(k) = MIN(DBLE(rg(k)*odts), prr_gml(k))
+ pnr_gml(k) = N0_g(k)*cgg(2)*ilamg(k)**cge(2) / rg(k) & ! RAIN2M
+ * prr_gml(k) * 10.0**(-0.5*tempc)
+ elseif (ssati(k).lt. 0.) then
+ prr_gml(k) = 0.0
+ prg_gde(k) = C_cube*t1_subl*diffu(k)*ssati(k)*rvs &
+ * N0_g(k) * (t1_qg_sd*ilamg(k)**cge(10) &
+ + t2_qg_sd*vsc2(k)*rhof2(k)*ilamg(k)**cge(11))
+ prg_gde(k) = MAX(DBLE(-rg(k)*odts), prg_gde(k))
endif
endif
diff --git a/physics/module_sf_mynn.F90 b/physics/module_sf_mynn.F90
index c60247cf6..eecc5493c 100644
--- a/physics/module_sf_mynn.F90
+++ b/physics/module_sf_mynn.F90
@@ -106,6 +106,7 @@ MODULE module_sf_mynn
REAL(kind_phys), DIMENSION(0:1000 ),SAVE :: psim_stab,psim_unstab, &
psih_stab,psih_unstab
+!$acc declare create(psim_stab, psim_unstab, psih_stab, psih_unstab)
CONTAINS
@@ -371,6 +372,20 @@ SUBROUTINE SFCLAY_mynn( &
errflg = 0
errmsg = ''
+!$acc enter data copyin( dz8w,U3D,V3D,QV3D,QC3D,P3D,T3D, &
+!$acc pattern_spp_sfc)
+
+!$acc enter data copyin( UST_WAT(:), UST_LND(:), UST_ICE(:), &
+!$acc MOL(:), QFLX(:), HFLX(:), &
+!$acc QSFC(:), QSFC_WAT(:), QSFC_LND(:), &
+!$acc QSFC_ICE(:))
+
+!$acc enter data create( dz8w1d(:), dz2w1d(:), U1D(:), &
+!$acc V1D(:), U1D2(:), V1D2(:), &
+!$acc QV1D(:), QC1D(:), P1D(:), &
+!$acc T1D(:), rstoch1D(:), qstar(:))
+
+
IF (debug_code >= 1) THEN
write(*,*)"======= printing of constants:"
write(*,*)"cp=", cp," g=", grav
@@ -382,6 +397,10 @@ SUBROUTINE SFCLAY_mynn( &
itf=ite !MIN0(ite,ide-1)
ktf=kte !MIN0(kte,kde-1)
+!$acc parallel loop present(dz8w,U3D,V3D,QV3D,QC3D,P3D,T3D, &
+!$acc pattern_spp_sfc,dz8w1d,dz2w1d,U1D, &
+!$acc V1D,U1D2,V1D2,QV1D,QC1D,P1D,T1D, &
+!$acc rstoch1D,qstar)
DO i=its,ite
dz8w1d(I) = dz8w(i,kts)
dz2w1d(I) = dz8w(i,kts+1)
@@ -403,6 +422,9 @@ SUBROUTINE SFCLAY_mynn( &
ENDDO
IF (itimestep==1 .AND. iter==1) THEN
+!$acc parallel loop present(U1D,V1D,UST_WAT,UST_LND,UST_ICE,MOL, &
+!$acc QFLX,HFLX,QV3D,QSFC,QSFC_WAT, &
+!$acc QSFC_LND,QSFC_ICE)
DO i=its,ite
IF (.not. flag_restart) THEN
!Everything here is used before calculated
@@ -432,6 +454,9 @@ SUBROUTINE SFCLAY_mynn( &
ENDDO
ENDIF
+!$acc exit data delete( dz8w,U3D,V3D,QV3D,QC3D,P3D,T3D, &
+!$acc pattern_spp_sfc, QC1D)
+
CALL SFCLAY1D_mynn(flag_iter, &
J,U1D,V1D,T1D,QV1D,P1D,dz8w1d, &
U1D2,V1D2,dz2w1d, &
@@ -471,6 +496,16 @@ SUBROUTINE SFCLAY_mynn( &
its,ite, jts,jte, kts,kte, &
errmsg, errflg )
+!$acc exit data copyout( UST_WAT(:), UST_LND(:), UST_ICE(:), &
+!$acc MOL(:), QFLX(:), HFLX(:), &
+!$acc QSFC(:), QSFC_WAT(:), QSFC_LND(:), &
+!$acc QSFC_ICE(:))
+
+!$acc exit data delete( dz8w1d(:), dz2w1d(:), U1D(:), &
+!$acc V1D(:), U1D2(:), V1D2(:), &
+!$acc QV1D(:), T1D(:), P1D(:), &
+!$acc rstoch1D(:), qstar(:))
+
END SUBROUTINE SFCLAY_MYNN
!-------------------------------------------------------------------
@@ -629,6 +664,22 @@ SUBROUTINE SFCLAY1D_mynn(flag_iter, &
character(len=*), intent(inout) :: errmsg
integer, intent(inout) :: errflg
+! Local fixed-size errmsg character array for error messages on accelerator
+! devices distinct from the host (e.g. GPUs). Necessary since OpenACC does
+! not support assumed-size (len=*) arrays like errmsg. Additional
+! device_errflg integer to denote when device_errmsg needs to be synced
+! with errmsg.
+ character(len=512) :: device_errmsg
+ integer :: device_errflg
+
+! Special versions of the fixed-size errmsg character array for error messages
+! on the device and it's errflag counterpart. These are necessary to ensure
+! the return statements at lines 1417 and 2030 are executed only for this
+! special case, and not any and all error messages set on the device.
+ character(len=512) :: device_special_errmsg
+ integer :: device_special_errflg
+
+
!----------------------------------------------------------------
! LOCAL VARS
!----------------------------------------------------------------
@@ -678,7 +729,65 @@ SUBROUTINE SFCLAY1D_mynn(flag_iter, &
! Initialize error-handling
errflg = 0
errmsg = ''
+ device_errflg = errflg
+ device_errmsg = errmsg
+ device_special_errflg = errflg
+ device_special_errmsg = errmsg
!-------------------------------------------------------------------
+!$acc update device(psim_stab, psim_unstab, psih_stab, psih_unstab)
+
+!$acc enter data create( ZA, ZA2, THV1D, TH1D, TC1D, TV1D, &
+!$acc RHO1D, QVSH, PSIH2, PSIM10, PSIH10, WSPDI, &
+!$acc GOVRTH, PSFC, THCON, &
+!$acc zratio_lnd, zratio_ice, zratio_wat, &
+!$acc TSK_lnd, TSK_ice, TSK_wat, &
+!$acc THSK_lnd, THSK_ice, THSK_wat, &
+!$acc THVSK_lnd, THVSK_ice, THVSK_wat, &
+!$acc GZ1OZ0_lnd, GZ1OZ0_ice, GZ1OZ0_wat, &
+!$acc GZ1OZt_lnd, GZ1OZt_ice, GZ1OZt_wat, &
+!$acc GZ2OZ0_lnd, GZ2OZ0_ice, GZ2OZ0_wat, &
+!$acc GZ2OZt_lnd, GZ2OZt_ice, GZ2OZt_wat, &
+!$acc GZ10OZ0_lnd, GZ10OZ0_ice, GZ10OZ0_wat, &
+!$acc GZ10OZt_lnd, GZ10OZt_ice, GZ10OZt_wat, &
+!$acc ZNTstoch_lnd, ZNTstoch_ice, ZNTstoch_wat, &
+!$acc ZT_lnd, ZT_ice, ZT_wat, &
+!$acc ZQ_lnd, ZQ_ice, ZQ_wat, &
+!$acc PSIQ_lnd, PSIQ_ice, PSIQ_wat, &
+!$acc PSIQ2_lnd, PSIQ2_ice, PSIQ2_wat, &
+!$acc QSFCMR_lnd, QSFCMR_ice, QSFCMR_wat )
+
+!$acc enter data copyin(flag_iter, dry, wet, icy, CPM, MAVAIL, &
+!$acc QFX, FLHC, FLQC, CHS, CH, CHS2, CQS2, USTM, &
+!$acc HFX, LH, wstar, qstar, PBLH, ZOL, MOL, RMOL, &
+!$acc T2, TH2, Q2, QV1D, PSFCPA, &
+!$acc WSPD, U10, V10, U1D, V1D, U1D2, V1D2, &
+!$acc T1D, P1D, rstoch1D, sigmaf, &
+!$acc shdmax, vegtype, z0pert, ztpert, dx, QGH, &
+!$acc dz2w1d, dz8w1d, &
+!$acc stress_wat, stress_lnd, stress_ice, &
+!$acc rb_wat, rb_lnd, rb_ice, &
+!$acc tskin_wat, tskin_lnd, tskin_ice, &
+!$acc tsurf_wat, tsurf_lnd, tsurf_ice, &
+!$acc psim, psih, &
+!$acc UST_wat, UST_lnd, UST_ice, &
+!$acc ZNT_wat, ZNT_lnd, ZNT_ice, &
+!$acc QSFC, QSFC_lnd, QSFC_wat, QSFC_ice, &
+!$acc QFLX, QFLX_lnd, QFLX_wat, QFLX_ice, &
+!$acc HFLX, HFLX_lnd, HFLX_wat, HFLX_ice, &
+!$acc PSIX_wat, PSIX_lnd, PSIX_ice, &
+!$acc PSIX10_wat, PSIX10_lnd, PSIX10_ice, &
+!$acc PSIT2_lnd, PSIT2_wat, PSIT2_ice, &
+!$acc PSIT_lnd, PSIT_wat, PSIT_ice, &
+!$acc ch_lnd, ch_wat, ch_ice, &
+!$acc cm_lnd, cm_wat, cm_ice, &
+!$acc snowh_lnd, snowh_wat, snowh_ice, &
+!$acc device_errmsg, device_errflg, &
+!$acc device_special_errmsg, device_special_errflg)
+
+!$acc parallel loop present(PSFCPA, PSFC, QSFC, T1D, flag_iter, tsurf_lnd, &
+!$acc QSFC_wat, QSFCMR_wat, wet, TSK_wat, tskin_wat, &
+!$acc QSFC_lnd, QSFCMR_lnd, dry, TSK_lnd, tskin_lnd, &
+!$acc QSFC_ice, QSFCMR_ice, icy, TSK_ice, tskin_ice)
DO I=its,ite
! PSFC ( in cmb) is used later in saturation checks
@@ -791,6 +900,10 @@ SUBROUTINE SFCLAY1D_mynn(flag_iter, &
endif ! flag_iter
ENDDO
+!$acc serial present(pblh, PSFCPA, dz8w1d, qflx, hflx, &
+!$acc dry, tskin_lnd, tsurf_lnd, qsfc_lnd, znt_lnd, ust_lnd, snowh_lnd, &
+!$acc icy, tskin_ice, tsurf_ice, qsfc_ice, znt_ice, ust_ice, snowh_ice, &
+!$acc wet, tskin_wat, tsurf_wat, qsfc_wat, znt_wat, ust_wat, snowh_wat)
IF (debug_code >= 1) THEN
write(0,*)"ITIMESTEP=",ITIMESTEP," iter=",iter
DO I=its,ite
@@ -815,7 +928,12 @@ SUBROUTINE SFCLAY1D_mynn(flag_iter, &
ENDIF
ENDDO
ENDIF
+!$acc end serial
+!$acc parallel loop present(PSFC, PSFCPA, QVSH, QV1D, THCON, flag_iter, &
+!$acc dry, tskin_lnd, TSK_lnd, tsurf_lnd, THSK_lnd, THVSK_lnd, qsfc_lnd, &
+!$acc icy, tskin_ice, TSK_ice, tsurf_ice, THSK_ice, THVSK_ice, qsfc_ice, &
+!$acc wet, tskin_wat, TSK_wat, tsurf_wat, THSK_wat, THVSK_wat, qsfc_wat)
DO I=its,ite
! PSFC ( in cmb) is used later in saturation checks
PSFC(I)=PSFCPA(I)/1000.
@@ -853,18 +971,21 @@ SUBROUTINE SFCLAY1D_mynn(flag_iter, &
endif ! flag_iter
ENDDO
+!$acc parallel loop present(TH1D, T1D, P1D, TC1D)
DO I=its,ite
! CONVERT LOWEST LAYER TEMPERATURE TO POTENTIAL TEMPERATURE:
TH1D(I)=T1D(I)*(100000./P1D(I))**ROVCP !(Theta, K)
TC1D(I)=T1D(I)-273.15 !(T, Celsius)
ENDDO
+!$acc parallel loop present(THV1D, TH1D, QVSH, TV1D, T1D)
DO I=its,ite
! CONVERT TO VIRTUAL TEMPERATURE
THV1D(I)=TH1D(I)*(1.+EP1*QVSH(I)) !(K)
TV1D(I)=T1D(I)*(1.+EP1*QVSH(I)) !(K)
ENDDO
+!$acc parallel loop present(RHO1D, P1D, TV1D, TH1D, ZA, ZA2, dz2w1d, dz8w1d, GOVRTH)
DO I=its,ite
RHO1D(I)=P1D(I)/(Rd*TV1D(I)) !now using value calculated in sfc driver
ZA(I)=0.5*dz8w1d(I) !height of first half-sigma level
@@ -873,11 +994,16 @@ SUBROUTINE SFCLAY1D_mynn(flag_iter, &
ENDDO
!tgs - should QFX and HFX be separate for land, ice and water?
+!$acc parallel loop present(QFX, QFLX, RHO1D, HFX, HFLX)
DO I=its,ite
QFX(i)=QFLX(i)*RHO1D(I)
HFX(i)=HFLX(i)*RHO1D(I)*cp
ENDDO
+!$acc serial present(THV1D, TV1D, RHO1D, GOVRTH, &
+!$acc dry, tsk_lnd, thvsk_lnd, &
+!$acc icy, tsk_ice, thvsk_ice, &
+!$acc wet, tsk_wat, thvsk_wat)
IF (debug_code ==2) THEN
!write(*,*)"ITIMESTEP=",ITIMESTEP
DO I=its,ite
@@ -890,7 +1016,9 @@ SUBROUTINE SFCLAY1D_mynn(flag_iter, &
write(*,*)"RHO1D=", RHO1D(i)," GOVRTH=",GOVRTH(i)
ENDDO
ENDIF
+!$acc end serial
+!$acc parallel loop present(T1D,P1D,QGH,QV1D,CPM)
DO I=its,ite
! QGH CHANGED TO USE LOWEST-LEVEL AIR TEMP
! Q2SAT = QGH IN LSM
@@ -908,6 +1036,10 @@ SUBROUTINE SFCLAY1D_mynn(flag_iter, &
CPM(I)=CP*(1.+0.84*QV1D(I))
ENDDO
+!$acc serial present(QGH, &
+!$acc wet, QSFC_wat, QSFCMR_wat, &
+!$acc dry, QSFC_lnd, QSFCMR_lnd, &
+!$acc icy, QSFC_ice, QSFCMR_ice)
IF (debug_code == 2) THEN
write(*,*)"ITIMESTEP=",ITIMESTEP
DO I=its,ite
@@ -925,7 +1057,13 @@ SUBROUTINE SFCLAY1D_mynn(flag_iter, &
endif
ENDDO
ENDIF
+!$acc end serial
+!$acc parallel loop present(flag_iter,U1D,V1D,WSPD,wet,dry,icy, &
+!$acc THV1D,THVSK_wat,THVSK_lnd,THVSK_ice, &
+!$acc hfx,RHO1D,qfx,WSTAR,pblh,dx,GOVRTH,ZA, &
+!$acc TSK_wat,TSK_lnd,TSK_ice, &
+!$acc rb_wat,rb_lnd,rb_ice)
DO I=its,ite
if( flag_iter(i) ) then
! DH* 20200401 - note. A weird bug in Intel 18 on hera prevents using the
@@ -1067,6 +1205,35 @@ SUBROUTINE SFCLAY1D_mynn(flag_iter, &
!--------------------------------------------------------------------
!--------------------------------------------------------------------
+!$acc parallel loop present(flag_iter, PSFCPA, dz8w1d, pblh, &
+!$acc device_errmsg, device_errflg, &
+!$acc device_special_errmsg, device_special_errflg, &
+!$acc wet, dry, icy, &
+!$acc ZT_wat, ZT_lnd, ZT_ice, &
+!$acc ZNT_wat, ZNT_lnd, ZNT_ice, &
+!$acc ZNTstoch_wat, ZNTstoch_lnd, ZNTstoch_ice, &
+!$acc UST_wat, UST_lnd, UST_ice, &
+!$acc ZQ_wat, ZQ_lnd, ZQ_ice, &
+!$acc snowh_wat, snowh_lnd, snowh_ice, &
+!$acc THVSK_wat, THVSK_lnd, THVSK_ice, &
+!$acc tskin_wat, tskin_lnd, tskin_ice, &
+!$acc tsurf_wat, tsurf_lnd, tsurf_ice, &
+!$acc qsfc_wat, qsfc_lnd, qsfc_ice, &
+!$acc GZ1OZ0_wat, GZ1OZt_wat, GZ2OZ0_wat, GZ2OZt_wat, GZ10OZ0_wat, GZ10OZt_wat, &
+!$acc GZ1OZ0_lnd, GZ1OZt_lnd, GZ2OZ0_lnd, GZ2OZt_lnd, GZ10OZ0_lnd, GZ10OZt_lnd, &
+!$acc GZ1OZ0_ice, GZ1OZt_ice, GZ2OZ0_ice, GZ2OZt_ice, GZ10OZ0_ice, GZ10OZt_ice, &
+!$acc zratio_wat, zratio_lnd, zratio_ice, &
+!$acc stress_wat, stress_lnd, stress_ice, &
+!$acc rb_wat, rb_lnd, rb_ice, &
+!$acc qflx, qflx_lnd, &
+!$acc hflx, hflx_lnd, &
+!$acc psim, psih, psim10, psih10, psih2, &
+!$acc psix_wat, psix10_wat, psit_wat, psit2_wat, psiq_wat, psiq2_wat, &
+!$acc psix_lnd, psix10_lnd, psit_lnd, psit2_lnd, psiq_lnd, psiq2_lnd, &
+!$acc psix_ice, psix10_ice, psit_ice, psit2_ice, psiq_ice, psiq2_ice, &
+!$acc WSPD, WSPDI, U1D, V1D, TC1D, THV1D, rstoch1D, USTM, ZA, ZOL, QVSH, &
+!$acc shdmax, vegtype, z0pert, ztpert, mol, rmol, wstar, qstar, sigmaf)
+
DO I=its,ite
if( flag_iter(i) ) then
@@ -1082,10 +1249,12 @@ SUBROUTINE SFCLAY1D_mynn(flag_iter, &
if (sfc_z0_type >= 0) then ! Avoid calculation is using wave model
! CALCULATE z0 (znt)
!--------------------------------------
+
IF (debug_code == 2) THEN
write(*,*)"=============Input to ZNT over water:"
write(*,*)"u*:",UST_wat(i)," wspd=",WSPD(i)," visc=",visc," za=",ZA(I)
ENDIF
+
IF ( PRESENT(ISFTCFLX) ) THEN
IF ( ISFTCFLX .EQ. 0 ) THEN
IF (COARE_OPT .EQ. 3.0) THEN
@@ -1170,7 +1339,7 @@ SUBROUTINE SFCLAY1D_mynn(flag_iter, &
ENDIF
ELSEIF ( ISFTCFLX .EQ. 4 ) THEN
!GFS zt formulation
- CALL GFS_zt_wat(ZT_wat(i),ZNTstoch_wat(i),restar,WSPD(i),ZA(i),sfc_z0_type,errmsg,errflg)
+ CALL GFS_zt_wat(ZT_wat(i),ZNTstoch_wat(i),restar,WSPD(i),ZA(i),sfc_z0_type,device_errmsg,device_errflg)
ZQ_wat(i)=ZT_wat(i)
ENDIF
ELSE
@@ -1183,6 +1352,7 @@ SUBROUTINE SFCLAY1D_mynn(flag_iter, &
rstoch1D(i),spp_sfc)
ENDIF
ENDIF
+
IF (debug_code > 1) THEN
write(*,*)"=============Output ZT & ZQ over water:"
write(*,*)"ZT:",ZT_wat(i)," ZQ:",ZQ_wat(i)
@@ -1230,9 +1400,16 @@ SUBROUTINE SFCLAY1D_mynn(flag_iter, &
ELSEIF ( IZ0TLND .EQ. 2 ) THEN
! DH note - at this point, qstar is either not initialized
! or initialized to zero, but certainly not set correctly
- errmsg = 'Logic error: qstar is not set correctly when calling Yang_2008'
- errflg = 1
+ device_special_errmsg = 'Logic error: qstar is not set correctly when calling Yang_2008'
+ device_special_errflg = 1
+#ifndef _OPENACC
+! Necessary since OpenACC does not support branching in parallel code
+! Must sync errmsg and errflg with device_errmsg and device_errflg, respectively
+! so that proper error message and error flag codes are returned.
+ errmsg = device_special_errmsg
+ errflg = device_special_errflg
return
+#endif
CALL Yang_2008(ZNTSTOCH_lnd(i),ZT_lnd(i),ZQ_lnd(i),UST_lnd(i),MOL(I),&
qstar(I),restar,visc)
ELSEIF ( IZ0TLND .EQ. 3 ) THEN
@@ -1249,6 +1426,7 @@ SUBROUTINE SFCLAY1D_mynn(flag_iter, &
UST_lnd(I),KARMAN,1.0_kind_phys,0,spp_sfc,rstoch1D(i))
ENDIF
ENDIF
+
IF (ZNTstoch_lnd(i) < 1E-8 .OR. Zt_lnd(i) < 1E-10) THEN
write(0,*)"===(land) capture bad input in mynn sfc layer, i=:",i
write(0,*)" ZNT=", ZNTstoch_lnd(i)," ZT=",Zt_lnd(i)
@@ -1258,7 +1436,6 @@ SUBROUTINE SFCLAY1D_mynn(flag_iter, &
" dz=",dz8w1d(i)," qflx=",qflx_lnd(i)," hflx=",hflx_lnd(i)," hpbl=",pblh(i)
ENDIF
-
GZ1OZ0_lnd(I)= LOG((ZA(I)+ZNTstoch_lnd(I))/ZNTstoch_lnd(I))
GZ1OZt_lnd(I)= LOG((ZA(I)+ZNTstoch_lnd(i))/ZT_lnd(i))
GZ2OZ0_lnd(I)= LOG((2.0+ZNTstoch_lnd(I))/ZNTstoch_lnd(I))
@@ -1821,6 +1998,26 @@ SUBROUTINE SFCLAY1D_mynn(flag_iter, &
endif ! flag_iter
ENDDO ! end i-loop
+#ifdef _OPENACC
+! Necessary since OpenACC does not support branching in parallel code.
+! Must sync host errflg, errmsg to determine if return must be triggered
+! and correct error message and error flag code returned.
+! This code is being executed on the HOST side only, pulling data from DEVICE.
+!$acc exit data copyout(device_special_errflg, device_special_errmsg)
+ IF (device_special_errflg /= 0) THEN
+ errflg = device_special_errflg
+ errmsg = device_special_errmsg
+ return
+ ENDIF
+#endif
+
+!$acc serial present(wet, dry, icy, &
+!$acc PSIM, PSIH, CPM, RHO1D, ZOL, wspd, MOL, &
+!$acc wstar, qstar, THV1D, HFX, MAVAIL, QVSH, &
+!$acc THVSK_wat, THVSK_lnd, THVSK_ice, &
+!$acc UST_wat, UST_lnd, UST_ice, &
+!$acc ZNTstoch_wat, ZNTstoch_lnd, ZNTstoch_ice, &
+!$acc zt_wat, zt_lnd, zt_ice)
IF (debug_code == 2) THEN
DO I=its,ite
IF(wet(i))write(*,*)"==== AT END OF MAIN LOOP, i=",i, "(wet)"
@@ -1841,10 +2038,29 @@ SUBROUTINE SFCLAY1D_mynn(flag_iter, &
write(*,*)"============================================="
ENDDO ! end i-loop
ENDIF
+!$acc end serial
!----------------------------------------------------------
! COMPUTE SURFACE HEAT AND MOISTURE FLUXES
!----------------------------------------------------------
+!$acc parallel loop present(flag_iter, dry, wet, icy, &
+!$acc QFX, HFX, FLHC, FLQC, LH, CHS, CH, CHS2, CQS2, &
+!$acc RHO1D, MAVAIL, USTM, &
+!$acc UST_lnd, UST_wat, UST_ice, &
+!$acc PSIQ_lnd, PSIT_lnd, PSIX_lnd, &
+!$acc PSIQ_wat, PSIT_wat, PSIX_wat, &
+!$acc PSIQ_ice, PSIT_ice, PSIX_ice, &
+!$acc PSIQ2_lnd, PSIT2_lnd, &
+!$acc PSIQ2_wat, PSIT2_wat, &
+!$acc PSIQ2_ice, PSIT2_ice, &
+!$acc QSFC, QSFC_lnd, QSFC_wat, QSFC_ice, &
+!$acc QFLX, QFLX_lnd, QFLX_wat, QFLX_ice, &
+!$acc HFLX, HFLX_lnd, HFLX_wat, HFLX_ice, &
+!$acc QSFCMR_lnd, QSFCMR_wat, QSFCMR_ice, &
+!$acc QV1D, WSPD, WSPDI, CPM, TH1D, &
+!$acc THSK_lnd, THSK_wat, THSK_ice, &
+!$acc ch_lnd, ch_wat, ch_ice, &
+!$acc cm_lnd, cm_wat, cm_ice)
DO I=its,ite
if( flag_iter(i) ) then
@@ -2040,6 +2256,18 @@ SUBROUTINE SFCLAY1D_mynn(flag_iter, &
ENDDO ! end i-loop
IF (compute_diag) then
+ !$acc parallel loop present(flag_iter, dry, wet, icy, &
+ !$acc ZA, ZA2, T2, TH2, TH1D, Q2, QV1D, PSFCPA, &
+ !$acc THSK_lnd, THSK_wat, THSK_ice, &
+ !$acc QSFC_lnd, QSFC_wat, QSFC_ice, &
+ !$acc U10, V10, U1D, V1D, U1D2, V1D2, &
+ !$acc ZNTstoch_lnd, ZNTstoch_lnd, ZNTstoch_ice, &
+ !$acc PSIX_lnd, PSIX_wat, PSIX_ice, &
+ !$acc PSIX10_lnd, PSIX10_wat, PSIX10_ice, &
+ !$acc PSIT2_lnd, PSIT2_wat, PSIT2_ice, &
+ !$acc PSIT_lnd, PSIT_wat, PSIT_ice, &
+ !$acc PSIQ2_lnd, PSIQ2_wat, PSIQ2_ice, &
+ !$acc PSIQ_lnd, PSIQ_wat, PSIQ_ice)
DO I=its,ite
if( flag_iter(i) ) then
!-----------------------------------------------------
@@ -2153,6 +2381,16 @@ SUBROUTINE SFCLAY1D_mynn(flag_iter, &
!-----------------------------------------------------
! DEBUG - SUSPICIOUS VALUES
!-----------------------------------------------------
+!$acc serial present(dry, wet, icy, CPM, MAVAIL, &
+!$acc HFX, LH, wstar, RHO1D, PBLH, ZOL, ZA, MOL, &
+!$acc PSIM, PSIH, WSTAR, T1D, TH1D, THV1D, QVSH, &
+!$acc UST_wat, UST_lnd, UST_ice, &
+!$acc THSK_wat, THSK_lnd, THSK_ice, &
+!$acc THVSK_wat, THVSK_lnd, THVSK_ice, &
+!$acc ZNTstoch_wat, ZNTstoch_lnd, ZNTstoch_ice, &
+!$acc ZT_wat, ZT_lnd, ZT_ice, &
+!$acc QSFC_wat, QSFC_lnd, QSFC_ice, &
+!$acc PSIX_wat, PSIX_lnd, PSIX_ice)
IF ( debug_code == 2) THEN
DO I=its,ite
yesno = 0
@@ -2257,6 +2495,62 @@ SUBROUTINE SFCLAY1D_mynn(flag_iter, &
ENDIF
ENDDO ! end i-loop
ENDIF ! end debug option
+!$acc end serial
+
+!$acc exit data copyout(CPM, FLHC, FLQC, CHS, CH, CHS2, CQS2,&
+!$acc USTM, wstar, qstar, ZOL, MOL, RMOL, &
+!$acc HFX, QFX, LH, QSFC, QFLX, HFLX, &
+!$acc T2, TH2, Q2, WSPD, U10, V10, &
+!$acc QGH, psim, psih, &
+!$acc stress_wat, stress_lnd, stress_ice, &
+!$acc rb_wat, rb_lnd, rb_ice, &
+!$acc UST_wat, UST_lnd, UST_ice, &
+!$acc ZNT_wat, ZNT_lnd, ZNT_ice, &
+!$acc QSFC_lnd, QSFC_wat, QSFC_ice, &
+!$acc QFLX_lnd, QFLX_wat, QFLX_ice, &
+!$acc HFLX_lnd, HFLX_wat, HFLX_ice, &
+!$acc PSIX_wat, PSIX_lnd, PSIX_ice, &
+!$acc PSIX10_wat, PSIX10_lnd, PSIX10_ice, &
+!$acc PSIT2_lnd, PSIT2_wat, PSIT2_ice, &
+!$acc PSIT_lnd, PSIT_wat, PSIT_ice, &
+!$acc ch_lnd, ch_wat, ch_ice, &
+!$acc cm_lnd, cm_wat, cm_ice, &
+!$acc device_errmsg, device_errflg)
+
+! Final sync of device and host error flags and messages
+IF (device_errflg /= 0) THEN
+ errflg = device_errflg
+ errmsg = device_errmsg
+ENDIF
+
+!$acc exit data delete( flag_iter, dry, wet, icy, dx, &
+!$acc MAVAIL, PBLH, PSFCPA, z0pert, ztpert, &
+!$acc QV1D, U1D, V1D, U1D2, V1D2, T1D, P1D, &
+!$acc rstoch1D, sigmaf, shdmax, vegtype, &
+!$acc dz2w1d, dz8w1d, &
+!$acc snowh_wat, snowh_lnd, snowh_ice, &
+!$acc tskin_wat, tskin_lnd, tskin_ice, &
+!$acc tsurf_wat, tsurf_lnd, tsurf_ice)
+
+!$acc exit data delete( ZA, ZA2, THV1D, TH1D, TC1D, TV1D, &
+!$acc RHO1D, QVSH, PSIH2, PSIM10, PSIH10, WSPDI, &
+!$acc GOVRTH, PSFC, THCON, &
+!$acc zratio_lnd, zratio_ice, zratio_wat, &
+!$acc TSK_lnd, TSK_ice, TSK_wat, &
+!$acc THSK_lnd, THSK_ice, THSK_wat, &
+!$acc THVSK_lnd, THVSK_ice, THVSK_wat, &
+!$acc GZ1OZ0_lnd, GZ1OZ0_ice, GZ1OZ0_wat, &
+!$acc GZ1OZt_lnd, GZ1OZt_ice, GZ1OZt_wat, &
+!$acc GZ2OZ0_lnd, GZ2OZ0_ice, GZ2OZ0_wat, &
+!$acc GZ2OZt_lnd, GZ2OZt_ice, GZ2OZt_wat, &
+!$acc GZ10OZ0_lnd, GZ10OZ0_ice, GZ10OZ0_wat, &
+!$acc GZ10OZt_lnd, GZ10OZt_ice, GZ10OZt_wat, &
+!$acc ZNTstoch_lnd, ZNTstoch_ice, ZNTstoch_wat, &
+!$acc ZT_lnd, ZT_ice, ZT_wat, &
+!$acc ZQ_lnd, ZQ_ice, ZQ_wat, &
+!$acc PSIQ_lnd, PSIQ_ice, PSIQ_wat, &
+!$acc PSIQ2_lnd, PSIQ2_ice, PSIQ2_wat, &
+!$acc QSFCMR_lnd, QSFCMR_ice, QSFCMR_wat )
END SUBROUTINE SFCLAY1D_mynn
!-------------------------------------------------------------------
@@ -2272,6 +2566,7 @@ END SUBROUTINE SFCLAY1D_mynn
SUBROUTINE zilitinkevich_1995(Z_0,Zt,Zq,restar,ustar,KARMAN,&
& landsea,IZ0TLND2,spp_sfc,rstoch)
+ !$acc routine seq
IMPLICIT NONE
REAL(kind_phys), INTENT(IN) :: Z_0,restar,ustar,KARMAN,landsea
INTEGER, OPTIONAL, INTENT(IN) :: IZ0TLND2
@@ -2341,6 +2636,7 @@ SUBROUTINE davis_etal_2008(Z_0,ustar)
!This is an update version from Davis et al. 2008, which
!corrects a small-bias in Z_0 (AHW real-time 2012).
+ !$acc routine seq
IMPLICIT NONE
REAL(kind_phys), INTENT(IN) :: ustar
REAL(kind_phys), INTENT(OUT) :: Z_0
@@ -2368,7 +2664,7 @@ END SUBROUTINE davis_etal_2008
!>This formulation for roughness length was designed account for.
!!wave steepness.
SUBROUTINE Taylor_Yelland_2001(Z_0,ustar,wsp10)
-
+ !$acc routine seq
IMPLICIT NONE
REAL(kind_phys), INTENT(IN) :: ustar,wsp10
REAL(kind_phys), INTENT(OUT) :: Z_0
@@ -2396,7 +2692,7 @@ END SUBROUTINE Taylor_Yelland_2001
!! The Charnock parameter CZC is varied from .011 to .018.
!! between 10-m wsp = 10 and 18..
SUBROUTINE charnock_1955(Z_0,ustar,wsp10,visc,zu)
-
+ !$acc routine seq
IMPLICIT NONE
REAL(kind_phys), INTENT(IN) :: ustar, visc, wsp10, zu
REAL(kind_phys), INTENT(OUT) :: Z_0
@@ -2421,7 +2717,7 @@ END SUBROUTINE charnock_1955
!!The Charnock parameter CZC is varied from about .005 to .028
!!between 10-m wind speeds of 6 and 19 m/s.
SUBROUTINE edson_etal_2013(Z_0,ustar,wsp10,visc,zu)
-
+ !$acc routine seq
IMPLICIT NONE
REAL(kind_phys), INTENT(IN) :: ustar, visc, wsp10, zu
REAL(kind_phys), INTENT(OUT) :: Z_0
@@ -2450,7 +2746,7 @@ END SUBROUTINE edson_etal_2013
!!data. The formula for land uses a constant ratio (Z_0/7.4) taken
!!from Garratt (1992).
SUBROUTINE garratt_1992(Zt,Zq,Z_0,Ren,landsea)
-
+ !$acc routine seq
IMPLICIT NONE
REAL(kind_phys), INTENT(IN) :: Ren, Z_0,landsea
REAL(kind_phys), INTENT(OUT) :: Zt,Zq
@@ -2486,7 +2782,7 @@ END SUBROUTINE garratt_1992
!!
!!This is for use over water only.
SUBROUTINE fairall_etal_2003(Zt,Zq,Ren,ustar,visc,rstoch,spp_sfc)
-
+ !$acc routine seq
IMPLICIT NONE
REAL(kind_phys), INTENT(IN) :: Ren,ustar,visc,rstoch
INTEGER, INTENT(IN) :: spp_sfc
@@ -2530,7 +2826,7 @@ END SUBROUTINE fairall_etal_2003
!! The actual reference is unknown. This was passed along by Jim Edson (personal communication).
!! This is for use over water only, preferably open ocean.
SUBROUTINE fairall_etal_2014(Zt,Zq,Ren,ustar,visc,rstoch,spp_sfc)
-
+ !$acc routine seq
IMPLICIT NONE
REAL(kind_phys), INTENT(IN) :: Ren,ustar,visc,rstoch
INTEGER, INTENT(IN) :: spp_sfc
@@ -2578,6 +2874,7 @@ END SUBROUTINE fairall_etal_2014
!!This should only be used over land!
SUBROUTINE Yang_2008(Z_0,Zt,Zq,ustar,tstar,qst,Ren,visc)
+ !$acc routine seq
IMPLICIT NONE
REAL(kind_phys), INTENT(IN) :: Z_0, Ren, ustar, tstar, qst, visc
REAL(kind_phys) :: ht, &! roughness height at critical Reynolds number
@@ -2613,6 +2910,7 @@ END SUBROUTINE Yang_2008
!>\ingroup mynn_sfc
SUBROUTINE GFS_z0_lnd(z0max,shdmax,z1,vegtype,ivegsrc,z0pert)
+ !$acc routine seq
REAL(kind_phys), INTENT(OUT) :: z0max
REAL(kind_phys), INTENT(IN) :: shdmax,z1,z0pert
INTEGER, INTENT(IN) :: vegtype,ivegsrc
@@ -2673,6 +2971,7 @@ END SUBROUTINE GFS_z0_lnd
!>\ingroup mynn_sfc
SUBROUTINE GFS_zt_lnd(ztmax,z0max,sigmaf,ztpert,ustar_lnd)
+ !$acc routine seq
REAL(kind_phys), INTENT(OUT) :: ztmax
REAL(kind_phys), INTENT(IN) :: z0max,sigmaf,ztpert,ustar_lnd
REAL(kind_phys) :: czilc, tem1, tem2
@@ -2701,6 +3000,7 @@ END SUBROUTINE GFS_zt_lnd
!>\ingroup mynn_sfc
SUBROUTINE GFS_z0_wat(z0rl_wat,ustar_wat,WSPD,z1,sfc_z0_type,redrag)
+ !$acc routine seq
REAL(kind_phys), INTENT(OUT) :: z0rl_wat
REAL(kind_phys), INTENT(INOUT):: ustar_wat
REAL(kind_phys), INTENT(IN) :: wspd,z1
@@ -2752,19 +3052,27 @@ SUBROUTINE GFS_z0_wat(z0rl_wat,ustar_wat,WSPD,z1,sfc_z0_type,redrag)
END SUBROUTINE GFS_z0_wat
!--------------------------------------------------------------------
!>\ingroup mynn_sfc
- SUBROUTINE GFS_zt_wat(ztmax,z0rl_wat,restar,WSPD,z1,sfc_z0_type,errmsg,errflg)
-
+ SUBROUTINE GFS_zt_wat(ztmax,z0rl_wat,restar,WSPD,z1,sfc_z0_type,device_errmsg,device_errflg)
+ !$acc routine seq
real(kind_phys), INTENT(OUT) :: ztmax
real(kind_phys), INTENT(IN) :: wspd,z1,z0rl_wat,restar
INTEGER, INTENT(IN) :: sfc_z0_type
- character(len=*), intent(out) :: errmsg
- integer, intent(out) :: errflg
+
+! Using device_errmsg and device_errflg rather than the CCPP errmsg and errflg
+! so that this subroutine can be run on an accelerator device with OpenACC.
+! character(len=*), intent(out) :: errmsg
+! integer, intent(out) :: errflg
+ character(len=512), intent(out) :: device_errmsg
+ integer, intent(out) :: device_errflg
+
real(kind_phys) :: z0,z0max,wind10m,rat,ustar_wat
real(kind_phys), PARAMETER :: charnock = 0.014, z0s_max=.317e-2
! Initialize error-handling
- errflg = 0
- errmsg = ''
+! errflg = 0
+! errmsg = ''
+ device_errflg = 0
+ device_errmsg = ''
! z0 = 0.01 * z0rl_wat
!Already converted to meters in the wrapper
@@ -2795,9 +3103,12 @@ SUBROUTINE GFS_zt_wat(ztmax,z0rl_wat,restar,WSPD,z1,sfc_z0_type,errmsg,errflg)
call znot_t_v7(wind10m, ztmax) ! 10-m wind,m/s, ztmax(m)
else if (sfc_z0_type > 0) then
write(0,*)'no option for sfc_z0_type=',sfc_z0_type
- errflg = 1
- errmsg = 'ERROR(GFS_zt_wat): sfc_z0_type not valid.'
+! errflg = 1
+! errmsg = 'ERROR(GFS_zt_wat): sfc_z0_type not valid.'
+ device_errflg = 1
+ device_errmsg = 'ERROR(GFS_zt_wat): sfc_z0_type not valid.'
return
+
endif
END SUBROUTINE GFS_zt_wat
@@ -2807,6 +3118,7 @@ END SUBROUTINE GFS_zt_wat
!! Weiguo Wang, 2019-0425
SUBROUTINE znot_m_v6(uref, znotm)
+ !$acc routine seq
use machine , only : kind_phys
IMPLICIT NONE
! Calculate areodynamical roughness over water with input 10-m wind
@@ -2856,6 +3168,7 @@ END SUBROUTINE znot_m_v6
!!
SUBROUTINE znot_t_v6(uref, znott)
+ !$acc routine seq
IMPLICIT NONE
!> Calculate scalar roughness over water with input 10-m wind
!! For low-to-moderate winds, try to match the Ck-U10 relationship from COARE algorithm
@@ -2922,6 +3235,7 @@ END SUBROUTINE znot_t_v6
!!
SUBROUTINE znot_m_v7(uref, znotm)
+ !$acc routine seq
IMPLICIT NONE
!> Calculate areodynamical roughness over water with input 10-m wind
!! For low-to-moderate winds, try to match the Cd-U10 relationship from COARE V3.5 (Edson et al. 2013)
@@ -2971,6 +3285,7 @@ END SUBROUTINE znot_m_v7
!!
SUBROUTINE znot_t_v7(uref, znott)
+ !$acc routine seq
IMPLICIT NONE
!> Calculate scalar roughness over water with input 10-m wind
!! For low-to-moderate winds, try to match the Ck-U10 relationship from COARE algorithm
@@ -3040,6 +3355,7 @@ END SUBROUTINE znot_t_v7
!! This should only be used over snow/ice!
SUBROUTINE Andreas_2002(Z_0,bvisc,ustar,Zt,Zq)
+ !$acc routine seq
IMPLICIT NONE
REAL(kind_phys), INTENT(IN) :: Z_0, bvisc, ustar
REAL(kind_phys), INTENT(OUT) :: Zt, Zq
@@ -3313,6 +3629,7 @@ END SUBROUTINE PSI_CB2005
!! and Holtslag (1991) for stable conditions.
SUBROUTINE Li_etal_2010(zL, Rib, zaz0, z0zt)
+ !$acc routine seq
IMPLICIT NONE
REAL(kind_phys), INTENT(OUT) :: zL
REAL(kind_phys), INTENT(IN) :: Rib, zaz0, z0zt
@@ -3471,6 +3788,7 @@ REAL(kind_phys) function zolri2(zol2,ri2,za,z0,zt,psi_opt)
REAL(kind_phys) function zolrib(ri,za,z0,zt,logz0,logzt,zol1,psi_opt)
+ !$acc routine seq
! This iterative algorithm to compute z/L from bulk-Ri
IMPLICIT NONE
@@ -3480,7 +3798,7 @@ REAL(kind_phys) function zolrib(ri,za,z0,zt,logz0,logzt,zol1,psi_opt)
REAL(kind_phys) :: zol20,zol3,zolt,zolold
INTEGER :: n
INTEGER, PARAMETER :: nmax = 20
- REAL(kind_phys), DIMENSION(nmax):: zLhux
+ !REAL(kind_phys), DIMENSION(nmax):: zLhux
REAL(kind_phys) :: psit2,psix2
!print*,"+++++++INCOMING: z/L=",zol1," ri=",ri
@@ -3522,7 +3840,7 @@ REAL(kind_phys) function zolrib(ri,za,z0,zt,logz0,logzt,zol1,psi_opt)
endif
!print*,"n=",n," psit2=",psit2," psix2=",psix2
zolrib=ri*psix2**2/psit2
- zLhux(n)=zolrib
+ !zLhux(n)=zolrib
n=n+1
enddo
@@ -3530,7 +3848,7 @@ REAL(kind_phys) function zolrib(ri,za,z0,zt,logz0,logzt,zol1,psi_opt)
!print*,"iter FAIL, n=",n," Ri=",ri," z/L=",zolri
!if convergence fails, use approximate values:
CALL Li_etal_2010(zolrib, ri, za/z0, z0/zt)
- zLhux(n)=zolrib
+ !zLhux(n)=zolrib
!print*,"FAILED, n=",n," Ri=",ri," z0=",z0
!print*,"z/L=",zLhux(1:nmax)
else
@@ -3595,6 +3913,7 @@ END SUBROUTINE psi_init
!
!>\ingroup mynn_sfc
real(kind_phys) function psim_stable_full(zolf)
+ !$acc routine seq
real(kind_phys) :: zolf
!psim_stable_full=-6.1*log(zolf+(1+zolf**2.5)**(1./2.5))
@@ -3605,6 +3924,7 @@ real(kind_phys) function psim_stable_full(zolf)
!>\ingroup mynn_sfc
real(kind_phys) function psih_stable_full(zolf)
+ !$acc routine seq
real(kind_phys) :: zolf
!psih_stable_full=-5.3*log(zolf+(1+zolf**1.1)**(1./1.1))
@@ -3615,6 +3935,7 @@ real(kind_phys) function psih_stable_full(zolf)
!>\ingroup mynn_sfc
real(kind_phys) function psim_unstable_full(zolf)
+ !$acc routine seq
real(kind_phys) :: zolf,x,ym,psimc,psimk
x=(1.-16.*zolf)**.25
@@ -3633,6 +3954,7 @@ real(kind_phys) function psim_unstable_full(zolf)
!>\ingroup mynn_sfc
real(kind_phys) function psih_unstable_full(zolf)
+ !$acc routine seq
real(kind_phys) :: zolf,y,yh,psihc,psihk
y=(1.-16.*zolf)**.5
@@ -3654,6 +3976,7 @@ real(kind_phys) function psih_unstable_full(zolf)
!>\ingroup mynn_sfc
!!
REAL(kind_phys) function psim_stable_full_gfs(zolf)
+ !$acc routine seq
REAL(kind_phys) :: zolf
REAL(kind_phys), PARAMETER :: alpha4 = 20.
REAL(kind_phys) :: aa
@@ -3667,6 +3990,7 @@ REAL(kind_phys) function psim_stable_full_gfs(zolf)
!>\ingroup mynn_sfc
!!
real(kind_phys) function psih_stable_full_gfs(zolf)
+ !$acc routine seq
real(kind_phys) :: zolf
real(kind_phys), PARAMETER :: alpha4 = 20.
real(kind_phys) :: bb
@@ -3680,6 +4004,7 @@ real(kind_phys) function psih_stable_full_gfs(zolf)
!>\ingroup mynn_sfc
!!
real(kind_phys) function psim_unstable_full_gfs(zolf)
+ !$acc routine seq
real(kind_phys) :: zolf
real(kind_phys) :: hl1,tem1
real(kind_phys), PARAMETER :: a0=-3.975, a1=12.32, &
@@ -3700,6 +4025,7 @@ real(kind_phys) function psim_unstable_full_gfs(zolf)
!>\ingroup mynn_sfc
!!
real(kind_phys) function psih_unstable_full_gfs(zolf)
+ !$acc routine seq
real(kind_phys) :: zolf
real(kind_phys) :: hl1,tem1
real(kind_phys), PARAMETER :: a0p=-7.941, a1p=24.75, &
@@ -3720,6 +4046,7 @@ real(kind_phys) function psih_unstable_full_gfs(zolf)
!>\ingroup mynn_sfc
!! look-up table functions - or, if beyond -10 < z/L < 10, recalculate
real(kind_phys) function psim_stable(zolf,psi_opt)
+ !$acc routine seq
integer :: nzol,psi_opt
real(kind_phys) :: rzol,zolf
@@ -3740,6 +4067,7 @@ real(kind_phys) function psim_stable(zolf,psi_opt)
!>\ingroup mynn_sfc
real(kind_phys) function psih_stable(zolf,psi_opt)
+ !$acc routine seq
integer :: nzol,psi_opt
real(kind_phys) :: rzol,zolf
@@ -3760,6 +4088,7 @@ real(kind_phys) function psih_stable(zolf,psi_opt)
!>\ingroup mynn_sfc
real(kind_phys) function psim_unstable(zolf,psi_opt)
+ !$acc routine seq
integer :: nzol,psi_opt
real(kind_phys) :: rzol,zolf
@@ -3780,6 +4109,7 @@ real(kind_phys) function psim_unstable(zolf,psi_opt)
!>\ingroup mynn_sfc
real(kind_phys) function psih_unstable(zolf,psi_opt)
+ !$acc routine seq
integer :: nzol,psi_opt
real(kind_phys) :: rzol,zolf
diff --git a/physics/module_sf_ruclsm.F90 b/physics/module_sf_ruclsm.F90
index 6294bc068..160127e43 100644
--- a/physics/module_sf_ruclsm.F90
+++ b/physics/module_sf_ruclsm.F90
@@ -991,10 +991,12 @@ SUBROUTINE LSMRUC(xlat,xlon, &
if(mosaic_lu == 1) then
! greenness factor: between 0 for min greenness and 1 for max greenness.
factor = max(zero,min(one,(vegfra(i,j)-shdmin(i,j))/max(one,(shdmax(i,j)-shdmin(i,j)))))
+ if (debug_print ) then
if (abs(xlat-testptlat).lt.0.1 .and. &
abs(xlon-testptlon).lt.0.1)then
print *,' lat,lon=',xlat,xlon,' factor=',factor
endif
+ endif
if((ivgtyp(i,j) == natural .or. ivgtyp(i,j) == crop) .and. factor > 0.75) then
! cropland or grassland, apply irrigation during the growing seaspon when fraction
diff --git a/physics/mp_thompson.F90 b/physics/mp_thompson.F90
index 6a95a706c..c456e87cd 100644
--- a/physics/mp_thompson.F90
+++ b/physics/mp_thompson.F90
@@ -409,7 +409,7 @@ subroutine mp_thompson_run(ncol, nlev, con_g, con_rd, &
integer, intent(in) :: n_var_spp
real(kind_phys), intent(in) :: spp_wts_mp(:,:)
real(kind_phys), intent(in) :: spp_prt_list(:)
- character(len=3), intent(in) :: spp_var_list(:)
+ character(len=10), intent(in) :: spp_var_list(:)
real(kind_phys), intent(in) :: spp_stddev_cutoff(:)
logical, intent (in) :: cplchm
diff --git a/physics/mp_thompson.meta b/physics/mp_thompson.meta
index 691698281..5918e4dd9 100644
--- a/physics/mp_thompson.meta
+++ b/physics/mp_thompson.meta
@@ -725,7 +725,7 @@
units = none
dimensions = (number_of_perturbed_spp_schemes)
type = character
- kind = len=3
+ kind = len=10
intent = in
[cplchm]
standard_name = flag_for_chemistry_coupling
diff --git a/physics/mynnsfc_wrapper.F90 b/physics/mynnsfc_wrapper.F90
index 1a970c9f4..3c033e65e 100644
--- a/physics/mynnsfc_wrapper.F90
+++ b/physics/mynnsfc_wrapper.F90
@@ -191,6 +191,16 @@ SUBROUTINE mynnsfc_wrapper_run( &
& IMS,IME,JMS,JME,KMS,KME, &
& ITS,ITE,JTS,JTE,KTS,KTE
+!$acc enter data create(hfx, znt, psim, psih, chs, &
+!$acc mavail, xland, GZ1OZ0, cpm, qgh, &
+!$acc qfx, snowh_wat)
+
+!$acc enter data create(dz, th, qv)
+
+!$acc enter data copyin(rmol, phii, t3d, exner, qvsh, slmsk, xland)
+
+!$acc enter data copyin(dry, wet, icy, znt_lnd, znt_wat, znt_ice, qsfc_lnd, qsfc_ice, qsfc_lnd_ruc, qsfc_ice_ruc)
+
! Initialize CCPP error handling variables
errmsg = ''
errflg = 0
@@ -203,6 +213,7 @@ SUBROUTINE mynnsfc_wrapper_run( &
! write(0,*)"iter=",iter
! endif
+!$acc kernels
! prep MYNN-only variables
dz(:,:) = 0
th(:,:) = 0
@@ -210,6 +221,9 @@ SUBROUTINE mynnsfc_wrapper_run( &
hfx(:) = 0
qfx(:) = 0
rmol(:) = 0
+!$acc end kernels
+
+!$acc parallel loop collapse(2) present(dz, phii, th, t3d, exner, qv, qvsh)
do k=1,2 !levs
do i=1,im
dz(i,k)=(phii(i,k+1) - phii(i,k))*g_inv
@@ -219,6 +233,7 @@ SUBROUTINE mynnsfc_wrapper_run( &
enddo
enddo
+!$acc parallel loop present(slmsk, xland, qgh, mavail, cpm, snowh_wat)
do i=1,im
if (slmsk(i)==1. .or. slmsk(i)==2.)then !sea/land/ice mask (=0/1/2) in FV3
xland(i)=1.0 !but land/water = (1/2) in SFCLAY_mynn
@@ -235,6 +250,7 @@ SUBROUTINE mynnsfc_wrapper_run( &
snowh_wat(i) = 0.0
enddo
+!$acc kernels
! cm -> m
where (dry) znt_lnd=znt_lnd*0.01
where (wet) znt_wat=znt_wat*0.01
@@ -245,6 +261,7 @@ SUBROUTINE mynnsfc_wrapper_run( &
where (dry) qsfc_lnd = qsfc_lnd_ruc/(1.+qsfc_lnd_ruc) ! spec. hum
where (icy) qsfc_ice = qsfc_ice_ruc/(1.+qsfc_ice_ruc) ! spec. hum.
end if
+!$acc end kernels
! if (lprnt) then
! write(0,*)"CALLING SFCLAY_mynn; input:"
@@ -274,6 +291,8 @@ SUBROUTINE mynnsfc_wrapper_run( &
! write(0,*)"PBLH=",pblh(1)," xland=",xland(1)
! endif
+!$acc exit data delete(qsfc_lnd_ruc, qsfc_ice_ruc)
+!$acc exit data delete(phii, qvsh, slmsk)
CALL SFCLAY_mynn( &
u3d=u,v3d=v,t3d=t3d,qv3d=qv,p3d=prsl,dz8w=dz, &
@@ -318,6 +337,13 @@ SUBROUTINE mynnsfc_wrapper_run( &
errmsg=errmsg, errflg=errflg )
if (errflg/=0) return
+!$acc exit data delete(hfx, znt, psim, psih, chs, &
+!$acc mavail, xland, GZ1OZ0, cpm, qgh, &
+!$acc qfx, snowh_wat, t3d, exner)
+!$acc exit data delete(dz, th, qv)
+!$acc exit data copyout(rmol)
+!$acc exit data copyout(qsfc_lnd, qsfc_ice)
+
!! POST MYNN SURFACE LAYER (INTERSTITIAL) WORK:
!do i = 1, im
! !* Taken from sfc_nst.f
@@ -336,10 +362,15 @@ SUBROUTINE mynnsfc_wrapper_run( &
! znt_ice(i)=znt_ice(i)*100.
!enddo
+!$acc kernels
! m -> cm
where (dry) znt_lnd=znt_lnd*100.
where (wet) znt_wat=znt_wat*100.
where (icy) znt_ice=znt_ice*100.
+!$acc end kernels
+
+!$acc exit data delete(dry, wet, icy)
+!$acc exit data copyout(znt_lnd, znt_wat, znt_ice)
! if (lprnt) then
! write(0,*)
diff --git a/physics/noahmp_tables.f90 b/physics/noahmp_tables.f90
index 3b06d7f53..753c8ff24 100644
--- a/physics/noahmp_tables.f90
+++ b/physics/noahmp_tables.f90
@@ -484,6 +484,9 @@ subroutine read_mp_table_parameters(errmsg, errflg)
sr2006_psi_e_a, sr2006_psi_e_b, sr2006_psi_e_c, sr2006_smcmax_a, &
sr2006_smcmax_b
+ errmsg = ''
+ errflg = 0
+
! initialize our variables to bad values, so that if the namelist read fails, we come to a screeching halt as soon as we try to use anything.
! vegetation parameters
isurban_table = -99999
@@ -783,7 +786,7 @@ subroutine read_mp_table_parameters(errmsg, errflg)
open(15, status='old', form='formatted', action='read', iostat=ierr)
end if
if ( ierr /= 0 ) then
- errmsg = 'warning: cannot find file noahmptable.tb'
+ errmsg = 'warning: cannot find file noahmptable.tbl'
errflg = 1
return
! write(*,'("warning: cannot find file noahmptable.tbl")')
@@ -914,7 +917,7 @@ subroutine read_mp_table_parameters(errmsg, errflg)
open(15, status='old', form='formatted', action='read', iostat=ierr)
end if
if ( ierr /= 0 ) then
- errmsg = 'warning: cannot find file noahmptable.tb'
+ errmsg = 'warning: cannot find file noahmptable.tbl'
errflg = 1
return
! write(*,'("warning: cannot find file noahmptable.tbl")')
@@ -957,7 +960,7 @@ subroutine read_mp_table_parameters(errmsg, errflg)
open(15, status='old', form='formatted', action='read', iostat=ierr)
end if
if ( ierr /= 0 ) then
- errmsg = 'warning: cannot find file noahmptable.tb'
+ errmsg = 'warning: cannot find file noahmptable.tbl'
errflg = 1
return
! write(*,'("warning: cannot find file noahmptable.tbl")')
@@ -982,7 +985,7 @@ subroutine read_mp_table_parameters(errmsg, errflg)
open(15, status='old', form='formatted', action='read', iostat=ierr)
end if
if (ierr /= 0) then
- errmsg = 'warning: cannot find file noahmptable.tb'
+ errmsg = 'warning: cannot find file noahmptable.tbl'
errflg = 1
return
! write(*,'("warning: cannot find file noahmptable.tbl")')
@@ -1011,7 +1014,7 @@ subroutine read_mp_table_parameters(errmsg, errflg)
open(15, status='old', form='formatted', action='read', iostat=ierr)
end if
if (ierr /= 0) then
- errmsg = 'warning: cannot find file noahmptable.tb'
+ errmsg = 'warning: cannot find file noahmptable.tbl'
errflg = 1
return
! write(*,'("warning: cannot find file noahmptable.tbl")')
@@ -1069,7 +1072,7 @@ subroutine read_mp_table_parameters(errmsg, errflg)
open(15, status='old', form='formatted', action='read', iostat=ierr)
end if
if (ierr /= 0) then
- errmsg = 'warning: cannot find file noahmptable.tb'
+ errmsg = 'warning: cannot find file noahmptable.tbl'
errflg = 1
return
! write(*,'("warning: cannot find file noahmptable.tbl")')
@@ -1096,7 +1099,7 @@ subroutine read_mp_table_parameters(errmsg, errflg)
open(15, status='old', form='formatted', action='read', iostat=ierr)
end if
if (ierr /= 0) then
- errmsg = 'warning: cannot find file noahmptable.tb'
+ errmsg = 'warning: cannot find file noahmptable.tbl'
errflg = 1
return
! write(*,'("warning: cannot find file noahmptable.tbl")')
@@ -1249,7 +1252,7 @@ subroutine read_mp_table_parameters(errmsg, errflg)
open(15, status='old', form='formatted', action='read', iostat=ierr)
end if
if (ierr /= 0) then
- errmsg = 'warning: cannot find file noahmptable.tb'
+ errmsg = 'warning: cannot find file noahmptable.tbl'
errflg = 1
return
! write(*,'("warning: cannot find file noahmptable.tbl")')
@@ -1278,7 +1281,7 @@ subroutine read_mp_table_parameters(errmsg, errflg)
open(15, status='old', form='formatted', action='read', iostat=ierr)
end if
if (ierr /= 0) then
- errmsg = 'warning: cannot find file noahmptable.tb'
+ errmsg = 'warning: cannot find file noahmptable.tbl'
errflg = 1
return
! write(*,'("warning: cannot find file noahmptable.tbl")')
diff --git a/physics/satmedmfvdifq.F b/physics/satmedmfvdifq.F
index 75925a5f3..73fc4aff8 100644
--- a/physics/satmedmfvdifq.F
+++ b/physics/satmedmfvdifq.F
@@ -22,7 +22,7 @@ module satmedmfvdifq
!! and to reduce the negative wind speed bias in upper troposphere
!!
!! Incorporate the LES-based changes for TC simulation
-!! (Chen et al.,2022, https://doi.org/10.1175/WAF-D-21-0168.1)
+!! (Chen et al.,2022 \cite Chen_2022)
!! with additional improvements on MF working with Cu schemes
!! Xiaomin Chen, 5/2/2022
!!
@@ -443,7 +443,9 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, &
endif
enddo
!
-! compute a function for green vegetation fraction and surface roughness
+!> - Compute a function for green vegetation fraction and surface roughness.
+!! Entrainment rate in updraft is a function of vegetation fraction and surface
+!! roughness length
!
do i = 1,im
tem = (sigmaf(i) - vegflo) / (vegfup - vegflo)
@@ -745,7 +747,7 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, &
endif
enddo
!
-! compute mean tke within pbl
+!> - Compute mean tke within pbl
!
do i = 1, im
sumx(i) = 0.
@@ -766,8 +768,8 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, &
endif
enddo
!
-! compute wind shear term as a sink term for updraft and downdraft
-! velocity
+!> - Compute wind shear term as a sink term for updraft and downdraft
+!! velocity
!
kps = max(kmpbl, kmscu)
do k = 2, kps
diff --git a/physics/set_soilveg.f b/physics/set_soilveg.f
index 37f2c2a73..35f4ace37 100644
--- a/physics/set_soilveg.f
+++ b/physics/set_soilveg.f
@@ -44,6 +44,9 @@ subroutine set_soilveg(me,isot,ivet,nlunit,errmsg,errflg)
& DEFINED_SLOPE, FXEXP_DATA, NROOT_DATA, REFKDT_DATA, Z0_DATA,
& CZIL_DATA, LAI_DATA, CSOIL_DATA
+ errmsg = ''
+ errflg = 0
+
cmy end locals
if(ivet.eq.2) then
diff --git a/physics/smoke_dust/dust_data_mod.F90 b/physics/smoke_dust/dust_data_mod.F90
index a710701f1..eb809378d 100755
--- a/physics/smoke_dust/dust_data_mod.F90
+++ b/physics/smoke_dust/dust_data_mod.F90
@@ -44,24 +44,10 @@ module dust_data_mod
! Never used:
! real(kind_phys), parameter :: fengsha_alpha = 0.3
! real(kind_phys), parameter :: fengsha_gamma = 1.3
+
! -- FENGSHA threshold velocities based on Dale A. Gillette's data
integer, parameter :: fengsha_maxstypes = 13
-! real(kind_phys), dimension(fengsha_maxstypes) :: dust_uthres = &
-! (/ 0.065, & ! Sand - 1
-! 0.20, & ! Loamy Sand - 2
-! 0.52, & ! Sandy Loam - 3
-! 0.50, & ! Silt Loam - 4
-! 0.50, & ! Silt - 5
-! 0.60, & ! Loam - 6
-! 0.73, & ! Sandy Clay Loam - 7
-! 0.73, & ! Silty Clay Loam - 8
-! 0.80, & ! Clay Loam - 9
-! 0.95, & ! Sandy Clay - 10
-! 0.95, & ! Silty Clay - 11
-! 1.00, & ! Clay - 12
-! 9.999 /) ! Other - 13
-! dust_uthres = 0.065, 0.18, 0.27, 0.30, 0.35, 0.38, 0.35, 0.41, 0.41,
-! 0.45,0.50,0.45,9999.0
+
real(kind_phys), dimension(fengsha_maxstypes), parameter :: dust_uthres = &
(/ 0.065, & ! Sand - 1
0.18, & ! Loamy Sand - 2
@@ -76,12 +62,16 @@ module dust_data_mod
0.50, & ! Silty Clay - 11
0.45, & ! Clay - 12
9999.0 /) ! Other - 13
- ! -- FENGSHA uses precalculated drag partition from ASCAT. See: Prigent et al. (2012,2015)
- integer, parameter :: dust_calcdrag = 1
- real(kind_phys) :: dust_alpha = 2.2
+ ! -- FENGSHA uses precalculated drag partition
+ integer, parameter :: dust_calcdrag = 1
+ ! -- FENGSHA dust moisture parameterization 1:fecan - 2:shao
+ integer :: dust_moist_opt = 1
+
+ real(kind_phys) :: dust_alpha = 1.0
real(kind_phys) :: dust_gamma = 1.0
-
+ real(kind_phys) :: dust_moist_correction = 1.0
+ real(kind_phys) :: dust_drylimit_factor = 1.0
! -- sea salt parameters
integer, dimension(nsalt), parameter :: spoint = (/ 1, 2, 2, 2, 2, 2, 3, 3, 3 /) ! 1 Clay, 2 Silt, 3 Sand
@@ -93,7 +83,7 @@ module dust_data_mod
(/ 1., 0.2, 0.2, 0.2, 0.2, 0.2, 0.333, 0.333, 0.333 /)
- ! -- soil vagatation parameters
+ ! -- soil vegatation parameters
integer, parameter :: max_soiltyp = 30
real(kind_phys), dimension(max_soiltyp), parameter :: &
maxsmc = (/ 0.421, 0.464, 0.468, 0.434, 0.406, 0.465, &
diff --git a/physics/smoke_dust/dust_fengsha_mod.F90 b/physics/smoke_dust/dust_fengsha_mod.F90
index 54a64239d..1e24c8947 100755
--- a/physics/smoke_dust/dust_fengsha_mod.F90
+++ b/physics/smoke_dust/dust_fengsha_mod.F90
@@ -61,6 +61,8 @@ subroutine gocart_dust_fengsha_driver(dt, &
REAL(kind_phys), INTENT(IN) :: dt ! time step
REAL(kind_phys), INTENT(IN) :: g ! gravity (m/s**2)
+
+
! Local variables
integer :: nmx,i,j,k,imx,jmx,lmx
integer :: ilwi
@@ -75,6 +77,7 @@ subroutine gocart_dust_fengsha_driver(dt, &
real(kind_phys), DIMENSION (num_emis_dust) :: distribution
real(kind_phys), dimension (3) :: massfrac
real(kind_phys) :: erodtot
+ real(kind_phys) :: moist_volumetric
! conversion values
conver=1.e-9
@@ -174,10 +177,13 @@ subroutine gocart_dust_fengsha_driver(dt, &
endif
endif
+ ! soil moisture correction factor
+ moist_volumetric = dust_moist_correction * smois(i,2,j)
+
! Call dust emission routine.
call source_dust(imx,jmx, lmx, nmx, dt, tc, ustar, massfrac, &
- erodtot, dxy, smois(i,1,j), airden, airmas, bems, g, dust_alpha, dust_gamma, &
+ erodtot, dxy, moist_volumetric, airden, airmas, bems, g, dust_alpha, dust_gamma, &
R, uthr(i,j))
! convert back to concentration
@@ -457,10 +463,16 @@ subroutine DustEmissionFENGSHA(slc, clay, sand, silt, &
! Now compute size-dependent total emission flux
! ----------------------------------------------
- ! Fecan moisture correction
- ! -------------------------
- h = moistureCorrectionFecan(slc, sand, clay, rhop)
-
+
+ if (dust_moist_opt .eq. 1) then
+
+ ! Fecan moisture correction
+ ! -------------------------
+ h = moistureCorrectionFecan(slc, sand, clay)
+ else
+ ! shao soil moisture correction
+ h = moistureCorrectionShao(slc)
+ end if
! Adjust threshold
! ----------------
u_thresh = uthrs * h
@@ -478,7 +490,7 @@ subroutine DustEmissionFENGSHA(slc, clay, sand, silt, &
end subroutine DustEmissionFENGSHA
!-----------------------------------------------------------------
- real function soilMoistureConvertVol2Grav(vsoil, sandfrac, rhop)
+ real function soilMoistureConvertVol2Grav(vsoil, sandfrac)
! !USES:
implicit NONE
@@ -486,7 +498,6 @@ real function soilMoistureConvertVol2Grav(vsoil, sandfrac, rhop)
! !INPUT PARAMETERS:
REAL(kind_phys), intent(in) :: vsoil ! volumetric soil moisture fraction [1]
REAL(kind_phys), intent(in) :: sandfrac ! fractional sand content [1]
- REAL(kind_phys), intent(in) :: rhop ! dry dust density [kg m-3]
! !DESCRIPTION: Convert soil moisture fraction from volumetric to gravimetric.
!
@@ -500,20 +511,21 @@ real function soilMoistureConvertVol2Grav(vsoil, sandfrac, rhop)
! !CONSTANTS:
REAL(kind_phys), parameter :: rhow = 1000. ! density of water [kg m-3]
-
+ REAL(kind_phys), parameter :: rhop = 1700. ! density of dry soil
!EOP
!-------------------------------------------------------------------------
! Begin...
! Saturated volumetric water content (sand-dependent) ! [m3 m-3]
- vsat = 0.489 - 0.00126 * ( 100. * sandfrac )
+ vsat = 0.489 - 0.126 * sandfrac
+
! Gravimetric soil content
- soilMoistureConvertVol2Grav = vsoil * rhow / (rhop * (1. - vsat))
+ soilMoistureConvertVol2Grav = 100.0 * (vsoil * rhow / rhop / ( 1. - vsat))
end function soilMoistureConvertVol2Grav
!----------------------------------------------------------------
- real function moistureCorrectionFecan(slc, sand, clay, rhop)
+ real function moistureCorrectionFecan(slc, sand, clay)
! !USES:
implicit NONE
@@ -522,7 +534,6 @@ real function moistureCorrectionFecan(slc, sand, clay, rhop)
REAL(kind_phys), intent(in) :: slc ! liquid water content of top soil layer, volumetric fraction [1]
REAL(kind_phys), intent(in) :: sand ! fractional sand content [1]
REAL(kind_phys), intent(in) :: clay ! fractional clay content [1]
- REAL(kind_phys), intent(in) :: rhop ! dry dust density [kg m-3]
! !DESCRIPTION: Compute correction factor to account for Fecal soil moisture
!
@@ -540,15 +551,46 @@ real function moistureCorrectionFecan(slc, sand, clay, rhop)
! Begin...
! Convert soil moisture from volumetric to gravimetric
- grvsoilm = soilMoistureConvertVol2Grav(slc, sand, 2650.)
+ grvsoilm = soilMoistureConvertVol2Grav(slc, sand)
! Compute fecan dry limit
- drylimit = clay * (14.0 * clay + 17.0)
+ drylimit = dust_drylimit_factor * clay * (14.0 * clay + 17.0)
! Compute soil moisture correction
moistureCorrectionFecan = sqrt(1.0 + 1.21 * max(0., grvsoilm - drylimit)**0.68)
end function moistureCorrectionFecan
+!----------------------------------------------------------------
+ real function moistureCorrectionShao(slc)
+
+! !USES:
+ implicit NONE
+
+! !INPUT PARAMETERS:
+ REAL(kind_phys), intent(in) :: slc ! liquid water content of top soil layer, volumetric fraction [1]
+
+! !DESCRIPTION: Compute correction factor to account for Fecal soil moisture
+!
+! !REVISION HISTORY:
+!
+! 02Apr2020, B.Baker/NOAA - Original implementation
+! 01Apr2020, R.Montuoro/NOAA - Adapted for GOCART process library
+
+! !Local Variables
+ real :: grvsoilm
+ real :: drylimit
+
+!EOP
+!---------------------------------------------------------------
+! Begin...
+
+ if (slc < 0.03) then
+ moistureCorrectionShao = exp(22.7 * slc)
+ else
+ moistureCorrectionShao = exp(95.3 * slc - 2.029)
+ end if
+
+ end function moistureCorrectionShao
!---------------------------------------------------------------
real function DustFluxV2HRatioMB95(clay, kvhmax)
diff --git a/physics/smoke_dust/rrfs_smoke_wrapper.F90 b/physics/smoke_dust/rrfs_smoke_wrapper.F90
index 1f9ef6340..7b69fc9e3 100755
--- a/physics/smoke_dust/rrfs_smoke_wrapper.F90
+++ b/physics/smoke_dust/rrfs_smoke_wrapper.F90
@@ -12,7 +12,8 @@ module rrfs_smoke_wrapper
num_moist, num_chem, num_emis_seas, num_emis_dust, &
DUST_OPT_FENGSHA, p_qv, p_atm_shum, p_atm_cldq, &
p_smoke, p_dust_1, p_coarse_pm, epsilc
- use dust_data_mod, only : dust_alpha, dust_gamma
+ use dust_data_mod, only : dust_alpha, dust_gamma, dust_moist_opt, &
+ dust_moist_correction, dust_drylimit_factor
use plume_data_mod, only : p_frp_std, p_frp_hr, num_frp_plume
use seas_mod, only : gocart_seasalt_driver
use dust_fengsha_mod, only : gocart_dust_fengsha_driver
@@ -49,6 +50,7 @@ subroutine rrfs_smoke_wrapper_run(im, kte, kme, ktau, dt, garea, land, jdate,
ebb_smoke_hr, frp_hr, frp_std_hr, &
coef_bb, ebu_smoke,fhist, min_fplume, max_fplume, hwp, wetness, &
smoke_ext, dust_ext, ndvel, ddvel_inout,rrfs_sd, &
+ dust_moist_opt_in, dust_moist_correction_in, dust_drylimit_factor_in, &
dust_alpha_in, dust_gamma_in, fire_in, &
seas_opt_in, dust_opt_in, drydep_opt_in, coarsepm_settling_in, &
do_plumerise_in, plumerisefire_frq_in, addsmoke_flag_in, &
@@ -91,12 +93,15 @@ subroutine rrfs_smoke_wrapper_run(im, kte, kme, ktau, dt, garea, land, jdate,
real(kind_phys), dimension(:,:), intent(out) :: smoke_ext, dust_ext
real(kind_phys), dimension(:,:), intent(inout) :: nwfa, nifa
real(kind_phys), dimension(:,:), intent(inout) :: ddvel_inout
- real (kind=kind_phys), dimension(:), intent(in) :: wetness
- integer, intent(in ) :: imp_physics, imp_physics_thompson
- real (kind=kind_phys), intent(in) :: dust_alpha_in, dust_gamma_in, wetdep_ls_alpha_in
- integer, intent(in) :: seas_opt_in, dust_opt_in, drydep_opt_in, &
- coarsepm_settling_in, plumerisefire_frq_in, &
- addsmoke_flag_in, wetdep_ls_opt_in
+ real(kind_phys), dimension(:), intent(in) :: wetness
+ real(kind_phys), intent(in) :: dust_alpha_in, dust_gamma_in, wetdep_ls_alpha_in
+ real(kind_phys), intent(in) :: dust_moist_correction_in
+ real(kind_phys), intent(in) :: dust_drylimit_factor_in
+ integer, intent(in) :: dust_moist_opt_in
+ integer, intent(in) :: imp_physics, imp_physics_thompson
+ integer, intent(in) :: seas_opt_in, dust_opt_in, drydep_opt_in, &
+ coarsepm_settling_in, plumerisefire_frq_in, &
+ addsmoke_flag_in, wetdep_ls_opt_in
logical, intent(in ) :: do_plumerise_in, rrfs_sd
character(len=*), intent(out) :: errmsg
integer, intent(out) :: errflg
@@ -314,6 +319,9 @@ subroutine rrfs_smoke_wrapper_run(im, kte, kme, ktau, dt, garea, land, jdate,
! Set at compile time in dust_data_mod:
dust_alpha = dust_alpha_in
dust_gamma = dust_gamma_in
+ dust_moist_opt = dust_moist_opt_in
+ dust_moist_correction = dust_moist_correction_in
+ dust_drylimit_factor = dust_drylimit_factor_in
call gocart_dust_fengsha_driver(dt,chem,rho_phy,smois,p8w,ssm, &
isltyp,vegfrac,snowh,xland,dxy,g,emis_dust,ust,znt, &
clayf,sandf,rdrag,uthr, &
diff --git a/physics/smoke_dust/rrfs_smoke_wrapper.meta b/physics/smoke_dust/rrfs_smoke_wrapper.meta
index bf2fddd60..cddc20fbc 100755
--- a/physics/smoke_dust/rrfs_smoke_wrapper.meta
+++ b/physics/smoke_dust/rrfs_smoke_wrapper.meta
@@ -210,17 +210,17 @@
kind = kind_phys
intent = in
[nsoil]
- standard_name = vertical_dimension_of_soil
- long_name = soil vertical layer dimension
+ standard_name = vertical_dimension_of_soil_internal_to_land_surface_scheme
+ long_name = number of soil layers internal to land surface model
units = count
dimensions = ()
type = integer
intent = in
[smc]
- standard_name = volume_fraction_of_condensed_water_in_soil
- long_name = volumetric fraction of soil moisture
+ standard_name = volume_fraction_of_soil_moisture_for_land_surface_model
+ long_name = volumetric fraction of soil moisture for lsm
units = frac
- dimensions = (horizontal_loop_extent,vertical_dimension_of_soil)
+ dimensions = (horizontal_loop_extent,vertical_dimension_of_soil_internal_to_land_surface_scheme)
type = real
kind = kind_phys
intent = inout
@@ -612,6 +612,32 @@
dimensions = ()
type = logical
intent = in
+[dust_moist_opt_in]
+ standard_name = control_for_dust_soil_moisture_option
+ long_name = smoke dust moisture parameterization 1 - fecan 2 - shao
+ units = index
+ dimensions = ()
+ type = integer
+ active = (do_smoke_coupling)
+ intent = in
+[dust_moist_correction_in]
+ standard_name = dust_moist_correction_fengsha_dust_scheme
+ long_name = moisture correction term for fengsha dust emission
+ units = none
+ dimensions = ()
+ type = real
+ kind = kind_phys
+ active = (do_smoke_coupling)
+ intent = in
+[dust_drylimit_factor_in]
+ standard_name = dust_drylimit_factor_fengsha_dust_scheme
+ long_name = moisture correction term for drylimit in fengsha dust emission
+ units = none
+ dimensions = ()
+ type = real
+ kind = kind_phys
+ active = (do_smoke_coupling)
+ intent = in
[dust_alpha_in]
standard_name = alpha_fengsha_dust_scheme
long_name = alpha paramter for fengsha dust scheme
| | |