phys/module_sf_oml.F

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