!-- XLV latent heat of vaporization for water (J/kg) ! MODULE MODULE_SF_GFDL !real, dimension(-100:2000,-100:2000), save :: z00 ! HISTORY LOG ! 2014-06-24 Weiguo Wang, move from HWRF to NMMB, ! change (i,k,j) to (i,j,k), nmmb's z-coordinate is from top to ! bottom, use 'kflip' ! 2015-10-23, Weiguo Wang, port subroutines calculating z0 from hwrf ! 2016-03-16 Weiguo Wang, use 10m wind to compute z0, add more formula for z0 ! USE MODULE_KINDS USE module_sf_exchcoef use ifport REAL,PRIVATE,SAVE :: randnum1, firstcall CONTAINS !------------------------------------------------------------------- SUBROUTINE SF_GFDL(U3D,V3D,T3D,QV3D,P3D, & CP,ROVCP,R,XLV,PSFC,CHS,CHS2,CQS2, CPM, & DT, SMOIS,num_soil_layers,ISLTYP,ZNT, & !#if (HWRF==1) MZNT, & !#endif UST,PSIM,PSIH, & XLAND,HFX,QFX,TAUX,TAUY,LH,GSW,GLW,TSK,FLHC,FLQC, & ! gopal's doing for Ocean coupling QGH,QSFC,U10,V10, & GZ1OZ0,WSPD,BR,ISFFLX, & EP1,EP2,KARMAN,NTSFLG,SFENTH,ICOEF_SF, & R_SEED, PERT_Z0, & ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & its,ite, jts,jte, kts,kte ) !------------------------------------------------------------------- USE MACHINE, ONLY : kind_phys USE FUNCPHYS , ONLY : gfuncphys,fpvs USE PHYSCONS, grav => con_g !------------------------------------------------------------------- IMPLICIT NONE !------------------------------------------------------------------- !-- U3D 3D u-velocity interpolated to theta points (m/s) !-- V3D 3D v-velocity interpolated to theta points (m/s) !-- T3D temperature (K) !-- QV3D 3D water vapor mixing ratio (Kg/Kg) !-- P3D 3D pressure (Pa) !-- DT time step (second) !-- CP heat capacity at constant pressure for dry air (J/kg/K) !-- ROVCP R/CP !-- R gas constant for dry air (J/kg/K) !-- XLV latent heat of vaporization for water (J/kg) !-- PSFC surface pressure (Pa) !-- ZNT thermal roughness length (m) !-- MZNT momentum roughness length (m) !-- MAVAIL surface moisture availability (between 0 and 1) !-- UST u* in similarity theory (m/s) !-- PSIM similarity stability function for momentum !-- PSIH similarity stability function for heat !-- XLAND land mask (1 for land, 2 for water) !-- HFX upward heat flux at the surface (W/m^2) !-- QFX upward moisture flux at the surface (kg/m^2/s) !-- TAUX RHO*U**2 (Kg/m/s^2) ! gopal's doing for Ocean coupling !-- TAUY RHO*U**2 (Kg/m/s^2) ! gopal's doing for Ocean coupling !-- LH net upward latent heat flux at surface (W/m^2) !-- GSW downward short wave flux at ground surface (W/m^2) !-- GLW downward long wave flux at ground surface (W/m^2) !-- TSK surface temperature (K) !-- FLHC exchange coefficient for heat (m/s) !-- FLQC exchange coefficient for moisture (m/s) !-- QGH lowest-level saturated mixing ratio !-- U10 diagnostic 10m u wind !-- V10 diagnostic 10m v wind !-- GZ1OZ0 log(z/z0) where z0 is roughness length !-- WSPD wind speed at lowest model level (m/s) !-- BR bulk Richardson number in surface layer !-- ISFFLX isfflx=1 for surface heat and moisture fluxes !-- EP1 constant for virtual temperature (R_v/R_d - 1) (dimensionless) !-- KARMAN Von Karman constant !-- SFENTH enthalpy flux factor 0 zot via charnock ..>0 zot enhanced>15m/s !-- ids start index for i in domain !-- ide end index for i in domain !-- jds start index for j in domain !-- jde end index for j in domain !-- kds start index for k in domain !-- kde end index for k in domain !-- ims start index for i in memory !-- ime end index for i in memory !-- jms start index for j in memory !-- jme end index for j in memory !-- kms start index for k in memory !-- kme end index for k in memory !-- its start index for i in tile !-- ite end index for i in tile !-- jts start index for j in tile !-- jte end index for j in tile !-- kts start index for k in tile !-- kte end index for k in tile !------------------------------------------------------------------- character*255 :: message INTEGER, INTENT(IN) :: ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & its,ite, jts,jte, kts,kte, & ISFFLX,NUM_SOIL_LAYERS,NTSFLG, & ICOEF_SF REAL, INTENT(IN) :: & CP, & EP1, & EP2, & KARMAN, & R, & ROVCP, & DT, & SFENTH, & XLV, R_SEED , PERT_Z0 !REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(IN) :: & REAL, DIMENSION(ims:ime, jms:jme, kms:kme), INTENT(IN) :: & P3D, & QV3D, & T3D, & U3D, & V3D INTEGER, DIMENSION( ims:ime , jms:jme ), INTENT(IN ):: ISLTYP !REAL, DIMENSION( ims:ime , 1:num_soil_layers, jms:jme ), INTENT(INOUT):: SMOIS REAL, DIMENSION( ims:ime , jms:jme, 1:num_soil_layers ), INTENT(INOUT):: SMOIS REAL, DIMENSION(ims:ime, jms:jme), INTENT(IN) :: & PSFC, & GLW, & GSW, & XLAND REAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: & TSK, & BR, & CHS, & CHS2, & CPM, & CQS2, & FLHC, & FLQC, & GZ1OZ0, & HFX, & LH, & PSIM, & PSIH, & QFX, & QGH, & QSFC, & UST, & ZNT, & !#if (HWRF==1) MZNT, & !KWON momentum zo !#endif WSPD, & TAUX, & ! gopal's doing for Ocean coupling TAUY REAL, DIMENSION(ims:ime, jms:jme), INTENT(OUT) :: & U10, & V10 !--------------------------- LOCAL VARS ------------------------------ REAL :: ESAT, & cpcgs, & smc, & smcdry, & smcmax REAL (kind=kind_phys) :: & !REAL :: & RHOX REAL, DIMENSION(1:30) :: MAXSMC, & DRYSMC REAL (kind=kind_phys), DIMENSION(its:ite) :: & !REAL , DIMENSION(its:ite) :: & ! CH, & ! CM, & DDVEL, & DRAIN, & EP, & EVAP, & ! FH, & ! FH2, & ! FM, & HFLX, & PH, & PM, & PRSL1, & PRSLKI, & PS, & Q1, & Q2M, & QSS, & QSURF, & RB, & RCL, & RHO1, & SLIMSK, & STRESS, & T1, & T2M, & THGB, & THX, & TSKIN, & SHELEG, & U1, & U10M, & USTAR, & V1, & V10M, & WIND, & Z0RL, & Z1 REAL , DIMENSION(its:ite) :: & CH, & CM, & FH, & FH2, & FM REAL, DIMENSION(kms:kme, ims:ime) :: & rpc, & tpc, & upc, & vpc REAL, DIMENSION(ims:ime) :: & pspc, & pkmax, & tstrc, & ustold, & zoc, & mzoc, & !ADDED BY KWON FOR momentum Zo wetc, & slwdc, & rib, & zkmax, & tkmax, & fxmx, & fxmy, & cdm, & fxh, & fxe, & xxfh, & xxfh2, & wind10, & tjloc INTEGER :: & I, & II, & IGPVS, & IM, & J, & K, & KM, & KFLIP LOGICAL :: lpr real :: zhalf ! wang, height of the first half level real :: tmp88, tmp99 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/ DATA IGPVS/0/ save igpvs !INTEGER, INTENT(IN) :: ICOEF_SF !LOGICAL, INTENT(IN) :: LCURR_SF ! INTEGER :: ICOEF_SF ! in the future, will input from namelist LOGICAL :: LCURR_SF LCURR_SF=.false. ! ICOEF_SF=2 !0- old cd, old ch, 1-new cd new ch, 2- new cd old ch !icoef_sf=4 !ICOEF_SF=3 !3 --- lower cd, than (2), 4-- lower cd for wind > 20 m/s, !icoef_sf=5 !5- similar to (4), ch same as (2) ! icoef_sf=6 ! use HWRF17 cd ch, wind10m ! icoef_sf=7 ! use HWRF17 cd ch, wind10m, but increase CD when w<20m/s,then reduce to 0.8 at w-45, ch*1.12 ! icoef_sf=8 ! use HWRF17 cd ch, wind10m, change cd,ch depending on wind, see code, ! icoef_sf=9 ! use HWRF17 cd ch, wind10m, change cd,ch depending on wind, see code, ! icoef_sf=10 ! use HWRF17 cd ch, wind10m ! icoef_sf=11 ! use HWRF17 cd ch, wind10m ! icoef_sf=15 ! use HWRF17 cd ch, wind10m ! write(0,*)'icoef_sf=',icoef_sf if( firstcall /= -999. ) then ! firstcall is not initialized, no chance it is equal to -999 randnum1=rand( jts*ide*2+its+int(r_seed*1000) ) firstcall=-999. endif ! write(0,*)'r_seed,pert_z0,randnum1=',r_seed,pert_z0,randnum1 lpr=.false. if(igpvs.eq.0) then ! call readzo(glat,glon,6,ims,ime,jms,jme,its,ite,jts,jte,z00) endif igpvs=1 IM=ITE-ITS+1 KM=KTE-KTS+1 WRITE(message,*)'WITHIN THE GFDL SCHEME, NTSFLG=1 FOR GFDL SLAB 2010 UPGRADS',NTSFLG !call wrf_debug(1,message) kflip=(kte+1)-kts !NMMB input z is from top to bottom, !wang DO J=jts,jte DO i=its,ite if( lpr .and. i == (its+ite)/2 .and. j == (jts+jte)/2 ) then write(0,*)'smois(i,j,1)=',smois(i,j,1:4) write(0,*)'p3d=',p3d(i,j,kflip) write(0,*)'u3d=',u3d(i,j,kflip) write(0,*)'v3d=',v3d(i,j,kflip) write(0,*)'t3d=',t3d(i,j,kflip) write(0,*)'qv3d=',qv3d(i,j,kflip) write(0,*)'psfc=',psfc(i,j) write(0,*)'GLW=',GLW(i,j) write(0,*)'GSW=',GSW(i,j) write(0,*)'xland=',xland(i,j) write(0,*)'isltyp=',isltyp(i,j) write(0,*)'ZNT=',ZNT(i,j) write(0,*)'HFX=',HFX(i,j) write(0,*)'QFX=',QFX(i,j) write(0,*)'CP=',cp write(0,*)'EP1=',ep1 write(0,*)'EP2=',ep2 write(0,*)'R=',R write(0,*)'Rovcp=',rovcp write(0,*)'DT=',DT write(0,*)'SFENTH=',sfenth write(0,*)'XLV=',XLV endif DDVEL(I)=0. RCL(i)=1. PRSL1(i)=P3D(i,j,kflip)*.001 wetc(i)=1.0 if(xland(i,j).lt.1.99) then smc=smois(i,j,1) smcdry=drysmc(isltyp(i,j)) smcmax=maxsmc(isltyp(i,j)) wetc(i)=(smc-smcdry)/(smcmax-smcdry) wetc(i)=amin1(1.,amax1(wetc(i),0.)) endif ! convert from Pa to cgs... pspc(i)=PSFC(i,j)*10. pkmax(i)=P3D(i,j,kflip)*10. PS(i)=PSFC(i,j)*.001 !Q1(I) = QV3D(i,j,kflip) !rpc(kts,i)=QV3D(i,j,kflip) Q1(I) = QV3D(i,j,kflip)/( 1.0-QV3D(i,j,kflip) ) !specific humidity (input) to Mixing ratio. rpc(kts,i)=QV3D(i,j,kflip)/( 1.0-QV3D(i,j,kflip) ) ! QSURF(I)=QSFC(I,J) QSURF(I)=0. SHELEG(I)=0. SLIMSK(i)=ABS(XLAND(i,j)-2.) !TSKIN(i)=TSK(i,j) !tstrc(i)=TSK(i,j) TSKIN(i)=amin1(TSK(i,j),320.0) tstrc(i)=amin1(TSK(i,j),320.0) T1(I) = T3D(i,j,kflip) tpc(kts,i)=T3D(i,j,kflip) U1(I) = U3D(i,j,kflip) upc(kts,i)=U3D(i,j,kflip) * 100. USTAR(I) = UST(i,j) ustold(I) = UST(i,j) V1(I) = V3D(i,j,kflip) vpc(kts,i)=v3D(i,j,kflip) * 100. Z0RL(I) = ZNT(i,j)*100. zoc(i)=ZNT(i,j)*100. if(XLAND(i,j).gt.1.99) zoc(i)=- zoc(i) if(XLAND(i,j).le.1.99) zoc(i)=amax1(zoc(i),2.0) ! Z0RL(I) = z00(i,j)*100. ! 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,j)+glw(i,j) slwdc(i)=0.239*60.*slwdc(i)*1.e-4 tjloc(i)=float(j) !Wang, calulate height of the first half level ! use previous u10 v10 to compute wind10, input to MFLUX to compute z0 wind10(i)=sqrt(u10(i,j)*u10(i,j)+v10(i,j)*v10(i,j)) !first time step, u10 and v10 are all zero zhalf=0.0 if (wind10(i) <= 1.0e-10 .or. wind10(i) > 150.0) then ! write(0,*)'maybe first timestep' zhalf = -R*tpc(kts,i)*alog(pkmax(i)/pspc(i))/grav wind10(i)=sqrt(u1(i)*u1(i)+v1(i)*v1(i))*alog(10.0/znt(i,j))/alog(zhalf/znt(i,j)) endif wind10(i)=wind10(i)*100.0 !! m/s to cm/s zhalf = -R*tpc(kts,i)*alog(pkmax(i)/pspc(i))/grav if (lpr .and. i == (its+ite)/2 .and. j == (jts+jte)/2 ) then write(0,*)'before call mflux,zh,w10,u10,v10,WIND1,znt' write(0,*)zhalf,wind10(i),u10(i,j),v10(i,j),sqrt(u1(i)**2+v1(i)**2),znt(i,j) write(0,*)'uselog',sqrt(u1(i)*u1(i)+v1(i)*v1(i))*alog(10.0/znt(i,j))/alog(zhalf/znt(i,j )) endif ! if( lpr .and. i == (its+ite)/2 .and. j == (jts+jte)/2 ) then write(0,*)'isltyp(i,j),wetc(i),PSFC(i,j),Q1(I),SLIMSK(i),TSKIN(i),T1(I),U1(I),USTAR(I),V1(I),gsw(i,j),glw(i,j)' write(0,*)isltyp(i,j),wetc(i),PSFC(i,j),Q1(I),SLIMSK(i),TSKIN(i),T1(I),U1(I),USTAR(I),V1(I),gsw(i,j),glw(i,j) WRITE(0,*)'i,j,upc(kts,i),vpc(kts,i),tpc(kts,i),rpc(kts,i), & &pkmax(i),pspc(i),wetc(i),tjloc(i),zoc(i),tstrc(i)' WRITE(0,*)i,j,upc(kts,i),vpc(kts,i),tpc(kts,i),rpc(kts,i), & pkmax(i),pspc(i),wetc(i),tjloc(i),zoc(i),tstrc(i) endif ENDDO DO i=its,ite PRSLKI(i)=(PS(I)/PRSL1(I))**ROVCP THGB(I)=TSKIN(i)*(100./PS(I))**ROVCP THX(I)=T1(i)*(100./PRSL1(I))**ROVCP RHO1(I)=PRSL1(I)*1000./(R*T1(I)*(1.+EP1*Q1(I))) Q1(I)=Q1(I)/(1.+Q1(I)) ENDDO ! if(j==2)then ! write(0,*)'--------------------------------------------' ! write(0,*) 'u, v, t, r, pkmax, pspc,wetc, tjloc,zoc,tstr' ! write(0,*)'--------------------------------------------' ! endif ! do i = its,ite ! WRITE(0,1010)i,j,upc(kts,i),vpc(kts,i),tpc(kts,i),rpc(kts,i), & ! pkmax(i),pspc(i),wetc(i),tjloc(i),zoc(i),tstrc(i) ! enddo !CALL MFLUX2( fxh,fxe,fxmx,fxmy,cdm,rib,xxfh,zoc,mzoc,tstrc, & !mzoc for momentum Zo KWON CALL MFLUX2( fxh,fxe,fxmx,fxmy,cdm,rib,xxfh,ustold,zoc,mzoc,tstrc, & !mzoc for momentum Zo KWON pspc,pkmax,wetc,slwdc,tjloc, & icoef_sf,lcurr_sf, & upc,vpc,tpc,rpc,dt,J,wind10,xxfh2,ntsflg,SFENTH, & r_seed, pert_z0, & ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & its,ite, jts,jte, kts,kte ) ! if(j==2)then ! write(0,*)'--------------------------------------------' ! write(0,*) 'fxh, fxe, fxmx, fxmy, cdm, xxfh zoc,tstrc' ! write(0,*)'--------------------------------------------' ! endif ! do i = its,ite ! WRITE(0,1010)i,j,fxh(i),fxe(i),fxmx(i),fxmy(i),cdm(i),rib(i),xxfh(i),zoc(i),tstrc(i) ! enddo 1010 format(2I4,9F11.6) !GFDL CALL PROGTM(IM,KM,PS,U1,V1,T1,Q1, & !GFDL SHELEG,TSKIN,QSURF, & !WRF SMC,STC,DM,SOILTYP,SIGMAF,VEGTYPE,CANOPY,DLWFLX, & !WRF SLRAD,SNOWMT,DELT, & !GFDL Z0RL, & !WRF TG3,GFLUX,F10M, & !GFDL U10M,V10M,T2M,Q2M, & !WRF ZSOIL, & !GFDL CM,CH,RB, & !WRF RHSCNPY,RHSMC,AIM,BIM,CIM, & !GFDL RCL,PRSL1,PRSLKI,SLIMSK, & !GFDL DRAIN,EVAP,HFLX,STRESS,EP, & !GFDL FM,FH,USTAR,WIND,DDVEL, & !GFDL PM,PH,FH2,QSS,Z1 ) DO i=its,ite ! update skin temp only when using GFDL slab... IF(NTSFLG==1) then tsk(i,j) = tstrc(i) ! gopal's doing !bugfix 4 ! bob's doing patch tsk with neigboring values... are grid boundaries if(j.eq.jde) then tsk(i,j) = tsk(i,j-1) endif if(j.eq.jds) then tsk(i,j) = tsk(i,j+1) endif if(i.eq.ide) tsk(i,j) = tsk(i-1,j) if(i.eq.ids) tsk(i,j) = tsk(i+1,j) endif znt(i,j)= 0.01*abs(zoc(i)) !#if (HWRF==1) mznt(i,j)= 0.01*abs(mzoc(i)) !#endif wspd(i,j) = SQRT(upc(kts,i)*upc(kts,i) + vpc(kts,i)*vpc(kts,i)) wspd(i,j) = amax1(wspd(i,j) ,100.)/100. u10m(i) = u1(i)*(wind10(i)/wspd(i,j))/100. v10m(i) = v1(i)*(wind10(i)/wspd(i,j))/100. ! br =0.0001*zfull(i,kmax)*dthv/ ! & (gmks*theta(i,kmax)*wspd *wspd ) ! zkmax = rgas*tpc(kmax,i)*qqlog(kmax)*og zkmax(i) = -R*tpc(kts,i)*alog(pkmax(i)/pspc(i))/grav !------------------------------------------------------------------------ gz1oz0(i,j)=alog(zkmax(i)/znt(i,j)) ustar (i)= 0.01*sqrt(cdm(i)* & (upc(kts,i)*upc(kts,i) + vpc(kts,i)*vpc(kts,i))) ! amax1((upc(kts,i)*upc(kts,i) + vpc(kts,i)*vpc(kts,i)),1.0e4) ) !cm^2 ! convert from g/(cm*cm*sec) to kg/(m*m*sec) qfx (i,j)=-10.*fxe(i) ! BOB: qfx (i,1)=-10.*fxe(i) ! cpcgs = 1.00464e7 ! convert from ergs/gram/K to J/kg/K cpmks=1004 ! hfx (i,1)=-0.001*cpcgs*fxh(i) hfx (i,j)= -10.*CP*fxh(i) ! Bob: hfx (i,1)= -10.*CP*fxh(i) taux (i,j)= fxmx(i)/10. ! gopal's doing for Ocean coupling tauy (i,j)= fxmy(i)/10. ! gopal's doing for Ocean coupling fm(i) = karman/sqrt(cdm(i)) fh(i) = karman*xxfh(i) PSIM(i,j)=GZ1OZ0(i,j)-FM(i) PSIH(i,j)=GZ1OZ0(i,j)-FH(i) fh2(i) = karman*xxfh2(i) ch(i) = karman*karman/(fm(i) * fh(i)) cm(i) = cdm(i) U10(i,j)=U10M(i) V10(i,j)=V10M(i) BR(i,j)=rib(i) CHS(I,J)=CH(I)*wspd (i,j) CHS2(I,J)=USTAR(I)*KARMAN/FH2(I) !wang tmp88=alog( (zkmax(i)+znt(i,j))/znt(i,j) ) CHS(I,J)=amax1(wspd (i,j),5.0)*karman*karman/ & ( amin1( amax1(abs(fm(i)),0.5),1.8*tmp88) * amin1( amax1(abs(fh(i)),0.5),1.8*tmp88) ) CHS2(I,J)=USTAR(I)*KARMAN/amax1(FH2(I),1.0) if (tsk(i,j) > 360) then write(0,*)'tsk,chs,chs2, fh2, fh,fm,wspd,rib' write(0,*)tsk(i,j),chs(i,j),chs2(i,j), fh2(i),fh(i),fm(i),wspd(i,j),rib(i) endif !!!2015-1020 cap CHS over land points, following 2014version HWRF, if (xland(i,j) .lt. 1.9) then !CHS2(I,J)=amin1(CHS2(I,J), 0.05) !CHS(I,J)=amin1(CHS(I,J), 1.0) if (chs2(i,j) < 0) chs2(i,j)=1.0e-6 endif !!! CPM(I,J)=CP*(1.+0.8*QV3D(i,j,kflip)) esat = fpvs(t1(i)) QGH(I,J)=ep2*esat/(1000.*ps(i)-esat) esat = fpvs(tskin(i)) qss(i) = ep2*esat/(1000.*ps(i)-esat) QSFC(I,J)=qss(i) !! mixng ratio, from noah_lsm, qsfc is mixing ratio ! PSIH(i,j)=PH(i) ! PSIM(i,j)=PM(i) UST(i,j)=ustar(i) ! wspd (i,j) = SQRT(upc(kts,i)*upc(kts,i) + vpc(kts,i)*vpc(kts,i)) ! wspd (i,j) = amax1(wspd (i,j) ,100.)/100. ! WSPD(i,j)=WIND(i) ! ZNT(i,j)=Z0RL(i)*.01 ENDDO ! write(0,*)'fm,fh,cm,ch(125)', fm(125),fh(125),cm(125),ch(125) DO i=its,ite FLHC(i,j)=CPM(I,J)*RHO1(I)*CHS(I,J) FLQC(i,j)=RHO1(I)*CHS(I,J) ! GZ1OZ0(i,j)=LOG(Z1(I)/(Z0RL(I)*.01)) CQS2(i,j)=CHS2(I,J) ENDDO IF (ISFFLX.EQ.0) THEN DO i=its,ite HFX(i,j)=0. LH(i,j)=0. QFX(i,j)=0. ENDDO ELSE DO i=its,ite IF(XLAND(I,J)-1.5.GT.0.)THEN ! HFX(I,J)=FLHC(I,J)*(THGB(I)-THX(I)) ! cpcgs = 1.00464e7 ! convert from ergs/gram/K to J/kg/K cpmks=1004 ! hfx (i,j)=-0.001*cpcgs*fxh(i) hfx (i,j)= -10.*CP*fxh(i) ! Bob: hfx (i,1)= -10.*CP*fxh(i) ELSEIF(XLAND(I,J)-1.5.LT.0.)THEN ! HFX(I,J)=FLHC(I,J)*(THGB(I)-THX(I)) ! cpcgs = 1.00464e7 ! convert from ergs/gram/K to J/kg/K cpmks=1004 ! hfx (i,j)=-0.001*cpcgs*fxh(i) hfx (i,j)= -10.*CP*fxh(i) ! Bob: hfx (i,j)= -10.*CP*fxh(i) HFX(I,J)=AMAX1(HFX(I,J),-250.) ENDIF ! QFX(I,J)=FLQC(I,J)*(QSFC(I,J)-Q1(I)) ! convert from g/(cm*cm*sec) to kg/(m*m*sec) qfx(i,j)=-10.*fxe(i) QFX(I,J)=AMAX1(QFX(I,J),0.) LH(I,J)=XLV*QFX(I,J) ENDDO ENDIF ! if(j.eq.2) write(0,*) 'u3d ,ustar,cdm at end of gfdlsfcmod' ! write(0,*) j,(u3d(ii,1,j),ii=70,90) ! write(0,*) j,(ustar(ii),ii=70,90) ! write(0,*) j,(cdm(ii),ii=70,90) if(j.eq.jds.or.j.eq.jde) then write(message,*) "TSFC in gfdl sf mod,dt, its,ite,jts,jts", dt,its,ite,jts,jte,ids,ide,jds,jde !call wrf_debug(1,message) write(message,*) "TSFC", (TSK(i,j),i=its,ite) !call wrf_debug(1,message) endif ENDDO ! FOR THE J LOOP I PRESUME ! if(100.ge.its.and.100.le.ite.and.100.ge.jts.and.100.le.jte) then ! write(0,*) 'output vars of sf_gfdl at i,j=100' ! write(0,*) 'TSK',TSK(100,100) ! write(0,*) 'PSFC',PSFC(100,100) ! write(0,*) 'GLW',GLW(100,100) ! write(0,*) 'GSW',GSW(100,100) ! write(0,*) 'XLAND',XLAND(100,100) ! write(0,*) 'BR',BR(100,100) ! write(0,*) 'CHS',CHS(100,100) ! write(0,*) 'CHS2',CHS2(100,100) ! write(0,*) 'CPM',CPM(100,100) ! write(0,*) 'FLHC',FLHC(100,100) ! write(0,*) 'FLQC',FLQC(100,100) ! write(0,*) 'GZ1OZ0',GZ1OZ0(100,100) ! write(0,*) 'HFX',HFX(100,100) ! write(0,*) 'LH',LH(100,100) ! write(0,*) 'PSIM',PSIM(100,100) ! write(0,*) 'PSIH',PSIH(100,100) ! write(0,*) 'QFX',QFX(100,100) ! write(0,*) 'QGH',QGH(100,100) ! write(0,*) 'QSFC',QSFC(100,100) ! write(0,*) 'UST',UST(100,100) ! write(0,*) 'ZNT',ZNT(100,100) ! write(0,*) 'wet',wet(100) ! write(0,*) 'smois',smois(100,1,100) ! write(0,*) 'WSPD',WSPD(100,100) ! write(0,*) 'U10',U10(100,100) ! write(0,*) 'V10',V10(100,100) ! endif if (lpr) then write(0,*)'in SF_GFDL,ust(i,j)=',ust((its+ite)/2,(jts+jte)/2) write(0,*)'in SF_GFDL,znt(i,j)=',znt((its+ite)/2,(jts+jte)/2) write(0,*)'in SF_GFDL,mznt(i,j)=',mznt((its+ite)/2,(jts+jte)/2) write(0,*)'in SF_GFDL,xland(i,j)=',xland((its+ite)/2,(jts+jte)/2) write(0,*)'in SF_GFDL,ISLTYP(i,j)=',ISLTYP((its+ite)/2,(jts+jte)/2) write(0,*)'hfx(i,j)=',hfx((its+ite)/2,jts:jte) write(0,*)'chs(i,j)=',chs((its+ite)/2,jts:jte) write(0,*)'netural chs(i,j)=',ust((its+ite)/2,jts:jte)*0.4/gz1oz0((its+ite)/2,jts:jte) write(0,*)'min,max hfx=',minval(hfx(its:ite,jts:jte)),maxval(hfx(its:ite,jts:jte)) write(0,*)'min,max qfx=',minval(qfx(its:ite,jts:jte)),maxval(qfx(its:ite,jts:jte)) write(0,*)'min,max chs=',minval(chs(its:ite,jts:jte)),maxval(chs(its:ite,jts:jte)) write(0,*)'min,max ust=',minval(ust(its:ite,jts:jte)),maxval(ust(its:ite,jts:jte)) write(0,*)'min,max znt=',minval(znt(its:ite,jts:jte)),maxval(znt(its:ite,jts:jte)) endif END SUBROUTINE SF_GFDL !------------------------------------------------------------------- !SUBROUTINE MFLUX2( fxh,fxe,fxmx,fxmy,cdm,rib,xxfh,zoc,mzoc,tstrc, & !mzoc KWON SUBROUTINE MFLUX2( fxh,fxe,fxmx,fxmy,cdm,rib,xxfh,ustold,zoc,mzoc,tstrc, & !mzoc KWON pspc,pkmax,wetc,slwdc,tjloc, & icoef_sf,lcurr_sf, & upc,vpc,tpc,rpc,dt,jfix,wind10,xxfh2,ntsflg,sfenth, & r_seed,pert_z0, & ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & its,ite, jts,jte, kts,kte ) !------------------------------------------------------------------------ ! ! 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. ! !------------------------------------------------------------------------ IMPLICIT NONE ! use allocate_mod ! use module_TLDATA , ONLY : tab,table,cp,g,rgas,og ! include 'RESOLUTION.h' ! include 'PARAMETERS.h' ! include 'STDUNITS.h' stdout ! include 'FLAGS.h' ! include 'BKINFO.h' nstep ! include 'CONSLEV.h' ! include 'CONMLEV.h' ! include 'ESTAB.h' ! include 'FILEC.h' ! include 'FILEPC.h' ! include 'GDINFO.h' ngd ! include 'LIMIT.h' ! include 'QLOGS.h' ! include 'TIME.h' dt(nnst) ! include 'WINDD.h' ! include 'ZLDATA.h' old MOBFLX? !----------------------------------------------------------------------- ! user interface variables !----------------------------------------------------------------------- integer,intent(in) :: ids,ide, jds,jde, kds,kde integer,intent(in) :: ims,ime, jms,jme, kms,kme integer,intent(in) :: its,ite, jts,jte, kts,kte integer,intent(in) :: jfix,ntsflg integer,intent(in) :: icoef_sf logical,intent(in) :: lcurr_sf real, intent (out), dimension (ims :ime ) :: fxh real, intent (out), dimension (ims :ime ) :: fxe real, intent (out), dimension (ims :ime ) :: fxmx real, intent (out), dimension (ims :ime ) :: fxmy real, intent (out), dimension (ims :ime ) :: cdm ! real, intent (out), dimension (ims :ime ) :: cdm2 real, intent (out), dimension (ims :ime ) :: rib real, intent (out), dimension (ims :ime ) :: xxfh real, intent (out), dimension (ims :ime ) :: xxfh2 real, intent (out), dimension (ims :ime ) :: wind10 real, intent ( inout), dimension (ims :ime ) :: zoc,mzoc !KWON real, intent ( inout), dimension (ims :ime ) :: tstrc real, intent ( inout), dimension (ims :ime ) :: ustold real, intent ( in) :: dt real, intent ( in) :: sfenth real, intent ( in), dimension (ims :ime ) :: pspc real, intent ( in), dimension (ims :ime ) :: pkmax real, intent ( in), dimension (ims :ime ) :: wetc real, intent ( in), dimension (ims :ime ) :: slwdc real, intent ( in), dimension (ims :ime ) :: tjloc real, intent ( in), dimension (kms:kme, ims :ime ) :: upc real, intent ( in), dimension (kms:kme, ims :ime ) :: vpc real, intent ( in), dimension (kms:kme, ims :ime ) :: tpc real, intent ( in), dimension (kms:kme, ims :ime ) :: rpc real, intent(in) :: r_seed, pert_z0 !----------------------------------------------------------------------- ! 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, dimension(1 :ime) :: aap real, dimension(1 :ime) :: bq1 real, dimension(1 :ime) :: bq1p real, dimension(1 :ime) :: delsrad real, dimension(1 :ime) :: ecof real, dimension(1 :ime) :: ecofp real, dimension(1 :ime) :: estso real, dimension(1 :ime) :: estsop real, dimension(1 :ime) :: fmz1 real, dimension(1 :ime) :: fmz10 real, dimension(1 :ime) :: fmz2 real, dimension(1 :ime) :: fmzo1 real, dimension(1 :ime) :: foft real, dimension(1 :ime) :: foftm real, dimension(1 :ime) :: frac real, dimension(1 :ime) :: land real, dimension(1 :ime) :: pssp real, dimension(1 :ime) :: qf real, dimension(1 :ime) :: rdiff real, dimension(1 :ime) :: rho real, dimension(1 :ime) :: rkmaxp real, dimension(1 :ime) :: rstso real, dimension(1 :ime) :: rstsop real, dimension(1 :ime) :: sf10 real, dimension(1 :ime) :: sf2 real, dimension(1 :ime) :: sfm real, dimension(1 :ime) :: sfzo real, dimension(1 :ime) :: sgzm real, dimension(1 :ime) :: slwa real, dimension(1 :ime) :: szeta real, dimension(1 :ime) :: szetam real, dimension(1 :ime) :: t1 real, dimension(1 :ime) :: t2 real, dimension(1 :ime) :: tab1 real, dimension(1 :ime) :: tab2 real, dimension(1 :ime) :: tempa1 real, dimension(1 :ime) :: tempa2 real, dimension(1 :ime) :: theta real, dimension(1 :ime) :: thetap real, dimension(1 :ime) :: tsg real, dimension(1 :ime) :: tsm real, dimension(1 :ime) :: tsp real, dimension(1 :ime) :: tss real, dimension(1 :ime) :: ucom real, dimension(1 :ime) :: uf10 real, dimension(1 :ime) :: uf2 real, dimension(1 :ime) :: ufh real, dimension(1 :ime) :: ufm real, dimension(1 :ime) :: ufzo real, dimension(1 :ime) :: ugzm real, dimension(1 :ime) :: uzeta real, dimension(1 :ime) :: uzetam real, dimension(1 :ime) :: vcom real, dimension(1 :ime) :: vrtkx real, dimension(1 :ime) :: vrts real, dimension(1 :ime) :: wind real, dimension(1 :ime) :: windp real, dimension(1 :ime) :: wind10p ! wang, save 10m wind at previous time step real, dimension(1 :ime) :: windx ! wang, used for 10m wind calculation ! real, dimension(1 :ime) :: xxfh real, dimension(1 :ime) :: xxfm real, dimension(1 :ime) :: xxsh real, dimension(1 :ime) :: z10 real, dimension(1 :ime) :: z2 real, dimension(1 :ime) :: zeta real, dimension(1 :ime) :: zkmax real, dimension(1 :ime) :: pss real, dimension(1 :ime) :: tstar real, dimension(1 :ime) :: ukmax real, dimension(1 :ime) :: vkmax real, dimension(1 :ime) :: tkmax real, dimension(1 :ime) :: rkmax real, dimension(1 :ime) :: zot real, dimension(1 :ime) :: fhzo1 real, dimension(1 :ime) :: sfh real :: ux13, yo, y,xo,x,ux21,ugzzo,ux11,ux12,uzetao,xnum,alll real :: ux1,ugz,x10,uzo,uq,ux2,ux3,xtan,xden,y10,uzet1o,ugz10 real :: szet2, zal2,ugz2 real :: rovcp,boycon,cmo2,psps1,zog,enrca,rca,cmo1,amask,en,ca,a,c real :: sgz,zal10,szet10,fmz,szo,sq,fmzo,rzeta1,zal1g,szetao,rzeta2,zal2g real :: hcap,xks,pith,teps,diffot,delten,alevp,psps2,alfus,nstep real :: shfx,sigt4,reflect real :: cor1,cor2,szetho,zal2gh,cons_p000001,cons_7,vis,ustar,restar,rat real :: wndm,ckg real :: yz,y1,y2,y3,y4,windmks,znott,znotm integer:: i,j,ii,iq,nnest,icnt,ngd,ip !----------------------------------------------------------------------- ! internal variables !----------------------------------------------------------------------- real, dimension (223) :: tab real, dimension (223) :: table real, dimension (101) :: tab11 real, dimension (41) :: table4 real, dimension (42) :: tab3 real, dimension (54) :: table2 real, dimension (54) :: table3 real, dimension (74) :: table1 real, dimension (80) :: tab22 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 ! real,parameter :: cp = 1.00464e7 real,parameter :: g = 980.6 real,parameter :: rgas = 2.87e6 real,parameter :: og = 1./g character*255 :: message real, parameter :: CZETMAX = 10. ,CZIL=0.1,SQVISC=258.2,EPSZT=1.0e-28 real, parameter :: ZILFC=-CZIL*0.4*SQVISC, RIC=0.505 real zzil ! ! character*10 routine ! routine = 'mflux2' ! !------------------------------------------------------------------------ ! set water availability constant "ecof" and land mask "land". ! limit minimum wind speed to 100 cm/s !------------------------------------------------------------------------ j = IFIX(tjloc(its)) ! 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. ! yz= -0.0001344 y1= 3.015e-05 y2= 1.517e-06 y3= -3.567e-08 y4= 2.046e-10 ifz(:)=0 !write(0,*)'ifz=',ifz(:) do i = its,ite z10(i) = 1000. z2 (i) = 200. pss(i) = pspc(i) tstar(i) = tstrc(i) ukmax(i) = upc(1,i) vkmax(i) = vpc(1,i) tkmax(i) = tpc(1,i) rkmax(i) = rpc(1,i) enddo ! write(0,*)'z10,pss,tstar,u...rkmax(ite)', & ! z10(ite), pss(ite),tstar(ite),ukmax(ite), & ! vkmax(ite),tkmax(ite),rkmax(ite) do i = its,ite windp(i) = SQRT(ukmax(i)*ukmax(i) + vkmax(i)*vkmax(i)) !wind (i) = amax1(windp(i),100.) !wind (i) = amax1(windp(i),100.) wind (i) = amax1(windp(i),500.) windx (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*wind(i)*wind(i)*og ! 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) !! wang, 2017, use a simple model to estimate zt,Zeng Dickinson,1998 !zot (i) = zoc(i)*amax1(1e-5,exp(-0.13*(amax1(ustold(i),0.1)*(zoc(i)/100)/1.5e-5)**0.45) ) ! zoc in cm ! zot (i) = zoc(i)/exp(2.0) ! zoc in cm zot (i) = zoc(i)*0.1 ! zoc in cm ! ZILFC=-0.1*0.4*258.2 !ZILFC=-CZIL*VKARMAN*SQVISC else ecof(i) = wetc(i) land(i) = 0.0 zot (i) = zoc(i) ! now use 2 regime fit for znot thermal !windmks=wind(i)*.01 windmks=wind10p(i)*.01 ! if ( icoef_sf .EQ. 1) then ! call znot_m_v1(windmks,znotm) ! call znot_t_v1(windmks,znott) ! else if ( icoef_sf .EQ. 0 ) then ! call znot_m_v0(windmks,znotm) ! call znot_t_v0(windmks,znott) ! else ! call znot_m_v1(windmks,znotm) ! call znot_t_v2(windmks,znott) ! endif call znot_wind10m(windmks,znott,znotm,icoef_sf) zoc(i) = -100.*znotm zot(i) = -100* znott ! write(0,*)'in random1=',r_seed,randnum1, zoc(i), pert_z0 if(r_seed > 0.0) then randnum1=rand( 10 + int(randnum1*1e6) ) zoc(i)=zoc(i)*( 1.0 + 2.0*(randnum1-0.5)*pert_z0 ) endif ! write(0,*)'out random1=',randnum1, zoc(i) endif !wang, foloowing MYJ z2(i)=200.0+ abs( zoc(i) ) z10(i)=1000.0 + abs( zoc(i) ) !wang !w znott=0.2375*exp(-0.5250*windmks) + 0.0025*exp(-0.0211*windmks) !w znott=0.01*znott ! go back to moon et al for below 7m/s !w if(windmks.le. 7.) & !w znott = (0.0185/9.8*(7.59e-8*wind(i)**2+ & !w 2.46e-4*wind(i))**2) ! back to cgs !w zot (i) = 100.*znott ! end of kwon correction.... ! in hwrf, thermal znot(zot) is passed as argument zoc ! in hwrf, momentum znot is recalculated internally !w zoc(i)=-(0.0185/9.8*(7.59e-8*wind(i)**2+ & !w 2.46e-4*wind(i))**2)*100. !w if(wind(i).ge.1250.0) & !w zoc(i)=-(.000739793 * wind(i) -0.58)/10 !w if(wind(i).ge.3000.) then !w windmks=wind(i)*.01 ! kwon znotm !w znotm = yz +windmks*y1 +windmks**2*y2 +windmks**3*y3 +windmks**4*y4 !powell 2003 ! back to cgs !w zoc(i) = 100.*znotm !w endif !w endif !------------------------------------------------------------------------ ! where necessary modify zo values over ocean. !------------------------------------------------------------------------ ! mzoc(i) = zoc(i) !FOR SAVE MOMENTUM Zo 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 ! write(0,*)'rgas,cp,rovcp ', rgas,cp,rovcp ! write(0,*)'--------------------------------------------------' ! write(0,*)'pkmax, pspc, theta, zkmax, zoc' ! write(0,*)'--------------------------------------------------' do i = its,ite ! theta(i) = tkmax(i)*rqc9 theta(i) = tkmax(i)/((pkmax(i)/pspc(i))**rovcp) vrtkx(i) = 1.0 + boycon*rkmax(i) ! zkmax(i) = rgas*tkmax(i)*qqlog(kmax)*og zkmax(i) = -rgas*tkmax(i)*alog(pkmax(i)/pspc(i))*og !wang zkmax(i) = zkmax(i) + abs( zoc(i) ) ! following MYJ's ! ! IF(I==78)write(0,*)I,JFIX,pkmax(i),pspc(i),theta(i),zkmax(i),zoc(i) enddo ! write(0,*)'pkmax,pspc ', pkmax,pspc ! write(0,*)'theta, zkmax, zoc ', theta, zkmax, zoc !------------------------------------------------------------------------ ! get saturation mixing ratios at surface !------------------------------------------------------------------------ do i = its,ite tsg (i) = tstar(i) !tab1 (i) = tstar(i) - 153.16 tab1 (i) = amin1(tstar(i),330.0) - 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))) ! t1 (i) = tab(it(i) + 1) ! t2 (i) = table(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) = 0.622*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) !wang ! if (rib (i) <= -50.0) rib(i)=-50.0 !wang ! calculate Zot for land point ! if (zoc(i) > 0.0 ) then ! stable ! IF(rib(i)>0.)THEN ! IF (rib(i)