!> \file gfdl_sfc_layer.f !! This file contains ... !> This module contains the CCPP-compliant GFDL surface layer scheme. module gfdl_sfc_layer use machine , only : kind_phys implicit none public :: gfdl_sfc_layer_init, gfdl_sfc_layer_run private contains !> \section arg_table_gfdl_sfc_layer_init Argument Table !! \htmlinclude gfdl_sfc_layer_init.html !! subroutine gfdl_sfc_layer_init (icoef_sf, cplwav, cplwav2atm, lcurr_sf, & pert_cd, ntsflg, errmsg, errflg) implicit none integer, intent(in) :: icoef_sf, ntsflg logical, intent(in) :: cplwav, cplwav2atm, lcurr_sf, pert_cd character(len=*), intent(out) :: errmsg integer, intent(out) :: errflg ! Initialize CCPP error handling variables errmsg = '' errflg = 0 #if HWRF==1 write(errmsg,'(*(a))') 'The GFDL surface layer scheme does not support '& //'use of the HWRF preprocessor flag in gfdl_sfc_layer.F90' errflg = 1 return #endif if (icoef_sf < 0 .or. icoef_sf > 8) then write(errmsg,'(*(a))') 'The value of icoef_sf is outside of the ' & //'supported range (0-8) in gfdl_sfc_layer.F90' errflg = 1 return end if if (cplwav .or. cplwav2atm) then write(errmsg,'(*(a))') 'The GFDL surface layer scheme is not set up ' & //'to be coupled to waves in gfdl_sfc_layer.F90' errflg = 1 return end if if (lcurr_sf) then write(errmsg,'(*(a))') 'The GFDL surface layer scheme is not set up ' & //'to be used with the lcurr_sf option in gfdl_sfc_layer.F90' errflg = 1 return end if if (pert_cd) then write(errmsg,'(*(a))') 'The GFDL surface layer scheme is not set up ' & //'to be used with the pert_cd option in gfdl_sfc_layer.F90' errflg = 1 return end if if (ntsflg > 0) then !GJF: In order to enable ntsflg > 0, the variable 'tstrc' passed into MFLUX2 should be set ! to the surface_skin_temperature_over_X_interstitial rather than the average of it and ! surface_skin_temperature_after_iteration_over_X write(errmsg,'(*(a))') 'Setting ntsflg > 0 is currently not supported'& //' in gfdl_sfc_layer.F90' errflg = 1 return end if !GJF: Initialization notes: In WRF, the subroutine module_sf_myjsfc/myjsfcinit ! is called for initialization of the GFDL surface layer scheme from ! the module_physics_init subroutine. It contains the following ! initializations which should already have been done by other ! code in UFS-related host models: ! IF(.NOT.RESTART)THEN ! DO J=JTS,JTE ! DO I=ITS,ITF ! USTAR(I,J)=0.1 ! ENDDO ! ENDDO ! ENDIF !also initialize surface roughness length end subroutine gfdl_sfc_layer_init !> \section arg_table_gfdl_sfc_layer_run Argument Table !! \htmlinclude gfdl_sfc_layer_run.html !! subroutine gfdl_sfc_layer_run (im, nsoil, km, xlat, xlon, flag_iter, lsm, & lsm_noah, lsm_noahmp, lsm_ruc, icoef_sf, cplwav, karman, & cplwav2atm, lcurr_sf, pert_Cd, ntsflg, sfenth, z1, shdmax, ivegsrc, & vegtype, sigmaf, dt, wet, dry, icy, isltyp, rd, grav, ep1, ep2, smois, & psfc, prsl1, q1, t1, u1, v1, wspd, u10, v10, gsw, glw, tsurf_wat, & tsurf_lnd, tsurf_ice, tskin_wat, tskin_lnd, tskin_ice, ustar_wat, & ustar_lnd, ustar_ice, znt_wat, znt_lnd, znt_ice, cdm_wat, cdm_lnd, & cdm_ice, stress_wat, stress_lnd, stress_ice, rib_wat, rib_lnd, rib_ice, & fm_wat, fm_lnd, fm_ice, fh_wat, fh_lnd, fh_ice, fh2_wat, fh2_lnd, & fh2_ice, ch_wat, ch_lnd, ch_ice, fm10_wat, fm10_lnd, fm10_ice, qss_wat, & qss_lnd, qss_ice, errmsg, errflg) use funcphys, only: fpvs !#### GJF: temporarily grab parameters from LSM-specific modules -- should go through CCPP #### ! (fixing this involves replacing the functionality of set_soilveg and namelist_soilveg) use namelist_soilveg, only: maxsmc_noah => maxsmc, drysmc_noah => drysmc use namelist_soilveg_ruc, only: maxsmc_ruc => maxsmc, drysmc_ruc => drysmc use noahmp_tables, only: maxsmc_noahmp => smcmax_table, drysmc_noahmp => smcdry_table !################################################################################################ implicit none integer, intent(in) :: im, nsoil, km, ivegsrc integer, intent(in) :: lsm, lsm_noah, lsm_noahmp, & lsm_ruc, icoef_sf, ntsflg logical, intent(in) :: cplwav, cplwav2atm !GJF: this scheme has not been tested with these on logical, intent(in) :: lcurr_sf !GJF: this scheme has not been tested with this option turned on; the variables scurx and scury need to be input in order to use this logical, intent(in) :: pert_Cd !GJF: this scheme has not been tested with this option turned on; the variables ens_random_seed and ens_Cdamp need to be input in order to use this logical, dimension(:), intent(in) :: flag_iter, wet, dry, icy integer, dimension(:), intent(in) :: isltyp, vegtype real(kind=kind_phys), intent(in) :: dt, sfenth real(kind=kind_phys), intent(in) :: rd,grav,ep1,ep2 real(kind=kind_phys), dimension(:,:), intent(in) :: smois real(kind=kind_phys), dimension(:), intent(in) :: psfc, prsl1, & q1, t1, u1, v1, wspd, u10, v10, gsw, glw, z1, shdmax, sigmaf, xlat, & xlon, tsurf_wat, tsurf_lnd, tsurf_ice real(kind=kind_phys), intent(inout), dimension(:) :: tskin_wat, & tskin_lnd, tskin_ice, ustar_wat, ustar_lnd, ustar_ice, & znt_wat, znt_lnd, znt_ice, cdm_wat, cdm_lnd, cdm_ice, & stress_wat, stress_lnd, stress_ice, rib_wat, rib_lnd, rib_ice, & fm_wat, fm_lnd, fm_ice, fh_wat, fh_lnd, fh_ice, fh2_wat, fh2_lnd, & fh2_ice, ch_wat, ch_lnd, ch_ice, fm10_wat, fm10_lnd, fm10_ice, & qss_wat, qss_lnd, qss_ice character(len=*), intent(out) :: errmsg integer, intent(out) :: errflg !local variables integer :: i, its, ite, ims, ime logical :: ch_bound_excursion !GJF: the vonKarman constant should come in through the CCPP and be defined by the host model real (kind=kind_phys), intent(in) :: karman real (kind=kind_phys), parameter :: log01=log(0.01), log05=log(0.05), & log07=log(0.07) !GJF: if the following variables will be used, they should be turned into intent(in) namelist options integer :: iwavecpl, ens_random_seed, issflx logical :: diag_wind10m, diag_qss real(kind=kind_phys) :: ens_Cdamp real(kind=kind_phys), dimension(im) :: wetc, pspc, pkmax, tstrc, upc, & vpc, mznt, slwdc, wind10, qfx, qgh, zkmax, z1_cm, z0max, ztmax real(kind=kind_phys), dimension(im) :: u10_lnd, u10_ocn, u10_ice, & v10_lnd, v10_ocn, v10_ice !GJF: the following variables are identified as: !"SCURX" "Surface Currents(X)" "m s-1" !"SCURY" "Surface Currents(Y)" "m s-1 !"CHARN" "Charnock Coeff" " " !"MSANG" "Wind/Stress Angle" "Radian" real(kind=kind_phys), dimension(im) :: charn, msang, scurx, scury real(kind=kind_phys), dimension(im) :: fxh, fxe, fxmx, fxmy, xxfh, & xxfh2, tzot real(kind=kind_phys), dimension(1:30) :: maxsmc, drysmc real(kind=kind_phys) :: smcmax, smcdry, zhalf, cd10, & esat, fm_lnd_old, fh_lnd_old, tem1, tem2, czilc, cd_low_limit, & cd_high_limit, ch_low_limit, ch_high_limit, fh2_fh_ratio !#### This block will become unnecessary when maxsmc and drysmc come through the CCPP #### if (lsm == lsm_noah) then maxsmc = maxsmc_noah drysmc = drysmc_noah else if (lsm == lsm_noahmp) then maxsmc = maxsmc_noahmp drysmc = drysmc_noahmp else if (lsm == lsm_ruc) then maxsmc = maxsmc_ruc drysmc = drysmc_ruc else !GJF: These data were from the original GFDL surface layer scheme, but ! rather than being hard-coded here, they should be shared with the ! LSM. These data are kept for legacy purposes. Note that these only ! have nonzero values for 16 soil types vs 19 for other STAS datasets data maxsmc/0.339, 0.421, 0.434, 0.476, 0.476, 0.439, & 0.404, 0.464, 0.465, 0.406, 0.468, 0.468, & 0.439, 1.000, 0.200, 0.421, 0.000, 0.000, & 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, & 0.000, 0.000, 0.000, 0.000, 0.000, 0.000/ data drysmc/0.010, 0.028, 0.047, 0.084, 0.084, 0.066, & 0.067, 0.120, 0.103, 0.100, 0.126, 0.138, & 0.066, 0.000, 0.006, 0.028, 0.000, 0.000, & 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, & 0.000, 0.000, 0.000, 0.000, 0.000, 0.000/ end if !######################################################################## !GJF: This code has not been tested with iwavecpl = 1; the variables 'charn' and 'msang' (and others?) need to be input in order to use this ! if (cplwav .or. cplwav2atm) then ! iwavecpl = 1 ! else ! iwavecpl = 0 ! end if iwavecpl = 0 !GJF: temporary setting of variables that should be moved to namelist is they are used ens_random_seed = 0 !used for HWRF ensemble? ens_Cdamp = 0.0 !used for HWRF ensemble? issflx = 0 !GJF: 1 = calculate surface fluxes, 0 = don't diag_wind10m = .false. !GJF: if one wants 10m wind speeds to come from this scheme, set this to True, ! put [u,v]10_[lnd/ocn/ice] in the scheme argument list (and metadata), and modify ! GFS_surface_compsites to receive the individual components and calculate an all-grid value diag_qss = .false. !GJF: saturation specific humidities are calculated by LSM, sea surface, and sea ice schemes in ! GFS-based suites ! Initialize CCPP error handling variables errmsg = '' errflg = 0 its = 1 ims = 1 ite = im ime = im do i=its, ite if (flag_iter(i)) then !GJF: Perform data preparation that is the same for all surface types pspc(i) = psfc(i)*10. ! convert from Pa to cgs pkmax(i) = prsl1(i)*10. ! convert from Pa to cgs upc(i) = u1(i)*100. ! convert from m s-1 to cm s-1 vpc(i) = v1(i)*100. ! convert from m s-1 to cm s-1 !Wang: use previous u10 v10 to compute wind10, input to MFLUX2 to compute z0 (for first time step, u10 and v10 may be zero) wind10(i)=sqrt(u10(i)*u10(i)+v10(i)*v10(i)) !m s-1 !Wang: calulate height of the first half level ! if (wind10(i) <= 1.0e-10 .or. wind10(i) > 150.0) then ! zhalf = -rd*t1(i)*alog(pkmax(i)/pspc(i))/grav !m ! endif !GJF: rather than calculate the height of the first half level, if it is precalculated ! in a different scheme, pass it in and use it; note that in FV3, calculating via the hypsometric equation ! occasionally produced values much shallower than those passed in !zkmax(i) = -rd*t1(i)*alog(pkmax(i)/pspc(i))/grav !m zkmax(i) = z1(i) z1_cm(i) = 100.0*z1(i) !GJF: these drag coefficient limits were suggested by Chunxi Zhang via his module_sf_sfclayrev.f90 cd_low_limit = 1.0e-5/zkmax(i) cd_high_limit = 0.1 !GJF: use the lower of 0.1 from Chunxi Zhang or 0.05/wspd from WRF's module_sf_gfdl.F ! (this will always be the latter if wspd has a minimum of 1.0 m s-1 from above) ch_low_limit = cd_low_limit ch_high_limit = min(0.1,0.05/wspd(i)) !slwdc... GFDL downward net flux in units of cal/(cm**2/min) !also divide by 10**4 to convert from /m**2 to /cm**2 slwdc(i)=gsw(i)+glw(i) slwdc(i)=0.239*60.*slwdc(i)*1.e-4 !GJF: these variables should be passed in if these options are used charn(i) = 0.0 !used with wave coupling (iwavecpl == 1) msang(i) = 0.0 !used with wave coupling (iwavecpl == 1) scurx(i) = 0.0 !used with ocean currents? (lcurr_sf == T) scury(i) = 0.0 !used with ocean currents? (lcurr_sf == T) if (diag_qss) then esat = fpvs(t1(i)) qgh(i) = ep2*esat/(psfc(i)-esat) end if !GJF: these vars are not needed in a GFS-based suite !rho1(i)=prsl1(i)/(rd*t1(i)*(1.+ep1*q1(i))) !cpm(i)=cp*(1.+0.8*q1(i)) !GJF: perform data preparation that depends on surface types and call the mflux2 subroutine for each surface type ! Note that this is different than the original WRF module_sf_gfdl.F where mflux2 is called once for all surface ! types, with negative roughness lengths denoting open ocean. if (dry(i)) then !GJF: from WRF's module_sf_gfdl.F smcdry=drysmc(isltyp(i)) smcmax=maxsmc(isltyp(i)) wetc(i)=(smois(i,1)-smcdry)/(smcmax-smcdry) wetc(i)=amin1(1.,amax1(wetc(i),0.)) !GJF: the lower boundary temperature passed in to MFLUX2 either follows GFS: tstrc(i) = 0.5*(tskin_lnd(i) + tsurf_lnd(i)) !averaging tskin_lnd and tsurf_lnd as in GFS surface layer breaks ntsflg functionality !GJF: or WRF module_sf_gfdl.F: !tstrc(i) = tskin_lnd(i) !GJF: Roughness Length Limitation section ! The WRF version of module_sf_gfdl.F has no checks on the roughness lengths prior to entering MFLUX2. ! The following limits were placed on roughness lengths from the GFS surface layer scheme at the suggestion ! of Chunxi Zhang. Using the GFDL surface layer without such checks can lead to instability in the UFS. !znt_lnd is in cm, z0max/ztmax are in m at this point z0max(i) = max(1.0e-6, min(0.01 * znt_lnd(i), zkmax(i))) tem1 = 1.0 - shdmax(i) tem2 = tem1 * tem1 tem1 = 1.0 - tem2 if( ivegsrc == 1 ) then if (vegtype(i) == 10) then z0max(i) = exp( tem2*log01 + tem1*log07 ) elseif (vegtype(i) == 6) then z0max(i) = exp( tem2*log01 + tem1*log05 ) elseif (vegtype(i) == 7) then ! z0max(i) = exp( tem2*log01 + tem1*log01 ) z0max(i) = 0.01 elseif (vegtype(i) == 16) then ! z0max(i) = exp( tem2*log01 + tem1*log01 ) z0max(i) = 0.01 else z0max(i) = exp( tem2*log01 + tem1*log(z0max(i)) ) endif elseif (ivegsrc == 2 ) then if (vegtype(i) == 7) then z0max(i) = exp( tem2*log01 + tem1*log07 ) elseif (vegtype(i) == 8) then z0max(i) = exp( tem2*log01 + tem1*log05 ) elseif (vegtype(i) == 9) then ! z0max(i) = exp( tem2*log01 + tem1*log01 ) z0max(i) = 0.01 elseif (vegtype(i) == 11) then ! z0max(i) = exp( tem2*log01 + tem1*log01 ) z0max(i) = 0.01 else z0max(i) = exp( tem2*log01 + tem1*log(z0max(i)) ) endif endif z0max(i) = max(z0max(i), 1.0e-6) ! czilc = 10.0 ** (- (0.40/0.07) * z0) ! fei's canopy height dependance of czil czilc = 0.8 tem1 = 1.0 - sigmaf(i) ztmax(i) = z0max(i)*exp( - tem1*tem1 & & * czilc*karman*sqrt(ustar_lnd(i)*(0.01/1.5e-05))) ztmax(i) = max(ztmax(i), 1.0e-6) !GJF: from WRF's module_sf_gfdl.F if (wind10(i) <= 1.0e-10 .or. wind10(i) > 150.0) then wind10(i)=wspd(i)*alog(10.0/z0max(i))/alog(z1(i)/z0max(i)) !m s-1 end if wind10(i)=wind10(i)*100.0 !convert from m/s to cm/s ztmax(i) = ztmax(i)*100.0 !convert from m to cm z0max(i) = z0max(i)*100.0 !convert from m to cm call mflux2 (fxh(i), fxe(i), fxmx(i), fxmy(i), cdm_lnd(i:i), rib_lnd(i:i), & xxfh(i), ztmax(i), z0max(i), tstrc(i), & pspc(i), pkmax(i), wetc(i), slwdc(i), z1_cm(i), icoef_sf, iwavecpl, lcurr_sf, charn(i), msang(i), & scurx(i), scury(i), pert_Cd, ens_random_seed, ens_Cdamp, upc(i), vpc(i), t1(i:i), q1(i:i), & dt, wind10(i), xxfh2(i), ntsflg, sfenth, tzot(i), ep2, errmsg, & errflg) if (errflg /= 0) return !GJF: this is broken when tstrc is set to an average of two variables if (ntsflg==1) then tskin_lnd(i) = tstrc(i) ! gopal's doing end if if (diag_wind10m) then u10_lnd(i) = u1(i)*(0.01*wind10(i)/wspd(i)) v10_lnd(i) = v1(i)*(0.01*wind10(i)/wspd(i)) end if !GJF: these variables are not needed in a GFS-based suite, but are found in WRF's module_sf_gfdl.F and kept in comments for legacy !gz1oz0(i) = alog(zkmax(i)/(0.01*znt_lnd(i))) !taux(i) = fxmx(i)/10. ! gopal's doing for Ocean coupling !tauy(i) = fxmy(i)/10. ! gopal's doing for Ocean coupling cdm_lnd(i) = max(cdm_lnd(i), cd_low_limit) cdm_lnd(i) = min(cdm_lnd(i), cd_high_limit) fm_lnd(i) = karman/sqrt(cdm_lnd(i)) !1) try fh_lnd from MFLUX2 fh_lnd(i) = karman*xxfh(i) !2) calc ch_lnd from fm_lnd and fh_lnd ch_lnd(i) = karman*karman/(fm_lnd(i) * fh_lnd(i)) !3) check if ch_lnd is out of bounds (if so, recalculate fh_lnd from bounded value) ch_bound_excursion = .false. if (ch_lnd(i) < ch_low_limit) then ch_bound_excursion = .true. ch_lnd(i) = ch_low_limit else if (ch_lnd(i) > ch_high_limit) then ch_bound_excursion = .true. ch_lnd(i) = ch_high_limit end if fh2_lnd(i) = karman*xxfh2(i) if (ch_bound_excursion) then fh2_fh_ratio = min(xxfh2(i)/xxfh(i), 1.0) fh_lnd(i) = karman*karman/(fm_lnd(i)*ch_lnd(i)) fh2_lnd(i) = fh2_fh_ratio*fh_lnd(i) end if !GJF: Other CCPP schemes (PBL) ask for fm/fh instead of psim/psih !psim_lnd(i)=gz1oz0(i)-fm_lnd(i) !psih_lnd(i)=gz1oz0(i)-fh_lnd(i) !GJF: from WRF's module_sf_gfdl.F ustar_lnd(i) = 0.01*sqrt(cdm_lnd(i)* & (upc(i)*upc(i) + vpc(i)*vpc(i))) !GJF: from Chunxi Zhang's module_sf_sfclayrev.f90 (I'm not sure it's necessary.) ustar_lnd(i) = amax1(ustar_lnd(i),0.001) stress_lnd(i) = cdm_lnd(i)*wspd(i)*wspd(i) !GJF: from WRF's module_sf_gfdl.F ! convert cd, ch to values at 10m, for output cd10 = cdm_lnd(i) if ( wind10(i) .ge. 0.1 ) then cd10=cdm_lnd(i)* (wspd(i)/(0.01*wind10(i)) )**2 !tmp9=0.01*abs(tzot(i)) !ch_out(i)=ch_lnd(i)*(wspd(i)/(0.01*wind10(i)) ) * & ! (alog(zkmax(i)/tmp9)/alog(10.0/tmp9)) end if fm10_lnd(i) = karman/sqrt(cd10) !GJF: conductances aren't used in other CCPP schemes, but this limit ! might be able to replace the limits on drag coefficients above !chs_lnd(i)=ch_lnd(i)*wspd (i) !conductance !chs2_lnd(i)=ustar_lnd(i)*karman/fh2_lnd(i) !2m conductance !!!2014-0922 cap CHS over land points ! chs_lnd(i)=amin1(chs_lnd(i), 0.05) ! chs2_lnd(i)=amin1(chs2_lnd(i), 0.05) ! if (chs2_lnd(i) < 0) chs2_lnd(i)=1.0e-6 if (diag_qss) then esat = fpvs(tskin_lnd(i)) qss_lnd(i) = ep2*esat/(psfc(i)-esat) end if !GJF: not used in CCPP !flhc_lnd(i)=cpm(i)*rho1(i)*chs_lnd(i) !flqc_lnd(i)=rho1(i)*chs_lnd(i) !cqs2_lnd(i)=chs2_lnd(i) end if !dry if (icy(i)) then !GJF: from WRF's module_sf_gfdl.F smcdry=drysmc(isltyp(i)) smcmax=maxsmc(isltyp(i)) wetc(i)=(smois(i,1)-smcdry)/(smcmax-smcdry) wetc(i)=amin1(1.,amax1(wetc(i),0.)) !GJF: the lower boundary temperature passed in to MFLUX2 either follows GFS: tstrc(i) = 0.5*(tskin_ice(i) + tsurf_ice(i)) !averaging tskin_ice and tsurf_ice as in GFS surface layer breaks ntsflg functionality !GJF: or WRF module_sf_gfdl.F: !tstrc(i) = tskin_ice(i) !averaging tskin_ice and tsurf_ice as in GFS surface layer breaks ntsflg functionality !GJF: Roughness Length Limitation section ! The WRF version of module_sf_gfdl.F has no checks on the roughness lengths prior to entering MFLUX2. ! The following limits were placed on roughness lengths from the GFS surface layer scheme at the suggestion ! of Chunxi Zhang. Using the GFDL surface layer without such checks can lead to instability in the UFS. !znt_ice is in cm, z0max/ztmax are in m at this point z0max(i) = max(1.0e-6, min(0.01 * znt_ice(i), zkmax(i))) !** xubin's new z0 over land and sea ice tem1 = 1.0 - shdmax(i) tem2 = tem1 * tem1 tem1 = 1.0 - tem2 if( ivegsrc == 1 ) then z0max(i) = exp( tem2*log01 + tem1*log(z0max(i)) ) elseif (ivegsrc == 2 ) then z0max(i) = exp( tem2*log01 + tem1*log(z0max(i)) ) endif z0max(i) = max(z0max(i), 1.0e-6) ! czilc = 10.0 ** (- (0.40/0.07) * z0) ! fei's canopy height ! dependance of czil czilc = 0.8 tem1 = 1.0 - sigmaf(i) ztmax(i) = z0max(i)*exp( - tem1*tem1 & & * czilc*karman*sqrt(ustar_ice(i)*(0.01/1.5e-05))) ztmax(i) = max(ztmax(i), 1.0e-6) !GJF: from WRF's module_sf_gfdl.F if (wind10(i) <= 1.0e-10 .or. wind10(i) > 150.0) then wind10(i)=wspd(i)*alog(10.0/z0max(i))/alog(z1(i)/z0max(i)) end if wind10(i)=wind10(i)*100.0 !! m/s to cm/s ztmax(i) = ztmax(i)*100.0 !m to cm z0max(i) = z0max(i)*100.0 !m to cm call mflux2 (fxh(i), fxe(i), fxmx(i), fxmy(i), cdm_ice(i:i), rib_ice(i:i), & xxfh(i), ztmax(i), z0max(i), tstrc(i), & pspc(i), pkmax(i), wetc(i), slwdc(i), z1_cm(i), icoef_sf, iwavecpl, lcurr_sf, charn(i), msang(i), & scurx(i), scury(i), pert_Cd, ens_random_seed, ens_Cdamp, upc(i), vpc(i), t1(i:i), q1(i:i), & dt, wind10(i), xxfh2(i), ntsflg, sfenth, tzot(i), ep2, errmsg, & errflg) if (errflg /= 0) return !GJF: this is broken when tstrc is set to an average of two variables if (ntsflg==1) then tskin_ice(i) = tstrc(i) ! gopal's doing end if if (diag_wind10m) then u10_ice(i) = u1(i)*(0.01*wind10(i)/wspd(i)) v10_ice(i) = v1(i)*(0.01*wind10(i)/wspd(i)) end if !GJF: these variables are not needed in a GFS-based suite, but are found in WRF's module_sf_gfdl.F and kept in comments for legacy !gz1oz0(i) = alog(zkmax(i)/znt_ice(i)) !taux(i) = fxmx(i)/10. ! gopal's doing for Ocean coupling !tauy(i) = fxmy(i)/10. ! gopal's doing for Ocean coupling cdm_ice(i) = max(cdm_ice(i), cd_low_limit) cdm_ice(i) = min(cdm_ice(i), cd_high_limit) fm_ice(i) = karman/sqrt(cdm_ice(i)) !1) try fh_ice from MFLUX2 fh_ice(i) = karman*xxfh(i) !2) calc ch_ice from fm_ice and fh_ice ch_ice(i) = karman*karman/(fm_ice(i) * fh_ice(i)) !3) check if ch_ice is out of bounds (if so, recalculate fh_ice from bounded value) ch_bound_excursion = .false. if (ch_ice(i) < ch_low_limit) then ch_bound_excursion = .true. ch_ice(i) = ch_low_limit else if (ch_ice(i) > ch_high_limit) then ch_bound_excursion = .true. ch_ice(i) = ch_high_limit end if fh2_ice(i) = karman*xxfh2(i) if (ch_bound_excursion) then fh2_fh_ratio = min(xxfh2(i)/xxfh(i), 1.0) fh_ice(i) = karman*karman/(fm_ice(i)*ch_ice(i)) fh2_ice(i) = fh2_fh_ratio*fh_ice(i) end if !Other CCPP schemes (PBL) ask for fm/fh instead of psim/psih !psim_ice(i)=gz1oz0(i)-fm_ice(i) !psih_ice(i)=gz1oz0(i)-fh_ice(i) ustar_ice(i) = 0.01*sqrt(cdm_ice(i)* & (upc(i)*upc(i) + vpc(i)*vpc(i))) !GJF: from Chunxi Zhang's module_sf_sfclayrev.f90 (I'm not sure it's necessary.) ustar_ice(i) = amax1(ustar_ice(i),0.001) stress_ice(i) = cdm_ice(i)*wspd(i)*wspd(i) !GJF: from WRF's module_sf_gfdl.F !!! convert cd, ch to values at 10m, for output cd10 = cdm_ice(i) if ( wind10(i) .ge. 0.1 ) then cd10=cdm_ice(i)* (wspd(i)/(0.01*wind10(i)) )**2 !tmp9=0.01*abs(tzot(i)) !ch_out(i)=ch_ice(i)*(wspd(i)/(0.01*wind10(i)) ) * & ! (alog(zkmax(i)/tmp9)/alog(10.0/tmp9)) end if fm10_ice(i) = karman/sqrt(cd10) !GJF: conductances aren't used in other CCPP schemes !chs_ice(i)=ch_ice(i)*wspd (i) !conductance !chs2_ice(i)=ustar_ice(i)*karman/fh2_ice(i) !2m conductance if (diag_qss) then esat = fpvs(tskin_ice(i)) qss_ice(i) = ep2*esat/(psfc(i)-esat) end if !flhc_ice(i)=cpm(i)*rho1(i)*chs_ice(i) !flqc_ice(i)=rho1(i)*chs_ice(i) !cqs2_ice(i)=chs2_ice(i) end if !ice if (wet(i)) then wetc(i) = 1.0 !GJF: the lower boundary temperature passed in to MFLUX2 either follows GFS: tstrc(i) = 0.5*(tskin_wat(i) + tsurf_wat(i)) !averaging tskin_wat and tsurf_wat as in GFS surface layer breaks ntsflg functionality !GJF: or WRF module_sf_gfdl.F: !tstrc(i) = tskin_wat(i) ! DH* 20201009: these bounds on ocean roughness lengths are from Chunxi Zhang's module_sf_sfclayrev.f90 (in cm) znt_wat(i)=min(2.85e-1,max(znt_wat(i),1.27e-5)) !GJF: from WRF's module_sf_gfdl.F if (wind10(i) <= 1.0e-10 .or. wind10(i) > 150.0) then wind10(i)=wspd(i)*alog(10.0/(0.01*znt_wat(i)))/alog(z1(i)/(0.01*znt_wat(i))) end if wind10(i)=wind10(i)*100.0 !! m/s to cm/s !GJF: mflux2 expects negative roughness length for ocean points znt_wat(i) = -znt_wat(i) call mflux2 (fxh(i), fxe(i), fxmx(i), fxmy(i), cdm_wat(i:i), rib_wat(i:i), & xxfh(i), znt_wat(i:i), mznt(i), tstrc(i), & pspc(i), pkmax(i), wetc(i), slwdc(i), z1_cm(i), icoef_sf, iwavecpl, lcurr_sf, charn(i), msang(i), & scurx(i), scury(i), pert_Cd, ens_random_seed, ens_Cdamp, upc(i), vpc(i), t1(i:i), q1(i:i), & dt, wind10(i), xxfh2(i), ntsflg, sfenth, tzot(i), ep2, errmsg, & errflg) if (errflg /= 0) return !GJF: this is broken when tstrc is set to an average of two variables if (ntsflg==1) then tskin_wat(i) = tstrc(i) ! gopal's doing end if znt_wat(i)= abs(znt_wat(i)) mznt(i)= abs(mznt(i)) !GJF: these bounds on ocean roughness lengths are from Chunxi Zhang's module_sf_sfclayrev.f90 (in cm) znt_wat(i)=min(2.85e-1,max(znt_wat(i),1.27e-5)) if (diag_wind10m) then u10_ocn(i) = u1(i)*(0.01*wind10(i)/wspd(i)) v10_ocn(i) = v1(i)*(0.01*wind10(i)/wspd(i)) end if !GJF: these variables are not needed in a GFS-based suite, but are found in WRF's module_sf_gfdl.F and kept in comments for legacy !gz1oz0(i) = alog(zkmax(i)/znt_wat(i)) !taux(i) = fxmx(i)/10. ! gopal's doing for Ocean coupling !tauy(i) = fxmy(i)/10. ! gopal's doing for Ocean coupling cdm_wat(i) = max(cdm_wat(i), cd_low_limit) cdm_wat(i) = min(cdm_wat(i), cd_high_limit) fm_wat(i) = karman/sqrt(cdm_wat(i)) !1) try fh_wat from MFLUX2 fh_wat(i) = karman*xxfh(i) !2) calc ch_wat from fm_wat and fh_wat ch_wat(i) = karman*karman/(fm_wat(i) * fh_wat(i)) !3) check if ch_lnd is out of bounds (if so, recalculate fh_lnd from bounded value) ch_bound_excursion = .false. if (ch_wat(i) < ch_low_limit) then ch_bound_excursion = .true. ch_wat(i) = ch_low_limit else if (ch_wat(i) > ch_high_limit) then ch_bound_excursion = .true. ch_wat(i) = ch_high_limit end if fh2_wat(i) = karman*xxfh2(i) if (ch_bound_excursion) then fh2_fh_ratio = min(xxfh2(i)/xxfh(i), 1.0) fh_wat(i) = karman*karman/(fm_wat(i)*ch_wat(i)) fh2_wat(i) = fh2_fh_ratio*fh_wat(i) end if !Other CCPP schemes (PBL) ask for fm/fh instead of psim/psih !psim_ocn(i)=gz1oz0(i)-fm_wat(i) !psih_ocn(i)=gz1oz0(i)-fh_wat(i) ustar_wat(i) = 0.01*sqrt(cdm_wat(i)* & (upc(i)*upc(i) + vpc(i)*vpc(i))) !GJF: from Chunxi Zhang's module_sf_sfclayrev.f90 (I'm not sure it's necessary.) ustar_wat(i) = amax1(ustar_wat(i),0.001) stress_wat(i) = cdm_wat(i)*wspd(i)*wspd(i) !GJF: from WRF's module_sf_gfdl.F !!! convert cd, ch to values at 10m, for output cd10 = cdm_wat(i) if ( wind10(i) .ge. 0.1 ) then cd10=cdm_wat(i)* (wspd(i)/(0.01*wind10(i)) )**2 !tmp9=0.01*abs(tzot(i)) !ch_out(i)=ch_wat(i)*(wspd(i)/(0.01*wind10(i)) ) * & ! (alog(zkmax(i)/tmp9)/alog(10.0/tmp9)) end if fm10_wat(i) = karman/sqrt(cd10) !GJF: conductances aren't used in other CCPP schemes !chs_ocn(i)=ch_wat(i)*wspd (i) !conductance !chs2_ocn(i)=ustar_wat(i)*karman/fh2_wat(i) !2m conductance if (diag_qss) then esat = fpvs(tskin_wat(i)) qss_wat(i) = ep2*esat/(psfc(i)-esat) end if end if !wet !flhc_ocn(i)=cpm(i)*rho1(i)*chs_ocn(i) !flqc_ocn(i)=rho1(i)*chs_ocn(i) !cqs2_ocn(i)=chs2_ocn(i) end if !flag_iter end do !GJF: this code has not been updated since GFS suites don't require this; one would need to have different values of hfx, qfx, lh for each surface type ! if (isfflx.eq.0) then ! do i=its,ite ! hfx(i)=0. ! lh(i)=0. ! qfx(i)=0. ! enddo ! else ! do i=its,ite ! if(islmsk == 0) then ! !water ! hfx(i)= -10.*cp*fxh(i) ! else if (islmsk == 1) then ! hfx(i)= -10.*cp*fxh(i) ! hfx(i)=amax1(hfx(i),-250.) ! end if ! qfx(j)=-10.*fxe(i) ! qfx(i)=amax1(qfx(i),0.) ! lh(i)=xlv*qfx(i) ! enddo ! endif end subroutine gfdl_sfc_layer_run !--------------------------------- !GJF (2020/04/21): The starting point for the MFLUX2 subroutine here was module_sf_gfdl.F in WRF SUBROUTINE MFLUX2( fxh,fxe,fxmx,fxmy,cdm,rib,xxfh,zoc,mzoc,tstrc, & !mzoc KWON pspc,pkmax,wetc,slwdc,z1, & icoef_sf,iwavecpl,lcurr_sf,alpha,gamma,xcur,ycur, & pert_Cd, ens_random_seed, ens_Cdamp, & upc,vpc,tpc,rpc,dt,wind10,xxfh2,ntsflg,sfenth, & tzot, ep2, errmsg, errflg) !------------------------------------------------------------------------ ! ! MFLUX2 computes surface fluxes of momentum, heat,and moisture ! using monin-obukhov. the roughness length "z0" is prescribed ! over land and over ocean "z0" is computed using charnocks formula. ! the universal functions (from similarity theory approach) are ! those of hicks. This is Bob's doing. ! !------------------------------------------------------------------------ USE module_sf_exchcoef IMPLICIT NONE !----------------------------------------------------------------------- ! user interface variables !----------------------------------------------------------------------- !GJF: This subroutine was converted to expect data from a single point instead of a horizontal array to accommodate a fractional landmask !integer,intent(in) :: ims,ime !integer,intent(in) :: its,ite integer, parameter :: ims = 1 integer, parameter :: ime = 1 integer, parameter :: its = 1 integer, parameter :: ite = 1 integer,intent(in) :: ntsflg integer,intent(in) :: icoef_sf integer,intent(in) :: iwavecpl logical,intent(in) :: lcurr_sf logical,intent(in) :: pert_Cd integer,intent(in) :: ens_random_seed real(kind=kind_phys),intent(in) :: ens_Cdamp real(kind=kind_phys), intent (out), dimension (ims :ime ) :: fxh real(kind=kind_phys), intent (out), dimension (ims :ime ) :: fxe real(kind=kind_phys), intent (out), dimension (ims :ime ) :: fxmx real(kind=kind_phys), intent (out), dimension (ims :ime ) :: fxmy real(kind=kind_phys), intent (inout), dimension (ims :ime ) :: cdm ! real, intent (out), dimension (ims :ime ) :: cdm2 real(kind=kind_phys), intent (out), dimension (ims :ime ) :: rib real(kind=kind_phys), intent (out), dimension (ims :ime ) :: xxfh real(kind=kind_phys), intent (out), dimension (ims :ime ) :: xxfh2 real(kind=kind_phys), intent (out), dimension (ims :ime ) :: wind10 real(kind=kind_phys), intent ( inout), dimension (ims :ime ) :: zoc,mzoc !KWON real(kind=kind_phys), intent ( inout), dimension (ims :ime ) :: tzot !WANG real(kind=kind_phys), intent ( inout), dimension (ims :ime ) :: tstrc real(kind=kind_phys), intent ( in) :: dt real(kind=kind_phys), intent ( in) :: sfenth real(kind=kind_phys), intent ( in), dimension (ims :ime ) :: pspc real(kind=kind_phys), intent ( in), dimension (ims :ime ) :: pkmax real(kind=kind_phys), intent ( in), dimension (ims :ime ) :: wetc real(kind=kind_phys), intent ( in), dimension (ims :ime ) :: slwdc real(kind=kind_phys), intent ( in), dimension (ims :ime ) :: alpha, gamma real(kind=kind_phys), intent ( in), dimension (ims :ime ) :: xcur, ycur real(kind=kind_phys), intent ( in), dimension (ims :ime ) :: z1 real(kind=kind_phys), intent ( in), dimension (ims :ime ) :: upc real(kind=kind_phys), intent ( in), dimension (ims :ime ) :: vpc real(kind=kind_phys), intent ( in), dimension (ims :ime ) :: tpc real(kind=kind_phys), intent ( in), dimension (ims :ime ) :: rpc real(kind=kind_phys), intent ( in) :: ep2 character(len=*), intent(out) :: errmsg integer, intent(out) :: errflg !----------------------------------------------------------------------- ! internal variables !----------------------------------------------------------------------- integer, parameter :: icntx = 30 integer, dimension(1 :ime) :: ifz integer, dimension(1 :ime) :: indx integer, dimension(1 :ime) :: istb integer, dimension(1 :ime) :: it integer, dimension(1 :ime) :: iutb real(kind=kind_phys), dimension(1 :ime) :: aap real(kind=kind_phys), dimension(1 :ime) :: bq1 real(kind=kind_phys), dimension(1 :ime) :: bq1p real(kind=kind_phys), dimension(1 :ime) :: delsrad real(kind=kind_phys), dimension(1 :ime) :: ecof real(kind=kind_phys), dimension(1 :ime) :: ecofp real(kind=kind_phys), dimension(1 :ime) :: estso real(kind=kind_phys), dimension(1 :ime) :: estsop real(kind=kind_phys), dimension(1 :ime) :: fmz1 real(kind=kind_phys), dimension(1 :ime) :: fmz10 real(kind=kind_phys), dimension(1 :ime) :: fmz2 real(kind=kind_phys), dimension(1 :ime) :: fmzo1 real(kind=kind_phys), dimension(1 :ime) :: foft real(kind=kind_phys), dimension(1 :ime) :: foftm real(kind=kind_phys), dimension(1 :ime) :: frac real(kind=kind_phys), dimension(1 :ime) :: land real(kind=kind_phys), dimension(1 :ime) :: pssp real(kind=kind_phys), dimension(1 :ime) :: qf real(kind=kind_phys), dimension(1 :ime) :: rdiff real(kind=kind_phys), dimension(1 :ime) :: rho real(kind=kind_phys), dimension(1 :ime) :: rkmaxp real(kind=kind_phys), dimension(1 :ime) :: rstso real(kind=kind_phys), dimension(1 :ime) :: rstsop real(kind=kind_phys), dimension(1 :ime) :: sf10 real(kind=kind_phys), dimension(1 :ime) :: sf2 real(kind=kind_phys), dimension(1 :ime) :: sfm real(kind=kind_phys), dimension(1 :ime) :: sfzo real(kind=kind_phys), dimension(1 :ime) :: sgzm real(kind=kind_phys), dimension(1 :ime) :: slwa real(kind=kind_phys), dimension(1 :ime) :: szeta real(kind=kind_phys), dimension(1 :ime) :: szetam real(kind=kind_phys), dimension(1 :ime) :: t1 real(kind=kind_phys), dimension(1 :ime) :: t2 real(kind=kind_phys), dimension(1 :ime) :: tab1 real(kind=kind_phys), dimension(1 :ime) :: tab2 real(kind=kind_phys), dimension(1 :ime) :: tempa1 real(kind=kind_phys), dimension(1 :ime) :: tempa2 real(kind=kind_phys), dimension(1 :ime) :: theta real(kind=kind_phys), dimension(1 :ime) :: thetap real(kind=kind_phys), dimension(1 :ime) :: tsg real(kind=kind_phys), dimension(1 :ime) :: tsm real(kind=kind_phys), dimension(1 :ime) :: tsp real(kind=kind_phys), dimension(1 :ime) :: tss real(kind=kind_phys), dimension(1 :ime) :: ucom real(kind=kind_phys), dimension(1 :ime) :: uf10 real(kind=kind_phys), dimension(1 :ime) :: uf2 real(kind=kind_phys), dimension(1 :ime) :: ufh real(kind=kind_phys), dimension(1 :ime) :: ufm real(kind=kind_phys), dimension(1 :ime) :: ufzo real(kind=kind_phys), dimension(1 :ime) :: ugzm real(kind=kind_phys), dimension(1 :ime) :: uzeta real(kind=kind_phys), dimension(1 :ime) :: uzetam real(kind=kind_phys), dimension(1 :ime) :: vcom real(kind=kind_phys), dimension(1 :ime) :: vrtkx real(kind=kind_phys), dimension(1 :ime) :: vrts real(kind=kind_phys), dimension(1 :ime) :: wind real(kind=kind_phys), dimension(1 :ime) :: windp real(kind=kind_phys), dimension(1 :ime) :: wind10p !WANG, 10m wind previous step real(kind=kind_phys), dimension(1 :ime) :: uvs1 ! real(kind=kind_phys), dimension(1 :ime) :: xxfh real(kind=kind_phys), dimension(1 :ime) :: xxfm real(kind=kind_phys), dimension(1 :ime) :: xxsh real(kind=kind_phys), dimension(1 :ime) :: z10 real(kind=kind_phys), dimension(1 :ime) :: z2 real(kind=kind_phys), dimension(1 :ime) :: zeta real(kind=kind_phys), dimension(1 :ime) :: zkmax real(kind=kind_phys), dimension(1 :ime) :: pss real(kind=kind_phys), dimension(1 :ime) :: tstar real(kind=kind_phys), dimension(1 :ime) :: ukmax real(kind=kind_phys), dimension(1 :ime) :: vkmax real(kind=kind_phys), dimension(1 :ime) :: tkmax real(kind=kind_phys), dimension(1 :ime) :: rkmax real(kind=kind_phys), dimension(1 :ime) :: zot real(kind=kind_phys), dimension(1 :ime) :: fhzo1 real(kind=kind_phys), dimension(1 :ime) :: sfh real(kind=kind_phys) :: ux13, yo, y,xo,x,ux21,ugzzo,ux11,ux12,uzetao,xnum,alll real(kind=kind_phys) :: ux1,ugz,x10,uzo,uq,ux2,ux3,xtan,xden,y10,uzet1o,ugz10 real(kind=kind_phys) :: szet2, zal2,ugz2 real(kind=kind_phys) :: rovcp,boycon,cmo2,psps1,zog,enrca,rca,cmo1,amask,en,ca,a,c real(kind=kind_phys) :: sgz,zal10,szet10,fmz,szo,sq,fmzo,rzeta1,zal1g,szetao,rzeta2,zal2g real(kind=kind_phys) :: hcap,xks,pith,teps,diffot,delten,alevp,psps2,alfus,nstep real(kind=kind_phys) :: shfx,sigt4,reflect real(kind=kind_phys) :: cor1,cor2,szetho,zal2gh,cons_p000001,cons_7,vis,ustar,restar,rat real(kind=kind_phys) :: wndm,ckg real(kind=kind_phys) :: windmks,znott,znotm real(kind=kind_phys) :: ubot, vbot integer:: i,j,ii,iq,nnest,icnt,ngd,ip !----------------------------------------------------------------------- ! internal variables !----------------------------------------------------------------------- real(kind=kind_phys), dimension (223) :: tab real(kind=kind_phys), dimension (223) :: table real(kind=kind_phys), dimension (101) :: tab11 real(kind=kind_phys), dimension (41) :: table4 real(kind=kind_phys), dimension (42) :: tab3 real(kind=kind_phys), dimension (54) :: table2 real(kind=kind_phys), dimension (54) :: table3 real(kind=kind_phys), dimension (74) :: table1 real(kind=kind_phys), dimension (80) :: tab22 character(len=255) :: message equivalence (tab(1),tab11(1)) equivalence (tab(102),tab22(1)) equivalence (tab(182),tab3(1)) equivalence (table(1),table1(1)) equivalence (table(75),table2(1)) equivalence (table(129),table3(1)) equivalence (table(183),table4(1)) data amask/ -98.0/ !----------------------------------------------------------------------- ! tables used to obtain the vapor pressures or saturated vapor ! pressure !----------------------------------------------------------------------- data tab11/21*0.01403,0.01719,0.02101,0.02561,0.03117,0.03784, & &.04584,.05542,.06685,.08049,.09672,.1160,.1388,.1658,.1977,.2353, & &.2796,.3316,.3925,.4638,.5472,.6444,.7577,.8894,1.042,1.220,1.425, & &1.662,1.936,2.252,2.615,3.032,3.511,4.060,4.688,5.406,6.225,7.159, & &8.223,9.432,10.80,12.36,14.13,16.12,18.38,20.92,23.80,27.03,30.67, & &34.76,39.35,44.49,50.26,56.71,63.93,71.98,80.97,90.98,102.1,114.5, & &128.3,143.6,160.6,179.4,200.2,223.3,248.8,276.9,307.9,342.1,379.8, & &421.3,466.9,517.0,572.0,632.3,698.5,770.9,850.2,937.0,1032./ data tab22/1146.6,1272.0,1408.1,1556.7,1716.9,1890.3,2077.6,2279.6 & &,2496.7,2729.8,2980.0,3247.8,3534.1,3839.8,4164.8,4510.5,4876.9, & &5265.1,5675.2,6107.8,6566.2,7054.7,7575.3,8129.4,8719.2,9346.5, & &10013.,10722.,11474.,12272.,13119.,14017.,14969.,15977.,17044., & &18173.,19367.,20630.,21964.,23373.,24861.,26430.,28086.,29831., & &31671.,33608.,35649.,37796.,40055.,42430.,44927.,47551.,50307., & &53200.,56236.,59422.,62762.,66264.,69934.,73777.,77802.,82015., & &86423.,91034.,95855.,100890.,106160.,111660.,117400.,123400., & &129650.,136170.,142980.,150070.,157460.,165160.,173180.,181530., & &190220.,199260./ data tab3/208670.,218450.,228610.,239180.,250160.,261560.,273400., & &285700.,298450.,311690.,325420.,339650.,354410.,369710.,385560., & &401980.,418980.,436590.,454810.,473670.,493170.,513350.,534220., & &555800.,578090.,601130.,624940.,649530.,674920.,701130.,728190., & &756110.,784920.,814630.,845280.,876880.,909450.,943020.,977610., & &1013250.,1049940.,1087740./ data table1/20*0.0,.3160e-02,.3820e-02,.4600e-02,.5560e-02,.6670e-02, & & .8000e-02,.9580e-02,.1143e-01,.1364e-01,.1623e-01,.1928e-01, & &.2280e-01,.2700e-01,.3190e-01,.3760e-01,.4430e-01,.5200e-01, & &.6090e-01,.7130e-01,.8340e-01,.9720e-01,.1133e+00,.1317e-00, & &.1526e-00,.1780e-00,.2050e-00,.2370e-00,.2740e-00,.3160e-00, & &.3630e-00,.4170e-00,.4790e-00,.5490e-00,.6280e-00,.7180e-00, & &.8190e-00,.9340e-00,.1064e+01,.1209e+01,.1368e+01,.1560e+01, & &.1770e+01,.1990e+01,.2260e+01,.2540e+01,.2880e+01,.3230e+01, & &.3640e+01,.4090e+01,.4590e+01,.5140e+01,.5770e+01,.6450e+01, & &.7220e+01/ data table2/.8050e+01,.8990e+01,.1001e+02,.1112e+02,.1240e+02, & &.1380e+02,.1530e+02,.1700e+02,.1880e+02,.2080e+02,.2310e+02, & &.2550e+02,.2810e+02,.3100e+02,.3420e+02,.3770e+02,.4150e+02, & &.4560e+02,.5010e+02,.5500e+02,.6030e+02,.6620e+02,.7240e+02, & &.7930e+02,.8680e+02,.9500e+02,.1146e+03,.1254e+03,.1361e+03, & &.1486e+03,.1602e+03,.1734e+03,.1873e+03,.2020e+03,.2171e+03, & &.2331e+03,.2502e+03,.2678e+03,.2863e+03,.3057e+03,.3250e+03, & &.3457e+03,.3664e+03,.3882e+03,.4101e+03,.4326e+03,.4584e+03, & &.4885e+03,.5206e+03,.5541e+03,.5898e+03,.6273e+03,.6665e+03, & &.7090e+03/ data table3/.7520e+03,.7980e+03,.8470e+03,.8980e+03,.9520e+03, & &.1008e+04,.1067e+04,.1129e+04,.1194e+04,.1263e+04,.1334e+04, & &.1409e+04,.1488e+04,.1569e+04,.1656e+04,.1745e+04,.1840e+04, & &.1937e+04,.2041e+04,.2147e+04,.2259e+04,.2375e+04,.2497e+04, & &.2624e+04,.2756e+04,.2893e+04,.3036e+04,.3186e+04,.3340e+04, & &.3502e+04,.3670e+04,.3843e+04,.4025e+04,.4213e+04,.4408e+04, & &.4611e+04,.4821e+04,.5035e+04,.5270e+04,.5500e+04,.5740e+04, & &.6000e+04,.6250e+04,.6520e+04,.6810e+04,.7090e+04,.7390e+04, & &.7700e+04,.8020e+04,.8350e+04,.8690e+04,.9040e+04,.9410e+04, & &.9780e+04/ data table4/.1016e+05,.1057e+05,.1098e+05,.1140e+05,.1184e+05, & &.1230e+05,.1275e+05,.1324e+05,.1373e+05,.1423e+05,.1476e+05, & &.1530e+05,.1585e+05,.1642e+05,.1700e+05,.1761e+05,.1822e+05, & &.1886e+05,.1950e+05,.2018e+05,.2087e+05,.2158e+05,.2229e+05, & &.2304e+05,.2381e+05,.2459e+05,.2539e+05,.2621e+05,.2706e+05, & &.2792e+05,.2881e+05,.2971e+05,.3065e+05,.3160e+05,.3257e+05, & &.3357e+05,.3459e+05,.3564e+05,.3669e+05,.3780e+05,.0000e+00/ ! ! spcify constants needed by MFLUX2 ! !GJF: should send through argument list, but these have nonstandard units real,parameter :: cp = 1.00464e7 real,parameter :: g = 980.6 real,parameter :: rgas = 2.87e6 real,parameter :: og = 1./g integer :: ntstep = 0 ! Initialize CCPP error handling variables errmsg = '' errflg = 0 ! #if HWRF==1 real*8 :: gasdev,ran1 !zhang real :: rr !zhang logical,save :: pert_Cd_local !zhang CHARACTER(len=3) :: env_memb,env_pp integer,save :: ens_random_seed_local,env_pp_local !zhang integer :: ensda_physics_pert !zhang real,save :: ens_Cdamp_local !zhang data ens_random_seed_local/0/ data env_pp_local/0/ if ( ens_random_seed_local .eq. 0 ) then CALL nl_get_ensda_physics_pert(1,ensda_physics_pert) ens_random_seed_local=ens_random_seed env_pp_local=ensda_physics_pert pert_Cd_local=.false. ens_Cdamp_local=0.0 ! env_pp=1: do physics perturbations for ensda members, ens_random_seed must be 99 if ( env_pp_local .eq. 1 ) then if ( ens_random_seed .ne. 99 ) then pert_Cd_local=.true. ens_Cdamp_local=ens_Cdamp else ! ens_random_seed=99 do physics perturbation for ensemble forecasts, env_pp must be zero ens_random_seed_local=ens_random_seed pert_Cd_local=pert_Cd ens_Cdamp_local=ens_Cdamp endif else ens_random_seed_local=ens_random_seed pert_Cd_local=pert_Cd ens_Cdamp_local=ens_Cdamp endif print*, "Cd ===", ens_random_seed_local,pert_Cd_local,ens_Cdamp_local,ensda_physics_pert endif #endif ! character*10 routine ! routine = 'mflux2' ! !------------------------------------------------------------------------ ! set water availability constant "ecof" and land mask "land". ! limit minimum wind speed to 100 cm/s !------------------------------------------------------------------------ ! constants for 10 m winds (correction for knots ! cor1 = .120 cor2 = 720. ! KWON : remove the artificial increase of 10m wind speed over 60kts ! which comes from GFDL hurricane model cor1 = 0. cor2 = 0. ! do i = its,ite z10(i) = 1000. z2 (i) = 200. pss(i) = pspc(i) tstar(i) = tstrc(i) if ( lcurr_sf .and. zoc(i) .le. 0.0 ) then ubot = upc(i) - xcur(i) * 100.0 vbot = vpc(i) - ycur(i) * 100.0 ! ubot = upc(i) ! vbot = vpc(i) else ubot = upc(i) vbot = vpc(i) endif uvs1(i)= amax1( SQRT(ubot*ubot + & vbot*vbot), 100.0) if ( iwavecpl .eq. 1 .and. zoc(i) .le. 0.0 ) then ukmax(i) = ( ubot * cos(gamma(i)) - & vbot * sin(gamma(i)) ) & * cos(gamma(i)) vkmax(i) = ( vbot * cos(gamma(i)) - & ubot * sin(gamma(i)) ) & * cos(gamma(i)) else ukmax(i) = ubot vkmax(i) = vbot endif ! ukmax(i) = upc(i) ! vkmax(i) = vpc(i) tkmax(i) = tpc(i) rkmax(i) = rpc(i) enddo do i = its,ite windp(i) = SQRT(ukmax(i)*ukmax(i) + vkmax(i)*vkmax(i)) wind (i) = amax1(windp(i),100.) !! use wind10 previous step wind10p(i) = wind10(i) !! cm/s wind10p(i) = amax1(wind10p(i),100.) !! if (zoc(i) .LT. amask) zoc(i) = -0.0185*0.001*wind10p(i)*wind10p(i)*og if (zoc(i) .GT. 0.0) then ecof(i) = wetc(i) land(i) = 1.0 zot (i) = zoc(i) else ecof(i) = wetc(i) land(i) = 0.0 windmks=wind10p(i)*.01 if ( iwavecpl .eq. 1 ) then call znot_wind10m(windmks,znott,znotm,icoef_sf,errmsg,errflg) !Check if Charnock parameter ratio is received in a proper range. if ( alpha(i) .ge. 0.2 .and. alpha(i) .le. 5. ) then znotm = znotm*alpha(i) endif zoc(i) = -100.*znotm zot(i) = -100* znott else call znot_wind10m(windmks,znott,znotm,icoef_sf,errmsg,errflg) zoc(i) = -100.*znotm zot(i) = -100* znott endif endif !------------------------------------------------------------------------ ! where necessary modify zo values over ocean. !------------------------------------------------------------------------ ! mzoc(i) = zoc(i) !FOR SAVE MOMENTUM Zo tzot(i) = zot(i) !output wang enddo !------------------------------------------------------------------------ ! define constants: ! a and c = constants used in evaluating universal function for ! stable case ! ca = karmen constant ! cm01 = constant part of vertical integral of universal ! function; stable case ( 0.5 < zeta < or = 10.0) ! cm02 = constant part of vertical integral of universal ! function; stable case ( zeta > 10.0) !------------------------------------------------------------------------ en = 2. c = .76 a = 5. ca = .4 cmo1 = .5*a - 1.648 cmo2 = 17.193 + .5*a - 10.*c boycon = .61 rovcp=rgas/cp do i = its,ite theta(i) = tkmax(i)/((pkmax(i)/pspc(i))**rovcp) vrtkx(i) = 1.0 + boycon*rkmax(i) !zkmax(i) = -rgas*tkmax(i)*alog(pkmax(i)/pspc(i))*og zkmax(i) = z1(i) !use precalculated height of first model layer center enddo !------------------------------------------------------------------------ ! get saturation mixing ratios at surface !------------------------------------------------------------------------ do i = its,ite tsg (i) = tstar(i) tab1 (i) = tstar(i) - 153.16 it (i) = IFIX(tab1(i)) tab2 (i) = tab1(i) - FLOAT(it(i)) t1 (i) = tab(min(223,max(1,it(i) + 1))) t2 (i) = table(min(223,max(1,it(i) + 1))) estso(i) = t1(i) + tab2(i)*t2(i) psps1 = (pss(i) - estso(i)) if(psps1 .EQ. 0.0)then psps1 = .1 endif rstso(i) = ep2*estso(i)/psps1 vrts (i) = 1. + boycon*ecof(i)*rstso(i) enddo !------------------------------------------------------------------------ ! check if consideration of virtual temperature changes stability. ! if so, set "dthetav" to near neutral value (1.0e-4). also check ! for very small lapse rates; if ABS(tempa1) <1.0e-4 then ! tempa1=1.0e-4 !------------------------------------------------------------------------ do i = its,ite tempa1(i) = theta(i)*vrtkx(i) - tstar(i)*vrts(i) tempa2(i) = tempa1(i)*(theta(i) - tstar(i)) if (tempa2(i) .LT. 0.) tempa1(i) = 1.0e-4 tab1(i) = ABS(tempa1(i)) if (tab1(i) .LT. 1.0e-4) tempa1(i) = 1.0e-4 !------------------------------------------------------------------------ ! compute bulk richardson number "rib" at each point. if "rib" ! exceeds 95% of critical richardson number "tab1" then "rib = tab1" !------------------------------------------------------------------------ rib (i) = g*zkmax(i)*tempa1(i)/ & (tkmax(i)*vrtkx(i)*wind(i)*wind(i)) tab2(i) = ABS(zoc(i)) tab1(i) = 0.95/(c*(1. - tab2(i)/zkmax(i))) if (rib(i) .GT. tab1(i)) rib(i) = tab1(i) enddo do i = its,ite zeta(i) = ca*rib(i)/0.03 enddo !------------------------------------------------------------------------ ! begin looping through points on line, solving wegsteins iteration ! for zeta at each point, and using hicks functions !------------------------------------------------------------------------ !------------------------------------------------------------------------ ! set initial guess of zeta=non - dimensional height "szeta" for ! stable points !------------------------------------------------------------------------ rca = 1./ca enrca = en*rca ! turn off interfacial layer by zeroing out enrca enrca = 0.0 zog = .0185*og !------------------------------------------------------------------------ ! stable points !------------------------------------------------------------------------ ip = 0 do i = its,ite if (zeta(i) .GE. 0.0) then ip = ip + 1 istb(ip) = i endif enddo if (ip .EQ. 0) go to 170 do i = 1,ip szetam(i) = 1.0e+30 sgzm(i) = 0.0e+00 szeta(i) = zeta(istb(i)) ifz(i) = 1 enddo !------------------------------------------------------------------------ ! begin wegstein iteration for "zeta" at stable points using ! hicks(1976) !------------------------------------------------------------------------ do icnt = 1,icntx do i = 1,ip if (ifz(i) .EQ. 0) go to 80 zal1g = ALOG(szeta(i)) if (szeta(i) .LE. 0.5) then fmz1(i) = (zal1g + a*szeta(i))*rca else if (szeta(i) .GT. 0.5 .AND. szeta(i) .LE. 10.) then rzeta1 = 1./szeta(i) fmz1(i) = (8.*zal1g + 4.25*rzeta1 - & 0.5*rzeta1*rzeta1 + cmo1)*rca else if (szeta(i) .GT. 10.) then fmz1(i) = (c*szeta(i) + cmo2)*rca endif szetao = ABS(zoc(istb(i)))/zkmax(istb(i))*szeta(i) zal2g = ALOG(szetao) if (szetao .LE. 0.5) then fmzo1(i) = (zal2g + a*szetao)*rca sfzo (i) = 1. + a*szetao else if (szetao .GT. 0.5 .AND. szetao .LE. 10.) then rzeta2 = 1./szetao fmzo1(i) = (8.*zal2g + 4.25*rzeta2 - & 0.5*rzeta2*rzeta2 + cmo1)*rca sfzo (i) = 8.0 - 4.25*rzeta2 + rzeta2*rzeta2 else if (szetao .GT. 10.) then fmzo1(i) = (c*szetao + cmo2)*rca sfzo (i) = c*szetao endif ! compute heat & moisture parts of zot.. for calculation of sfh szetho = ABS(zot(istb(i)))/zkmax(istb(i))*szeta(i) zal2gh = ALOG(szetho) if (szetho .LE. 0.5) then fhzo1(i) = (zal2gh + a*szetho)*rca sfzo (i) = 1. + a*szetho else if (szetho .GT. 0.5 .AND. szetho .LE. 10.) then rzeta2 = 1./szetho fhzo1(i) = (8.*zal2gh + 4.25*rzeta2 - & 0.5*rzeta2*rzeta2 + cmo1)*rca sfzo (i) = 8.0 - 4.25*rzeta2 + rzeta2*rzeta2 else if (szetho .GT. 10.) then fhzo1(i) = (c*szetho + cmo2)*rca sfzo (i) = c*szetho endif !------------------------------------------------------------------------ ! compute universal function at 10 meters for diagnostic purposes !------------------------------------------------------------------------ szet10 = ABS(z10(istb(i)))/zkmax(istb(i))*szeta(i) zal10 = ALOG(szet10) if (szet10 .LE. 0.5) then fmz10(i) = (zal10 + a*szet10)*rca else if (szet10 .GT. 0.5 .AND. szet10 .LE. 10.) then rzeta2 = 1./szet10 fmz10(i) = (8.*zal10 + 4.25*rzeta2 - & 0.5*rzeta2*rzeta2 + cmo1)*rca else if (szet10 .GT. 10.) then fmz10(i) = (c*szet10 + cmo2)*rca endif sf10(i) = fmz10(i) - fmzo1(i) ! compute 2m values for diagnostics in HWRF szet2 = ABS(z2 (istb(i)))/zkmax(istb(i))*szeta(i) zal2 = ALOG(szet2 ) if (szet2 .LE. 0.5) then fmz2 (i) = (zal2 + a*szet2 )*rca else if (szet2 .GT. 0.5 .AND. szet2 .LE. 2.) then rzeta2 = 1./szet2 fmz2 (i) = (8.*zal2 + 4.25*rzeta2 - & 0.5*rzeta2*rzeta2 + cmo1)*rca else if (szet2 .GT. 2.) then fmz2 (i) = (c*szet2 + cmo2)*rca endif sf2 (i) = fmz2 (i) - fmzo1(i) sfm(i) = fmz1(i) - fmzo1(i) sfh(i) = fmz1(i) - fhzo1(i) sgz = ca*rib(istb(i))*sfm(i)*sfm(i)/ & (sfh(i) + enrca*sfzo(i)) fmz = (sgz - szeta(i))/szeta(i) fmzo = ABS(fmz) if (fmzo .GE. 5.0e-5) then sq = (sgz - sgzm(i))/(szeta(i) - szetam(i)) if(sq .EQ. 1) then write(errmsg,'(*(a))') 'NCO ERROR DIVIDE BY ZERO IN gfdl_sfc_layer.F90/MFLUX2 (STABLE CASE)'// & 'sq is 1 ',fmzo,sgz,sgzm(i),szeta(i),szetam(i) errflg = 1 return endif szetam(i) = szeta(i) szeta (i) = (sgz - szeta(i)*sq)/(1.0 - sq) sgzm (i) = sgz else ifz(i) = 0 endif 80 continue enddo enddo do i = 1,ip if (ifz(i) .GE. 1) go to 110 enddo go to 130 110 continue write(errmsg,'(*(a))') 'NON-CONVERGENCE FOR STABLE ZETA IN gfdl_sfc_layer.F90/MFLUX2' errflg = 1 return ! call MPI_CLOSE(1,routine) !------------------------------------------------------------------------ ! update "zo" for ocean points. "zo"cannot be updated within the ! wegsteins iteration as the scheme (for the near neutral case) ! can become unstable !------------------------------------------------------------------------ 130 continue do i = 1,ip szo = zoc(istb(i)) if (szo .LT. 0.0) then wndm=wind(istb(i))*0.01 if(wndm.lt.15.0) then ckg=0.0185*og else ckg=(sfenth*(4*0.000308*wndm) + (1.-sfenth)*0.0185 )*og endif szo = - ckg*wind(istb(i))*wind(istb(i))/ & (sfm(i)*sfm(i)) cons_p000001 = .000001 cons_7 = 7. vis = 1.4E-1 ustar = sqrt( -szo / zog) restar = -ustar * szo / vis restar = max(restar,cons_p000001) ! Rat taken from Zeng, Zhao and Dickinson 1997 rat = 2.67 * restar ** .25 - 2.57 rat = min(rat ,cons_7) !constant rat=0. zot(istb(i)) = szo * exp(-rat) else zot(istb(i)) = zoc(istb(i)) endif ! in hwrf thermal znot is loaded back into the zoc array for next step zoc(istb(i)) = szo enddo do i = 1,ip xxfm(istb(i)) = sfm(i) xxfh(istb(i)) = sfh(i) xxfh2(istb(i)) = sf2 (i) xxsh(istb(i)) = sfzo(i) enddo !------------------------------------------------------------------------ ! obtain wind at 10 meters for diagnostic purposes !------------------------------------------------------------------------ do i = 1,ip wind10(istb(i)) = sf10(i)*uvs1(istb(i))/sfm(i) wind10(istb(i)) = wind10(istb(i)) * 1.944 if(wind10(istb(i)) .GT. 6000.0) then wind10(istb(i))=wind10(istb(i))+wind10(istb(i))*cor1 & - cor2 endif ! the above correction done by GFDL in centi-kts!!!-change back wind10(istb(i)) = wind10(istb(i)) / 1.944 enddo !------------------------------------------------------------------------ ! unstable points !------------------------------------------------------------------------ 170 continue iq = 0 do i = its,ite if (zeta(i) .LT. 0.0) then iq = iq + 1 iutb(iq) = i endif enddo if (iq .EQ. 0) go to 290 do i = 1,iq uzeta (i) = zeta(iutb(i)) ifz (i) = 1 uzetam(i) = 1.0e+30 ugzm (i) = 0.0e+00 enddo !------------------------------------------------------------------------ ! begin wegstein iteration for "zeta" at unstable points using ! hicks functions !------------------------------------------------------------------------ do icnt = 1,icntx do i = 1,iq if (ifz(i) .EQ. 0) go to 200 ugzzo = ALOG(zkmax(iutb(i))/ABS(zot(iutb(i)))) uzetao = ABS(zot(iutb(i)))/zkmax(iutb(i))*uzeta(i) ux11 = 1. - 16.*uzeta(i) ux12 = 1. - 16.*uzetao y = SQRT(ux11) yo = SQRT(ux12) ufzo(i) = 1./yo ux13 = (1. + y)/(1. + yo) ux21 = ALOG(ux13) ufh(i) = (ugzzo - 2.*ux21)*rca ! recompute scalers for ufm in terms of mom znot... zoc ugzzo = ALOG(zkmax(iutb(i))/ABS(zoc(iutb(i)))) uzetao = ABS(zoc(iutb(i)))/zkmax(iutb(i))*uzeta(i) ux11 = 1. - 16.*uzeta(i) ux12 = 1. - 16.*uzetao y = SQRT(ux11) yo = SQRT(ux12) ux13 = (1. + y)/(1. + yo) ux21 = ALOG(ux13) ! ufzo(i) = 1./yo x = SQRT(y) xo = SQRT(yo) xnum = (x**2 + 1.)*((x + 1.)**2) xden = (xo**2 + 1.)*((xo + 1.)**2) xtan = ATAN(x) - ATAN(xo) ux3 = ALOG(xnum/xden) ufm(i) = (ugzzo - ux3 + 2.*xtan)*rca !------------------------------------------------------------------------ ! obtain ten meter winds for diagnostic purposes !------------------------------------------------------------------------ ugz10 = ALOG(z10(iutb(i))/ABS(zoc(iutb(i)))) uzet1o = ABS(z10(iutb(i)))/zkmax(iutb(i))*uzeta(i) uzetao = ABS(zoc(iutb(i)))/zkmax(iutb(i))*uzeta(i) ux11 = 1. - 16.*uzet1o ux12 = 1. - 16.*uzetao y = SQRT(ux11) y10 = SQRT(ux12) ux13 = (1. + y)/(1. + y10) ux21 = ALOG(ux13) x = SQRT(y) x10 = SQRT(y10) xnum = (x**2 + 1.)*((x + 1.)**2) xden = (x10**2 + 1.)*((x10 + 1.)**2) xtan = ATAN(x) - ATAN(x10) ux3 = ALOG(xnum/xden) uf10(i) = (ugz10 - ux3 + 2.*xtan)*rca ! obtain 2m values for diagnostics... ugz2 = ALOG(z2 (iutb(i))/ABS(zoc(iutb(i)))) uzet1o = ABS(z2 (iutb(i)))/zkmax(iutb(i))*uzeta(i) uzetao = ABS(zoc(iutb(i)))/zkmax(iutb(i))*uzeta(i) ux11 = 1. - 16.*uzet1o ux12 = 1. - 16.*uzetao y = SQRT(ux11) yo = SQRT(ux12) ux13 = (1. + y)/(1. + yo) ux21 = ALOG(ux13) uf2 (i) = (ugzzo - 2.*ux21)*rca ugz = ca*rib(iutb(i))*ufm(i)*ufm(i)/(ufh(i) + enrca*ufzo(i)) ux1 = (ugz - uzeta(i))/uzeta(i) ux2 = ABS(ux1) if (ux2 .GE. 5.0e-5) then uq = (ugz - ugzm(i))/(uzeta(i) - uzetam(i)) uzetam(i) = uzeta(i) if(uq .EQ. 1) then write(errmsg,'(*(a))') 'NCO ERROR DIVIDE BY ZERO IN gfdl_sfc_layer.F90/MFLUX2 (UNSTABLE CASE)'// & 'uq is 1 ',ux2,ugz,ugzm(i),uzeta(i),uzetam(i) errflg = 1 return endif uzeta (i) = (ugz - uzeta(i)*uq)/(1.0 - uq) ugzm (i) = ugz else ifz(i) = 0 endif 200 continue enddo enddo do i = 1,iq if (ifz(i) .GE. 1) go to 230 enddo go to 250 230 continue write(errmsg,'(*(a))') 'NON-CONVERGENCE FOR UNSTABLE ZETA IN ROW'// & 'uq is 1 ',ux2,ugz,ugzm(i),uzeta(i),uzetam(i) errflg = 1 return ! call MPI_CLOSE(1,routine) !------------------------------------------------------------------------ ! gather unstable values !------------------------------------------------------------------------ 250 continue !------------------------------------------------------------------------ ! update "zo" for ocean points. zo cannot be updated within the ! wegsteins iteration as the scheme (for the near neutral case) ! can become unstable. !------------------------------------------------------------------------ do i = 1,iq uzo = zoc(iutb(i)) if (zoc(iutb(i)) .LT. 0.0) then wndm=wind(iutb(i))*0.01 if(wndm.lt.15.0) then ckg=0.0185*og else ckg=(4*0.000308*wndm)*og ckg=(sfenth*(4*0.000308*wndm) + (1.-sfenth)*0.0185 )*og endif uzo =-ckg*wind(iutb(i))*wind(iutb(i))/(ufm(i)*ufm(i)) cons_p000001 = .000001 cons_7 = 7. vis = 1.4E-1 ustar = sqrt( -uzo / zog) restar = -ustar * uzo / vis restar = max(restar,cons_p000001) ! Rat taken from Zeng, Zhao and Dickinson 1997 rat = 2.67 * restar ** .25 - 2.57 rat = min(rat ,cons_7) !constant rat=0.0 zot(iutb(i)) = uzo * exp(-rat) else zot(iutb(i)) = zoc(iutb(i)) endif ! in hwrf thermal znot is loaded back into the zoc array for next step zoc(iutb(i)) = uzo enddo !------------------------------------------------------------------------ ! obtain wind at ten meters for diagnostic purposes !------------------------------------------------------------------------ do i = 1,iq wind10(iutb(i)) = uf10(i)*uvs1(iutb(i))/ufm(i) wind10(iutb(i)) = wind10(iutb(i)) * 1.944 if(wind10(iutb(i)) .GT. 6000.0) then wind10(iutb(i))=wind10(iutb(i))+wind10(iutb(i))*cor1 & - cor2 endif ! the above correction done by GFDL in centi-kts!!!-change back wind10(iutb(i)) = wind10(iutb(i)) / 1.944 enddo do i = 1,iq xxfm(iutb(i)) = ufm(i) xxfh(iutb(i)) = ufh(i) xxfh2(iutb(i)) = uf2 (i) xxsh(iutb(i)) = ufzo(i) enddo 290 continue do i = its,ite ucom(i) = ukmax(i) vcom(i) = vkmax(i) if (windp(i) .EQ. 0.0) then windp(i) = 100.0 ucom (i) = 100.0/SQRT(2.0) vcom (i) = 100.0/SQRT(2.0) endif rho(i) = pss(i)/(rgas*(tsg(i) + enrca*(theta(i) - & tsg(i))*xxsh(i)/(xxfh(i) + enrca*xxsh(i)))) bq1(i) = wind(i)*rho(i)/(xxfm(i)*(xxfh(i) + enrca*xxsh(i))) enddo ! do land sfc temperature prediction if ntsflg=1 ! ntsflg = 1 ! gopal's doing if (ntsflg .EQ. 0) go to 370 alll = 600. xks = 0.01 hcap = .5/2.39e-8 pith = SQRT(4.*ATAN(1.0)) alfus = alll/2.39e-8 teps = 0.1 ! slwdc... in units of cal/min ???? ! slwa... in units of ergs/sec/cm*2 ! 1 erg=2.39e-8 cal !------------------------------------------------------------------------ ! pack land and sea ice points !------------------------------------------------------------------------ ip = 0 do i = its,ite if (land(i) .EQ. 1) then ip = ip + 1 indx (ip) = i ! slwa is defined as positive down.... slwa (ip) = slwdc(i)/(2.39e-8*60.) tss (ip) = tstar(i) thetap (ip) = theta(i) rkmaxp (ip) = rkmax(i) aap (ip) = 5.673e-5 pssp (ip) = pss(i) ecofp (ip) = ecof(i) estsop (ip) = estso(i) rstsop (ip) = rstso(i) bq1p (ip) = bq1(i) bq1p (ip) = amax1(bq1p(ip),0.1e-3) delsrad(ip) = dt *pith/(hcap*SQRT(3600.*24.*xks)) endif enddo !------------------------------------------------------------------------ ! initialize variables for first pass of iteration !------------------------------------------------------------------------ do i = 1,ip ifz (i) = 1 tsm (i) = tss(i) rdiff(i) = amin1(0.0,(rkmaxp(i) - rstsop(i))) 300 format(2X, ' SURFACE EQUILIBRIUM CALCULATION ') foftm(i) = tss(i) + delsrad(i)*(slwa(i) - aap(i)*tsm(i)**4 - & cp*bq1p(i)*(tsm(i) - thetap(i)) + ecofp(i)*alfus*bq1p(i)* & rdiff(i)) tsp(i) = foftm(i) enddo !------------------------------------------------------------------------ ! do iteration to determine "tstar" at new time level !------------------------------------------------------------------------ do icnt = 1,icntx do i = 1,ip if (ifz(i) .EQ. 0) go to 330 tab1 (i) = tsp(i) - 153.16 it (i) = IFIX(tab1(i)) tab2 (i) = tab1(i) - FLOAT(it(i)) t1 (i) = tab(min(223,max(1,it(i) + 1))) t2 (i) = table(min(223,max(1,it(i) + 1))) estsop(i) = t1(i) + tab2(i)*t2(i) psps2 = (pssp(i) - estsop(i)) if(psps2 .EQ. 0.0)then psps2 = .1 endif rstsop(i) = ep2*estsop(i)/psps2 rdiff (i) = amin1(0.0,(rkmaxp(i) - rstsop(i))) foft(i) = tss(i) + delsrad(i)*(slwa(i) - aap(i)*tsp(i)**4 - & cp*bq1p(i)*(tsp(i) - thetap(i)) + ecofp(i)*alfus*bq1p(i)* & rdiff(i)) frac(i) = ABS((foft(i) - tsp(i))/tsp(i)) !------------------------------------------------------------------------ ! check for convergence of all points use wegstein iteration !------------------------------------------------------------------------ if (frac(i) .GE. teps) then qf (i) = (foft(i) - foftm(i))/(tsp(i) - tsm(i)) tsm (i) = tsp(i) tsp (i) = (foft(i) - tsp(i)*qf(i))/(1. - qf(i)) foftm(i) = foft(i) else ifz(i) = 0 endif 330 continue enddo enddo !------------------------------------------------------------------------ ! check for convergence of "t star" prediction !------------------------------------------------------------------------ do i = 1,ip if (ifz(i) .EQ. 1) then write(errmsg,'(*(a))') 'NON-CONVERGENCE OF T* PREDICTED (T*,I) = ', & tsp(i), i errflg = 1 return ! call MPI_CLOSE(1,routine) endif enddo do i = 1,ip ii = indx(i) tstrc(ii) = tsp (i) enddo !------------------------------------------------------------------------ ! compute fluxes and momentum drag coef !------------------------------------------------------------------------ 370 continue do i = its,ite !!! if ( iwavecpl .eq. 1 .and. zoc(i) .le. 0.0 ) then windmks = wind10(i) * 0.01 call znot_wind10m(windmks,znott,znotm,icoef_sf,errmsg,errflg) !Check if Charnock parameter ratio is received in a proper range. if ( alpha(i) .ge. 0.2 .and. alpha(i) .le. 5. ) then znotm = znotm*alpha(i) endif zoc(i) = -100.*znotm zot(i) = -100* znott endif !!!! fxh(i) = bq1(i)*(theta(i) - tsg(i)) fxe(i) = ecof(i)*bq1(i)*(rkmax(i) - rstso(i)) if (fxe(i) .GT. 0.0) fxe(i) = 0.0 fxmx(i) = rho(i)/(xxfm(i)*xxfm(i))*wind(i)*wind(i)*ucom(i)/ & windp(i) fxmy(i) = rho(i)/(xxfm(i)*xxfm(i))*wind(i)*wind(i)*vcom(i)/ & windp(i) cdm(i) = 1./(xxfm(i)*xxfm(i)) #if HWRF==1 ! randomly perturb the Cd !zzz if( pert_Cd_local .and. ens_random_seed_local .gt. 0 ) then if( pert_Cd_local ) then ens_random_seed_local=ran1(-ens_random_seed_local)*1000 rr=2.0*ens_Cdamp_local*ran1(-ens_random_seed_local)-ens_Cdamp_local cdm(i) = cdm(i) *(1.0+rr) endif #endif enddo ntstep = ntstep + 1 return end subroutine MFLUX2 end module gfdl_sfc_layer