1 !-- XLV latent heat of vaporization for water (J/kg)
5 !real, dimension(-100:2000,-100:2000), save :: z00
9 !-------------------------------------------------------------------
10 SUBROUTINE SF_GFDL(U3D,V3D,T3D,QV3D,P3D, &
11 CP,ROVCP,R,XLV,PSFC,CHS,CHS2,CQS2, CPM, &
12 DT, SMOIS,num_soil_layers,ISLTYP,ZNT, &
17 XLAND,HFX,QFX,TAUX,TAUY,LH,GSW,GLW,TSK,FLHC,FLQC, & ! gopal's doing for Ocean coupling
19 ICOEF_SF,IWAVECPL,LCURR_SF,CHARN,MSANG,SCURX, SCURY,&
20 pert_Cd, ens_random_seed, ens_Cdamp, &
21 GZ1OZ0,WSPD,BR,ZKMAX, ISFFLX, &
22 EP1,EP2,KARMAN,NTSFLG,SFENTH, &
24 ids,ide, jds,jde, kds,kde, &
25 ims,ime, jms,jme, kms,kme, &
26 its,ite, jts,jte, kts,kte )
27 !-------------------------------------------------------------------
28 USE MODULE_GFS_MACHINE, ONLY : kind_phys
29 USE MODULE_GFS_FUNCPHYS , ONLY : gfuncphys,fpvs
30 USE MODULE_GFS_PHYSCONS, grav => con_g
31 !-------------------------------------------------------------------
33 !-------------------------------------------------------------------
34 !-- U3D 3D u-velocity interpolated to theta points (m/s)
35 !-- V3D 3D v-velocity interpolated to theta points (m/s)
36 !-- T3D temperature (K)
37 !-- QV3D 3D water vapor mixing ratio (Kg/Kg)
38 !-- P3D 3D pressure (Pa)
39 !-- DT time step (second)
40 !-- CP heat capacity at constant pressure for dry air (J/kg/K)
42 !-- R gas constant for dry air (J/kg/K)
43 !-- XLV latent heat of vaporization for water (J/kg)
44 !-- PSFC surface pressure (Pa)
45 !-- ZNT thermal roughness length (m)
46 !-- MZNT momentum roughness length (m)
47 !-- MAVAIL surface moisture availability (between 0 and 1)
48 !-- UST u* in similarity theory (m/s)
49 !-- PSIM similarity stability function for momentum
50 !-- PSIH similarity stability function for heat
51 !-- XLAND land mask (1 for land, 2 for water)
52 !-- HFX upward heat flux at the surface (W/m^2)
53 !-- QFX upward moisture flux at the surface (kg/m^2/s)
54 !-- TAUX RHO*U**2 (Kg/m/s^2) ! gopal's doing for Ocean coupling
55 !-- TAUY RHO*U**2 (Kg/m/s^2) ! gopal's doing for Ocean coupling
56 !-- LH net upward latent heat flux at surface (W/m^2)
57 !-- GSW downward short wave flux at ground surface (W/m^2)
58 !-- GLW downward long wave flux at ground surface (W/m^2)
59 !-- TSK surface temperature (K)
60 !-- FLHC exchange coefficient for heat (m/s)
61 !-- FLQC exchange coefficient for moisture (m/s)
62 !-- QGH lowest-level saturated mixing ratio
63 !-- U10 diagnostic 10m u wind
64 !-- V10 diagnostic 10m v wind
65 !-- GZ1OZ0 log(z/z0) where z0 is roughness length
66 !-- WSPD wind speed at lowest model level (m/s)
67 !-- BR bulk Richardson number in surface layer
68 !-- ISFFLX isfflx=1 for surface heat and moisture fluxes
69 !-- EP1 constant for virtual temperature (R_v/R_d - 1) (dimensionless)
70 !-- KARMAN Von Karman constant
71 !-- SFENTH enthalpy flux factor 0 zot via charnock ..>0 zot enhanced>15m/s
72 !-- ids start index for i in domain
73 !-- ide end index for i in domain
74 !-- jds start index for j in domain
75 !-- jde end index for j in domain
76 !-- kds start index for k in domain
77 !-- kde end index for k in domain
78 !-- ims start index for i in memory
79 !-- ime end index for i in memory
80 !-- jms start index for j in memory
81 !-- jme end index for j in memory
82 !-- kms start index for k in memory
83 !-- kme end index for k in memory
84 !-- its start index for i in tile
85 !-- ite end index for i in tile
86 !-- jts start index for j in tile
87 !-- jte end index for j in tile
88 !-- kts start index for k in tile
89 !-- kte end index for k in tile
90 !-------------------------------------------------------------------
91 character*255 :: message
92 INTEGER, INTENT(IN) :: ids,ide, jds,jde, kds,kde, &
93 ims,ime, jms,jme, kms,kme, &
94 its,ite, jts,jte, kts,kte, &
95 ISFFLX,NUM_SOIL_LAYERS,NTSFLG
96 INTEGER, INTENT(IN) :: ICOEF_SF
97 INTEGER, INTENT(IN) :: IWAVECPL
98 LOGICAL, INTENT(IN) :: LCURR_SF
99 logical,intent(in) :: pert_Cd
100 integer,intent(in) :: ens_random_seed
101 real,intent(in) :: ens_Cdamp
103 REAL, INTENT(IN) :: &
114 REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(IN) :: &
120 INTEGER, DIMENSION( ims:ime , jms:jme ), INTENT(IN ):: ISLTYP
121 REAL, DIMENSION( ims:ime , 1:num_soil_layers, jms:jme ), INTENT(INOUT):: SMOIS
123 REAL, DIMENSION(ims:ime, jms:jme), INTENT(IN) :: &
128 REAL, DIMENSION(ims:ime, jms:jme), INTENT(IN) :: &
133 REAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: &
156 MZNT, & !KWON momentum zo
159 TAUX, & ! gopal's doing for Ocean coupling
162 REAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: &
166 !--------------------------- LOCAL VARS ------------------------------
174 REAL (kind=kind_phys) :: &
176 REAL, DIMENSION(1:30) :: MAXSMC, &
178 REAL (kind=kind_phys), DIMENSION(its:ite) :: &
218 REAL, DIMENSION(kms:kme, ims:ime) :: &
224 REAL, DIMENSION(ims:ime) :: &
229 mzoc, & !ADDED BY KWON FOR momentum Zo
259 real :: tmp9,zhalf ! wang, height of the first half level
261 DATA MAXSMC/0.339, 0.421, 0.434, 0.476, 0.476, 0.439, &
262 0.404, 0.464, 0.465, 0.406, 0.468, 0.468, &
263 0.439, 1.000, 0.200, 0.421, 0.000, 0.000, &
264 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, &
265 0.000, 0.000, 0.000, 0.000, 0.000, 0.000/
267 DATA DRYSMC/0.010, 0.028, 0.047, 0.084, 0.084, 0.066, &
268 0.067, 0.120, 0.103, 0.100, 0.126, 0.138, &
269 0.066, 0.000, 0.006, 0.028, 0.000, 0.000, &
270 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, &
271 0.000, 0.000, 0.000, 0.000, 0.000, 0.000/
277 ! call readzo(glat,glon,6,ims,ime,jms,jme,its,ite,jts,jte,z00)
284 WRITE(message,*)'WITHIN THE GFDL SCHEME, NTSFLG=1 FOR GFDL SLAB 2010 UPGRADS',NTSFLG
285 call wrf_debug(2,message)
287 !! write(0,*)'icoef_sf,lcurr_sf=',icoef_sf,lcurr_sf
294 PRSL1(i)=P3D(i,kts,j)*.001
296 if(xland(i,j).lt.1.99) then
298 smcdry=drysmc(isltyp(i,j))
299 smcmax=maxsmc(isltyp(i,j))
300 wetc(i)=(smc-smcdry)/(smcmax-smcdry)
301 wetc(i)=amin1(1.,amax1(wetc(i),0.))
303 ! convert from Pa to cgs...
304 pspc(i)=PSFC(i,j)*10.
305 pkmax(i)=P3D(i,kts,j)*10.
307 Q1(I) = QV3D(i,kts,j)
308 rpc(kts,i)=QV3D(i,kts,j)
312 SLIMSK(i)=ABS(XLAND(i,j)-2.)
316 tpc(kts,i)=T3D(i,kts,j)
318 upc(kts,i)=U3D(i,kts,j) * 100.
321 vpc(kts,i)=v3D(i,kts,j) * 100.
322 Z0RL(I) = ZNT(i,j)*100.
325 if(XLAND(i,j).gt.1.99) zoc(i)=- zoc(i)
326 ! Z0RL(I) = z00(i,j)*100.
327 ! slwdc... GFDL downward net flux in units of cal/(cm**2/min)
328 ! also divide by 10**4 to convert from /m**2 to /cm**2
329 slwdc(i)=gsw(i,j)+glw(i,j)
330 slwdc(i)=0.239*60.*slwdc(i)*1.e-4
332 alpha(i) = charn(i,j)
333 gamma(i) = msang(i,j)
338 !Wang, calulate height of the first half level
339 ! use previous u10 v10 to compute wind10, input to MFLUX to compute z0
340 wind10(i)=sqrt(u10(i,j)*u10(i,j)+v10(i,j)*v10(i,j))
341 !first time step, u10 and v10 may be zero
343 if (wind10(i) <= 1.0e-10 .or. wind10(i) > 150.0) then
344 ! write(0,*)'maybe first timestep'
345 zhalf = -R*tpc(kts,i)*alog(pkmax(i)/pspc(i))/grav
346 wind10(i)=sqrt(u1(i)*u1(i)+v1(i)*v1(i))*alog(10.0/znt(i,j))/alog(zhalf/znt(i,j))
348 wind10(i)=wind10(i)*100.0 !! m/s to cm/s
349 zhalf = -R*tpc(kts,i)*alog(pkmax(i)/pspc(i))/grav
350 ! if (i == (its+ite)/2 .and. j == (jts+jte)/2 ) then
351 ! write(0,*)'before call mflux,zh,w10,u10,v10,WIND1,znt'
352 ! write(0,*)zhalf,wind10(i),u10(i,j),v10(i,j),sqrt(u1(i)**2+v1(i)**2),znt(i,j)
353 ! write(0,*)'uselog', sqrt(u1(i)*u1(i)+v1(i)*v1(i))*alog(10.0/znt(i,j))/alog(zhalf/znt(i,j))
360 PRSLKI(i)=(PS(I)/PRSL1(I))**ROVCP
361 THGB(I)=TSKIN(i)*(100./PS(I))**ROVCP
362 THX(I)=T1(i)*(100./PRSL1(I))**ROVCP
363 RHO1(I)=PRSL1(I)*1000./(R*T1(I)*(1.+EP1*Q1(I)))
364 Q1(I)=Q1(I)/(1.+Q1(I))
368 ! write(0,*)'--------------------------------------------'
369 ! write(0,*) 'u, v, t, r, pkmax, pspc,wetc, tjloc,zoc,tstr'
370 ! write(0,*)'--------------------------------------------'
374 ! WRITE(0,1010)i,j,upc(kts,i),vpc(kts,i),tpc(kts,i),rpc(kts,i), &
375 ! pkmax(i),pspc(i),wetc(i),tjloc(i),zoc(i),tstrc(i)
378 CALL MFLUX2( fxh,fxe,fxmx,fxmy,cdm,rib,xxfh,zoc,mzoc,tstrc, & !mzoc for momentum Zo KWON
379 pspc,pkmax,wetc,slwdc,tjloc, &
380 icoef_sf,iwavecpl,lcurr_sf,alpha,gamma,xcur,ycur, &
381 pert_Cd, ens_random_seed, ens_Cdamp, &
382 upc,vpc,tpc,rpc,dt,J,wind10,xxfh2,ntsflg,SFENTH, &
383 tzot, & ! tzot , zot for thermal z0 Wang
384 ids,ide, jds,jde, kds,kde, &
385 ims,ime, jms,jme, kms,kme, &
386 its,ite, jts,jte, kts,kte )
387 !!! if(j == (jts+jte)/2) write(0,*)'after call mflux,w10',wind10( (its+ite)/2 )
389 ! write(0,*)'--------------------------------------------'
390 ! write(0,*) 'fxh, fxe, fxmx, fxmy, cdm, xxfh zoc,tstrc'
391 ! write(0,*)'--------------------------------------------'
395 ! WRITE(0,1010)i,j,fxh(i),fxe(i),fxmx(i),fxmy(i),cdm(i),rib(i),xxfh(i),zoc(i),tstrc(i)
398 1010 format(2I4,9F11.6)
402 !GFDL CALL PROGTM(IM,KM,PS,U1,V1,T1,Q1, &
403 !GFDL SHELEG,TSKIN,QSURF, &
404 !WRF SMC,STC,DM,SOILTYP,SIGMAF,VEGTYPE,CANOPY,DLWFLX, &
405 !WRF SLRAD,SNOWMT,DELT, &
407 !WRF TG3,GFLUX,F10M, &
408 !GFDL U10M,V10M,T2M,Q2M, &
411 !WRF RHSCNPY,RHSMC,AIM,BIM,CIM, &
412 !GFDL RCL,PRSL1,PRSLKI,SLIMSK, &
413 !GFDL DRAIN,EVAP,HFLX,STRESS,EP, &
414 !GFDL FM,FH,USTAR,WIND,DDVEL, &
415 !GFDL PM,PH,FH2,QSS,Z1 )
419 ! update skin temp only when using GFDL slab...
422 tsk(i,j) = tstrc(i) ! gopal's doing
424 ! bob's doing patch tsk with neigboring values... are grid boundaries
426 tsk(i,j) = tsk(i,j-1)
430 tsk(i,j) = tsk(i,j+1)
433 if(i.eq.ide) tsk(i,j) = tsk(i-1,j)
434 if(i.eq.ids) tsk(i,j) = tsk(i+1,j)
438 znt(i,j)= 0.01*abs(zoc(i))
440 mznt(i,j)= 0.01*abs(mzoc(i))
442 wspd(i,j) = SQRT(upc(kts,i)*upc(kts,i) + vpc(kts,i)*vpc(kts,i))
443 wspd(i,j) = amax1(wspd(i,j) ,100.)/100.
444 u10m(i) = u1(i)*(wind10(i)/wspd(i,j))/100.
445 v10m(i) = v1(i)*(wind10(i)/wspd(i,j))/100.
446 ! br =0.0001*zfull(i,kmax)*dthv/
447 ! & (gmks*theta(i,kmax)*wspd *wspd )
448 ! zkmax = rgas*tpc(kmax,i)*qqlog(kmax)*og
449 zkmax(i,j) = -R*tpc(kts,i)*alog(pkmax(i)/pspc(i))/grav
450 !------------------------------------------------------------------------
452 gz1oz0(i,j)=alog(zkmax(i,j)/znt(i,j))
453 ustar (i)= 0.01*sqrt(cdm(i)* &
454 (upc(kts,i)*upc(kts,i) + vpc(kts,i)*vpc(kts,i)))
455 ! convert from g/(cm*cm*sec) to kg/(m*m*sec)
456 qfx (i,j)=-10.*fxe(i) ! BOB: qfx (i,1)=-10.*fxe(i)
458 ! convert from ergs/gram/K to J/kg/K cpmks=1004
459 ! hfx (i,1)=-0.001*cpcgs*fxh(i)
460 hfx (i,j)= -10.*CP*fxh(i) ! Bob: hfx (i,1)= -10.*CP*fxh(i)
461 taux (i,j)= fxmx(i)/10. ! gopal's doing for Ocean coupling
462 tauy (i,j)= fxmy(i)/10. ! gopal's doing for Ocean coupling
463 fm(i) = karman/sqrt(cdm(i))
464 fh(i) = karman*xxfh(i)
465 PSIM(i,j)=GZ1OZ0(i,j)-FM(i)
466 PSIH(i,j)=GZ1OZ0(i,j)-FH(i)
467 fh2(i) = karman*xxfh2(i)
468 ch(i) = karman*karman/(fm(i) * fh(i))
473 !!! convert cd, ch to values at 10m, for output
474 if ( wind10(i) .ge. 0.1 ) then
475 cd_out(i,j)=cm(i)* (wspd(i,j)/(0.01*wind10(i)) )**2
476 tmp9=0.01*abs(tzot(i))
477 ch_out(i,j)=ch(i)*(wspd(i,j)/(0.01*wind10(i)) ) * &
478 (alog(zkmax(i,j)/tmp9)/alog(10.0/tmp9))
479 ! if(j == (jts+jte)/2 .and. i == (its+ite)/2 ) then
480 ! write(0,*)'cm,cd10=',cm(i),cd_out(i,j)
481 ! write(0,*)'ch,ch10=',ch(i),ch_out(i,j)
482 ! write(0,*)'w/w10,zkmax(i),tzot=',wspd(i,j)/(0.01*wind10(i)),zkmax(i,j),tmp9
489 CHS(I,J)=CH(I)*wspd (i,j)
490 CHS2(I,J)=USTAR(I)*KARMAN/FH2(I)
492 !!!2014-0922 cap CHS over land points
493 if (xland(i,j) .lt. 1.9) then
494 CHS(I,J)=amin1(CHS(I,J), 0.05)
495 CHS2(I,J)=amin1(CHS2(I,J), 0.05)
496 if (chs2(i,j) < 0) chs2(i,j)=1.0e-6
501 CPM(I,J)=CP*(1.+0.8*QV3D(i,kts,j))
503 QGH(I,J)=ep2*esat/(1000.*ps(i)-esat)
504 esat = fpvs(tskin(i))
505 qss(i) = ep2*esat/(1000.*ps(i)-esat)
510 ! wspd (i,j) = SQRT(upc(kts,i)*upc(kts,i) + vpc(kts,i)*vpc(kts,i))
511 ! wspd (i,j) = amax1(wspd (i,j) ,100.)/100.
513 ! ZNT(i,j)=Z0RL(i)*.01
516 ! write(0,*)'fm,fh,cm,ch(125)', fm(125),fh(125),cm(125),ch(125)
519 FLHC(i,j)=CPM(I,J)*RHO1(I)*CHS(I,J)
520 FLQC(i,j)=RHO1(I)*CHS(I,J)
521 ! GZ1OZ0(i,j)=LOG(Z1(I)/(Z0RL(I)*.01))
525 IF (ISFFLX.EQ.0) THEN
533 IF(XLAND(I,J)-1.5.GT.0.)THEN
534 ! HFX(I,J)=FLHC(I,J)*(THGB(I)-THX(I))
536 ! convert from ergs/gram/K to J/kg/K cpmks=1004
537 ! hfx (i,j)=-0.001*cpcgs*fxh(i)
538 hfx (i,j)= -10.*CP*fxh(i) ! Bob: hfx (i,1)= -10.*CP*fxh(i)
539 ELSEIF(XLAND(I,J)-1.5.LT.0.)THEN
540 ! HFX(I,J)=FLHC(I,J)*(THGB(I)-THX(I))
542 ! convert from ergs/gram/K to J/kg/K cpmks=1004
543 ! hfx (i,j)=-0.001*cpcgs*fxh(i)
544 hfx (i,j)= -10.*CP*fxh(i) ! Bob: hfx (i,j)= -10.*CP*fxh(i)
545 HFX(I,J)=AMAX1(HFX(I,J),-250.)
547 ! QFX(I,J)=FLQC(I,J)*(QSFC(I,J)-Q1(I))
548 ! convert from g/(cm*cm*sec) to kg/(m*m*sec)
550 QFX(I,J)=AMAX1(QFX(I,J),0.)
554 ! if(j.eq.2) write(0,*) 'u3d ,ustar,cdm at end of gfdlsfcmod'
555 ! write(0,*) j,(u3d(ii,1,j),ii=70,90)
556 ! write(0,*) j,(ustar(ii),ii=70,90)
557 ! write(0,*) j,(cdm(ii),ii=70,90)
558 if(j.eq.jds.or.j.eq.jde) then
560 write(message,*) "TSFC in gfdl sf mod,dt, its,ite,jts,jts", dt,its,ite,jts,jte,ids,ide,jds,jde
561 call wrf_debug(1,message)
562 write(message,*) "TSFC", (TSK(i,j),i=its,ite)
563 call wrf_debug(1,message)
566 ENDDO ! FOR THE J LOOP I PRESUME
567 ! if(100.ge.its.and.100.le.ite.and.100.ge.jts.and.100.le.jte) then
568 ! write(0,*) 'output vars of sf_gfdl at i,j=100'
569 ! write(0,*) 'TSK',TSK(100,100)
570 ! write(0,*) 'PSFC',PSFC(100,100)
571 ! write(0,*) 'GLW',GLW(100,100)
572 ! write(0,*) 'GSW',GSW(100,100)
573 ! write(0,*) 'XLAND',XLAND(100,100)
574 ! write(0,*) 'BR',BR(100,100)
575 ! write(0,*) 'CHS',CHS(100,100)
576 ! write(0,*) 'CHS2',CHS2(100,100)
577 ! write(0,*) 'CPM',CPM(100,100)
578 ! write(0,*) 'FLHC',FLHC(100,100)
579 ! write(0,*) 'FLQC',FLQC(100,100)
580 ! write(0,*) 'GZ1OZ0',GZ1OZ0(100,100)
581 ! write(0,*) 'HFX',HFX(100,100)
582 ! write(0,*) 'LH',LH(100,100)
583 ! write(0,*) 'PSIM',PSIM(100,100)
584 ! write(0,*) 'PSIH',PSIH(100,100)
585 ! write(0,*) 'QFX',QFX(100,100)
586 ! write(0,*) 'QGH',QGH(100,100)
587 ! write(0,*) 'QSFC',QSFC(100,100)
588 ! write(0,*) 'UST',UST(100,100)
589 ! write(0,*) 'ZNT',ZNT(100,100)
590 ! write(0,*) 'wet',wet(100)
591 ! write(0,*) 'smois',smois(100,1,100)
592 ! write(0,*) 'WSPD',WSPD(100,100)
593 ! write(0,*) 'U10',U10(100,100)
594 ! write(0,*) 'V10',V10(100,100)
598 END SUBROUTINE SF_GFDL
600 !-------------------------------------------------------------------
601 SUBROUTINE MFLUX2( fxh,fxe,fxmx,fxmy,cdm,rib,xxfh,zoc,mzoc,tstrc, & !mzoc KWON
602 pspc,pkmax,wetc,slwdc,tjloc, &
603 icoef_sf,iwavecpl,lcurr_sf,alpha,gamma,xcur,ycur, &
604 pert_Cd, ens_random_seed, ens_Cdamp, &
605 upc,vpc,tpc,rpc,dt,jfix,wind10,xxfh2,ntsflg,sfenth, &
607 ids,ide, jds,jde, kds,kde, &
608 ims,ime, jms,jme, kms,kme, &
609 its,ite, jts,jte, kts,kte )
611 !------------------------------------------------------------------------
613 ! MFLUX2 computes surface fluxes of momentum, heat,and moisture
614 ! using monin-obukhov. the roughness length "z0" is prescribed
615 ! over land and over ocean "z0" is computed using charnocks formula.
616 ! the universal functions (from similarity theory approach) are
617 ! those of hicks. This is Bob's doing.
619 !------------------------------------------------------------------------
621 USE module_sf_exchcoef
624 ! use module_TLDATA , ONLY : tab,table,cp,g,rgas,og
626 ! include 'RESOLUTION.h'
627 ! include 'PARAMETERS.h'
628 ! include 'STDUNITS.h' stdout
631 ! include 'BKINFO.h' nstep
632 ! include 'CONSLEV.h'
633 ! include 'CONMLEV.h'
637 ! include 'GDINFO.h' ngd
640 ! include 'TIME.h' dt(nnst)
642 ! include 'ZLDATA.h' old MOBFLX?
644 !-----------------------------------------------------------------------
645 ! user interface variables
646 !-----------------------------------------------------------------------
647 integer,intent(in) :: ids,ide, jds,jde, kds,kde
648 integer,intent(in) :: ims,ime, jms,jme, kms,kme
649 integer,intent(in) :: its,ite, jts,jte, kts,kte
650 integer,intent(in) :: jfix,ntsflg
651 integer,intent(in) :: icoef_sf
652 integer,intent(in) :: iwavecpl
653 logical,intent(in) :: lcurr_sf
654 logical,intent(in) :: pert_Cd
655 integer,intent(in) :: ens_random_seed
656 real,intent(in) :: ens_Cdamp
658 real, intent (out), dimension (ims :ime ) :: fxh
659 real, intent (out), dimension (ims :ime ) :: fxe
660 real, intent (out), dimension (ims :ime ) :: fxmx
661 real, intent (out), dimension (ims :ime ) :: fxmy
662 real, intent (inout), dimension (ims :ime ) :: cdm
663 ! real, intent (out), dimension (ims :ime ) :: cdm2
664 real, intent (out), dimension (ims :ime ) :: rib
665 real, intent (out), dimension (ims :ime ) :: xxfh
666 real, intent (out), dimension (ims :ime ) :: xxfh2
667 real, intent (out), dimension (ims :ime ) :: wind10
669 real, intent ( inout), dimension (ims :ime ) :: zoc,mzoc !KWON
670 real, intent ( inout), dimension (ims :ime ) :: tzot !WANG
671 real, intent ( inout), dimension (ims :ime ) :: tstrc
673 real, intent ( in) :: dt
674 real, intent ( in) :: sfenth
675 real, intent ( in), dimension (ims :ime ) :: pspc
676 real, intent ( in), dimension (ims :ime ) :: pkmax
677 real, intent ( in), dimension (ims :ime ) :: wetc
678 real, intent ( in), dimension (ims :ime ) :: slwdc
679 real, intent ( in), dimension (ims :ime ) :: tjloc
680 real, intent ( in), dimension (ims :ime ) :: alpha, gamma
681 real, intent ( in), dimension (ims :ime ) :: xcur, ycur
683 real, intent ( in), dimension (kms:kme, ims :ime ) :: upc
684 real, intent ( in), dimension (kms:kme, ims :ime ) :: vpc
685 real, intent ( in), dimension (kms:kme, ims :ime ) :: tpc
686 real, intent ( in), dimension (kms:kme, ims :ime ) :: rpc
688 !-----------------------------------------------------------------------
690 !-----------------------------------------------------------------------
692 integer, parameter :: icntx = 30
694 integer, dimension(1 :ime) :: ifz
695 integer, dimension(1 :ime) :: indx
696 integer, dimension(1 :ime) :: istb
697 integer, dimension(1 :ime) :: it
698 integer, dimension(1 :ime) :: iutb
700 real, dimension(1 :ime) :: aap
701 real, dimension(1 :ime) :: bq1
702 real, dimension(1 :ime) :: bq1p
703 real, dimension(1 :ime) :: delsrad
704 real, dimension(1 :ime) :: ecof
705 real, dimension(1 :ime) :: ecofp
706 real, dimension(1 :ime) :: estso
707 real, dimension(1 :ime) :: estsop
708 real, dimension(1 :ime) :: fmz1
709 real, dimension(1 :ime) :: fmz10
710 real, dimension(1 :ime) :: fmz2
711 real, dimension(1 :ime) :: fmzo1
712 real, dimension(1 :ime) :: foft
713 real, dimension(1 :ime) :: foftm
714 real, dimension(1 :ime) :: frac
715 real, dimension(1 :ime) :: land
716 real, dimension(1 :ime) :: pssp
717 real, dimension(1 :ime) :: qf
718 real, dimension(1 :ime) :: rdiff
719 real, dimension(1 :ime) :: rho
720 real, dimension(1 :ime) :: rkmaxp
721 real, dimension(1 :ime) :: rstso
722 real, dimension(1 :ime) :: rstsop
723 real, dimension(1 :ime) :: sf10
724 real, dimension(1 :ime) :: sf2
725 real, dimension(1 :ime) :: sfm
726 real, dimension(1 :ime) :: sfzo
727 real, dimension(1 :ime) :: sgzm
728 real, dimension(1 :ime) :: slwa
729 real, dimension(1 :ime) :: szeta
730 real, dimension(1 :ime) :: szetam
731 real, dimension(1 :ime) :: t1
732 real, dimension(1 :ime) :: t2
733 real, dimension(1 :ime) :: tab1
734 real, dimension(1 :ime) :: tab2
735 real, dimension(1 :ime) :: tempa1
736 real, dimension(1 :ime) :: tempa2
737 real, dimension(1 :ime) :: theta
738 real, dimension(1 :ime) :: thetap
739 real, dimension(1 :ime) :: tsg
740 real, dimension(1 :ime) :: tsm
741 real, dimension(1 :ime) :: tsp
742 real, dimension(1 :ime) :: tss
743 real, dimension(1 :ime) :: ucom
744 real, dimension(1 :ime) :: uf10
745 real, dimension(1 :ime) :: uf2
746 real, dimension(1 :ime) :: ufh
747 real, dimension(1 :ime) :: ufm
748 real, dimension(1 :ime) :: ufzo
749 real, dimension(1 :ime) :: ugzm
750 real, dimension(1 :ime) :: uzeta
751 real, dimension(1 :ime) :: uzetam
752 real, dimension(1 :ime) :: vcom
753 real, dimension(1 :ime) :: vrtkx
754 real, dimension(1 :ime) :: vrts
755 real, dimension(1 :ime) :: wind
756 real, dimension(1 :ime) :: windp
757 real, dimension(1 :ime) :: wind10p !WANG, 10m wind previous step
758 real, dimension(1 :ime) :: uvs1
759 ! real, dimension(1 :ime) :: xxfh
760 real, dimension(1 :ime) :: xxfm
761 real, dimension(1 :ime) :: xxsh
762 real, dimension(1 :ime) :: z10
763 real, dimension(1 :ime) :: z2
764 real, dimension(1 :ime) :: zeta
765 real, dimension(1 :ime) :: zkmax
767 real, dimension(1 :ime) :: pss
768 real, dimension(1 :ime) :: tstar
769 real, dimension(1 :ime) :: ukmax
770 real, dimension(1 :ime) :: vkmax
771 real, dimension(1 :ime) :: tkmax
772 real, dimension(1 :ime) :: rkmax
773 real, dimension(1 :ime) :: zot
774 real, dimension(1 :ime) :: fhzo1
775 real, dimension(1 :ime) :: sfh
777 real :: ux13, yo, y,xo,x,ux21,ugzzo,ux11,ux12,uzetao,xnum,alll
778 real :: ux1,ugz,x10,uzo,uq,ux2,ux3,xtan,xden,y10,uzet1o,ugz10
779 real :: szet2, zal2,ugz2
780 real :: rovcp,boycon,cmo2,psps1,zog,enrca,rca,cmo1,amask,en,ca,a,c
781 real :: sgz,zal10,szet10,fmz,szo,sq,fmzo,rzeta1,zal1g,szetao,rzeta2,zal2g
782 real :: hcap,xks,pith,teps,diffot,delten,alevp,psps2,alfus,nstep
783 real :: shfx,sigt4,reflect
784 real :: cor1,cor2,szetho,zal2gh,cons_p000001,cons_7,vis,ustar,restar,rat
786 real :: windmks,znott,znotm
788 integer:: i,j,ii,iq,nnest,icnt,ngd,ip
790 !-----------------------------------------------------------------------
792 !-----------------------------------------------------------------------
794 real, dimension (223) :: tab
795 real, dimension (223) :: table
796 real, dimension (101) :: tab11
797 real, dimension (41) :: table4
798 real, dimension (42) :: tab3
799 real, dimension (54) :: table2
800 real, dimension (54) :: table3
801 real, dimension (74) :: table1
802 real, dimension (80) :: tab22
804 character(len=255) :: message
806 equivalence (tab(1),tab11(1))
807 equivalence (tab(102),tab22(1))
808 equivalence (tab(182),tab3(1))
809 equivalence (table(1),table1(1))
810 equivalence (table(75),table2(1))
811 equivalence (table(129),table3(1))
812 equivalence (table(183),table4(1))
815 !-----------------------------------------------------------------------
816 ! tables used to obtain the vapor pressures or saturated vapor
818 !-----------------------------------------------------------------------
820 data tab11/21*0.01403,0.01719,0.02101,0.02561,0.03117,0.03784, &
821 &.04584,.05542,.06685,.08049,.09672,.1160,.1388,.1658,.1977,.2353, &
822 &.2796,.3316,.3925,.4638,.5472,.6444,.7577,.8894,1.042,1.220,1.425, &
823 &1.662,1.936,2.252,2.615,3.032,3.511,4.060,4.688,5.406,6.225,7.159, &
824 &8.223,9.432,10.80,12.36,14.13,16.12,18.38,20.92,23.80,27.03,30.67, &
825 &34.76,39.35,44.49,50.26,56.71,63.93,71.98,80.97,90.98,102.1,114.5, &
826 &128.3,143.6,160.6,179.4,200.2,223.3,248.8,276.9,307.9,342.1,379.8, &
827 &421.3,466.9,517.0,572.0,632.3,698.5,770.9,850.2,937.0,1032./
829 data tab22/1146.6,1272.0,1408.1,1556.7,1716.9,1890.3,2077.6,2279.6 &
830 &,2496.7,2729.8,2980.0,3247.8,3534.1,3839.8,4164.8,4510.5,4876.9, &
831 &5265.1,5675.2,6107.8,6566.2,7054.7,7575.3,8129.4,8719.2,9346.5, &
832 &10013.,10722.,11474.,12272.,13119.,14017.,14969.,15977.,17044., &
833 &18173.,19367.,20630.,21964.,23373.,24861.,26430.,28086.,29831., &
834 &31671.,33608.,35649.,37796.,40055.,42430.,44927.,47551.,50307., &
835 &53200.,56236.,59422.,62762.,66264.,69934.,73777.,77802.,82015., &
836 &86423.,91034.,95855.,100890.,106160.,111660.,117400.,123400., &
837 &129650.,136170.,142980.,150070.,157460.,165160.,173180.,181530., &
840 data tab3/208670.,218450.,228610.,239180.,250160.,261560.,273400., &
841 &285700.,298450.,311690.,325420.,339650.,354410.,369710.,385560., &
842 &401980.,418980.,436590.,454810.,473670.,493170.,513350.,534220., &
843 &555800.,578090.,601130.,624940.,649530.,674920.,701130.,728190., &
844 &756110.,784920.,814630.,845280.,876880.,909450.,943020.,977610., &
845 &1013250.,1049940.,1087740./
847 data table1/20*0.0,.3160e-02,.3820e-02,.4600e-02,.5560e-02,.6670e-02, &
848 & .8000e-02,.9580e-02,.1143e-01,.1364e-01,.1623e-01,.1928e-01, &
849 &.2280e-01,.2700e-01,.3190e-01,.3760e-01,.4430e-01,.5200e-01, &
850 &.6090e-01,.7130e-01,.8340e-01,.9720e-01,.1133e+00,.1317e-00, &
851 &.1526e-00,.1780e-00,.2050e-00,.2370e-00,.2740e-00,.3160e-00, &
852 &.3630e-00,.4170e-00,.4790e-00,.5490e-00,.6280e-00,.7180e-00, &
853 &.8190e-00,.9340e-00,.1064e+01,.1209e+01,.1368e+01,.1560e+01, &
854 &.1770e+01,.1990e+01,.2260e+01,.2540e+01,.2880e+01,.3230e+01, &
855 &.3640e+01,.4090e+01,.4590e+01,.5140e+01,.5770e+01,.6450e+01, &
858 data table2/.8050e+01,.8990e+01,.1001e+02,.1112e+02,.1240e+02, &
859 &.1380e+02,.1530e+02,.1700e+02,.1880e+02,.2080e+02,.2310e+02, &
860 &.2550e+02,.2810e+02,.3100e+02,.3420e+02,.3770e+02,.4150e+02, &
861 &.4560e+02,.5010e+02,.5500e+02,.6030e+02,.6620e+02,.7240e+02, &
862 &.7930e+02,.8680e+02,.9500e+02,.1146e+03,.1254e+03,.1361e+03, &
863 &.1486e+03,.1602e+03,.1734e+03,.1873e+03,.2020e+03,.2171e+03, &
864 &.2331e+03,.2502e+03,.2678e+03,.2863e+03,.3057e+03,.3250e+03, &
865 &.3457e+03,.3664e+03,.3882e+03,.4101e+03,.4326e+03,.4584e+03, &
866 &.4885e+03,.5206e+03,.5541e+03,.5898e+03,.6273e+03,.6665e+03, &
869 data table3/.7520e+03,.7980e+03,.8470e+03,.8980e+03,.9520e+03, &
870 &.1008e+04,.1067e+04,.1129e+04,.1194e+04,.1263e+04,.1334e+04, &
871 &.1409e+04,.1488e+04,.1569e+04,.1656e+04,.1745e+04,.1840e+04, &
872 &.1937e+04,.2041e+04,.2147e+04,.2259e+04,.2375e+04,.2497e+04, &
873 &.2624e+04,.2756e+04,.2893e+04,.3036e+04,.3186e+04,.3340e+04, &
874 &.3502e+04,.3670e+04,.3843e+04,.4025e+04,.4213e+04,.4408e+04, &
875 &.4611e+04,.4821e+04,.5035e+04,.5270e+04,.5500e+04,.5740e+04, &
876 &.6000e+04,.6250e+04,.6520e+04,.6810e+04,.7090e+04,.7390e+04, &
877 &.7700e+04,.8020e+04,.8350e+04,.8690e+04,.9040e+04,.9410e+04, &
880 data table4/.1016e+05,.1057e+05,.1098e+05,.1140e+05,.1184e+05, &
881 &.1230e+05,.1275e+05,.1324e+05,.1373e+05,.1423e+05,.1476e+05, &
882 &.1530e+05,.1585e+05,.1642e+05,.1700e+05,.1761e+05,.1822e+05, &
883 &.1886e+05,.1950e+05,.2018e+05,.2087e+05,.2158e+05,.2229e+05, &
884 &.2304e+05,.2381e+05,.2459e+05,.2539e+05,.2621e+05,.2706e+05, &
885 &.2792e+05,.2881e+05,.2971e+05,.3065e+05,.3160e+05,.3257e+05, &
886 &.3357e+05,.3459e+05,.3564e+05,.3669e+05,.3780e+05,.0000e+00/
888 ! spcify constants needed by MFLUX2
890 real,parameter :: cp = 1.00464e7
891 real,parameter :: g = 980.6
892 real,parameter :: rgas = 2.87e6
893 real,parameter :: og = 1./g
894 integer :: ntstep = 0
897 real*8 :: gasdev,ran1 !zhang
899 logical,save :: pert_Cd_local !zhang
900 CHARACTER(len=3) :: env_memb,env_pp
901 integer,save :: ens_random_seed_local,env_pp_local !zhang
902 integer :: ensda_physics_pert !zhang
903 real,save :: ens_Cdamp_local !zhang
904 data ens_random_seed_local/0/
906 if ( ens_random_seed_local .eq. 0 ) then
907 CALL nl_get_ensda_physics_pert(1,ensda_physics_pert)
908 ens_random_seed_local=ens_random_seed
909 env_pp_local=ensda_physics_pert
910 pert_Cd_local=.false.
912 ! env_pp=1: do physics perturbations for ensda members, ens_random_seed must be 99
913 if ( env_pp_local .eq. 1 ) then
914 if ( ens_random_seed .ne. 99 ) then
916 ens_Cdamp_local=ens_Cdamp
918 ! ens_random_seed=99 do physics perturbation for ensemble forecasts, env_pp must be zero
919 ens_random_seed_local=ens_random_seed
920 pert_Cd_local=pert_Cd
921 ens_Cdamp_local=ens_Cdamp
924 ens_random_seed_local=ens_random_seed
925 pert_Cd_local=pert_Cd
926 ens_Cdamp_local=ens_Cdamp
928 print*, "Cd ===", ens_random_seed_local,pert_Cd_local,ens_Cdamp_local,ensda_physics_pert
932 ! character*10 routine
935 !------------------------------------------------------------------------
936 ! set water availability constant "ecof" and land mask "land".
937 ! limit minimum wind speed to 100 cm/s
938 !------------------------------------------------------------------------
940 ! constants for 10 m winds (correction for knots
944 ! KWON : remove the artificial increase of 10m wind speed over 60kts
945 ! which comes from GFDL hurricane model
956 if ( lcurr_sf .and. zoc(i) .le. 0.0 ) then
957 ubot = upc(1,i) - xcur(i) * 100.0
958 vbot = vpc(1,i) - ycur(i) * 100.0
965 uvs1(i)= amax1( SQRT(ubot*ubot + &
967 if ( iwavecpl .eq. 1 .and. zoc(i) .le. 0.0 ) then
968 ukmax(i) = ( ubot * cos(gamma(i)) - &
969 vbot * sin(gamma(i)) ) &
971 vkmax(i) = ( vbot * cos(gamma(i)) - &
972 ubot * sin(gamma(i)) ) &
980 ! ukmax(i) = upc(1,i)
981 ! vkmax(i) = vpc(1,i)
986 ! write(0,*)'z10,pss,tstar,u...rkmax(ite)', &
987 ! z10(ite), pss(ite),tstar(ite),ukmax(ite), &
988 ! vkmax(ite),tkmax(ite),rkmax(ite)
991 windp(i) = SQRT(ukmax(i)*ukmax(i) + vkmax(i)*vkmax(i))
992 wind (i) = amax1(windp(i),100.)
994 !! use wind10 previous step
995 wind10p(i) = wind10(i) !! cm/s
996 wind10p(i) = amax1(wind10p(i),100.)
999 ! if (zoc(i) .LT. amask) zoc(i) = -0.0185*0.001*wind(i)*wind(i)*og
1000 if (zoc(i) .LT. amask) zoc(i) = -0.0185*0.001*wind10p(i)*wind10p(i)*og
1001 if (zoc(i) .GT. 0.0) then
1008 !windmks=wind(i)*.01
1009 windmks=wind10p(i)*.01
1010 if ( iwavecpl .eq. 1 ) then
1011 call znot_wind10m(windmks,znott,znotm,icoef_sf)
1012 !Check if Charnock parameter ratio is received in a proper range.
1013 if ( alpha(i) .ge. 0.2 .and. alpha(i) .le. 5. ) then
1014 znotm = znotm*alpha(i)
1016 zoc(i) = -100.*znotm
1017 zot(i) = -100* znott
1019 ! if ( icoef_sf .EQ. 1) then
1020 ! call znot_m_v1(windmks,znotm)
1021 ! call znot_t_v1(windmks,znott)
1022 ! else if ( icoef_sf .EQ. 0 ) then
1023 ! call znot_m_v0(windmks,znotm)
1024 ! call znot_t_v0(windmks,znott)
1026 ! call znot_m_v1(windmks,znotm)
1027 ! call znot_t_v2(windmks,znott)
1029 call znot_wind10m(windmks,znott,znotm,icoef_sf)
1030 zoc(i) = -100.*znotm
1031 zot(i) = -100* znott
1035 ! if(i == (its+ite)/2 .and. j == (jts+jte)/2) then
1036 ! write(0,*)'wind10p(i),windp(i),windmks,zoc(i),zot, land'
1037 ! write(0,*)wind10p(i),windp(i),windmks,zoc(i),zot(i),land(i)
1041 !------------------------------------------------------------------------
1042 ! where necessary modify zo values over ocean.
1043 !------------------------------------------------------------------------
1045 mzoc(i) = zoc(i) !FOR SAVE MOMENTUM Zo
1046 tzot(i) = zot(i) !output wang
1049 !------------------------------------------------------------------------
1051 ! a and c = constants used in evaluating universal function for
1053 ! ca = karmen constant
1054 ! cm01 = constant part of vertical integral of universal
1055 ! function; stable case ( 0.5 < zeta < or = 10.0)
1056 ! cm02 = constant part of vertical integral of universal
1057 ! function; stable case ( zeta > 10.0)
1058 !------------------------------------------------------------------------
1065 cmo2 = 17.193 + .5*a - 10.*c
1068 ! write(0,*)'rgas,cp,rovcp ', rgas,cp,rovcp
1070 ! write(0,*)'--------------------------------------------------'
1071 ! write(0,*)'pkmax, pspc, theta, zkmax, zoc'
1072 ! write(0,*)'--------------------------------------------------'
1075 ! theta(i) = tkmax(i)*rqc9
1076 theta(i) = tkmax(i)/((pkmax(i)/pspc(i))**rovcp)
1077 vrtkx(i) = 1.0 + boycon*rkmax(i)
1078 ! zkmax(i) = rgas*tkmax(i)*qqlog(kmax)*og
1079 zkmax(i) = -rgas*tkmax(i)*alog(pkmax(i)/pspc(i))*og
1080 ! IF(I==78)write(0,*)I,JFIX,pkmax(i),pspc(i),theta(i),zkmax(i),zoc(i)
1083 ! write(0,*)'pkmax,pspc ', pkmax,pspc
1084 ! write(0,*)'theta, zkmax, zoc ', theta, zkmax, zoc
1086 !------------------------------------------------------------------------
1087 ! get saturation mixing ratios at surface
1088 !------------------------------------------------------------------------
1092 tab1 (i) = tstar(i) - 153.16
1093 it (i) = IFIX(tab1(i))
1094 tab2 (i) = tab1(i) - FLOAT(it(i))
1095 t1 (i) = tab(min(223,max(1,it(i) + 1)))
1096 t2 (i) = table(min(223,max(1,it(i) + 1)))
1097 estso(i) = t1(i) + tab2(i)*t2(i)
1098 psps1 = (pss(i) - estso(i))
1099 if(psps1 .EQ. 0.0)then
1102 rstso(i) = 0.622*estso(i)/psps1
1103 vrts (i) = 1. + boycon*ecof(i)*rstso(i)
1106 !------------------------------------------------------------------------
1107 ! check if consideration of virtual temperature changes stability.
1108 ! if so, set "dthetav" to near neutral value (1.0e-4). also check
1109 ! for very small lapse rates; if ABS(tempa1) <1.0e-4 then
1111 !------------------------------------------------------------------------
1114 tempa1(i) = theta(i)*vrtkx(i) - tstar(i)*vrts(i)
1115 tempa2(i) = tempa1(i)*(theta(i) - tstar(i))
1116 if (tempa2(i) .LT. 0.) tempa1(i) = 1.0e-4
1117 tab1(i) = ABS(tempa1(i))
1118 if (tab1(i) .LT. 1.0e-4) tempa1(i) = 1.0e-4
1119 !------------------------------------------------------------------------
1120 ! compute bulk richardson number "rib" at each point. if "rib"
1121 ! exceeds 95% of critical richardson number "tab1" then "rib = tab1"
1122 !------------------------------------------------------------------------
1124 rib (i) = g*zkmax(i)*tempa1(i)/ &
1125 (tkmax(i)*vrtkx(i)*wind(i)*wind(i))
1126 tab2(i) = ABS(zoc(i))
1127 tab1(i) = 0.95/(c*(1. - tab2(i)/zkmax(i)))
1128 if (rib(i) .GT. tab1(i)) rib(i) = tab1(i)
1132 zeta(i) = ca*rib(i)/0.03
1135 ! write(0,*)'rib,zeta,vrtkx,vrts(ite) ', rib(ite),zeta(ite), &
1136 ! vrtkx(ite),vrts(ite)
1137 !------------------------------------------------------------------------
1138 ! begin looping through points on line, solving wegsteins iteration
1139 ! for zeta at each point, and using hicks functions
1140 !------------------------------------------------------------------------
1142 !------------------------------------------------------------------------
1143 ! set initial guess of zeta=non - dimensional height "szeta" for
1145 !------------------------------------------------------------------------
1149 ! turn off interfacial layer by zeroing out enrca
1153 !------------------------------------------------------------------------
1155 !------------------------------------------------------------------------
1159 if (zeta(i) .GE. 0.0) then
1165 if (ip .EQ. 0) go to 170
1169 szeta(i) = zeta(istb(i))
1173 !------------------------------------------------------------------------
1174 ! begin wegstein iteration for "zeta" at stable points using
1176 !------------------------------------------------------------------------
1180 if (ifz(i) .EQ. 0) go to 80
1181 zal1g = ALOG(szeta(i))
1182 if (szeta(i) .LE. 0.5) then
1183 fmz1(i) = (zal1g + a*szeta(i))*rca
1184 else if (szeta(i) .GT. 0.5 .AND. szeta(i) .LE. 10.) then
1185 rzeta1 = 1./szeta(i)
1186 fmz1(i) = (8.*zal1g + 4.25*rzeta1 - &
1187 0.5*rzeta1*rzeta1 + cmo1)*rca
1188 else if (szeta(i) .GT. 10.) then
1189 fmz1(i) = (c*szeta(i) + cmo2)*rca
1191 szetao = ABS(zoc(istb(i)))/zkmax(istb(i))*szeta(i)
1192 zal2g = ALOG(szetao)
1193 if (szetao .LE. 0.5) then
1194 fmzo1(i) = (zal2g + a*szetao)*rca
1195 sfzo (i) = 1. + a*szetao
1196 else if (szetao .GT. 0.5 .AND. szetao .LE. 10.) then
1198 fmzo1(i) = (8.*zal2g + 4.25*rzeta2 - &
1199 0.5*rzeta2*rzeta2 + cmo1)*rca
1200 sfzo (i) = 8.0 - 4.25*rzeta2 + rzeta2*rzeta2
1201 else if (szetao .GT. 10.) then
1202 fmzo1(i) = (c*szetao + cmo2)*rca
1207 ! compute heat & moisture parts of zot.. for calculation of sfh
1209 szetho = ABS(zot(istb(i)))/zkmax(istb(i))*szeta(i)
1210 zal2gh = ALOG(szetho)
1211 if (szetho .LE. 0.5) then
1212 fhzo1(i) = (zal2gh + a*szetho)*rca
1213 sfzo (i) = 1. + a*szetho
1214 else if (szetho .GT. 0.5 .AND. szetho .LE. 10.) then
1216 fhzo1(i) = (8.*zal2gh + 4.25*rzeta2 - &
1217 0.5*rzeta2*rzeta2 + cmo1)*rca
1218 sfzo (i) = 8.0 - 4.25*rzeta2 + rzeta2*rzeta2
1219 else if (szetho .GT. 10.) then
1220 fhzo1(i) = (c*szetho + cmo2)*rca
1224 !------------------------------------------------------------------------
1225 ! compute universal function at 10 meters for diagnostic purposes
1226 !------------------------------------------------------------------------
1228 !!!! if (ngd .EQ. nNEST) then
1229 szet10 = ABS(z10(istb(i)))/zkmax(istb(i))*szeta(i)
1230 zal10 = ALOG(szet10)
1231 if (szet10 .LE. 0.5) then
1232 fmz10(i) = (zal10 + a*szet10)*rca
1233 else if (szet10 .GT. 0.5 .AND. szet10 .LE. 10.) then
1235 fmz10(i) = (8.*zal10 + 4.25*rzeta2 - &
1236 0.5*rzeta2*rzeta2 + cmo1)*rca
1237 else if (szet10 .GT. 10.) then
1238 fmz10(i) = (c*szet10 + cmo2)*rca
1240 sf10(i) = fmz10(i) - fmzo1(i)
1241 ! compute 2m values for diagnostics in HWRF
1242 szet2 = ABS(z2 (istb(i)))/zkmax(istb(i))*szeta(i)
1244 if (szet2 .LE. 0.5) then
1245 fmz2 (i) = (zal2 + a*szet2 )*rca
1246 else if (szet2 .GT. 0.5 .AND. szet2 .LE. 2.) then
1248 fmz2 (i) = (8.*zal2 + 4.25*rzeta2 - &
1249 0.5*rzeta2*rzeta2 + cmo1)*rca
1250 else if (szet2 .GT. 2.) then
1251 fmz2 (i) = (c*szet2 + cmo2)*rca
1253 sf2 (i) = fmz2 (i) - fmzo1(i)
1256 sfm(i) = fmz1(i) - fmzo1(i)
1257 sfh(i) = fmz1(i) - fhzo1(i)
1258 sgz = ca*rib(istb(i))*sfm(i)*sfm(i)/ &
1259 (sfh(i) + enrca*sfzo(i))
1260 fmz = (sgz - szeta(i))/szeta(i)
1262 if (fmzo .GE. 5.0e-5) then
1263 sq = (sgz - sgzm(i))/(szeta(i) - szetam(i))
1265 write(message,*)'NCO ERROR DIVIDE BY ZERO IN MFLUX2 (STABLE CASE)'
1266 call wrf_message(message)
1267 write(message,*)'sq is 1 ',fmzo,sgz,sgzm(i),szeta(i),szetam(i)
1268 call wrf_message(message)
1270 szetam(i) = szeta(i)
1271 szeta (i) = (sgz - szeta(i)*sq)/(1.0 - sq)
1281 if (ifz(i) .GE. 1) go to 110
1289 120 format(2X, ' NON-CONVERGENCE FOR STABLE ZETA IN ROW ')
1290 ! call MPI_CLOSE(1,routine)
1292 !------------------------------------------------------------------------
1293 ! update "zo" for ocean points. "zo"cannot be updated within the
1294 ! wegsteins iteration as the scheme (for the near neutral case)
1295 ! can become unstable
1296 !------------------------------------------------------------------------
1301 if (szo .LT. 0.0) then
1302 wndm=wind(istb(i))*0.01
1303 if(wndm.lt.15.0) then
1306 !! ckg=(0.000308*wndm+0.00925)*og
1307 !! ckg=(0.000616*wndm)*og
1308 ckg=(sfenth*(4*0.000308*wndm) + (1.-sfenth)*0.0185 )*og
1311 szo = - ckg*wind(istb(i))*wind(istb(i))/ &
1313 cons_p000001 = .000001
1317 ustar = sqrt( -szo / zog)
1318 restar = -ustar * szo / vis
1319 restar = max(restar,cons_p000001)
1320 ! Rat taken from Zeng, Zhao and Dickinson 1997
1321 rat = 2.67 * restar ** .25 - 2.57
1322 rat = min(rat ,cons_7) !constant
1324 zot(istb(i)) = szo * exp(-rat)
1326 zot(istb(i)) = zoc(istb(i))
1329 ! in hwrf thermal znot is loaded back into the zoc array for next step
1334 xxfm(istb(i)) = sfm(i)
1335 xxfh(istb(i)) = sfh(i)
1336 xxfh2(istb(i)) = sf2 (i)
1337 xxsh(istb(i)) = sfzo(i)
1340 !------------------------------------------------------------------------
1341 ! obtain wind at 10 meters for diagnostic purposes
1342 !------------------------------------------------------------------------
1344 !!! if (ngd .EQ. nNEST) then
1346 wind10(istb(i)) = sf10(i)*uvs1(istb(i))/sfm(i)
1347 wind10(istb(i)) = wind10(istb(i)) * 1.944
1348 if(wind10(istb(i)) .GT. 6000.0) then
1349 wind10(istb(i))=wind10(istb(i))+wind10(istb(i))*cor1 &
1352 ! the above correction done by GFDL in centi-kts!!!-change back
1353 wind10(istb(i)) = wind10(istb(i)) / 1.944
1356 !!! if (ngd .EQ. nNEST-1 .AND. llwe .EQ. 1 ) then
1358 !!! wind10c(istb(i),j) = sf10(i)*wind(istb(i))/sfm(i)
1359 !!! wind10c(istb(i),j) = wind10c(istb(i),j) * 1.944
1360 !!! if(wind10c(istb(i),j) .GT. 6000.0) then
1361 !!! wind10c(istb(i),j)=wind10c(istb(i),j)+wind10c(istb(i),j)*cor1
1367 !------------------------------------------------------------------------
1369 !------------------------------------------------------------------------
1375 if (zeta(i) .LT. 0.0) then
1381 if (iq .EQ. 0) go to 290
1383 uzeta (i) = zeta(iutb(i))
1389 !------------------------------------------------------------------------
1390 ! begin wegstein iteration for "zeta" at unstable points using
1392 !------------------------------------------------------------------------
1396 if (ifz(i) .EQ. 0) go to 200
1397 ugzzo = ALOG(zkmax(iutb(i))/ABS(zot(iutb(i))))
1398 uzetao = ABS(zot(iutb(i)))/zkmax(iutb(i))*uzeta(i)
1399 ux11 = 1. - 16.*uzeta(i)
1400 ux12 = 1. - 16.*uzetao
1404 ux13 = (1. + y)/(1. + yo)
1406 ufh(i) = (ugzzo - 2.*ux21)*rca
1407 ! recompute scalers for ufm in terms of mom znot... zoc
1408 ugzzo = ALOG(zkmax(iutb(i))/ABS(zoc(iutb(i))))
1409 uzetao = ABS(zoc(iutb(i)))/zkmax(iutb(i))*uzeta(i)
1410 ux11 = 1. - 16.*uzeta(i)
1411 ux12 = 1. - 16.*uzetao
1414 ux13 = (1. + y)/(1. + yo)
1419 xnum = (x**2 + 1.)*((x + 1.)**2)
1420 xden = (xo**2 + 1.)*((xo + 1.)**2)
1421 xtan = ATAN(x) - ATAN(xo)
1422 ux3 = ALOG(xnum/xden)
1423 ufm(i) = (ugzzo - ux3 + 2.*xtan)*rca
1424 !!!! if (ngd .EQ. nNEST) then
1426 !------------------------------------------------------------------------
1427 ! obtain ten meter winds for diagnostic purposes
1428 !------------------------------------------------------------------------
1430 ugz10 = ALOG(z10(iutb(i))/ABS(zoc(iutb(i))))
1431 uzet1o = ABS(z10(iutb(i)))/zkmax(iutb(i))*uzeta(i)
1432 uzetao = ABS(zoc(iutb(i)))/zkmax(iutb(i))*uzeta(i)
1433 ux11 = 1. - 16.*uzet1o
1434 ux12 = 1. - 16.*uzetao
1437 ux13 = (1. + y)/(1. + y10)
1441 xnum = (x**2 + 1.)*((x + 1.)**2)
1442 xden = (x10**2 + 1.)*((x10 + 1.)**2)
1443 xtan = ATAN(x) - ATAN(x10)
1444 ux3 = ALOG(xnum/xden)
1445 uf10(i) = (ugz10 - ux3 + 2.*xtan)*rca
1447 ! obtain 2m values for diagnostics...
1450 ugz2 = ALOG(z2 (iutb(i))/ABS(zoc(iutb(i))))
1451 uzet1o = ABS(z2 (iutb(i)))/zkmax(iutb(i))*uzeta(i)
1452 uzetao = ABS(zoc(iutb(i)))/zkmax(iutb(i))*uzeta(i)
1453 ux11 = 1. - 16.*uzet1o
1454 ux12 = 1. - 16.*uzetao
1457 ux13 = (1. + y)/(1. + yo)
1459 uf2 (i) = (ugzzo - 2.*ux21)*rca
1463 ugz = ca*rib(iutb(i))*ufm(i)*ufm(i)/(ufh(i) + enrca*ufzo(i))
1464 ux1 = (ugz - uzeta(i))/uzeta(i)
1466 if (ux2 .GE. 5.0e-5) then
1467 uq = (ugz - ugzm(i))/(uzeta(i) - uzetam(i))
1468 uzetam(i) = uzeta(i)
1470 call wrf_message('NCO ERROR DIVIDE BY ZERO IN MFLUX2 (UNSTABLE CASE)')
1471 write(message,*)'uq is 1 ',ux2,ugz,ugzm(i),uzeta(i),uzetam(i)
1472 call wrf_message(message)
1474 uzeta (i) = (ugz - uzeta(i)*uq)/(1.0 - uq)
1485 if (ifz(i) .GE. 1) go to 230
1492 call wrf_message(message)
1493 240 format(2X, ' NON-CONVERGENCE FOR UNSTABLE ZETA IN ROW ')
1494 ! call MPI_CLOSE(1,routine)
1496 !------------------------------------------------------------------------
1497 ! gather unstable values
1498 !------------------------------------------------------------------------
1502 !------------------------------------------------------------------------
1503 ! update "zo" for ocean points. zo cannot be updated within the
1504 ! wegsteins iteration as the scheme (for the near neutral case)
1505 ! can become unstable.
1506 !------------------------------------------------------------------------
1510 if (zoc(iutb(i)) .LT. 0.0) then
1511 wndm=wind(iutb(i))*0.01
1512 if(wndm.lt.15.0) then
1515 !! ckg=(0.000308*wndm+0.00925)*og <
1516 !! ckg=(0.000616*wndm)*og <
1517 ckg=(4*0.000308*wndm)*og
1518 ckg=(sfenth*(4*0.000308*wndm) + (1.-sfenth)*0.0185 )*og
1520 uzo =-ckg*wind(iutb(i))*wind(iutb(i))/(ufm(i)*ufm(i))
1521 cons_p000001 = .000001
1525 ustar = sqrt( -uzo / zog)
1526 restar = -ustar * uzo / vis
1527 restar = max(restar,cons_p000001)
1528 ! Rat taken from Zeng, Zhao and Dickinson 1997
1529 rat = 2.67 * restar ** .25 - 2.57
1530 rat = min(rat ,cons_7) !constant
1532 zot(iutb(i)) = uzo * exp(-rat)
1534 zot(iutb(i)) = zoc(iutb(i))
1536 ! in hwrf thermal znot is loaded back into the zoc array for next step
1540 !------------------------------------------------------------------------
1541 ! obtain wind at ten meters for diagnostic purposes
1542 !------------------------------------------------------------------------
1544 !!! if (ngd .EQ. nNEST) then
1546 wind10(iutb(i)) = uf10(i)*uvs1(iutb(i))/ufm(i)
1547 wind10(iutb(i)) = wind10(iutb(i)) * 1.944
1548 if(wind10(iutb(i)) .GT. 6000.0) then
1549 wind10(iutb(i))=wind10(iutb(i))+wind10(iutb(i))*cor1 &
1552 ! the above correction done by GFDL in centi-kts!!!-change back
1553 wind10(iutb(i)) = wind10(iutb(i)) / 1.944
1556 !!! if (ngd .EQ. nNEST-1) then
1558 !!! wind10c(iutb(i),j) = uf10(i)*wind(iutb(i))/ufm(i)
1559 !!! wind10c(iutb(i),j) = wind10c(iutb(i),j) * 1.944
1560 !!! if(wind10c(iutb(i),j) .GT. 6000.0) then
1561 !!! wind10c(iutb(i),j)=wind10c(iutb(i),j)+wind10c(iutb(i),j)*cor1
1568 xxfm(iutb(i)) = ufm(i)
1569 xxfh(iutb(i)) = ufh(i)
1570 xxfh2(iutb(i)) = uf2 (i)
1571 xxsh(iutb(i)) = ufzo(i)
1579 if (windp(i) .EQ. 0.0) then
1581 ucom (i) = 100.0/SQRT(2.0)
1582 vcom (i) = 100.0/SQRT(2.0)
1584 rho(i) = pss(i)/(rgas*(tsg(i) + enrca*(theta(i) - &
1585 tsg(i))*xxsh(i)/(xxfh(i) + enrca*xxsh(i))))
1586 bq1(i) = wind(i)*rho(i)/(xxfm(i)*(xxfh(i) + enrca*xxsh(i)))
1589 ! do land sfc temperature prediction if ntsflg=1
1590 ! ntsflg = 1 ! gopal's doing
1592 if (ntsflg .EQ. 0) go to 370
1596 pith = SQRT(4.*ATAN(1.0))
1597 alfus = alll/2.39e-8
1599 ! slwdc... in units of cal/min ????
1600 ! slwa... in units of ergs/sec/cm*2
1602 !------------------------------------------------------------------------
1603 ! pack land and sea ice points
1604 !------------------------------------------------------------------------
1608 if (land(i) .EQ. 1) then
1611 ! slwa is defined as positive down....
1612 slwa (ip) = slwdc(i)/(2.39e-8*60.)
1614 thetap (ip) = theta(i)
1615 rkmaxp (ip) = rkmax(i)
1618 ecofp (ip) = ecof(i)
1619 estsop (ip) = estso(i)
1620 rstsop (ip) = rstso(i)
1622 bq1p (ip) = amax1(bq1p(ip),0.1e-3)
1623 delsrad(ip) = dt *pith/(hcap*SQRT(3600.*24.*xks))
1627 !------------------------------------------------------------------------
1628 ! initialize variables for first pass of iteration
1629 !------------------------------------------------------------------------
1634 rdiff(i) = amin1(0.0,(rkmaxp(i) - rstsop(i)))
1635 !!! if (nstep .EQ. -99 .AND. ngd .GT. 1 .OR. &
1636 !!! nstep .EQ. -99 .AND. ngd .EQ. 1) then
1638 !!! if (j .EQ. 1 .AND. i .EQ. 1) write(6,300)
1639 300 format(2X, ' SURFACE EQUILIBRIUM CALCULATION ')
1641 !! foftm(i) = thetap(i) + 1./(cp*bq1p(i))*(slwa(i) - aap(i)* &
1642 !! tsm(i)**4 + ecofp(i)*alfus*bq1p(i)*rdiff(i))
1645 foftm(i) = tss(i) + delsrad(i)*(slwa(i) - aap(i)*tsm(i)**4 - &
1646 cp*bq1p(i)*(tsm(i) - thetap(i)) + ecofp(i)*alfus*bq1p(i)* &
1652 !------------------------------------------------------------------------
1653 ! do iteration to determine "tstar" at new time level
1654 !------------------------------------------------------------------------
1658 if (ifz(i) .EQ. 0) go to 330
1659 tab1 (i) = tsp(i) - 153.16
1660 it (i) = IFIX(tab1(i))
1661 tab2 (i) = tab1(i) - FLOAT(it(i))
1662 t1 (i) = tab(min(223,max(1,it(i) + 1)))
1663 t2 (i) = table(min(223,max(1,it(i) + 1)))
1664 estsop(i) = t1(i) + tab2(i)*t2(i)
1665 psps2 = (pssp(i) - estsop(i))
1666 if(psps2 .EQ. 0.0)then
1669 rstsop(i) = 0.622*estsop(i)/psps2
1670 rdiff (i) = amin1(0.0,(rkmaxp(i) - rstsop(i)))
1671 !!! if (nstep .EQ. -99 .AND. ngd .GT. 1 .OR. &
1672 !!! nstep .EQ. -99 .AND. ngd .EQ. 1) then
1673 !!! foft(i) = thetap(i) + (1./(cp*bq1p(i)))*(slwa(i) - aap(i)* &
1674 !!! tsp(i)**4 + ecofp(i)*alfus*bq1p(i)*rdiff(i))
1676 foft(i) = tss(i) + delsrad(i)*(slwa(i) - aap(i)*tsp(i)**4 - &
1677 cp*bq1p(i)*(tsp(i) - thetap(i)) + ecofp(i)*alfus*bq1p(i)* &
1680 !!! if (ngd .EQ. 1 .AND. j .EQ. 48 .AND. i .EQ. 19) then
1681 !!! reflect = slwa(i)
1682 !!! sigt4 = -aap(i)*tsp(i)**4
1683 !!! shfx = -cp*bq1p(i)*(tsp(i) - thetap(i))
1684 !!! alevp = ecofp(i)*alfus*bq1p(i)*rdiff(i)
1685 !!! delten = delsrad(i)
1686 !!! diffot = foft(i) - tss(i)
1688 frac(i) = ABS((foft(i) - tsp(i))/tsp(i))
1690 !------------------------------------------------------------------------
1691 ! check for convergence of all points use wegstein iteration
1692 !------------------------------------------------------------------------
1694 if (frac(i) .GE. teps) then
1695 qf (i) = (foft(i) - foftm(i))/(tsp(i) - tsm(i))
1697 tsp (i) = (foft(i) - tsp(i)*qf(i))/(1. - qf(i))
1706 !------------------------------------------------------------------------
1707 ! check for convergence of "t star" prediction
1708 !------------------------------------------------------------------------
1711 if (ifz(i) .EQ. 1) then
1712 write(message, 340) tsp(i), i, j
1713 340 format(2X, ' NON-CONVERGENCE OF T* PREDICTED (T*,I,J) = ', E14.8, &
1715 call wrf_message(message)
1717 write(message,345) indx(i), j, tstar(indx(i)), tsp(i), ip
1718 345 format(2X, ' I, J, OLD T*, NEW T*, NPTS ', 2I5, 2E14.8, I5)
1719 call wrf_message(message)
1721 ! write(6,350) reflect, sigt4, shfx, alevp, delten, diffot
1722 350 format(2X, ' REFLECT, SIGT4, SHFX, ALEVP, DELTEN, DIFFOT ', &
1725 ! call MPI_CLOSE(1,routine)
1734 !------------------------------------------------------------------------
1735 ! compute fluxes and momentum drag coef
1736 !------------------------------------------------------------------------
1741 if ( iwavecpl .eq. 1 .and. zoc(i) .le. 0.0 ) then
1742 windmks = wind10(i) * 0.01
1743 call znot_wind10m(windmks,znott,znotm,icoef_sf)
1744 !Check if Charnock parameter ratio is received in a proper range.
1745 if ( alpha(i) .ge. 0.2 .and. alpha(i) .le. 5. ) then
1746 znotm = znotm*alpha(i)
1748 zoc(i) = -100.*znotm
1749 zot(i) = -100* znott
1752 fxh(i) = bq1(i)*(theta(i) - tsg(i))
1753 fxe(i) = ecof(i)*bq1(i)*(rkmax(i) - rstso(i))
1754 if (fxe(i) .GT. 0.0) fxe(i) = 0.0
1755 fxmx(i) = rho(i)/(xxfm(i)*xxfm(i))*wind(i)*wind(i)*ucom(i)/ &
1757 fxmy(i) = rho(i)/(xxfm(i)*xxfm(i))*wind(i)*wind(i)*vcom(i)/ &
1759 cdm(i) = 1./(xxfm(i)*xxfm(i))
1760 ! print *, 'i,zot,zoc,cdm,cdm2,tsg,wind', &
1761 ! i, zot(i),zoc(i), cdm(i),cdm2(i), tsg(i),wind(i)
1763 ! randomly perturb the Cd
1764 !zzz if( pert_Cd_local .and. ens_random_seed_local .gt. 0 ) then
1765 if( pert_Cd_local ) then
1766 ! print*, "zhang cd ens_random_seed==",ens_random_seed,ens_random_seed_local
1767 ens_random_seed_local=ran1(-ens_random_seed_local)*1000
1768 rr=2.0*ens_Cdamp_local*ran1(-ens_random_seed_local)-ens_Cdamp_local
1769 ! print*, "zhang cd a=", rr,cdm(i),ens_Cdamp_local,ens_random_seed_local
1770 cdm(i) = cdm(i) *(1.0+rr)
1771 ! print*, "zhang cd b=", rr,cdm(i),ens_Cdamp_local,ens_random_seed_local
1778 end subroutine MFLUX2
1780 SUBROUTINE hwrfsfcinit(isn,XICE,VEGFRA,SNOW,SNOWC,CANWAT,SMSTAV, &
1781 SMSTOT, SFCRUNOFF,UDRUNOFF,GRDFLX,ACSNOW, &
1782 ACSNOM,IVGTYP,ISLTYP,TSLB,SMOIS,DZS,SFCEVP, & ! STEMP
1786 ids,ide, jds,jde, kds,kde, &
1787 ims,ime, jms,jme, kms,kme, &
1788 its,ite, jts,jte, kts,kte )
1793 INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, &
1794 ims,ime, jms,jme, kms,kme, &
1795 its,ite, jts,jte, kts,kte
1797 INTEGER, INTENT(IN) :: num_soil_layers
1799 REAL, DIMENSION( num_soil_layers), INTENT(IN) :: DZS
1801 REAL, DIMENSION( ims:ime, num_soil_layers, jms:jme ) , &
1802 INTENT(INOUT) :: SMOIS, &
1805 REAL, DIMENSION( ims:ime, jms:jme ) , &
1806 INTENT(INOUT) :: SNOW, &
1821 INTEGER, DIMENSION( ims:ime, jms:jme ) , &
1822 INTENT(INOUT) :: IVGTYP, &
1827 INTEGER, INTENT(IN) :: isn
1828 LOGICAL, INTENT(IN) :: allowed_to_read
1831 INTEGER :: icm,jcm,itf,jtf
1857 END SUBROUTINE hwrfsfcinit
1859 END MODULE module_sf_gfdl