wrfv2_fire/phys/module_sf_oml.F

   1 !WRF:MODEL_LAYER:PHYSICS
   2 !
   3 MODULE module_sf_oml
   4
   5 CONTAINS
   6
   7 !----------------------------------------------------------------
   8    SUBROUTINE OCEANML(tml,t0ml,hml,h0ml,huml,hvml,ust,u_phy,v_phy, &
   9                       tmoml,f,g,oml_gamma,                         &
  10                       XLAND,HFX,LH,TSK,GSW,GLW,EMISS,              &
  11                       DELTSM,STBOLT,                               &
  12                       ids,ide, jds,jde, kds,kde,                   &
  13                       ims,ime, jms,jme, kms,kme,                   &
  14                       its,ite, jts,jte, kts,kte                    )
  15
  16 !----------------------------------------------------------------
  17    IMPLICIT NONE
  18 !----------------------------------------------------------------
  19 !
  20 !  SUBROUTINE OCEANML CALCULATES THE SEA SURFACE TEMPERATURE (TSK)
  21 !  FROM A SIMPLE OCEAN MIXED LAYER MODEL BASED ON
  22 !  (Pollard, Rhines and Thompson (1973).
  23 !
  24 !-- TML         ocean mixed layer temperature (K)
  25 !-- T0ML        ocean mixed layer temperature (K) at initial time
  26 !-- TMOML       top 200 m ocean mean temperature (K) at initial time
  27 !-- HML         ocean mixed layer depth (m)
  28 !-- H0ML        ocean mixed layer depth (m) at initial time
  29 !-- HUML        ocean mixed layer u component of wind
  30 !-- HVML        ocean mixed layer v component of wind
  31 !-- OML_GAMMA   deep water lapse rate (K m-1)
  32 !-- UAIR,VAIR   lowest model level wind component
  33 !-- UST         frictional velocity
  34 !-- HFX         upward heat flux at the surface (W/m^2)
  35 !-- LH          latent heat flux at the surface (W/m^2)
  36 !-- TSK         surface temperature (K)
  37 !-- GSW         downward short wave flux at ground surface (W/m^2)
  38 !-- GLW         downward long wave flux at ground surface (W/m^2)
  39 !-- EMISS       emissivity of the surface
  40 !-- XLAND       land mask (1 for land, 2 for water)
  41 !-- STBOLT      Stefan-Boltzmann constant (W/m^2/K^4)
  42 !-- F           Coriolis parameter
  43 !-- DT          time step (second)
  44 !-- G           acceleration due to gravity
  45
  46    INTEGER,  INTENT(IN   )   ::     ids,ide, jds,jde, kds,kde,  &
  47                                     ims,ime, jms,jme, kms,kme,  &
  48                                     its,ite, jts,jte, kts,kte
  49
  50    REAL,     INTENT(IN   )   ::     DELTSM, STBOLT
  51
  52    REAL,    DIMENSION( ims:ime, jms:jme )                     , &
  53             INTENT(IN   )    ::                          EMISS, &
  54                                                          XLAND, &
  55                                                            GSW, &
  56                                                            GLW, &
  57                                                            HFX, &
  58                                                             LH
  59
  60    REAL,    DIMENSION( ims:ime, jms:jme )                     , &
  61             INTENT(INOUT)    ::                            TSK
  62
  63    REAL,    DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) ::     &
  64                                     TML,T0ML,HML,H0ML,HUML,HVML
  65
  66    REAL,    DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN   ) ::     &
  67                                              U_PHY,V_PHY
  68
  69    REAL,    DIMENSION( ims:ime, jms:jme ), INTENT(IN   ) ::     &
  70                                              UST, F, TMOML
  71
  72    REAL,    INTENT(IN   )   ::     G
  73    REAL,    INTENT(IN   )   ::     OML_GAMMA
  74
  75 ! LOCAL VARS
  76
  77    INTEGER ::  I,J
  78
  79    DO J=jts,jte
  80
  81          DO i=its,ite
  82             IF (XLAND(I,J).GT.1.5) THEN
  83                CALL OML1D(I,J,TML(i,j),T0ML(i,j),HML(i,j),H0ML(i,j),           &
  84                           HUML(i,j),HVML(i,j),TSK(i,j),HFX(i,j),               &
  85                           LH(i,j),GSW(i,j),GLW(i,j),TMOML(i,j),                &
  86                           U_PHY(i,kts,j),V_PHY(i,kts,j),UST(i,j),F(i,j),       &
  87                           EMISS(i,j),STBOLT,G,DELTSM,OML_GAMMA,                &
  88                           ids,ide, jds,jde, kds,kde,                           &
  89                           ims,ime, jms,jme, kms,kme,                           &
  90                           its,ite, jts,jte, kts,kte                            )
  91             ENDIF
  92          ENDDO
  93
  94    ENDDO
  95
  96    END SUBROUTINE OCEANML
  97
  98 !----------------------------------------------------------------
  99    SUBROUTINE OML1D(I,J,TML,T0ML,H,H0,HUML,                              &
 100                     HVML,TSK,HFX,                                        &
 101                     LH,GSW,GLW,TMOML,                                    &
 102                     UAIR,VAIR,UST,F,EMISS,STBOLT,G,DT,OML_GAMMA,         &
 103                     ids,ide, jds,jde, kds,kde,                           &
 104                     ims,ime, jms,jme, kms,kme,                           &
 105                     its,ite, jts,jte, kts,kte                            )
 106
 107 !----------------------------------------------------------------
 108    IMPLICIT NONE
 109 !----------------------------------------------------------------
 110 !
 111 !  SUBROUTINE OCEANML CALCULATES THE SEA SURFACE TEMPERATURE (TSK)
 112 !  FROM A SIMPLE OCEAN MIXED LAYER MODEL BASED ON
 113 !  (Pollard, Rhines and Thompson (1973).
 114 !
 115 !-- TML         ocean mixed layer temperature (K)
 116 !-- T0ML        ocean mixed layer temperature (K) at initial time
 117 !-- TMOML       top 200 m ocean mean temperature (K) at initial time
 118 !-- H           ocean mixed layer depth (m)
 119 !-- H0          ocean mixed layer depth (m) at initial time
 120 !-- HUML        ocean mixed layer u component of wind
 121 !-- HVML        ocean mixed layer v component of wind
 122 !-- OML_GAMMA   deep water lapse rate (K m-1)
 123 !-- OMLCALL     whether to call oml model
 124 !-- UAIR,VAIR   lowest model level wind component
 125 !-- UST         frictional velocity
 126 !-- HFX         upward heat flux at the surface (W/m^2)
 127 !-- LH          latent heat flux at the surface (W/m^2)
 128 !-- TSK         surface temperature (K)
 129 !-- GSW         downward short wave flux at ground surface (W/m^2)
 130 !-- GLW         downward long wave flux at ground surface (W/m^2)
 131 !-- EMISS       emissivity of the surface
 132 !-- STBOLT      Stefan-Boltzmann constant (W/m^2/K^4)
 133 !-- F           Coriolis parameter
 134 !-- DT          time step (second)
 135 !-- G           acceleration due to gravity
 136 !
 137 !----------------------------------------------------------------
 138    INTEGER, INTENT(IN   )    ::      I, J
 139    INTEGER, INTENT(IN   )    ::      ids,ide, jds,jde, kds,kde, &
 140                                      ims,ime, jms,jme, kms,kme, &
 141                                      its,ite, jts,jte, kts,kte
 142
 143    REAL,    INTENT(INOUT)    :: TML, H, H0, HUML, HVML, TSK
 144
 145    REAL,    INTENT(IN   )    :: T0ML, HFX, LH, GSW, GLW,        &
 146                                 UAIR, VAIR, UST, F, EMISS, TMOML
 147
 148    REAL,    INTENT(IN) :: STBOLT, G, DT, OML_GAMMA
 149
 150 ! Local
 151    REAL :: rhoair, rhowater, Gam, alp, BV2, A1, A2, B2, u, v, wspd, &
 152            hu1, hv1, hu2, hv2, taux, tauy, tauxair, tauyair, q, hold, &
 153            hsqrd, thp, cwater, ust2
 154    CHARACTER(LEN=120) :: time_series
 155
 156       hu1=huml
 157       hv1=hvml
 158       rhoair=1.
 159       rhowater=1000.
 160       cwater=4200.
 161 ! Deep ocean lapse rate (K/m) - from Rich
 162       Gam=oml_gamma
 163 !     if(i.eq.1 .and. j.eq.1 .or. i.eq.105.and.j.eq.105) print *, 'gamma = ', gam
 164 !     Gam=0.14
 165 !     Gam=5.6/40.
 166 !     Gam=5./100.
 167 ! Thermal expansion coeff (/K)
 168 !     alp=.0002
 169 !     temp dependence (/K)
 170       alp=max((tml-273.15)*1.e-5, 1.e-6)
 171       BV2=alp*g*Gam
 172       thp=t0ml-Gam*(h-h0)
 173       A1=(tml-thp)*h - 0.5*Gam*h*h
 174       if(h.ne.0.)then
 175         u=hu1/h
 176         v=hv1/h
 177       else
 178         u=0.
 179         v=0.
 180       endif
 181
 182 !  time step
 183
 184         q=(-hfx-lh+gsw+glw-stbolt*emiss*tml*tml*tml*tml)/(rhowater*cwater)
 185 !       wspd=max(sqrt(uair*uair+vair*vair),0.1)
 186         wspd=sqrt(uair*uair+vair*vair)
 187         if (wspd .lt. 1.e-10 ) then
 188 !          print *, 'i,j,wspd are ', i,j,wspd
 189            wspd = 1.e-10
 190         endif
 191 ! limit ust to 1.6 to give a value of ust for water of 0.05
 192 !       ust2=min(ust, 1.6)
 193 ! new limit for ust: reduce atmospheric ust by half for ocean
 194         ust2=0.5*ust
 195         tauxair=ust2*ust2*uair/wspd
 196         taux=rhoair/rhowater*tauxair
 197         tauyair=ust2*ust2*vair/wspd
 198         tauy=rhoair/rhowater*tauyair
 199 ! note: forward-backward coriolis force for effective time-centering
 200         hu2=hu1+dt*( f*hv1 + taux)
 201         hv2=hv1+dt*(-f*hu2 + tauy)
 202 ! consider the flux effect
 203         A2=A1+q*dt
 204
 205         huml=hu2
 206         hvml=hv2
 207
 208         hold=h
 209         B2=hu2*hu2+hv2*hv2
 210         hsqrd=-A2/Gam + sqrt(A2*A2/(Gam*Gam) + 2.*B2/BV2)
 211         h=sqrt(max(hsqrd,0.0))
 212 ! limit to positive h change
 213         if(h.lt.hold)h=hold
 214
 215         if(h.ne.0.)then
 216 ! no change unless tml is warmer than layer mean temp tmoml or tsk-5 (see omlinit)
 217           if(tml.ge.tmoml)then
 218             tml=max(t0ml - Gam*(h-h0) + 0.5*Gam*h + A2/h, tmoml)
 219           else
 220             tml=tmoml
 221           endif
 222           u=hu2/h
 223           v=hv2/h
 224         else
 225           tml=t0ml
 226           u=0.
 227           v=0.
 228         endif
 229         tsk=tml
 230 !        if(h.gt.100.)print *,i,j,h,tml,' h,tml'
 231
 232 ! ww: output point data
 233 !     if( (i.eq.190 .and. j.eq.115) .or. (i.eq.170 .and. j.eq.125) ) then
 234 !        write(jtime,fmt='("TS ",f10.0)') float(itimestep)
 235 !        CALL wrf_message ( TRIM(jtime) )
 236 !        write(time_series,fmt='("OML",2I4,2F9.5,2F8.2,2E15.5,F8.3)') &
 237 !              i,j,u,v,tml,h,taux,tauy,a2
 238 !        CALL wrf_message ( TRIM(time_series) )
 239 !     end if
 240
 241    END SUBROUTINE OML1D
 242
 243 !================================================================
 244    SUBROUTINE omlinit(oml_hml0, tsk,                           &
 245                       tml,t0ml,hml,h0ml,huml,hvml,tmoml,       &
 246                       allowed_to_read, start_of_simulation,    &
 247                       ids,ide, jds,jde, kds,kde,               &
 248                       ims,ime, jms,jme, kms,kme,               &
 249                       its,ite, jts,jte, kts,kte                )
 250 !----------------------------------------------------------------
 251    IMPLICIT NONE
 252 !----------------------------------------------------------------
 253    LOGICAL , INTENT(IN)      ::      allowed_to_read
 254    LOGICAL , INTENT(IN)      ::      start_of_simulation
 255    INTEGER, INTENT(IN   )    ::      ids,ide, jds,jde, kds,kde, &
 256                                      ims,ime, jms,jme, kms,kme, &
 257                                      its,ite, jts,jte, kts,kte
 258
 259    REAL,    DIMENSION( ims:ime, jms:jme )                     , &
 260             INTENT(IN)    ::                               TSK
 261
 262    REAL,    DIMENSION( ims:ime, jms:jme )                     , &
 263             INTENT(INOUT)    ::     TML, T0ML, HML, H0ML, HUML, HVML, TMOML
 264    REAL   , INTENT(IN   )    ::     oml_hml0
 265
 266 !  LOCAR VAR
 267
 268    INTEGER                   ::      L,J,I,itf,jtf
 269    CHARACTER*1024 message
 270
 271 !----------------------------------------------------------------
 272
 273    itf=min0(ite,ide-1)
 274    jtf=min0(jte,jde-1)
 275
 276    IF(start_of_simulation) THEN
 277      DO J=jts,jtf
 278      DO I=its,itf
 279        TML(I,J)=TSK(I,J)
 280        T0ML(I,J)=TSK(I,J)
 281      ENDDO
 282      ENDDO
 283      IF (oml_hml0 .gt. 0.) THEN
 284         WRITE(message,*)'Initializing OML with HML0 = ', oml_hml0
 285         CALL wrf_debug (0, TRIM(message))
 286         DO J=jts,jtf
 287         DO I=its,itf
 288           HML(I,J)=oml_hml0
 289           H0ML(I,J)=HML(I,J)
 290           HUML(I,J)=0.
 291           HVML(I,J)=0.
 292           TMOML(I,J)=TSK(I,J)-5.
 293         ENDDO
 294         ENDDO
 295      ELSE
 296         WRITE(message,*)'Initializing OML with real HML0, h(1,1) = ', h0ml(1,1)
 297         CALL wrf_debug (0, TRIM(message))
 298         DO J=jts,jtf
 299         DO I=its,itf
 300           HML(I,J)=H0ML(I,J)
 301 ! fill in near coast area with SST: 200 K was set as missing value in ocean pre-processing code
 302           IF(TMOML(I,J).GT.200. .and. TMOML(I,J).LE.201.) TMOML(I,J)=TSK(I,J)
 303         ENDDO
 304         ENDDO
 305      ENDIF
 306    ENDIF
 307
 308    END SUBROUTINE omlinit
 309 !-------------------------------------------------------------------
 310 END MODULE module_sf_oml