99 subroutine gfdl_sfc_layer_run (im, nsoil, km, xlat, xlon, flag_iter, lsm, &
100 lsm_noah, lsm_noahmp, lsm_ruc, icoef_sf, cplwav, karman, &
101 cplwav2atm, lcurr_sf, pert_Cd, ntsflg, sfenth, z1, shdmax, ivegsrc, &
102 vegtype, sigmaf, dt, wet, dry, icy, isltyp, rd, grav, ep1, ep2, smois, &
103 psfc, prsl1, q1, t1, u1, v1, wspd, u10, v10, gsw, glw, tsurf_wat, &
104 tsurf_lnd, tsurf_ice, tskin_wat, tskin_lnd, tskin_ice, ustar_wat, &
105 ustar_lnd, ustar_ice, znt_wat, znt_lnd, znt_ice, cdm_wat, cdm_lnd, &
106 cdm_ice, stress_wat, stress_lnd, stress_ice, rib_wat, rib_lnd, rib_ice, &
107 fm_wat, fm_lnd, fm_ice, fh_wat, fh_lnd, fh_ice, fh2_wat, fh2_lnd, &
108 fh2_ice, ch_wat, ch_lnd, ch_ice, fm10_wat, fm10_lnd, fm10_ice, qss_wat, &
109 qss_lnd, qss_ice, errmsg, errflg)
117 use noahmp_tables,
only: maxsmc_noahmp => smcmax_table, drysmc_noahmp => smcdry_table
122 integer,
intent(in) :: im, nsoil, km, ivegsrc
123 integer,
intent(in) :: lsm,
lsm_noah, lsm_noahmp, &
125 logical,
intent(in) :: cplwav, cplwav2atm
126 logical,
intent(in) :: lcurr_sf
127 logical,
intent(in) :: pert_cd
128 logical,
dimension(:),
intent(in) :: flag_iter, wet, dry, icy
129 integer,
dimension(:),
intent(in) :: isltyp, vegtype
130 real(kind=kind_phys),
intent(in) :: dt, sfenth
131 real(kind=kind_phys),
intent(in) :: rd,grav,ep1,ep2
132 real(kind=kind_phys),
dimension(:,:),
intent(in) :: smois
133 real(kind=kind_phys),
dimension(:),
intent(in) :: psfc, prsl1, &
134 q1, t1, u1, v1, wspd, u10, v10, gsw, glw, z1, shdmax, sigmaf, xlat, &
135 xlon, tsurf_wat, tsurf_lnd, tsurf_ice
137 real(kind=kind_phys),
intent(inout),
dimension(:) :: tskin_wat, &
138 tskin_lnd, tskin_ice, ustar_wat, ustar_lnd, ustar_ice, &
139 znt_wat, znt_lnd, znt_ice, cdm_wat, cdm_lnd, cdm_ice, &
140 stress_wat, stress_lnd, stress_ice, rib_wat, rib_lnd, rib_ice, &
141 fm_wat, fm_lnd, fm_ice, fh_wat, fh_lnd, fh_ice, fh2_wat, fh2_lnd, &
142 fh2_ice, ch_wat, ch_lnd, ch_ice, fm10_wat, fm10_lnd, fm10_ice, &
143 qss_wat, qss_lnd, qss_ice
145 character(len=*),
intent(out) :: errmsg
146 integer,
intent(out) :: errflg
150 integer :: i, its, ite, ims, ime
152 logical :: ch_bound_excursion
155 real (kind=kind_phys),
intent(in) :: karman
156 real (kind=kind_phys),
parameter :: log01=log(0.01), log05=log(0.05), &
160 integer :: iwavecpl, ens_random_seed, issflx
161 logical :: diag_wind10m, diag_qss
162 real(kind=kind_phys) :: ens_cdamp
164 real(kind=kind_phys),
dimension(im) :: wetc, pspc, pkmax, tstrc, upc, &
165 vpc, mznt, slwdc, wind10, qfx, qgh, zkmax, z1_cm, z0max, ztmax
166 real(kind=kind_phys),
dimension(im) :: u10_lnd, u10_ocn, u10_ice, &
167 v10_lnd, v10_ocn, v10_ice
174 real(kind=kind_phys),
dimension(im) :: charn, msang, scurx, scury
176 real(kind=kind_phys),
dimension(im) :: fxh, fxe, fxmx, fxmy, xxfh, &
178 real(kind=kind_phys),
dimension(1:30) :: maxsmc, drysmc
179 real(kind=kind_phys) :: smcmax, smcdry, zhalf, cd10, &
180 esat, fm_lnd_old, fh_lnd_old, tem1, tem2, czilc, cd_low_limit, &
181 cd_high_limit, ch_low_limit, ch_high_limit, fh2_fh_ratio
187 else if (lsm == lsm_noahmp)
then
188 maxsmc = maxsmc_noahmp
189 drysmc = drysmc_noahmp
198 data maxsmc/0.339, 0.421, 0.434, 0.476, 0.476, 0.439, &
199 0.404, 0.464, 0.465, 0.406, 0.468, 0.468, &
200 0.439, 1.000, 0.200, 0.421, 0.000, 0.000, &
201 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, &
202 0.000, 0.000, 0.000, 0.000, 0.000, 0.000/
203 data drysmc/0.010, 0.028, 0.047, 0.084, 0.084, 0.066, &
204 0.067, 0.120, 0.103, 0.100, 0.126, 0.138, &
205 0.066, 0.000, 0.006, 0.028, 0.000, 0.000, &
206 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, &
207 0.000, 0.000, 0.000, 0.000, 0.000, 0.000/
224 diag_wind10m = .false.
240 if (flag_iter(i))
then
243 pspc(i) = psfc(i)*10.
244 pkmax(i) = prsl1(i)*10.
250 wind10(i)=sqrt(u10(i)*u10(i)+v10(i)*v10(i))
262 z1_cm(i) = 100.0*z1(i)
265 cd_low_limit = 1.0e-5/zkmax(i)
269 ch_low_limit = cd_low_limit
270 ch_high_limit = min(0.1,0.05/wspd(i))
274 slwdc(i)=gsw(i)+glw(i)
275 slwdc(i)=0.239*60.*slwdc(i)*1.e-4
297 smcdry=drysmc(isltyp(i))
298 smcmax=maxsmc(isltyp(i))
299 wetc(i)=(smois(i,1)-smcdry)/(smcmax-smcdry)
300 wetc(i)=amin1(1.,amax1(wetc(i),0.))
303 tstrc(i) = 0.5*(tskin_lnd(i) + tsurf_lnd(i))
313 z0max(i) = max(1.0e-6, min(0.01 * znt_lnd(i), zkmax(i)))
315 tem1 = 1.0 - shdmax(i)
319 if( ivegsrc == 1 )
then
320 if (vegtype(i) == 10)
then
321 z0max(i) = exp( tem2*log01 + tem1*log07 )
322 elseif (vegtype(i) == 6)
then
323 z0max(i) = exp( tem2*log01 + tem1*log05 )
324 elseif (vegtype(i) == 7)
then
327 elseif (vegtype(i) == 16)
then
331 z0max(i) = exp( tem2*log01 + tem1*log(z0max(i)) )
333 elseif (ivegsrc == 2 )
then
334 if (vegtype(i) == 7)
then
335 z0max(i) = exp( tem2*log01 + tem1*log07 )
336 elseif (vegtype(i) == 8)
then
337 z0max(i) = exp( tem2*log01 + tem1*log05 )
338 elseif (vegtype(i) == 9)
then
341 elseif (vegtype(i) == 11)
then
345 z0max(i) = exp( tem2*log01 + tem1*log(z0max(i)) )
349 z0max(i) = max(z0max(i), 1.0e-6)
354 tem1 = 1.0 - sigmaf(i)
355 ztmax(i) = z0max(i)*exp( - tem1*tem1 &
356 & * czilc*karman*sqrt(ustar_lnd(i)*(0.01/1.5e-05)))
357 ztmax(i) = max(ztmax(i), 1.0e-6)
360 if (wind10(i) <= 1.0e-10 .or. wind10(i) > 150.0)
then
361 wind10(i)=wspd(i)*alog(10.0/z0max(i))/alog(z1(i)/z0max(i))
363 wind10(i)=wind10(i)*100.0
365 ztmax(i) = ztmax(i)*100.0
366 z0max(i) = z0max(i)*100.0
368 call mflux2 (fxh(i), fxe(i), fxmx(i), fxmy(i), cdm_lnd(i:i), rib_lnd(i:i), &
369 xxfh(i), ztmax(i), z0max(i), tstrc(i), &
370 pspc(i), pkmax(i), wetc(i), slwdc(i), z1_cm(i), icoef_sf, iwavecpl, lcurr_sf, charn(i), msang(i), &
371 scurx(i), scury(i), pert_cd, ens_random_seed, ens_cdamp, upc(i), vpc(i), t1(i:i), q1(i:i), &
372 dt, wind10(i), xxfh2(i), ntsflg, sfenth, tzot(i), ep2, errmsg, &
374 if (errflg /= 0)
return
378 tskin_lnd(i) = tstrc(i)
381 if (diag_wind10m)
then
382 u10_lnd(i) = u1(i)*(0.01*wind10(i)/wspd(i))
383 v10_lnd(i) = v1(i)*(0.01*wind10(i)/wspd(i))
391 cdm_lnd(i) = max(cdm_lnd(i), cd_low_limit)
392 cdm_lnd(i) = min(cdm_lnd(i), cd_high_limit)
393 fm_lnd(i) = karman/sqrt(cdm_lnd(i))
396 fh_lnd(i) = karman*xxfh(i)
399 ch_lnd(i) = karman*karman/(fm_lnd(i) * fh_lnd(i))
402 ch_bound_excursion = .false.
403 if (ch_lnd(i) < ch_low_limit)
then
404 ch_bound_excursion = .true.
405 ch_lnd(i) = ch_low_limit
406 else if (ch_lnd(i) > ch_high_limit)
then
407 ch_bound_excursion = .true.
408 ch_lnd(i) = ch_high_limit
411 fh2_lnd(i) = karman*xxfh2(i)
413 if (ch_bound_excursion)
then
414 fh2_fh_ratio = min(xxfh2(i)/xxfh(i), 1.0)
415 fh_lnd(i) = karman*karman/(fm_lnd(i)*ch_lnd(i))
416 fh2_lnd(i) = fh2_fh_ratio*fh_lnd(i)
424 ustar_lnd(i) = 0.01*sqrt(cdm_lnd(i)* &
425 (upc(i)*upc(i) + vpc(i)*vpc(i)))
427 ustar_lnd(i) = amax1(ustar_lnd(i),0.001)
429 stress_lnd(i) = cdm_lnd(i)*wspd(i)*wspd(i)
434 if ( wind10(i) .ge. 0.1 )
then
435 cd10=cdm_lnd(i)* (wspd(i)/(0.01*wind10(i)) )**2
440 fm10_lnd(i) = karman/sqrt(cd10)
454 esat = fpvs(tskin_lnd(i))
455 qss_lnd(i) = ep2*
esat/(psfc(i)-
esat)
466 smcdry=drysmc(isltyp(i))
467 smcmax=maxsmc(isltyp(i))
468 wetc(i)=(smois(i,1)-smcdry)/(smcmax-smcdry)
469 wetc(i)=amin1(1.,amax1(wetc(i),0.))
473 tstrc(i) = 0.5*(tskin_ice(i) + tsurf_ice(i))
484 z0max(i) = max(1.0e-6, min(0.01 * znt_ice(i), zkmax(i)))
486 tem1 = 1.0 - shdmax(i)
490 if( ivegsrc == 1 )
then
491 z0max(i) = exp( tem2*log01 + tem1*log(z0max(i)) )
492 elseif (ivegsrc == 2 )
then
493 z0max(i) = exp( tem2*log01 + tem1*log(z0max(i)) )
496 z0max(i) = max(z0max(i), 1.0e-6)
502 tem1 = 1.0 - sigmaf(i)
503 ztmax(i) = z0max(i)*exp( - tem1*tem1 &
504 & * czilc*karman*sqrt(ustar_ice(i)*(0.01/1.5e-05)))
505 ztmax(i) = max(ztmax(i), 1.0e-6)
509 if (wind10(i) <= 1.0e-10 .or. wind10(i) > 150.0)
then
510 wind10(i)=wspd(i)*alog(10.0/z0max(i))/alog(z1(i)/z0max(i))
512 wind10(i)=wind10(i)*100.0
514 ztmax(i) = ztmax(i)*100.0
515 z0max(i) = z0max(i)*100.0
517 call mflux2 (fxh(i), fxe(i), fxmx(i), fxmy(i), cdm_ice(i:i), rib_ice(i:i), &
518 xxfh(i), ztmax(i), z0max(i), tstrc(i), &
519 pspc(i), pkmax(i), wetc(i), slwdc(i), z1_cm(i), icoef_sf, iwavecpl, lcurr_sf, charn(i), msang(i), &
520 scurx(i), scury(i), pert_cd, ens_random_seed, ens_cdamp, upc(i), vpc(i), t1(i:i), q1(i:i), &
521 dt, wind10(i), xxfh2(i), ntsflg, sfenth, tzot(i), ep2, errmsg, &
523 if (errflg /= 0)
return
527 tskin_ice(i) = tstrc(i)
530 if (diag_wind10m)
then
531 u10_ice(i) = u1(i)*(0.01*wind10(i)/wspd(i))
532 v10_ice(i) = v1(i)*(0.01*wind10(i)/wspd(i))
540 cdm_ice(i) = max(cdm_ice(i), cd_low_limit)
541 cdm_ice(i) = min(cdm_ice(i), cd_high_limit)
542 fm_ice(i) = karman/sqrt(cdm_ice(i))
545 fh_ice(i) = karman*xxfh(i)
548 ch_ice(i) = karman*karman/(fm_ice(i) * fh_ice(i))
551 ch_bound_excursion = .false.
552 if (ch_ice(i) < ch_low_limit)
then
553 ch_bound_excursion = .true.
554 ch_ice(i) = ch_low_limit
555 else if (ch_ice(i) > ch_high_limit)
then
556 ch_bound_excursion = .true.
557 ch_ice(i) = ch_high_limit
560 fh2_ice(i) = karman*xxfh2(i)
562 if (ch_bound_excursion)
then
563 fh2_fh_ratio = min(xxfh2(i)/xxfh(i), 1.0)
564 fh_ice(i) = karman*karman/(fm_ice(i)*ch_ice(i))
565 fh2_ice(i) = fh2_fh_ratio*fh_ice(i)
572 ustar_ice(i) = 0.01*sqrt(cdm_ice(i)* &
573 (upc(i)*upc(i) + vpc(i)*vpc(i)))
575 ustar_ice(i) = amax1(ustar_ice(i),0.001)
577 stress_ice(i) = cdm_ice(i)*wspd(i)*wspd(i)
582 if ( wind10(i) .ge. 0.1 )
then
583 cd10=cdm_ice(i)* (wspd(i)/(0.01*wind10(i)) )**2
588 fm10_ice(i) = karman/sqrt(cd10)
595 esat = fpvs(tskin_ice(i))
596 qss_ice(i) = ep2*
esat/(psfc(i)-
esat)
608 tstrc(i) = 0.5*(tskin_wat(i) + tsurf_wat(i))
613 znt_wat(i)=min(2.85e-1,max(znt_wat(i),1.27e-5))
616 if (wind10(i) <= 1.0e-10 .or. wind10(i) > 150.0)
then
617 wind10(i)=wspd(i)*alog(10.0/(0.01*znt_wat(i)))/alog(z1(i)/(0.01*znt_wat(i)))
619 wind10(i)=wind10(i)*100.0
622 znt_wat(i) = -znt_wat(i)
624 call mflux2 (fxh(i), fxe(i), fxmx(i), fxmy(i), cdm_wat(i:i), rib_wat(i:i), &
625 xxfh(i), znt_wat(i:i), mznt(i), tstrc(i), &
626 pspc(i), pkmax(i), wetc(i), slwdc(i), z1_cm(i), icoef_sf, iwavecpl, lcurr_sf, charn(i), msang(i), &
627 scurx(i), scury(i), pert_cd, ens_random_seed, ens_cdamp, upc(i), vpc(i), t1(i:i), q1(i:i), &
628 dt, wind10(i), xxfh2(i), ntsflg, sfenth, tzot(i), ep2, errmsg, &
630 if (errflg /= 0)
return
634 tskin_wat(i) = tstrc(i)
637 znt_wat(i)= abs(znt_wat(i))
638 mznt(i)= abs(mznt(i))
641 znt_wat(i)=min(2.85e-1,max(znt_wat(i),1.27e-5))
643 if (diag_wind10m)
then
644 u10_ocn(i) = u1(i)*(0.01*wind10(i)/wspd(i))
645 v10_ocn(i) = v1(i)*(0.01*wind10(i)/wspd(i))
653 cdm_wat(i) = max(cdm_wat(i), cd_low_limit)
654 cdm_wat(i) = min(cdm_wat(i), cd_high_limit)
655 fm_wat(i) = karman/sqrt(cdm_wat(i))
658 fh_wat(i) = karman*xxfh(i)
661 ch_wat(i) = karman*karman/(fm_wat(i) * fh_wat(i))
664 ch_bound_excursion = .false.
665 if (ch_wat(i) < ch_low_limit)
then
666 ch_bound_excursion = .true.
667 ch_wat(i) = ch_low_limit
668 else if (ch_wat(i) > ch_high_limit)
then
669 ch_bound_excursion = .true.
670 ch_wat(i) = ch_high_limit
673 fh2_wat(i) = karman*xxfh2(i)
675 if (ch_bound_excursion)
then
676 fh2_fh_ratio = min(xxfh2(i)/xxfh(i), 1.0)
677 fh_wat(i) = karman*karman/(fm_wat(i)*ch_wat(i))
678 fh2_wat(i) = fh2_fh_ratio*fh_wat(i)
685 ustar_wat(i) = 0.01*sqrt(cdm_wat(i)* &
686 (upc(i)*upc(i) + vpc(i)*vpc(i)))
688 ustar_wat(i) = amax1(ustar_wat(i),0.001)
690 stress_wat(i) = cdm_wat(i)*wspd(i)*wspd(i)
695 if ( wind10(i) .ge. 0.1 )
then
696 cd10=cdm_wat(i)* (wspd(i)/(0.01*wind10(i)) )**2
701 fm10_wat(i) = karman/sqrt(cd10)
708 esat = fpvs(tskin_wat(i))
709 qss_wat(i) = ep2*
esat/(psfc(i)-
esat)
746 SUBROUTINE mflux2( fxh,fxe,fxmx,fxmy,cdm,rib,xxfh,zoc,mzoc,tstrc, & !mzoc KWON
747 pspc,pkmax,wetc,slwdc,z1, &
748 icoef_sf,iwavecpl,lcurr_sf,alpha,gamma,xcur,ycur, &
749 pert_Cd, ens_random_seed, ens_Cdamp, &
750 upc,vpc,tpc,rpc,dt,wind10,xxfh2,ntsflg,sfenth, &
751 tzot, ep2, errmsg, errflg)
772 integer,
parameter :: ims = 1
773 integer,
parameter :: ime = 1
774 integer,
parameter :: its = 1
775 integer,
parameter :: ite = 1
776 integer,
intent(in) :: ntsflg
777 integer,
intent(in) :: icoef_sf
778 integer,
intent(in) :: iwavecpl
779 logical,
intent(in) :: lcurr_sf
780 logical,
intent(in) :: pert_Cd
781 integer,
intent(in) :: ens_random_seed
782 real(kind=kind_phys),
intent(in) :: ens_cdamp
784 real(kind=kind_phys),
intent (out),
dimension (ims :ime ) :: fxh
785 real(kind=kind_phys),
intent (out),
dimension (ims :ime ) :: fxe
786 real(kind=kind_phys),
intent (out),
dimension (ims :ime ) :: fxmx
787 real(kind=kind_phys),
intent (out),
dimension (ims :ime ) :: fxmy
788 real(kind=kind_phys),
intent (inout),
dimension (ims :ime ) :: cdm
790 real(kind=kind_phys),
intent (out),
dimension (ims :ime ) :: rib
791 real(kind=kind_phys),
intent (out),
dimension (ims :ime ) :: xxfh
792 real(kind=kind_phys),
intent (out),
dimension (ims :ime ) :: xxfh2
793 real(kind=kind_phys),
intent (out),
dimension (ims :ime ) :: wind10
795 real(kind=kind_phys),
intent ( inout),
dimension (ims :ime ) :: zoc,mzoc
796 real(kind=kind_phys),
intent ( inout),
dimension (ims :ime ) :: tzot
797 real(kind=kind_phys),
intent ( inout),
dimension (ims :ime ) :: tstrc
799 real(kind=kind_phys),
intent ( in) :: dt
800 real(kind=kind_phys),
intent ( in) :: sfenth
801 real(kind=kind_phys),
intent ( in),
dimension (ims :ime ) :: pspc
802 real(kind=kind_phys),
intent ( in),
dimension (ims :ime ) :: pkmax
803 real(kind=kind_phys),
intent ( in),
dimension (ims :ime ) :: wetc
804 real(kind=kind_phys),
intent ( in),
dimension (ims :ime ) :: slwdc
805 real(kind=kind_phys),
intent ( in),
dimension (ims :ime ) :: alpha, gamma
806 real(kind=kind_phys),
intent ( in),
dimension (ims :ime ) :: xcur, ycur
807 real(kind=kind_phys),
intent ( in),
dimension (ims :ime ) :: z1
809 real(kind=kind_phys),
intent ( in),
dimension (ims :ime ) :: upc
810 real(kind=kind_phys),
intent ( in),
dimension (ims :ime ) :: vpc
811 real(kind=kind_phys),
intent ( in),
dimension (ims :ime ) :: tpc
812 real(kind=kind_phys),
intent ( in),
dimension (ims :ime ) :: rpc
814 real(kind=kind_phys),
intent ( in) :: ep2
816 character(len=*),
intent(out) :: errmsg
817 integer,
intent(out) :: errflg
823 integer,
parameter :: icntx = 30
825 integer,
dimension(1 :ime) :: ifz
826 integer,
dimension(1 :ime) :: indx
827 integer,
dimension(1 :ime) :: istb
828 integer,
dimension(1 :ime) :: it
829 integer,
dimension(1 :ime) :: iutb
831 real(kind=kind_phys),
dimension(1 :ime) :: aap
832 real(kind=kind_phys),
dimension(1 :ime) :: bq1
833 real(kind=kind_phys),
dimension(1 :ime) :: bq1p
834 real(kind=kind_phys),
dimension(1 :ime) :: delsrad
835 real(kind=kind_phys),
dimension(1 :ime) :: ecof
836 real(kind=kind_phys),
dimension(1 :ime) :: ecofp
837 real(kind=kind_phys),
dimension(1 :ime) :: estso
838 real(kind=kind_phys),
dimension(1 :ime) :: estsop
839 real(kind=kind_phys),
dimension(1 :ime) :: fmz1
840 real(kind=kind_phys),
dimension(1 :ime) :: fmz10
841 real(kind=kind_phys),
dimension(1 :ime) :: fmz2
842 real(kind=kind_phys),
dimension(1 :ime) :: fmzo1
843 real(kind=kind_phys),
dimension(1 :ime) :: foft
844 real(kind=kind_phys),
dimension(1 :ime) :: foftm
845 real(kind=kind_phys),
dimension(1 :ime) :: frac
846 real(kind=kind_phys),
dimension(1 :ime) :: land
847 real(kind=kind_phys),
dimension(1 :ime) :: pssp
848 real(kind=kind_phys),
dimension(1 :ime) :: qf
849 real(kind=kind_phys),
dimension(1 :ime) :: rdiff
850 real(kind=kind_phys),
dimension(1 :ime) :: rho
851 real(kind=kind_phys),
dimension(1 :ime) :: rkmaxp
852 real(kind=kind_phys),
dimension(1 :ime) :: rstso
853 real(kind=kind_phys),
dimension(1 :ime) :: rstsop
854 real(kind=kind_phys),
dimension(1 :ime) :: sf10
855 real(kind=kind_phys),
dimension(1 :ime) :: sf2
856 real(kind=kind_phys),
dimension(1 :ime) :: sfm
857 real(kind=kind_phys),
dimension(1 :ime) :: sfzo
858 real(kind=kind_phys),
dimension(1 :ime) :: sgzm
859 real(kind=kind_phys),
dimension(1 :ime) :: slwa
860 real(kind=kind_phys),
dimension(1 :ime) :: szeta
861 real(kind=kind_phys),
dimension(1 :ime) :: szetam
862 real(kind=kind_phys),
dimension(1 :ime) :: t1
863 real(kind=kind_phys),
dimension(1 :ime) :: t2
864 real(kind=kind_phys),
dimension(1 :ime) :: tab1
865 real(kind=kind_phys),
dimension(1 :ime) :: tab2
866 real(kind=kind_phys),
dimension(1 :ime) :: tempa1
867 real(kind=kind_phys),
dimension(1 :ime) :: tempa2
868 real(kind=kind_phys),
dimension(1 :ime) :: theta
869 real(kind=kind_phys),
dimension(1 :ime) :: thetap
870 real(kind=kind_phys),
dimension(1 :ime) :: tsg
871 real(kind=kind_phys),
dimension(1 :ime) :: tsm
872 real(kind=kind_phys),
dimension(1 :ime) :: tsp
873 real(kind=kind_phys),
dimension(1 :ime) :: tss
874 real(kind=kind_phys),
dimension(1 :ime) :: ucom
875 real(kind=kind_phys),
dimension(1 :ime) :: uf10
876 real(kind=kind_phys),
dimension(1 :ime) :: uf2
877 real(kind=kind_phys),
dimension(1 :ime) :: ufh
878 real(kind=kind_phys),
dimension(1 :ime) :: ufm
879 real(kind=kind_phys),
dimension(1 :ime) :: ufzo
880 real(kind=kind_phys),
dimension(1 :ime) :: ugzm
881 real(kind=kind_phys),
dimension(1 :ime) :: uzeta
882 real(kind=kind_phys),
dimension(1 :ime) :: uzetam
883 real(kind=kind_phys),
dimension(1 :ime) :: vcom
884 real(kind=kind_phys),
dimension(1 :ime) :: vrtkx
885 real(kind=kind_phys),
dimension(1 :ime) :: vrts
886 real(kind=kind_phys),
dimension(1 :ime) :: wind
887 real(kind=kind_phys),
dimension(1 :ime) :: windp
888 real(kind=kind_phys),
dimension(1 :ime) :: wind10p
889 real(kind=kind_phys),
dimension(1 :ime) :: uvs1
891 real(kind=kind_phys),
dimension(1 :ime) :: xxfm
892 real(kind=kind_phys),
dimension(1 :ime) :: xxsh
893 real(kind=kind_phys),
dimension(1 :ime) :: z10
894 real(kind=kind_phys),
dimension(1 :ime) :: z2
895 real(kind=kind_phys),
dimension(1 :ime) :: zeta
896 real(kind=kind_phys),
dimension(1 :ime) :: zkmax
898 real(kind=kind_phys),
dimension(1 :ime) :: pss
899 real(kind=kind_phys),
dimension(1 :ime) :: tstar
900 real(kind=kind_phys),
dimension(1 :ime) :: ukmax
901 real(kind=kind_phys),
dimension(1 :ime) :: vkmax
902 real(kind=kind_phys),
dimension(1 :ime) :: tkmax
903 real(kind=kind_phys),
dimension(1 :ime) :: rkmax
904 real(kind=kind_phys),
dimension(1 :ime) :: zot
905 real(kind=kind_phys),
dimension(1 :ime) :: fhzo1
906 real(kind=kind_phys),
dimension(1 :ime) :: sfh
908 real(kind=kind_phys) :: ux13, yo, y,xo,x,ux21,ugzzo,ux11,ux12,uzetao,xnum,alll
909 real(kind=kind_phys) :: ux1,ugz,x10,uzo,uq,ux2,ux3,xtan,xden,y10,uzet1o,ugz10
910 real(kind=kind_phys) :: szet2, zal2,ugz2
911 real(kind=kind_phys) :: rovcp,boycon,cmo2,psps1,zog,enrca,rca,cmo1,amask,en,ca,a,c
912 real(kind=kind_phys) :: sgz,zal10,szet10,fmz,szo,sq,fmzo,rzeta1,zal1g,szetao,rzeta2,zal2g
913 real(kind=kind_phys) :: hcap,xks,pith,teps,diffot,delten,alevp,psps2,alfus,nstep
914 real(kind=kind_phys) :: shfx,sigt4,reflect
915 real(kind=kind_phys) :: cor1,cor2,szetho,zal2gh,cons_p000001,cons_7,vis,ustar,restar,rat
916 real(kind=kind_phys) :: wndm,ckg
917 real(kind=kind_phys) :: windmks,znott,znotm
918 real(kind=kind_phys) :: ubot, vbot
919 integer:: i,j,ii,iq,nnest,icnt,ngd,ip
925 real(kind=kind_phys),
dimension (223) :: tab
926 real(kind=kind_phys),
dimension (223) :: table
927 real(kind=kind_phys),
dimension (101) :: tab11
928 real(kind=kind_phys),
dimension (41) :: table4
929 real(kind=kind_phys),
dimension (42) :: tab3
930 real(kind=kind_phys),
dimension (54) :: table2
931 real(kind=kind_phys),
dimension (54) :: table3
932 real(kind=kind_phys),
dimension (74) :: table1
933 real(kind=kind_phys),
dimension (80) :: tab22
935 character(len=255) :: message
937 equivalence(tab(1),tab11(1))
938 equivalence(tab(102),tab22(1))
939 equivalence(tab(182),tab3(1))
940 equivalence(table(1),table1(1))
941 equivalence(table(75),table2(1))
942 equivalence(table(129),table3(1))
943 equivalence(table(183),table4(1))
951 data tab11/21*0.01403,0.01719,0.02101,0.02561,0.03117,0.03784, &
952 &.04584,.05542,.06685,.08049,.09672,.1160,.1388,.1658,.1977,.2353, &
953 &.2796,.3316,.3925,.4638,.5472,.6444,.7577,.8894,1.042,1.220,1.425, &
954 &1.662,1.936,2.252,2.615,3.032,3.511,4.060,4.688,5.406,6.225,7.159, &
955 &8.223,9.432,10.80,12.36,14.13,16.12,18.38,20.92,23.80,27.03,30.67, &
956 &34.76,39.35,44.49,50.26,56.71,63.93,71.98,80.97,90.98,102.1,114.5, &
957 &128.3,143.6,160.6,179.4,200.2,223.3,248.8,276.9,307.9,342.1,379.8, &
958 &421.3,466.9,517.0,572.0,632.3,698.5,770.9,850.2,937.0,1032./
960 data tab22/1146.6,1272.0,1408.1,1556.7,1716.9,1890.3,2077.6,2279.6 &
961 &,2496.7,2729.8,2980.0,3247.8,3534.1,3839.8,4164.8,4510.5,4876.9, &
962 &5265.1,5675.2,6107.8,6566.2,7054.7,7575.3,8129.4,8719.2,9346.5, &
963 &10013.,10722.,11474.,12272.,13119.,14017.,14969.,15977.,17044., &
964 &18173.,19367.,20630.,21964.,23373.,24861.,26430.,28086.,29831., &
965 &31671.,33608.,35649.,37796.,40055.,42430.,44927.,47551.,50307., &
966 &53200.,56236.,59422.,62762.,66264.,69934.,73777.,77802.,82015., &
967 &86423.,91034.,95855.,100890.,106160.,111660.,117400.,123400., &
968 &129650.,136170.,142980.,150070.,157460.,165160.,173180.,181530., &
971 data tab3/208670.,218450.,228610.,239180.,250160.,261560.,273400., &
972 &285700.,298450.,311690.,325420.,339650.,354410.,369710.,385560., &
973 &401980.,418980.,436590.,454810.,473670.,493170.,513350.,534220., &
974 &555800.,578090.,601130.,624940.,649530.,674920.,701130.,728190., &
975 &756110.,784920.,814630.,845280.,876880.,909450.,943020.,977610., &
976 &1013250.,1049940.,1087740./
978 data table1/20*0.0,.3160e-02,.3820e-02,.4600e-02,.5560e-02,.6670e-02, &
979 & .8000e-02,.9580e-02,.1143e-01,.1364e-01,.1623e-01,.1928e-01, &
980 &.2280e-01,.2700e-01,.3190e-01,.3760e-01,.4430e-01,.5200e-01, &
981 &.6090e-01,.7130e-01,.8340e-01,.9720e-01,.1133e+00,.1317e-00, &
982 &.1526e-00,.1780e-00,.2050e-00,.2370e-00,.2740e-00,.3160e-00, &
983 &.3630e-00,.4170e-00,.4790e-00,.5490e-00,.6280e-00,.7180e-00, &
984 &.8190e-00,.9340e-00,.1064e+01,.1209e+01,.1368e+01,.1560e+01, &
985 &.1770e+01,.1990e+01,.2260e+01,.2540e+01,.2880e+01,.3230e+01, &
986 &.3640e+01,.4090e+01,.4590e+01,.5140e+01,.5770e+01,.6450e+01, &
989 data table2/.8050e+01,.8990e+01,.1001e+02,.1112e+02,.1240e+02, &
990 &.1380e+02,.1530e+02,.1700e+02,.1880e+02,.2080e+02,.2310e+02, &
991 &.2550e+02,.2810e+02,.3100e+02,.3420e+02,.3770e+02,.4150e+02, &
992 &.4560e+02,.5010e+02,.5500e+02,.6030e+02,.6620e+02,.7240e+02, &
993 &.7930e+02,.8680e+02,.9500e+02,.1146e+03,.1254e+03,.1361e+03, &
994 &.1486e+03,.1602e+03,.1734e+03,.1873e+03,.2020e+03,.2171e+03, &
995 &.2331e+03,.2502e+03,.2678e+03,.2863e+03,.3057e+03,.3250e+03, &
996 &.3457e+03,.3664e+03,.3882e+03,.4101e+03,.4326e+03,.4584e+03, &
997 &.4885e+03,.5206e+03,.5541e+03,.5898e+03,.6273e+03,.6665e+03, &
1000 data table3/.7520e+03,.7980e+03,.8470e+03,.8980e+03,.9520e+03, &
1001 &.1008e+04,.1067e+04,.1129e+04,.1194e+04,.1263e+04,.1334e+04, &
1002 &.1409e+04,.1488e+04,.1569e+04,.1656e+04,.1745e+04,.1840e+04, &
1003 &.1937e+04,.2041e+04,.2147e+04,.2259e+04,.2375e+04,.2497e+04, &
1004 &.2624e+04,.2756e+04,.2893e+04,.3036e+04,.3186e+04,.3340e+04, &
1005 &.3502e+04,.3670e+04,.3843e+04,.4025e+04,.4213e+04,.4408e+04, &
1006 &.4611e+04,.4821e+04,.5035e+04,.5270e+04,.5500e+04,.5740e+04, &
1007 &.6000e+04,.6250e+04,.6520e+04,.6810e+04,.7090e+04,.7390e+04, &
1008 &.7700e+04,.8020e+04,.8350e+04,.8690e+04,.9040e+04,.9410e+04, &
1011 data table4/.1016e+05,.1057e+05,.1098e+05,.1140e+05,.1184e+05, &
1012 &.1230e+05,.1275e+05,.1324e+05,.1373e+05,.1423e+05,.1476e+05, &
1013 &.1530e+05,.1585e+05,.1642e+05,.1700e+05,.1761e+05,.1822e+05, &
1014 &.1886e+05,.1950e+05,.2018e+05,.2087e+05,.2158e+05,.2229e+05, &
1015 &.2304e+05,.2381e+05,.2459e+05,.2539e+05,.2621e+05,.2706e+05, &
1016 &.2792e+05,.2881e+05,.2971e+05,.3065e+05,.3160e+05,.3257e+05, &
1017 &.3357e+05,.3459e+05,.3564e+05,.3669e+05,.3780e+05,.0000e+00/
1022 real,
parameter :: cp = 1.00464e7
1023 real,
parameter :: g = 980.6
1024 real,
parameter :: rgas = 2.87e6
1025 real,
parameter :: og = 1./g
1026 integer :: ntstep = 0
1033 real*8 :: gasdev,ran1
1035 logical,
save :: pert_Cd_local
1036 CHARACTER(len=3) :: env_memb,env_pp
1037 integer,
save :: ens_random_seed_local,env_pp_local
1038 integer :: ensda_physics_pert
1039 real,
save :: ens_Cdamp_local
1040 data ens_random_seed_local/0/
1041 data env_pp_local/0/
1042 if ( ens_random_seed_local .eq. 0 )
then
1043 CALL nl_get_ensda_physics_pert(1,ensda_physics_pert)
1044 ens_random_seed_local=ens_random_seed
1045 env_pp_local=ensda_physics_pert
1046 pert_cd_local=.false.
1049 if ( env_pp_local .eq. 1 )
then
1050 if ( ens_random_seed .ne. 99 )
then
1051 pert_cd_local=.true.
1052 ens_cdamp_local=ens_cdamp
1055 ens_random_seed_local=ens_random_seed
1056 pert_cd_local=pert_cd
1057 ens_cdamp_local=ens_cdamp
1060 ens_random_seed_local=ens_random_seed
1061 pert_cd_local=pert_cd
1062 ens_cdamp_local=ens_cdamp
1064 print*,
"Cd ===", ens_random_seed_local,pert_cd_local,ens_cdamp_local,ensda_physics_pert
1091 if ( lcurr_sf .and. zoc(i) .le. 0.0 )
then
1092 ubot = upc(i) - xcur(i) * 100.0
1093 vbot = vpc(i) - ycur(i) * 100.0
1100 uvs1(i)= amax1( sqrt(ubot*ubot + &
1102 if ( iwavecpl .eq. 1 .and. zoc(i) .le. 0.0 )
then
1103 ukmax(i) = ( ubot * cos(gamma(i)) - &
1104 vbot * sin(gamma(i)) ) &
1106 vkmax(i) = ( vbot * cos(gamma(i)) - &
1107 ubot * sin(gamma(i)) ) &
1122 windp(i) = sqrt(ukmax(i)*ukmax(i) + vkmax(i)*vkmax(i))
1123 wind(i) = amax1(windp(i),100.)
1126 wind10p(i) = wind10(i)
1127 wind10p(i) = amax1(wind10p(i),100.)
1130 if (zoc(i) .LT. amask) zoc(i) = -0.0185*0.001*wind10p(i)*wind10p(i)*og
1131 if (zoc(i) .GT. 0.0)
then
1138 windmks=wind10p(i)*.01
1139 if ( iwavecpl .eq. 1 )
then
1140 call znot_wind10m(windmks,znott,znotm,icoef_sf,errmsg,errflg)
1142 if ( alpha(i) .ge. 0.2 .and. alpha(i) .le. 5. )
then
1143 znotm = znotm*alpha(i)
1145 zoc(i) = -100.*znotm
1146 zot(i) = -100* znott
1148 call znot_wind10m(windmks,znott,znotm,icoef_sf,errmsg,errflg)
1149 zoc(i) = -100.*znotm
1150 zot(i) = -100* znott
1177 cmo2 = 17.193 + .5*a - 10.*c
1182 theta(i) = tkmax(i)/((pkmax(i)/pspc(i))**rovcp)
1183 vrtkx(i) = 1.0 + boycon*rkmax(i)
1194 tab1(i) = tstar(i) - 153.16
1195 it(i) = ifix(tab1(i))
1196 tab2(i) = tab1(i) - float(it(i))
1197 t1(i) = tab(min(223,max(1,it(i) + 1)))
1198 t2(i) = table(min(223,max(1,it(i) + 1)))
1199 estso(i) = t1(i) + tab2(i)*t2(i)
1200 psps1 = (pss(i) - estso(i))
1201 if(psps1 .EQ. 0.0)
then
1204 rstso(i) = ep2*estso(i)/psps1
1205 vrts(i) = 1. + boycon*ecof(i)*rstso(i)
1216 tempa1(i) = theta(i)*vrtkx(i) - tstar(i)*vrts(i)
1217 tempa2(i) = tempa1(i)*(theta(i) - tstar(i))
1218 if (tempa2(i) .LT. 0.) tempa1(i) = 1.0e-4
1219 tab1(i) = abs(tempa1(i))
1220 if (tab1(i) .LT. 1.0e-4) tempa1(i) = 1.0e-4
1226 rib(i) = g*zkmax(i)*tempa1(i)/ &
1227 (tkmax(i)*vrtkx(i)*wind(i)*wind(i))
1228 tab2(i) = abs(zoc(i))
1229 tab1(i) = 0.95/(c*(1. - tab2(i)/zkmax(i)))
1230 if (rib(i) .GT. tab1(i)) rib(i) = tab1(i)
1234 zeta(i) = ca*rib(i)/0.03
1259 if (zeta(i) .GE. 0.0)
then
1265 if (ip .EQ. 0)
go to 170
1269 szeta(i) = zeta(istb(i))
1280 if (ifz(i) .EQ. 0)
go to 80
1281 zal1g = alog(szeta(i))
1282 if (szeta(i) .LE. 0.5)
then
1283 fmz1(i) = (zal1g + a*szeta(i))*rca
1284 else if (szeta(i) .GT. 0.5 .AND. szeta(i) .LE. 10.)
then
1285 rzeta1 = 1./szeta(i)
1286 fmz1(i) = (8.*zal1g + 4.25*rzeta1 - &
1287 0.5*rzeta1*rzeta1 + cmo1)*rca
1288 else if (szeta(i) .GT. 10.)
then
1289 fmz1(i) = (c*szeta(i) + cmo2)*rca
1291 szetao = abs(zoc(istb(i)))/zkmax(istb(i))*szeta(i)
1292 zal2g = alog(szetao)
1293 if (szetao .LE. 0.5)
then
1294 fmzo1(i) = (zal2g + a*szetao)*rca
1295 sfzo(i) = 1. + a*szetao
1296 else if (szetao .GT. 0.5 .AND. szetao .LE. 10.)
then
1298 fmzo1(i) = (8.*zal2g + 4.25*rzeta2 - &
1299 0.5*rzeta2*rzeta2 + cmo1)*rca
1300 sfzo(i) = 8.0 - 4.25*rzeta2 + rzeta2*rzeta2
1301 else if (szetao .GT. 10.)
then
1302 fmzo1(i) = (c*szetao + cmo2)*rca
1309 szetho = abs(zot(istb(i)))/zkmax(istb(i))*szeta(i)
1310 zal2gh = alog(szetho)
1311 if (szetho .LE. 0.5)
then
1312 fhzo1(i) = (zal2gh + a*szetho)*rca
1313 sfzo(i) = 1. + a*szetho
1314 else if (szetho .GT. 0.5 .AND. szetho .LE. 10.)
then
1316 fhzo1(i) = (8.*zal2gh + 4.25*rzeta2 - &
1317 0.5*rzeta2*rzeta2 + cmo1)*rca
1318 sfzo(i) = 8.0 - 4.25*rzeta2 + rzeta2*rzeta2
1319 else if (szetho .GT. 10.)
then
1320 fhzo1(i) = (c*szetho + cmo2)*rca
1328 szet10 = abs(z10(istb(i)))/zkmax(istb(i))*szeta(i)
1329 zal10 = alog(szet10)
1330 if (szet10 .LE. 0.5)
then
1331 fmz10(i) = (zal10 + a*szet10)*rca
1332 else if (szet10 .GT. 0.5 .AND. szet10 .LE. 10.)
then
1334 fmz10(i) = (8.*zal10 + 4.25*rzeta2 - &
1335 0.5*rzeta2*rzeta2 + cmo1)*rca
1336 else if (szet10 .GT. 10.)
then
1337 fmz10(i) = (c*szet10 + cmo2)*rca
1339 sf10(i) = fmz10(i) - fmzo1(i)
1341 szet2 = abs(z2(istb(i)))/zkmax(istb(i))*szeta(i)
1343 if (szet2 .LE. 0.5)
then
1344 fmz2(i) = (zal2 + a*szet2 )*rca
1345 else if (szet2 .GT. 0.5 .AND. szet2 .LE. 2.)
then
1347 fmz2(i) = (8.*zal2 + 4.25*rzeta2 - &
1348 0.5*rzeta2*rzeta2 + cmo1)*rca
1349 else if (szet2 .GT. 2.)
then
1350 fmz2(i) = (c*szet2 + cmo2)*rca
1352 sf2(i) = fmz2(i) - fmzo1(i)
1354 sfm(i) = fmz1(i) - fmzo1(i)
1355 sfh(i) = fmz1(i) - fhzo1(i)
1356 sgz = ca*rib(istb(i))*sfm(i)*sfm(i)/ &
1357 (sfh(i) + enrca*sfzo(i))
1358 fmz = (sgz - szeta(i))/szeta(i)
1360 if (fmzo .GE. 5.0e-5)
then
1361 sq = (sgz - sgzm(i))/(szeta(i) - szetam(i))
1363 write(errmsg,
'(*(a))')
'NCO ERROR DIVIDE BY ZERO IN gfdl_sfc_layer.F90/MFLUX2 (STABLE CASE)'// &
1364 'sq is 1 ',fmzo,sgz,sgzm(i),szeta(i),szetam(i)
1368 szetam(i) = szeta(i)
1369 szeta(i) = (sgz - szeta(i)*sq)/(1.0 - sq)
1379 if (ifz(i) .GE. 1)
go to 110
1386 write(errmsg,
'(*(a))')
'NON-CONVERGENCE FOR STABLE ZETA IN gfdl_sfc_layer.F90/MFLUX2'
1400 if (szo .LT. 0.0)
then
1401 wndm=wind(istb(i))*0.01
1402 if(wndm.lt.15.0)
then
1405 ckg=(sfenth*(4*0.000308*wndm) + (1.-sfenth)*0.0185 )*og
1408 szo = - ckg*wind(istb(i))*wind(istb(i))/ &
1410 cons_p000001 = .000001
1414 ustar = sqrt( -szo / zog)
1415 restar = -ustar * szo / vis
1416 restar = max(restar,cons_p000001)
1418 rat = 2.67 * restar ** .25 - 2.57
1419 rat = min(rat ,cons_7)
1421 zot(istb(i)) = szo * exp(-rat)
1423 zot(istb(i)) = zoc(istb(i))
1431 xxfm(istb(i)) = sfm(i)
1432 xxfh(istb(i)) = sfh(i)
1433 xxfh2(istb(i)) = sf2(i)
1434 xxsh(istb(i)) = sfzo(i)
1442 wind10(istb(i)) = sf10(i)*uvs1(istb(i))/sfm(i)
1443 wind10(istb(i)) = wind10(istb(i)) * 1.944
1444 if(wind10(istb(i)) .GT. 6000.0)
then
1445 wind10(istb(i))=wind10(istb(i))+wind10(istb(i))*cor1 &
1449 wind10(istb(i)) = wind10(istb(i)) / 1.944
1460 if (zeta(i) .LT. 0.0)
then
1466 if (iq .EQ. 0)
go to 290
1468 uzeta(i) = zeta(iutb(i))
1481 if (ifz(i) .EQ. 0)
go to 200
1482 ugzzo = alog(zkmax(iutb(i))/abs(zot(iutb(i))))
1483 uzetao = abs(zot(iutb(i)))/zkmax(iutb(i))*uzeta(i)
1484 ux11 = 1. - 16.*uzeta(i)
1485 ux12 = 1. - 16.*uzetao
1489 ux13 = (1. + y)/(1. + yo)
1491 ufh(i) = (ugzzo - 2.*ux21)*rca
1493 ugzzo = alog(zkmax(iutb(i))/abs(zoc(iutb(i))))
1494 uzetao = abs(zoc(iutb(i)))/zkmax(iutb(i))*uzeta(i)
1495 ux11 = 1. - 16.*uzeta(i)
1496 ux12 = 1. - 16.*uzetao
1499 ux13 = (1. + y)/(1. + yo)
1504 xnum = (x**2 + 1.)*((x + 1.)**2)
1505 xden = (xo**2 + 1.)*((xo + 1.)**2)
1506 xtan = atan(x) - atan(xo)
1507 ux3 = alog(xnum/xden)
1508 ufm(i) = (ugzzo - ux3 + 2.*xtan)*rca
1514 ugz10 = alog(z10(iutb(i))/abs(zoc(iutb(i))))
1515 uzet1o = abs(z10(iutb(i)))/zkmax(iutb(i))*uzeta(i)
1516 uzetao = abs(zoc(iutb(i)))/zkmax(iutb(i))*uzeta(i)
1517 ux11 = 1. - 16.*uzet1o
1518 ux12 = 1. - 16.*uzetao
1521 ux13 = (1. + y)/(1. + y10)
1525 xnum = (x**2 + 1.)*((x + 1.)**2)
1526 xden = (x10**2 + 1.)*((x10 + 1.)**2)
1527 xtan = atan(x) - atan(x10)
1528 ux3 = alog(xnum/xden)
1529 uf10(i) = (ugz10 - ux3 + 2.*xtan)*rca
1534 ugz2 = alog(z2(iutb(i))/abs(zoc(iutb(i))))
1535 uzet1o = abs(z2(iutb(i)))/zkmax(iutb(i))*uzeta(i)
1536 uzetao = abs(zoc(iutb(i)))/zkmax(iutb(i))*uzeta(i)
1537 ux11 = 1. - 16.*uzet1o
1538 ux12 = 1. - 16.*uzetao
1541 ux13 = (1. + y)/(1. + yo)
1543 uf2(i) = (ugzzo - 2.*ux21)*rca
1546 ugz = ca*rib(iutb(i))*ufm(i)*ufm(i)/(ufh(i) + enrca*ufzo(i))
1547 ux1 = (ugz - uzeta(i))/uzeta(i)
1549 if (ux2 .GE. 5.0e-5)
then
1550 uq = (ugz - ugzm(i))/(uzeta(i) - uzetam(i))
1551 uzetam(i) = uzeta(i)
1553 write(errmsg,
'(*(a))')
'NCO ERROR DIVIDE BY ZERO IN gfdl_sfc_layer.F90/MFLUX2 (UNSTABLE CASE)'// &
1554 'uq is 1 ',ux2,ugz,ugzm(i),uzeta(i),uzetam(i)
1558 uzeta(i) = (ugz - uzeta(i)*uq)/(1.0 - uq)
1569 if (ifz(i) .GE. 1)
go to 230
1575 write(errmsg,
'(*(a))')
'NON-CONVERGENCE FOR UNSTABLE ZETA IN ROW'// &
1576 'uq is 1 ',ux2,ugz,ugzm(i),uzeta(i),uzetam(i)
1596 if (zoc(iutb(i)) .LT. 0.0)
then
1597 wndm=wind(iutb(i))*0.01
1598 if(wndm.lt.15.0)
then
1601 ckg=(4*0.000308*wndm)*og
1602 ckg=(sfenth*(4*0.000308*wndm) + (1.-sfenth)*0.0185 )*og
1604 uzo =-ckg*wind(iutb(i))*wind(iutb(i))/(ufm(i)*ufm(i))
1605 cons_p000001 = .000001
1609 ustar = sqrt( -uzo / zog)
1610 restar = -ustar * uzo / vis
1611 restar = max(restar,cons_p000001)
1613 rat = 2.67 * restar ** .25 - 2.57
1614 rat = min(rat ,cons_7)
1616 zot(iutb(i)) = uzo * exp(-rat)
1618 zot(iutb(i)) = zoc(iutb(i))
1628 wind10(iutb(i)) = uf10(i)*uvs1(iutb(i))/ufm(i)
1629 wind10(iutb(i)) = wind10(iutb(i)) * 1.944
1630 if(wind10(iutb(i)) .GT. 6000.0)
then
1631 wind10(iutb(i))=wind10(iutb(i))+wind10(iutb(i))*cor1 &
1635 wind10(iutb(i)) = wind10(iutb(i)) / 1.944
1639 xxfm(iutb(i)) = ufm(i)
1640 xxfh(iutb(i)) = ufh(i)
1641 xxfh2(iutb(i)) = uf2(i)
1642 xxsh(iutb(i)) = ufzo(i)
1650 if (windp(i) .EQ. 0.0)
then
1652 ucom(i) = 100.0/sqrt(2.0)
1653 vcom(i) = 100.0/sqrt(2.0)
1655 rho(i) = pss(i)/(rgas*(tsg(i) + enrca*(theta(i) - &
1656 tsg(i))*xxsh(i)/(xxfh(i) + enrca*xxsh(i))))
1657 bq1(i) = wind(i)*rho(i)/(xxfm(i)*(xxfh(i) + enrca*xxsh(i)))
1663 if (ntsflg .EQ. 0)
go to 370
1667 pith = sqrt(4.*atan(1.0))
1668 alfus = alll/2.39e-8
1679 if (land(i) .EQ. 1)
then
1683 slwa(ip) = slwdc(i)/(2.39e-8*60.)
1685 thetap(ip) = theta(i)
1686 rkmaxp(ip) = rkmax(i)
1690 estsop(ip) = estso(i)
1691 rstsop(ip) = rstso(i)
1693 bq1p(ip) = amax1(bq1p(ip),0.1e-3)
1694 delsrad(ip) = dt *pith/(hcap*sqrt(3600.*24.*xks))
1705 rdiff(i) = amin1(0.0,(rkmaxp(i) - rstsop(i)))
1707300
format(2x,
' SURFACE EQUILIBRIUM CALCULATION ')
1709 foftm(i) = tss(i) + delsrad(i)*(slwa(i) - aap(i)*tsm(i)**4 - &
1710 cp*bq1p(i)*(tsm(i) - thetap(i)) + ecofp(i)*alfus*bq1p(i)* &
1721 if (ifz(i) .EQ. 0)
go to 330
1722 tab1(i) = tsp(i) - 153.16
1723 it(i) = ifix(tab1(i))
1724 tab2(i) = tab1(i) - float(it(i))
1725 t1(i) = tab(min(223,max(1,it(i) + 1)))
1726 t2(i) = table(min(223,max(1,it(i) + 1)))
1727 estsop(i) = t1(i) + tab2(i)*t2(i)
1728 psps2 = (pssp(i) - estsop(i))
1729 if(psps2 .EQ. 0.0)
then
1732 rstsop(i) = ep2*estsop(i)/psps2
1733 rdiff(i) = amin1(0.0,(rkmaxp(i) - rstsop(i)))
1735 foft(i) = tss(i) + delsrad(i)*(slwa(i) - aap(i)*tsp(i)**4 - &
1736 cp*bq1p(i)*(tsp(i) - thetap(i)) + ecofp(i)*alfus*bq1p(i)* &
1739 frac(i) = abs((foft(i) - tsp(i))/tsp(i))
1745 if (frac(i) .GE. teps)
then
1746 qf(i) = (foft(i) - foftm(i))/(tsp(i) - tsm(i))
1748 tsp(i) = (foft(i) - tsp(i)*qf(i))/(1. - qf(i))
1762 if (ifz(i) .EQ. 1)
then
1763 write(errmsg,
'(*(a))')
'NON-CONVERGENCE OF T* PREDICTED (T*,I) = ', &
1783 if ( iwavecpl .eq. 1 .and. zoc(i) .le. 0.0 )
then
1784 windmks = wind10(i) * 0.01
1785 call znot_wind10m(windmks,znott,znotm,icoef_sf,errmsg,errflg)
1787 if ( alpha(i) .ge. 0.2 .and. alpha(i) .le. 5. )
then
1788 znotm = znotm*alpha(i)
1790 zoc(i) = -100.*znotm
1791 zot(i) = -100* znott
1794 fxh(i) = bq1(i)*(theta(i) - tsg(i))
1795 fxe(i) = ecof(i)*bq1(i)*(rkmax(i) - rstso(i))
1796 if (fxe(i) .GT. 0.0) fxe(i) = 0.0
1797 fxmx(i) = rho(i)/(xxfm(i)*xxfm(i))*wind(i)*wind(i)*ucom(i)/ &
1799 fxmy(i) = rho(i)/(xxfm(i)*xxfm(i))*wind(i)*wind(i)*vcom(i)/ &
1801 cdm(i) = 1./(xxfm(i)*xxfm(i))
1805 if( pert_cd_local )
then
1806 ens_random_seed_local=ran1(-ens_random_seed_local)*1000
1807 rr=2.0*ens_cdamp_local*ran1(-ens_random_seed_local)-ens_cdamp_local
1808 cdm(i) = cdm(i) *(1.0+rr)