phys/module_sf_gfdl.F

   1 !-- XLV         latent heat of vaporization for water (J/kg)
   2 !
   3 MODULE module_sf_gfdl
   4
   5 !real, dimension(-100:2000,-100:2000), save :: z00
   6
   7 CONTAINS
   8
   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,      &
  13 #if (HWRF==1)
  14                      MZNT,                                      &
  15 #endif
  16                      UST,PSIM,PSIH,                             &
  17                      XLAND,HFX,QFX,TAUX,TAUY,LH,GSW,GLW,TSK,FLHC,FLQC,    & ! gopal's doing for Ocean coupling
  18                      QGH,QSFC,U10,V10,                          &
  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,              &
  23                      Cd_out,Ch_out,                             &
  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 !-------------------------------------------------------------------
  32       IMPLICIT NONE
  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)
  41 !-- ROVCP       R/CP
  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
 102
 103       REAL,    INTENT(IN) ::                                            &
 104                                         CP,                             &
 105                                         EP1,                            &
 106                                         EP2,                            &
 107                                         KARMAN,                         &
 108                                         R,                              &
 109                                         ROVCP,                          &
 110                                         DT,                             &
 111                                         SFENTH,                         &
 112                                         XLV
 113
 114       REAL,    DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(IN) ::      &
 115                                         P3D,                            &
 116                                         QV3D,                           &
 117                                         T3D,                            &
 118                                         U3D,                            &
 119                                         V3D
 120       INTEGER, DIMENSION( ims:ime , jms:jme ), INTENT(IN )::   ISLTYP
 121       REAL, DIMENSION( ims:ime , 1:num_soil_layers, jms:jme ), INTENT(INOUT)::   SMOIS
 122
 123       REAL,    DIMENSION(ims:ime, jms:jme), INTENT(IN) ::               &
 124                                         PSFC,                           &
 125                                         GLW,                            &
 126                                         GSW,                            &
 127                                         XLAND
 128       REAL,    DIMENSION(ims:ime, jms:jme), INTENT(IN) ::               &
 129                                         CHARN,                          &
 130                                         MSANG,                          &
 131                                         SCURX,                          &
 132                                         SCURY
 133       REAL,    DIMENSION(ims:ime, jms:jme), INTENT(INOUT) ::            &
 134                                         TSK,                            &
 135                                         BR,                             &
 136                                         ZKMAX,                          &
 137                                         CHS,                            &
 138                                         Cd_out,                         &
 139                                         Ch_out,                         &
 140                                         CHS2,                           &
 141                                         CPM,                            &
 142                                         CQS2,                           &
 143                                         FLHC,                           &
 144                                         FLQC,                           &
 145                                         GZ1OZ0,                         &
 146                                         HFX,                            &
 147                                         LH,                             &
 148                                         PSIM,                           &
 149                                         PSIH,                           &
 150                                         QFX,                            &
 151                                         QGH,                            &
 152                                         QSFC,                           &
 153                                         UST,                            &
 154                                         ZNT,                            &
 155 #if (HWRF==1)
 156                                         MZNT,                           &   !KWON momentum zo
 157 #endif
 158                                         WSPD,                           &
 159                                         TAUX,                           & ! gopal's doing for Ocean coupling
 160                                         TAUY
 161
 162       REAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) ::               &
 163                                         U10,                            &
 164                                         V10
 165
 166 !--------------------------- LOCAL VARS ------------------------------
 167
 168       REAL ::                           ESAT,                           &
 169                                         cpcgs,                          &
 170                                         smc,                            &
 171                                         smcdry,                         &
 172                                         smcmax
 173
 174       REAL     (kind=kind_phys) ::                                      &
 175                                         RHOX
 176       REAL, DIMENSION(1:30)    ::       MAXSMC, &
 177                                         DRYSMC
 178       REAL     (kind=kind_phys), DIMENSION(its:ite) ::                  &
 179                                         CH,                             &
 180                                         CM,                             &
 181                                         DDVEL,                          &
 182                                         DRAIN,                          &
 183                                         EP,                             &
 184                                         EVAP,                           &
 185                                         FH,                             &
 186                                         FH2,                            &
 187                                         FM,                             &
 188                                         HFLX,                           &
 189                                         PH,                             &
 190                                         PM,                             &
 191                                         PRSL1,                          &
 192                                         PRSLKI,                         &
 193                                         PS,                             &
 194                                         Q1,                             &
 195                                         Q2M,                            &
 196                                         QSS,                            &
 197                                         QSURF,                          &
 198                                         RB,                             &
 199                                         RCL,                            &
 200                                         RHO1,                           &
 201                                         SLIMSK,                         &
 202                                         STRESS,                         &
 203                                         T1,                             &
 204                                         T2M,                            &
 205                                         THGB,                           &
 206                                         THX,                            &
 207                                         TSKIN,                          &
 208                                         SHELEG,                         &
 209                                         U1,                             &
 210                                         U10M,                           &
 211                                         USTAR,                          &
 212                                         V1,                             &
 213                                         V10M,                           &
 214                                         WIND,                           &
 215                                         Z0RL,                           &
 216                                         Z1
 217
 218       REAL,    DIMENSION(kms:kme, ims:ime) ::                           &
 219                                         rpc,                            &
 220                                         tpc,                            &
 221                                         upc,                            &
 222                                         vpc
 223
 224      REAL,    DIMENSION(ims:ime) ::                                     &
 225                                         pspc,                           &
 226                                         pkmax,                          &
 227                                         tstrc,                          &
 228                                         zoc,                            &
 229                                         mzoc,                           &  !ADDED BY KWON FOR momentum Zo
 230                                         tzot,                           &  !Wang
 231                                         wetc,                           &
 232                                         slwdc,                          &
 233                                         rib,                            &
 234 !!!                                        zkmax,                          &
 235                                         tkmax,                          &
 236                                         fxmx,                           &
 237                                         fxmy,                           &
 238                                         cdm,                            &
 239                                         fxh,                            &
 240                                         fxe,                            &
 241                                         xxfh,                           &
 242                                         xxfh2,                          &
 243                                         wind10,                         &
 244                                         tjloc,                          &
 245                                         alpha,                          &
 246                                         gamma,                          &
 247                                         xcur,                           &
 248                                         ycur
 249
 250       INTEGER ::                                                        &
 251                                         I,                              &
 252                                         II,                             &
 253                                         IGPVS,                          &
 254                                         IM,                             &
 255                                         J,                              &
 256                                         K,                              &
 257                                         KM
 258
 259       real :: tmp9,zhalf  ! wang, height of the first half level
 260
 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/
 266
 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/
 272       DATA IGPVS/0/
 273       save igpvs
 274
 275
 276    if(igpvs.eq.0) then
 277 !     call readzo(glat,glon,6,ims,ime,jms,jme,its,ite,jts,jte,z00)
 278    endif
 279    igpvs=1
 280
 281    IM=ITE-ITS+1
 282    KM=KTE-KTS+1
 283
 284    WRITE(message,*)'WITHIN THE GFDL SCHEME, NTSFLG=1 FOR GFDL SLAB  2010 UPGRADS',NTSFLG
 285    call wrf_debug(2,message)
 286
 287 !!     write(0,*)'icoef_sf,lcurr_sf=',icoef_sf,lcurr_sf
 288
 289    DO J=jts,jte
 290
 291       DO i=its,ite
 292         DDVEL(I)=0.
 293         RCL(i)=1.
 294         PRSL1(i)=P3D(i,kts,j)*.001
 295          wetc(i)=1.0
 296          if(xland(i,j).lt.1.99) then
 297          smc=smois(i,1,j)
 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.))
 302          endif
 303 !       convert from Pa to cgs...
 304         pspc(i)=PSFC(i,j)*10.
 305         pkmax(i)=P3D(i,kts,j)*10.
 306         PS(i)=PSFC(i,j)*.001
 307         Q1(I) = QV3D(i,kts,j)
 308         rpc(kts,i)=QV3D(i,kts,j)
 309 !        QSURF(I)=QSFC(I,J)
 310         QSURF(I)=0.
 311         SHELEG(I)=0.
 312         SLIMSK(i)=ABS(XLAND(i,j)-2.)
 313         TSKIN(i)=TSK(i,j)
 314         tstrc(i)=TSK(i,j)
 315         T1(I) = T3D(i,kts,j)
 316         tpc(kts,i)=T3D(i,kts,j)
 317         U1(I) = U3D(i,kts,j)
 318         upc(kts,i)=U3D(i,kts,j) * 100.
 319         USTAR(I) = UST(i,j)
 320         V1(I) = V3D(i,kts,j)
 321         vpc(kts,i)=v3D(i,kts,j) * 100.
 322         Z0RL(I) = ZNT(i,j)*100.
 323         zoc(i)=ZNT(i,j)*100.
 324         cdm(i) = Cd_out(i,j)
 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
 331          tjloc(i)=float(j)
 332         alpha(i) = charn(i,j)
 333         gamma(i) = msang(i,j)
 334         xcur(i) = scurx(i,j)
 335         ycur(i) = scury(i,j)
 336
 337
 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
 342            zhalf=0.0
 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))
 347           endif
 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))
 354 !         endif
 355         !
 356
 357       ENDDO
 358
 359       DO i=its,ite
 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))
 365       ENDDO
 366
 367 !     if(j==2)then
 368 !       write(0,*)'--------------------------------------------'
 369 !       write(0,*) 'u, v, t, r, pkmax, pspc,wetc, tjloc,zoc,tstr'
 370 !       write(0,*)'--------------------------------------------'
 371 !     endif
 372
 373 !     do i = its,ite
 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)
 376 !     enddo
 377
 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  )
 388 !     if(j==2)then
 389 !       write(0,*)'--------------------------------------------'
 390 !       write(0,*) 'fxh, fxe, fxmx, fxmy, cdm, xxfh zoc,tstrc'
 391 !       write(0,*)'--------------------------------------------'
 392 !     endif
 393
 394 !     do i = its,ite
 395 !       WRITE(0,1010)i,j,fxh(i),fxe(i),fxmx(i),fxmy(i),cdm(i),rib(i),xxfh(i),zoc(i),tstrc(i)
 396 !     enddo
 397
 398 1010  format(2I4,9F11.6)
 399
 400
 401
 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,                                    &
 406 !GFDL             Z0RL,                                                 &
 407 !WRF              TG3,GFLUX,F10M,                                       &
 408 !GFDL             U10M,V10M,T2M,Q2M,                                    &
 409 !WRF              ZSOIL,                                                &
 410 !GFDL             CM,CH,RB,                                             &
 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                                      )
 416
 417
 418       DO i=its,ite
 419 !       update skin temp only when using GFDL slab...
 420
 421         IF(NTSFLG==1)    then
 422         tsk(i,j) = tstrc(i)      ! gopal's doing
 423 !bugfix 4
 424 !       bob's doing patch tsk with neigboring values... are grid boundaries
 425         if(j.eq.jde) then
 426          tsk(i,j) = tsk(i,j-1)
 427          endif
 428
 429         if(j.eq.jds) then
 430          tsk(i,j) = tsk(i,j+1)
 431          endif
 432
 433         if(i.eq.ide) tsk(i,j) = tsk(i-1,j)
 434         if(i.eq.ids) tsk(i,j) = tsk(i+1,j)
 435         endif
 436
 437
 438         znt(i,j)= 0.01*abs(zoc(i))
 439 #if (HWRF==1)
 440         mznt(i,j)= 0.01*abs(mzoc(i))
 441 #endif
 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 !------------------------------------------------------------------------
 451
 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)
 457 !       cpcgs  = 1.00464e7
 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))
 469         cm(i)      = cdm(i)
 470         Cd_out(i,j) = cm(i)
 471         Ch_out(i,j) = ch(i)
 472
 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
 483 !           endif
 484          endif
 485 !!!
 486         U10(i,j)=U10M(i)
 487         V10(i,j)=V10M(i)
 488         BR(i,j)=rib(i)
 489         CHS(I,J)=CH(I)*wspd (i,j)
 490         CHS2(I,J)=USTAR(I)*KARMAN/FH2(I)
 491
 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
 497           endif
 498 !!!
 499
 500
 501         CPM(I,J)=CP*(1.+0.8*QV3D(i,kts,j))
 502         esat = fpvs(t1(i))
 503         QGH(I,J)=ep2*esat/(1000.*ps(i)-esat)
 504         esat = fpvs(tskin(i))
 505         qss(i) = ep2*esat/(1000.*ps(i)-esat)
 506         QSFC(I,J)=qss(i)
 507 !       PSIH(i,j)=PH(i)
 508 !       PSIM(i,j)=PM(i)
 509         UST(i,j)=ustar(i)
 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.
 512 !       WSPD(i,j)=WIND(i)
 513 !       ZNT(i,j)=Z0RL(i)*.01
 514       ENDDO
 515
 516 !       write(0,*)'fm,fh,cm,ch(125)', fm(125),fh(125),cm(125),ch(125)
 517
 518       DO i=its,ite
 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))
 522         CQS2(i,j)=CHS2(I,J)
 523       ENDDO
 524
 525       IF (ISFFLX.EQ.0) THEN
 526         DO i=its,ite
 527           HFX(i,j)=0.
 528           LH(i,j)=0.
 529           QFX(i,j)=0.
 530         ENDDO
 531       ELSE
 532         DO i=its,ite
 533           IF(XLAND(I,J)-1.5.GT.0.)THEN
 534 !           HFX(I,J)=FLHC(I,J)*(THGB(I)-THX(I))
 535 !           cpcgs  = 1.00464e7
 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))
 541 !           cpcgs  = 1.00464e7
 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.)
 546           ENDIF
 547 !         QFX(I,J)=FLQC(I,J)*(QSFC(I,J)-Q1(I))
 548 !       convert from g/(cm*cm*sec) to kg/(m*m*sec)
 549           qfx(i,j)=-10.*fxe(i)
 550           QFX(I,J)=AMAX1(QFX(I,J),0.)
 551           LH(I,J)=XLV*QFX(I,J)
 552         ENDDO
 553       ENDIF
 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
 559
 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)
 564     endif
 565
 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)
 595 !      endif
 596
 597
 598    END SUBROUTINE SF_GFDL
 599
 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,    &
 606                            tzot, &
 607                           ids,ide, jds,jde, kds,kde,                      &
 608                           ims,ime, jms,jme, kms,kme,                      &
 609                           its,ite, jts,jte, kts,kte                       )
 610
 611 !------------------------------------------------------------------------
 612 !
 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.
 618 !
 619 !------------------------------------------------------------------------
 620
 621       USE module_sf_exchcoef
 622       IMPLICIT NONE
 623 !     use allocate_mod
 624 !     use module_TLDATA , ONLY : tab,table,cp,g,rgas,og
 625
 626 !     include 'RESOLUTION.h'
 627 !     include 'PARAMETERS.h'
 628 !     include 'STDUNITS.h'     stdout
 629
 630 !     include 'FLAGS.h'
 631 !     include 'BKINFO.h'        nstep
 632 !     include 'CONSLEV.h'
 633 !     include 'CONMLEV.h'
 634 !     include 'ESTAB.h'
 635 !     include 'FILEC.h'
 636 !     include 'FILEPC.h'
 637 !     include 'GDINFO.h'        ngd
 638 !     include 'LIMIT.h'
 639 !     include 'QLOGS.h'
 640 !     include 'TIME.h'          dt(nnst)
 641 !     include 'WINDD.h'
 642 !     include 'ZLDATA.h'        old MOBFLX?
 643
 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
 657
 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
 668
 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
 672
 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
 682
 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
 687
 688 !-----------------------------------------------------------------------
 689 !     internal variables
 690 !-----------------------------------------------------------------------
 691
 692       integer, parameter :: icntx = 30
 693
 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
 699
 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
 766
 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
 776
 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
 785       real :: wndm,ckg
 786       real :: windmks,znott,znotm
 787       real :: ubot, vbot
 788       integer:: i,j,ii,iq,nnest,icnt,ngd,ip
 789
 790 !-----------------------------------------------------------------------
 791 !     internal variables
 792 !-----------------------------------------------------------------------
 793
 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
 803
 804       character(len=255) :: message
 805
 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))
 813
 814       data amask/ -98.0/
 815 !-----------------------------------------------------------------------
 816 !     tables used to obtain the vapor pressures or saturated vapor
 817 !     pressure
 818 !-----------------------------------------------------------------------
 819
 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./
 828
 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.,    &
 838      &190220.,199260./
 839
 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./
 846
 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,          &
 856      &.7220e+01/
 857
 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,        &
 867      &.7090e+03/
 868
 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,        &
 878      &.9780e+04/
 879
 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/
 887 !
 888 ! spcify constants needed by MFLUX2
 889 !
 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
 895 !
 896 #if HWRF==1
 897       real*8 :: gasdev,ran1          !zhang
 898       real :: rr                     !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/
 905       data env_pp_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.
 911          ens_Cdamp_local=0.0
 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
 915                pert_Cd_local=.true.
 916                ens_Cdamp_local=ens_Cdamp
 917             else
 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
 922             endif
 923          else
 924             ens_random_seed_local=ens_random_seed
 925             pert_Cd_local=pert_Cd
 926             ens_Cdamp_local=ens_Cdamp
 927          endif
 928       print*, "Cd ===", ens_random_seed_local,pert_Cd_local,ens_Cdamp_local,ensda_physics_pert
 929       endif
 930 #endif
 931
 932 !     character*10 routine
 933 !     routine = 'mflux2'
 934 !
 935 !------------------------------------------------------------------------
 936 !     set water availability constant "ecof" and land mask "land".
 937 !     limit minimum wind speed to 100 cm/s
 938 !------------------------------------------------------------------------
 939       j = IFIX(tjloc(its))
 940 !  constants for 10 m winds  (correction for knots
 941 !
 942        cor1 = .120
 943        cor2 = 720.
 944 ! KWON : remove the artificial increase of 10m wind speed over 60kts
 945 !        which comes from GFDL hurricane model
 946         cor1 = 0.
 947         cor2 = 0.
 948 !
 949
 950       do i = its,ite
 951         z10(i) = 1000.
 952         z2 (i) =  200.
 953         pss(i) = pspc(i)
 954         tstar(i) = tstrc(i)
 955
 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
 959 !          ubot = upc(1,i)
 960 !          vbot = vpc(1,i)
 961         else
 962           ubot = upc(1,i)
 963           vbot = vpc(1,i)
 964         endif
 965         uvs1(i)= amax1( SQRT(ubot*ubot +    &
 966                              vbot*vbot), 100.0)
 967         if ( iwavecpl .eq. 1 .and. zoc(i) .le. 0.0 ) then
 968           ukmax(i) = ( ubot * cos(gamma(i))  -          &
 969                       vbot * sin(gamma(i)) )            &
 970                                   * cos(gamma(i))
 971           vkmax(i) = ( vbot * cos(gamma(i))  -          &
 972                       ubot * sin(gamma(i)) )            &
 973                                   * cos(gamma(i))
 974
 975         else
 976           ukmax(i) = ubot
 977           vkmax(i) = vbot
 978         endif
 979
 980 !       ukmax(i) = upc(1,i)
 981 !        vkmax(i) = vpc(1,i)
 982         tkmax(i) = tpc(1,i)
 983         rkmax(i) = rpc(1,i)
 984       enddo
 985
 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)
 989
 990       do i = its,ite
 991         windp(i) = SQRT(ukmax(i)*ukmax(i) + vkmax(i)*vkmax(i))
 992         wind (i) = amax1(windp(i),100.)
 993
 994 !! use wind10 previous step
 995          wind10p(i) = wind10(i)  !! cm/s
 996         wind10p(i) = amax1(wind10p(i),100.)
 997 !!
 998
 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
1002           ecof(i) = wetc(i)
1003           land(i) = 1.0
1004           zot (i) = zoc(i)
1005         else
1006           ecof(i) = wetc(i)
1007           land(i) = 0.0
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)
1015             endif
1016             zoc(i)  = -100.*znotm
1017             zot(i) =  -100* znott
1018           else
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)
1025 !            else
1026 !              call  znot_m_v1(windmks,znotm)
1027 !              call  znot_t_v2(windmks,znott)
1028 !            endif
1029             call znot_wind10m(windmks,znott,znotm,icoef_sf)
1030             zoc(i)  = -100.*znotm
1031             zot(i) =  -100* znott
1032           endif
1033         endif
1034
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)
1038 !          endif
1039
1040
1041 !------------------------------------------------------------------------
1042 !     where necessary modify zo values over ocean.
1043 !------------------------------------------------------------------------
1044 !
1045       mzoc(i) = zoc(i)                !FOR SAVE MOMENTUM Zo
1046       tzot(i) = zot(i)                 !output wang
1047       enddo
1048
1049 !------------------------------------------------------------------------
1050 !     define constants:
1051 !     a and c = constants used in evaluating universal function for
1052 !               stable case
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 !------------------------------------------------------------------------
1059
1060       en     = 2.
1061       c      = .76
1062       a      = 5.
1063       ca     = .4
1064       cmo1   = .5*a - 1.648
1065       cmo2   = 17.193 + .5*a - 10.*c
1066       boycon = .61
1067         rovcp=rgas/cp
1068 !       write(0,*)'rgas,cp,rovcp ', rgas,cp,rovcp
1069
1070 !     write(0,*)'--------------------------------------------------'
1071 !     write(0,*)'pkmax,    pspc,    theta,    zkmax,        zoc'
1072 !     write(0,*)'--------------------------------------------------'
1073
1074       do i = its,ite
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)
1081       enddo
1082
1083 !        write(0,*)'pkmax,pspc ',   pkmax,pspc
1084 !        write(0,*)'theta, zkmax, zoc ', theta, zkmax, zoc
1085
1086 !------------------------------------------------------------------------
1087 !     get saturation mixing ratios at surface
1088 !------------------------------------------------------------------------
1089
1090       do i = its,ite
1091         tsg  (i) = tstar(i)
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
1100            psps1 = .1
1101           endif
1102         rstso(i) = 0.622*estso(i)/psps1
1103         vrts (i) = 1. + boycon*ecof(i)*rstso(i)
1104       enddo
1105
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
1110 !     tempa1=1.0e-4
1111 !------------------------------------------------------------------------
1112
1113       do i = its,ite
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 !------------------------------------------------------------------------
1123
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)
1129       enddo
1130
1131       do i = its,ite
1132         zeta(i) = ca*rib(i)/0.03
1133       enddo
1134
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 !------------------------------------------------------------------------
1141
1142 !------------------------------------------------------------------------
1143 !     set initial guess of zeta=non - dimensional height "szeta" for
1144 !     stable points
1145 !------------------------------------------------------------------------
1146
1147       rca   = 1./ca
1148       enrca = en*rca
1149 !     turn off interfacial layer by zeroing out enrca
1150       enrca = 0.0
1151       zog   = .0185*og
1152
1153 !------------------------------------------------------------------------
1154 !     stable points
1155 !------------------------------------------------------------------------
1156
1157       ip    = 0
1158       do i = its,ite
1159         if (zeta(i) .GE. 0.0) then
1160           ip = ip + 1
1161           istb(ip) = i
1162         endif
1163       enddo
1164
1165       if (ip .EQ. 0) go to 170
1166       do i = 1,ip
1167         szetam(i) = 1.0e+30
1168         sgzm(i)   = 0.0e+00
1169         szeta(i)  = zeta(istb(i))
1170         ifz(i)    = 1
1171       enddo
1172
1173 !------------------------------------------------------------------------
1174 !     begin wegstein iteration for "zeta" at stable points using
1175 !     hicks(1976)
1176 !------------------------------------------------------------------------
1177
1178       do icnt = 1,icntx
1179         do i = 1,ip
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
1190           endif
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
1197             rzeta2   = 1./szetao
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
1203             sfzo (i) = c*szetao
1204           endif
1205
1206
1207 !        compute heat & moisture parts of zot.. for calculation of sfh
1208
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
1215             rzeta2   = 1./szetho
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
1221             sfzo (i) = c*szetho
1222           endif
1223
1224 !------------------------------------------------------------------------
1225 !      compute universal function at 10 meters for diagnostic purposes
1226 !------------------------------------------------------------------------
1227
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
1234               rzeta2   = 1./szet10
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
1239             endif
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)
1243             zal2   = ALOG(szet2 )
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
1247               rzeta2   = 1./szet2
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
1252             endif
1253             sf2 (i) = fmz2 (i) - fmzo1(i)
1254
1255 !!!!      endif
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)
1261           fmzo   = ABS(fmz)
1262           if (fmzo .GE. 5.0e-5) then
1263             sq        = (sgz - sgzm(i))/(szeta(i) - szetam(i))
1264             if(sq .EQ. 1) then
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)
1269             endif
1270             szetam(i) = szeta(i)
1271             szeta (i) = (sgz - szeta(i)*sq)/(1.0 - sq)
1272             sgzm  (i) = sgz
1273           else
1274             ifz(i) = 0
1275           endif
1276 80        continue
1277         enddo
1278       enddo
1279
1280       do i = 1,ip
1281         if (ifz(i) .GE. 1) go to 110
1282       enddo
1283
1284       go to 130
1285
1286 110   continue
1287
1288       write(6,120)
1289 120   format(2X, ' NON-CONVERGENCE FOR STABLE ZETA IN ROW ')
1290 !     call MPI_CLOSE(1,routine)
1291
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 !------------------------------------------------------------------------
1297
1298 130   continue
1299       do i = 1,ip
1300         szo = zoc(istb(i))
1301         if (szo .LT. 0.0)  then
1302         wndm=wind(istb(i))*0.01
1303           if(wndm.lt.15.0) then
1304           ckg=0.0185*og
1305           else
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
1309           endif
1310
1311        szo =  - ckg*wind(istb(i))*wind(istb(i))/  &
1312                              (sfm(i)*sfm(i))
1313         cons_p000001    =    .000001
1314         cons_7          =         7.
1315         vis             =     1.4E-1
1316
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
1323         rat=0.
1324         zot(istb(i)) = szo   * exp(-rat)
1325          else
1326          zot(istb(i)) = zoc(istb(i))
1327         endif
1328
1329 !      in hwrf thermal znot is loaded back into the zoc array for next step
1330              zoc(istb(i)) = szo
1331       enddo
1332
1333       do i = 1,ip
1334         xxfm(istb(i)) = sfm(i)
1335         xxfh(istb(i)) = sfh(i)
1336         xxfh2(istb(i)) = sf2 (i)
1337         xxsh(istb(i)) = sfzo(i)
1338       enddo
1339
1340 !------------------------------------------------------------------------
1341 !     obtain wind at 10 meters for diagnostic purposes
1342 !------------------------------------------------------------------------
1343
1344 !!!   if (ngd .EQ. nNEST) then
1345         do i = 1,ip
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 &
1350                         - cor2
1351            endif
1352 !     the above correction done by GFDL in centi-kts!!!-change back
1353          wind10(istb(i)) = wind10(istb(i)) / 1.944
1354         enddo
1355 !!!   endif
1356 !!!   if (ngd .EQ. nNEST-1 .AND. llwe .EQ. 1 ) then
1357 !!!     do i = 1,ip
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
1362 !!!  *                   - cor2
1363 !!!        endif
1364 !!!     enddo
1365 !!!   endif
1366
1367 !------------------------------------------------------------------------
1368 !     unstable points
1369 !------------------------------------------------------------------------
1370
1371 170   continue
1372
1373       iq = 0
1374       do i = its,ite
1375         if (zeta(i) .LT. 0.0) then
1376           iq       = iq + 1
1377           iutb(iq) = i
1378         endif
1379       enddo
1380
1381       if (iq .EQ. 0) go to 290
1382       do i = 1,iq
1383         uzeta (i) = zeta(iutb(i))
1384         ifz   (i) = 1
1385         uzetam(i) = 1.0e+30
1386         ugzm  (i) = 0.0e+00
1387       enddo
1388
1389 !------------------------------------------------------------------------
1390 !     begin wegstein iteration for "zeta" at unstable points using
1391 !     hicks functions
1392 !------------------------------------------------------------------------
1393
1394       do icnt = 1,icntx
1395         do i = 1,iq
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
1401           y       = SQRT(ux11)
1402           yo      = SQRT(ux12)
1403           ufzo(i) = 1./yo
1404           ux13    = (1. + y)/(1. + yo)
1405           ux21    = ALOG(ux13)
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
1412           y       = SQRT(ux11)
1413           yo      = SQRT(ux12)
1414           ux13    = (1. + y)/(1. + yo)
1415           ux21    = ALOG(ux13)
1416 !          ufzo(i) = 1./yo
1417           x       = SQRT(y)
1418           xo      = SQRT(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
1425
1426 !------------------------------------------------------------------------
1427 !     obtain ten meter winds for diagnostic purposes
1428 !------------------------------------------------------------------------
1429
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
1435             y       = SQRT(ux11)
1436             y10     = SQRT(ux12)
1437             ux13    = (1. + y)/(1. + y10)
1438             ux21    = ALOG(ux13)
1439             x       = SQRT(y)
1440             x10     = SQRT(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
1446
1447 !   obtain 2m values for diagnostics...
1448
1449
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
1455           y       = SQRT(ux11)
1456           yo      = SQRT(ux12)
1457           ux13    = (1. + y)/(1. + yo)
1458           ux21    = ALOG(ux13)
1459           uf2 (i)  = (ugzzo - 2.*ux21)*rca
1460
1461
1462 !!!       endif
1463           ugz = ca*rib(iutb(i))*ufm(i)*ufm(i)/(ufh(i) + enrca*ufzo(i))
1464           ux1 = (ugz - uzeta(i))/uzeta(i)
1465           ux2 = ABS(ux1)
1466           if (ux2 .GE. 5.0e-5) then
1467             uq        = (ugz - ugzm(i))/(uzeta(i) - uzetam(i))
1468             uzetam(i) = uzeta(i)
1469             if(uq .EQ. 1) then
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)
1473             endif
1474             uzeta (i) = (ugz - uzeta(i)*uq)/(1.0 - uq)
1475             ugzm  (i) = ugz
1476           else
1477             ifz(i) = 0
1478           endif
1479 200       continue
1480         enddo
1481       enddo
1482
1483
1484       do i = 1,iq
1485         if (ifz(i) .GE. 1) go to 230
1486       enddo
1487
1488       go to 250
1489
1490 230   continue
1491       write(message,240)
1492       call wrf_message(message)
1493 240   format(2X, ' NON-CONVERGENCE FOR UNSTABLE ZETA IN ROW ')
1494 !     call MPI_CLOSE(1,routine)
1495
1496 !------------------------------------------------------------------------
1497 !     gather unstable values
1498 !------------------------------------------------------------------------
1499
1500 250   continue
1501
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 !------------------------------------------------------------------------
1507
1508       do i = 1,iq
1509         uzo = zoc(iutb(i))
1510         if (zoc(iutb(i)) .LT. 0.0)   then
1511         wndm=wind(iutb(i))*0.01
1512          if(wndm.lt.15.0) then
1513          ckg=0.0185*og
1514          else
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
1519          endif
1520         uzo         =-ckg*wind(iutb(i))*wind(iutb(i))/(ufm(i)*ufm(i))
1521         cons_p000001    =    .000001
1522         cons_7          =         7.
1523         vis             =     1.4E-1
1524
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
1531         rat=0.0
1532         zot(iutb(i)) =  uzo   * exp(-rat)
1533          else
1534          zot(iutb(i)) = zoc(iutb(i))
1535         endif
1536 !      in hwrf thermal znot is loaded back into the zoc array for next step
1537             zoc(iutb(i)) = uzo
1538       enddo
1539
1540 !------------------------------------------------------------------------
1541 !     obtain wind at ten meters for diagnostic purposes
1542 !------------------------------------------------------------------------
1543
1544 !!!   if (ngd .EQ. nNEST) then
1545         do i = 1,iq
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 &
1550                          - cor2
1551            endif
1552  !     the above correction done by GFDL in centi-kts!!!-change back
1553           wind10(iutb(i)) = wind10(iutb(i)) / 1.944
1554         enddo
1555 !!!   endif
1556 !!!   if (ngd .EQ. nNEST-1) then
1557 !!!     do i = 1,iq
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
1562 !!!  *                    - cor2
1563 !!!        endif
1564 !!!     enddo
1565 !!!   endif
1566
1567       do i = 1,iq
1568         xxfm(iutb(i)) = ufm(i)
1569         xxfh(iutb(i)) = ufh(i)
1570         xxfh2(iutb(i)) = uf2 (i)
1571         xxsh(iutb(i)) = ufzo(i)
1572       enddo
1573
1574 290   continue
1575
1576       do i = its,ite
1577         ucom(i) = ukmax(i)
1578         vcom(i) = vkmax(i)
1579         if (windp(i) .EQ. 0.0) then
1580           windp(i) = 100.0
1581           ucom (i) = 100.0/SQRT(2.0)
1582           vcom (i) = 100.0/SQRT(2.0)
1583         endif
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)))
1587       enddo
1588
1589 !     do land sfc temperature prediction if ntsflg=1
1590 !     ntsflg = 1                                    ! gopal's doing
1591
1592       if (ntsflg .EQ. 0) go to 370
1593       alll = 600.
1594       xks   = 0.01
1595       hcap  = .5/2.39e-8
1596       pith  = SQRT(4.*ATAN(1.0))
1597       alfus = alll/2.39e-8
1598       teps  = 0.1
1599 !     slwdc... in units of cal/min ????
1600 !     slwa...  in units of ergs/sec/cm*2
1601 !     1 erg=2.39e-8 cal
1602 !------------------------------------------------------------------------
1603 !     pack land and sea ice points
1604 !------------------------------------------------------------------------
1605
1606       ip    = 0
1607       do i = its,ite
1608         if (land(i) .EQ. 1) then
1609           ip = ip + 1
1610           indx   (ip) = i
1611 !         slwa is defined as positive down....
1612           slwa   (ip) =    slwdc(i)/(2.39e-8*60.)
1613           tss    (ip) = tstar(i)
1614           thetap (ip) = theta(i)
1615           rkmaxp (ip) = rkmax(i)
1616           aap    (ip) = 5.673e-5
1617           pssp   (ip) = pss(i)
1618           ecofp  (ip) = ecof(i)
1619           estsop (ip) = estso(i)
1620           rstsop (ip) = rstso(i)
1621           bq1p   (ip) = bq1(i)
1622           bq1p   (ip) = amax1(bq1p(ip),0.1e-3)
1623           delsrad(ip) = dt   *pith/(hcap*SQRT(3600.*24.*xks))
1624         endif
1625       enddo
1626
1627 !------------------------------------------------------------------------
1628 !     initialize variables for first pass of iteration
1629 !------------------------------------------------------------------------
1630
1631       do i = 1,ip
1632         ifz  (i) = 1
1633         tsm  (i) = tss(i)
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
1637
1638 !!!      if (j .EQ. 1 .AND. i .EQ. 1) write(6,300)
1639 300   format(2X, ' SURFACE EQUILIBRIUM CALCULATION ')
1640
1641 !!        foftm(i) = thetap(i) + 1./(cp*bq1p(i))*(slwa(i) - aap(i)* &
1642 !!                  tsm(i)**4 + ecofp(i)*alfus*bq1p(i)*rdiff(i))
1643 !!      else
1644
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)* &
1647            rdiff(i))
1648 !!      endif
1649         tsp(i) = foftm(i)
1650       enddo
1651
1652 !------------------------------------------------------------------------
1653 !     do iteration to determine "tstar" at new time level
1654 !------------------------------------------------------------------------
1655
1656       do icnt = 1,icntx
1657         do i = 1,ip
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
1667                psps2 = .1
1668              endif
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))
1675 !!!       else
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)* &
1678              rdiff(i))
1679 !!!       endif
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)
1687 !!!       endif
1688           frac(i) = ABS((foft(i) - tsp(i))/tsp(i))
1689
1690 !------------------------------------------------------------------------
1691 !      check for convergence of all points use wegstein iteration
1692 !------------------------------------------------------------------------
1693
1694           if (frac(i) .GE. teps) then
1695             qf   (i) = (foft(i) - foftm(i))/(tsp(i) - tsm(i))
1696             tsm  (i) = tsp(i)
1697             tsp  (i) = (foft(i) - tsp(i)*qf(i))/(1. - qf(i))
1698             foftm(i) = foft(i)
1699           else
1700             ifz(i) = 0
1701           endif
1702 330       continue
1703         enddo
1704       enddo
1705
1706 !------------------------------------------------------------------------
1707 !     check for convergence of "t star" prediction
1708 !------------------------------------------------------------------------
1709
1710       do i = 1,ip
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, &
1714             2I5)
1715         call wrf_message(message)
1716
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)
1720
1721 !       write(6,350) reflect, sigt4, shfx, alevp, delten, diffot
1722 350   format(2X, ' REFLECT, SIGT4, SHFX, ALEVP, DELTEN, DIFFOT ', &
1723             6E14.8)
1724
1725 !         call MPI_CLOSE(1,routine)
1726         endif
1727       enddo
1728
1729       do i = 1,ip
1730         ii        = indx(i)
1731         tstrc(ii) = tsp (i)
1732       enddo
1733
1734 !------------------------------------------------------------------------
1735 !     compute fluxes  and momentum drag coef
1736 !------------------------------------------------------------------------
1737
1738 370   continue
1739       do i = its,ite
1740 !!!
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)
1747           endif
1748           zoc(i)  = -100.*znotm
1749           zot(i) =  -100* znott
1750         endif
1751 !!!!
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)/ &
1756                  windp(i)
1757         fxmy(i) = rho(i)/(xxfm(i)*xxfm(i))*wind(i)*wind(i)*vcom(i)/ &
1758         windp(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)
1762 #if HWRF==1
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
1772         endif
1773 #endif
1774
1775       enddo
1776       ntstep = ntstep + 1
1777       return
1778       end subroutine MFLUX2
1779
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
1783                         TMN,                                            &
1784                         num_soil_layers,                                &
1785                         allowed_to_read,                                &
1786                         ids,ide, jds,jde, kds,kde,                      &
1787                         ims,ime, jms,jme, kms,kme,                      &
1788                         its,ite, jts,jte, kts,kte                     )
1789
1790    IMPLICIT NONE
1791
1792 ! Arguments
1793    INTEGER,  INTENT(IN   )   ::     ids,ide, jds,jde, kds,kde,  &
1794                                     ims,ime, jms,jme, kms,kme,  &
1795                                     its,ite, jts,jte, kts,kte
1796
1797    INTEGER, INTENT(IN)       ::     num_soil_layers
1798
1799    REAL,    DIMENSION( num_soil_layers), INTENT(IN) :: DZS
1800
1801    REAL,    DIMENSION( ims:ime, num_soil_layers, jms:jme )    , &
1802             INTENT(INOUT)    ::                          SMOIS, &
1803                                                          TSLB      !STEMP
1804
1805    REAL,    DIMENSION( ims:ime, jms:jme )                     , &
1806             INTENT(INOUT)    ::                           SNOW, &
1807                                                          SNOWC, &
1808                                                         CANWAT, &
1809                                                         SMSTAV, &
1810                                                         SMSTOT, &
1811                                                      SFCRUNOFF, &
1812                                                       UDRUNOFF, &
1813                                                         SFCEVP, &
1814                                                         GRDFLX, &
1815                                                         ACSNOW, &
1816                                                           XICE, &
1817                                                         VEGFRA, &
1818                                                         TMN, &
1819                                                         ACSNOM
1820
1821    INTEGER, DIMENSION( ims:ime, jms:jme )                     , &
1822             INTENT(INOUT)    ::                         IVGTYP, &
1823                                                         ISLTYP
1824
1825 !
1826
1827   INTEGER, INTENT(IN) :: isn
1828   LOGICAL, INTENT(IN) :: allowed_to_read
1829 ! Local
1830   INTEGER             :: iseason
1831   INTEGER :: icm,jcm,itf,jtf
1832   INTEGER ::  I,J,L
1833
1834
1835    itf=min0(ite,ide-1)
1836    jtf=min0(jte,jde-1)
1837
1838    icm = ide/2
1839    jcm = jde/2
1840
1841    iseason=isn
1842
1843    DO J=jts,jtf
1844        DO I=its,itf
1845 !      SNOW(i,j)=0.
1846        SNOWC(i,j)=0.
1847 !      SMSTAV(i,j)=
1848 !      SMSTOT(i,j)=
1849 !      SFCRUNOFF(i,j)=
1850 !      UDRUNOFF(i,j)=
1851 !      GRDFLX(i,j)=
1852 !      ACSNOW(i,j)=
1853 !      ACSNOM(i,j)=
1854     ENDDO
1855    ENDDO
1856
1857   END SUBROUTINE hwrfsfcinit
1858
1859 END MODULE module_sf_gfdl