| 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 tmol or tsk-5 (see omlinit) |
|---|
| 217 | if(tml.ge.tmoml .and. h.ne.0.)then |
|---|
| 218 | tml=max(t0ml - Gam*(h-h0) + 0.5*Gam*h + A2/h, tmoml) |
|---|
| 219 | u=hu2/h |
|---|
| 220 | v=hv2/h |
|---|
| 221 | else |
|---|
| 222 | tml=t0ml |
|---|
| 223 | u=0. |
|---|
| 224 | v=0. |
|---|
| 225 | endif |
|---|
| 226 | tsk=tml |
|---|
| 227 | ! if(h.gt.100.)print *,i,j,h,tml,' h,tml' |
|---|
| 228 | |
|---|
| 229 | ! ww: output point data |
|---|
| 230 | ! if( (i.eq.190 .and. j.eq.115) .or. (i.eq.170 .and. j.eq.125) ) then |
|---|
| 231 | ! write(jtime,fmt='("TS ",f10.0)') float(itimestep) |
|---|
| 232 | ! CALL wrf_message ( TRIM(jtime) ) |
|---|
| 233 | ! write(time_series,fmt='("OML",2I4,2F9.5,2F8.2,2E15.5,F8.3)') & |
|---|
| 234 | ! i,j,u,v,tml,h,taux,tauy,a2 |
|---|
| 235 | ! CALL wrf_message ( TRIM(time_series) ) |
|---|
| 236 | ! end if |
|---|
| 237 | |
|---|
| 238 | END SUBROUTINE OML1D |
|---|
| 239 | |
|---|
| 240 | !================================================================ |
|---|
| 241 | SUBROUTINE omlinit(oml_hml0, tsk, & |
|---|
| 242 | tml,t0ml,hml,h0ml,huml,hvml,tmoml, & |
|---|
| 243 | allowed_to_read, start_of_simulation, & |
|---|
| 244 | ids,ide, jds,jde, kds,kde, & |
|---|
| 245 | ims,ime, jms,jme, kms,kme, & |
|---|
| 246 | its,ite, jts,jte, kts,kte ) |
|---|
| 247 | !---------------------------------------------------------------- |
|---|
| 248 | IMPLICIT NONE |
|---|
| 249 | !---------------------------------------------------------------- |
|---|
| 250 | LOGICAL , INTENT(IN) :: allowed_to_read |
|---|
| 251 | LOGICAL , INTENT(IN) :: start_of_simulation |
|---|
| 252 | INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, & |
|---|
| 253 | ims,ime, jms,jme, kms,kme, & |
|---|
| 254 | its,ite, jts,jte, kts,kte |
|---|
| 255 | |
|---|
| 256 | REAL, DIMENSION( ims:ime, jms:jme ) , & |
|---|
| 257 | INTENT(IN) :: TSK |
|---|
| 258 | |
|---|
| 259 | REAL, DIMENSION( ims:ime, jms:jme ) , & |
|---|
| 260 | INTENT(INOUT) :: TML, T0ML, HML, H0ML, HUML, HVML, TMOML |
|---|
| 261 | REAL , INTENT(IN ) :: oml_hml0 |
|---|
| 262 | |
|---|
| 263 | ! LOCAR VAR |
|---|
| 264 | |
|---|
| 265 | INTEGER :: L,J,I,itf,jtf |
|---|
| 266 | CHARACTER*1024 message |
|---|
| 267 | |
|---|
| 268 | !---------------------------------------------------------------- |
|---|
| 269 | |
|---|
| 270 | itf=min0(ite,ide-1) |
|---|
| 271 | jtf=min0(jte,jde-1) |
|---|
| 272 | |
|---|
| 273 | IF(start_of_simulation) THEN |
|---|
| 274 | DO J=jts,jtf |
|---|
| 275 | DO I=its,itf |
|---|
| 276 | TML(I,J)=TSK(I,J) |
|---|
| 277 | T0ML(I,J)=TSK(I,J) |
|---|
| 278 | ENDDO |
|---|
| 279 | ENDDO |
|---|
| 280 | IF (oml_hml0 .gt. 0.) THEN |
|---|
| 281 | WRITE(message,*)'Initializing OML with HML0 = ', oml_hml0 |
|---|
| 282 | CALL wrf_debug (0, TRIM(message)) |
|---|
| 283 | DO J=jts,jtf |
|---|
| 284 | DO I=its,itf |
|---|
| 285 | HML(I,J)=oml_hml0 |
|---|
| 286 | H0ML(I,J)=HML(I,J) |
|---|
| 287 | HUML(I,J)=0. |
|---|
| 288 | HVML(I,J)=0. |
|---|
| 289 | TMOML(I,J)=TSK(I,J)-5. |
|---|
| 290 | ENDDO |
|---|
| 291 | ENDDO |
|---|
| 292 | ELSE |
|---|
| 293 | WRITE(message,*)'Initializing OML with real HML0, h(1,1) = ', h0ml(1,1) |
|---|
| 294 | CALL wrf_debug (0, TRIM(message)) |
|---|
| 295 | DO J=jts,jtf |
|---|
| 296 | DO I=its,itf |
|---|
| 297 | HML(I,J)=H0ML(I,J) |
|---|
| 298 | ! fill in near coast area with SST: 200 K was set as missing value in ocean pre-processing code |
|---|
| 299 | IF(TMOML(I,J).GT.200. .and. TMOML(I,J).LE.201.) TMOML(I,J)=TSK(I,J) |
|---|
| 300 | ENDDO |
|---|
| 301 | ENDDO |
|---|
| 302 | ENDIF |
|---|
| 303 | ENDIF |
|---|
| 304 | |
|---|
| 305 | END SUBROUTINE omlinit |
|---|
| 306 | !------------------------------------------------------------------- |
|---|
| 307 | END MODULE module_sf_oml |
|---|