| blob | 82b3f8442b1d885ec1450d8c39d169e67d8517d9 |
1 !WRF:MODEL_LAYER:PHYSICS
2 !
3 MODULE module_sf_3dpwp
5 CONTAINS
7 !-------------------------------------------------------------------
8 ! cyl: Implemented a three-dimensional ocean model (DPWP) Chiaying Lee RSMAS/UM
9 !---------------------------------------------------------------------------
10 SUBROUTINE DPWP (ims,ime, jms,jme, kms,kme,its,ite, jts,jte, kts,kte, &
11 ids,ide, jds,jde, kds,kde,okms, okme, &
12 om_tmp,om_s,om_u, om_v, om_density, om_depth, om_ml, &
13 om_lat, om_lon, &
14 HFX, QFX, GSW, GLW, UST, U_PHY, V_PHY, &
15 STBOLT, DELTSM, TSK, LH, XLAND, &
16 rdx, rdy, msfu, msfv, msft,xtime,om_tini,om_sini,id,omdt)
17 !----------------------------------------------------------------------------
18 implicit none
19 !----------------------------------------------------------------------------
20 ! Subroutine DPWP is a three-dimensional full physics ocean circulation model based
21 ! on Price et al. (1994, JPO). This ocean model is coupled with WRF
22 ! at every ocean time step. It is developed as a part of the University of Miami Coupled Model
23 ! (UMCM, Chen et al. 2012, JAS). A detailed description of the coupling process can
24 ! be found in Lee and Chen (2013, MWR).
25 !
26 ! This ocean model can be used with multi-nested grids the same manner as in WRF. It is recommended
27 ! to use grid spacing > 3 km for the ocean model while WRF can be configured at smaller
28 ! grid spacing then the ocean model.
29 !
30 ! The ocean model can be initialized with observed satellite SST or global model analyzed SST.
31 ! The and upper ocean temperature and sanility profiles can be initialized from climatology,
32 ! in-situ observations, and global ocean model analysis fields.
33 !
34 !-- om_tmp 3-dimensional ocean temperature (k)
35 !-- om_s 3-dimensional ocean salinity (psu)
36 !-- om_u 3-dimensional ocean current in x-direction (ms-1)
37 !-- om_v 3-dimensional ocean current in y-direction (ms-1)
38 !-- om_depth 3-dimensional ocean depth (m)
39 !-- om_ml ocean mixed layer depth (m)
40 !-- om_tini temperature reference profile based on ocean initial condition (k)
41 !-- om_sini salinity reference profile based on ocean initial condition (psu)
42 !-- om_dt ocean timestep (s)
43 !-------------------------------------------------------------------------------
45 !logical number to control the dynamics of ocean model
46 integer :: noasym, nonlcl, nodeep, nopg, idia
47 ! spatial variables
48 INTEGER , INTENT(IN) :: id
49 INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, &
50 ims,ime, jms,jme, kms,kme, &
51 its,ite, jts,jte, kts,kte
52 REAL , INTENT(IN ) :: rdx, &
53 rdy,xtime
54 integer :: i_start, i_end, j_start, j_end
55 !physical variables from WRF
56 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: msfu, &
57 msfv, &
58 msft
59 real :: mrdx, mrdy
60 real :: pscale, hascale, vascale
61 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) ::HFX, QFX, LH, GSW, GLW, UST,&
62 XLAND
63 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(INOUT) :: TSK
64 REAL, INTENT(IN ) :: STBOLT, DELTSM
65 REAL , DIMENSION( ims:ime ,kms:kme, jms:jme ),INTENT(IN ) :: U_PHY, V_PHY
66 !constant variables
67 real :: cp, crbc, cadv, ddx, g
68 real :: d2r, ro, rhoair, cpw, rb, rg, udrag, ucon
69 real , dimension (ims:ime, jms:jme) :: f, dml, tr, sr
70 real :: beta1,beta2
71 real , dimension (ims:ime, jms:jme) :: dr, alpha, beta
72 integer :: i, j, k, im, ip,jm,jp,km,kp,kk
73 integer :: ii, ntss
74 !ocean variables ::
75 integer, intent(in) :: okms, okme
76 real, dimension(ims:ime, okms:okme, jms:jme), INTENT(INOUT):: OM_TMP,OM_S,OM_U,OM_V,OM_DEPTH
77 real, dimension(ims:ime, okms:okme, jms:jme), INTENT(IN):: OM_TINI,OM_SINI
78 real, dimension(ims:ime, okms:okme, jms:jme), INTENT(INOUT):: OM_DENSITY
79 real, dimension(ims:ime, okms:okme, jms:jme):: OM_W
80 real, dimension(ims:ime, jms:jme),INTENT(INOUT):: OM_ML, OM_LAT, OM_LON
81 real, dimension(ims:ime, okms:okme, jms:jme):: tt, st, ut, vt, t9, s9, u9,v9
82 real, dimension(ims:ime, okms:okme, jms:jme):: pp, tentm, t0, s0, u0,v0
83 real, dimension(ims:ime, okms:okme, jms:jme):: dz, tham, tvam
84 real :: divu, divv, tha, sha, uha, vha, tva, vva, sva, uva, pgx, pgy
85 real, dimension(okms:okme) :: tb5, sb5, ub5, vb5, tb4, sb4, ub4, vb4
86 real, dimension(ims:ime,okms:okme, jms:jme) :: trr, srr
87 real, dimension(ims:ime, okms:okme, jms:jme)::bbeta, aalpha,sg
88 ! other variables
89 real :: dt, dti,dt7,vs, us,omdt
90 real :: dtdy, dtdx, dsdy, dsdx, dudy, dudx, dvdy, dvdx, wavg, &
91 dtdz, dsdz, dudz, dvdz
92 integer :: iupw, kvs, kus, ia, ja
94 i_start = max(its,ids+1)
95 i_end = min(ite,ide-2)
96 j_start = max(jts,jds+1)
97 j_end = min(jte,jde-2)
100 print*,'3DPWP domain',id,' run ocean'
102 dtdy = 0.0
103 dtdx = 0.0
104 dsdy = 0.0
105 dudy = 0.0
106 dudx = 0.0
107 dsdx = 0.0
108 dvdy = 0.0
109 dvdx = 0.0
110 wavg = 0.0
111 dtdz = 0.0
112 dsdz = 0.0
113 dudz = 0.0
114 dvdz = 0.0
115 divu = 0.0
116 divv = 0.0
117 tha = 0.0
118 sha = 0.0
119 uha = 0.0
120 vha = 0.0
121 tva = 0.0
122 vva = 0.0
123 sva = 0.0
124 uva = 0.0
125 pgx = 0.0
126 pgy = 0.0
127 !----------------------------------------------------------------
128 ! these flags control the dynamics of the model:
129 !
130 ! noasym = 1 shuts off the asymmetry of the storm winds
131 ! = 2 also shuts off radial wind component of storm winds
132 ! nonlcl = 1 shuts off horizontal advection,
133 ! = 2 shuts off horizontal and vertical advection,
134 ! = 3 above plus applies spatially constant wind stress
135 ! (useful for testing 1-d dynamics)
136 ! nodeep = 1 shuts off ml deepening by entrainment
137 ! nopg = 1 shuts off pressure gradients
138 ! idia = 1 turns on diagnostics
139 ! to calculate 1D pwp just choose nonlcl = 2, and nopg = 1
140 !
141 ! set flags to default values
142 noasym = 0
143 nonlcl = 0
144 nodeep = 0
145 nopg = 0
146 idia = 0
147 ! set some scales that are used for diagnostics - be careful with these !!
148 pscale = 1.0
149 if(nopg.eq.1) pscale = 0.
150 hascale = 1.0
151 if(nonlcl.ge.1) hascale = 0.
152 vascale = 1.0
153 if(nonlcl.ge.2) then
154 vascale = 0.
155 hascale = 0.
156 end if
157 !------------------------------------------------------------------
159 ! calculate the expansion coefficient
160 d2r = 2*3.14159/360.
161 g = 9.8
162 ro = 1.024e3 ! kg/m3
163 rhoair = 1
164 cpw = 4183.3
165 beta1 = .6
166 beta2 = 20.0
167 rb = 0.65
168 rg = 0.25
169 cp = 1.5
171 do i = i_start-1, i_end+1
172 do j = j_start-1, j_end+1
173 f (i,j)= 2.0*7.29e-5*sin(om_lat(i,j)*d2r)
174 if (om_tmp(i,1,j) .gt. 0.0) then
175 do k = 1,okme
176 trr(i,k,j) = om_tini(i,k,j)
177 srr(i,k,j) = om_sini(i,k,j)
178 aalpha(i,k,j) = sgt(trr(i,k,j)+0.5, srr(i,k,j), sg(i,k,j))- &
179 sgt(trr(i,k,j)-0.5,srr(i,k,j), sg(i,k,j))
180 bbeta(i,k,j) = sgt(trr(i,k,j), srr(i,k,j)+0.5, sg(i,k,j) )-sgt(trr(i,k,j),srr(i,k,j)-0.5, sg(i,k,j))
181 enddo
182 endif
183 enddo
184 enddo
186 udrag = 9999.
187 ucon = 0.
188 !set time step
189 dt =omdt*60.0
190 dti = 2.*dt
191 ntss =2
192 do 1111 ii = 1, ntss
193 if (ii .eq. 2) go to 88
194 do i= i_start-1, i_end+1
195 do j = j_start-1, j_end+1
196 do k = 1, okme
197 tt(i,k,j) = 0.0
198 ut(i,k,j) = 0.0
199 vt(i,k,j) = 0.0
200 st(i,k,j) = 0.0
201 t9(i,k,j) = om_tmp(i,k,j)
202 u9(i,k,j) = om_u(i,k,j)
203 v9(i,k,j) = om_v(i,k,j)
204 s9(i,k,j) = om_s(i,k,j)
205 t0(i,k,j) = om_tmp(i,k,j)
206 u0(i,k,j) = om_u(i,k,j)
207 v0(i,k,j) = om_v(i,k,j)
208 s0(i,k,j) = om_s(i,k,j)
209 om_w(i,k,j) = 0.0
210 pp(i,k,j) = 0.0
211 tham(i,k,j) = 0.0
212 tentm(i,k,j) = 0.0
213 enddo
214 enddo
215 enddo
217 go to 89
218 !second iteration
219 88 continue
220 dti = dt/2.
221 do i = i_start-1, i_end+1
222 do j = j_start-1, j_end+1
223 do k = 1, okme
224 t0(i,k,j) = t9(i,k,j)
225 u0(i,k,j) = u9(i,k,j)
226 v0(i,k,j) = v9(i,k,j)
227 s0(i,k,j) = s9(i,k,j)
228 om_w(i,k,j) = 0.0
229 pp(i,k,j) = 0.0
230 tham(i,k,j) = 0.0
231 tentm(i,k,j) = 0.0
232 enddo
233 enddo
234 enddo
236 89 continue
237 do 635 j = j_start, j_end
238 do 635 i = i_start, i_end
240 jp = j + 1
241 jm = j - 1
242 ip = i + 1
243 im = i - 1
245 if (om_tmp(ip,1,j) .le. 0 .or. om_tmp(im,1,j) .le. 0 .or. om_tmp(i,1,jp) .le. 0 .or. om_tmp(i,1,jm) .le.0 ) then
247 do k = 1, okme
248 om_w(i,k,j) = 0.0
249 enddo
250 else
252 do 6381 k = 1, okme
253 mrdx = rdx* msfu(i,j)
254 mrdy = rdy* msfv(i,j)
255 if (k .eq. 1) then
256 dz(i,k,j) = om_depth(i,k,j)
257 else
258 dz(i,k,j) = om_depth(i,k,j)-om_depth(i,k-1,j)
259 endif
260 divu = dz(i,k,j)*(om_u(ip,k,j)-om_u(im,k,j))*mrdx/2.0
261 divv = dz(i,k,j)*(om_v(i,k,jp)-om_v(i,k,jm))*mrdy/2.0
262 om_w(i,k,j) = divu+divv
263 6381 continue
264 do k = 2, okme
265 om_w(i,k,j) = om_w(i,k,j)+om_w(i,k-1,j)
266 enddo
268 endif
269 635 continue
270 do 6328 k = 1, okme
271 if (jts .eq. jds) then
272 do i = its, ite
273 om_w(i,k,jts) = om_w(i,k,jts+1)
274 enddo
275 endif
276 if (jte .eq. jde-1) then
277 do i = its, ite
278 om_w(i,k,jte) = om_w(i,k,jte-1)
279 enddo
280 endif
281 if (its .eq. ids) then
282 do j = jts, jte
283 om_w(its,k,j) = om_w(its+1,k,j)
284 enddo
285 endif
286 if (ite.eq. ide-1) then
287 do j = jts, jte
288 om_w(ite,k,j) = om_w(ite-1,k,j)
289 enddo
290 endif
291 if(its.eq. ids.and.jts .eq. jds) then
292 om_w(its,k,jts) = om_w(its+1,k,jts+1)
293 endif
294 if(its.eq. ids.and.jte .eq. jde-1) then
295 om_w(its,k,jte) = om_w(its+1,k,jte-1)
296 endif
297 if(ite.eq. ide-1.and.jte .eq. jde-1) then
298 om_w(ite,k,jte) = om_w(ite-1,k,jte-1)
299 endif
300 if(ite.eq. ide-1.and.jts .eq. jds) then
301 om_w(ite,k,jts) = om_w(ite-1,k,jts+1)
302 endif
305 6328 continue
306 ! compute the pressure perturbation by integrating
307 do 6327 i = i_start-1, i_end+1
308 do 6327 j = j_start-1, j_end+1
309 do k = 1, okme
310 if (k .eq. 1) then
311 dz(i,k,j) = om_depth(i,k,j)
312 else
313 dz(i,k,j) = om_depth(i,k,j)-om_depth(i,k-1,j)
314 endif
315 enddo
317 if (om_tmp(i,1,j) .gt. 0) then
318 pp (i,okme,j) = -g*dz(i,okme,j)*(aalpha(i,okme,j)*(om_tmp(i,okme,j)-trr(i,okme,j))+ &
319 bbeta(i,okme,j)*(om_s(i,okme,j)-srr(i,okme,j)))
321 do 632 k = 1,okme-1
322 kk = okme - k
323 pp(i,kk,j) = -g*dz(i,kk,j)*(aalpha(i,kk,j)*(om_tmp(i,kk,j)-trr(i,kk,j))+ &
324 bbeta(i,kk,j)*(om_s(i,kk,j)-srr(i,kk,j)))+pp(i,kk+1,j)
326 632 continue
327 endif
330 6327 continue
331 do 410 j = j_start, j_end
332 jp = j+1
333 jm = j-1
334 do 410 i = i_start, i_end
335 ip = i + 1
336 im = i - 1
337 do 420 k = 1 ,okme
338 mrdx = rdx* msfu(i,j)
339 mrdy = rdy* msfv(i,j)
340 tha = 0.
341 sha = 0.
342 uha = 0.
343 vha = 0.
345 tva = 0.
346 sva = 0.
347 vva = 0.
348 uva = 0.
350 pgx = 0.
351 pgy = 0.
352 ! estimate the tendencies due to horizontal advection
353 IUPW = 1
354 if (IUPW .eq. 1) then
355 vs = sign(1.,om_v(i,k,j))
356 kvs = -1* nint(vs)
357 us = sign(1.,om_u(i,k,j))
358 kus = -1* nint(us)
359 dtdy = vs*(om_tmp(i,k,j)-om_tmp(i,k,j+kvs))*mrdy
360 dtdx = us*(om_tmp(i,k,j)-om_tmp(i+kus,k,j))*mrdx
363 dsdy = vs*(om_s(i,k,j)-om_s(i,k,j+kvs))*mrdy
364 dsdx = us*(om_s(i,k,j)-om_s(i+kus,k,j))*mrdx
367 dudy = vs*(om_u(i,k,j)-om_u(i,k,j+kvs))*mrdy
368 dudx = us*(om_u(i,k,j)-om_u(i+kus,k,j))*mrdx
371 dvdy = vs*(om_v(i,k,j)-om_v(i,k,j+kvs))*mrdy
372 dvdx = us*(om_v(i,k,j)-om_v(i+kus,k,j))*mrdx
374 tha =-(om_v(i,k,j)*dtdy + om_u(i,k,j)*dtdx)
375 sha =-(om_v(i,k,j)*dsdy + om_u(i,k,j)*dsdx)
376 uha =-(om_v(i,k,j)*dudy + om_u(i,k,j)*dudx)
377 vha =-(om_v(i,k,j)*dvdy + om_u(i,k,j)*dvdx)
378 endif
379 ! estimate the tendencies due to vertical advection
380 ! skip if k = 1, assume that there is no vertical advection at surface
382 if ( k .eq. 1) go to 3323
383 wavg = 0.5*(om_w(i,k-1,j) + om_w(i,k,j))
385 km = k - 1
386 kp = k + 1
387 if (k .eq. okme) kp = k
388 dtdz = (om_tmp(i,km,j)-om_tmp(i,kp,j))/(om_depth(i,km,j)-om_depth(i,kp,j))
389 dsdz = (om_s(i,km,j)-om_s(i,kp,j))/(om_depth(i,km,j)-om_depth(i,kp,j))
390 dudz = (om_u(i,km,j)-om_u(i,kp,j))/(om_depth(i,km,j)-om_depth(i,kp,j))
391 dvdz = (om_v(i,km,j)-om_v(i,kp,j))/(om_depth(i,km,j)-om_depth(i,kp,j))
394 tva =wavg*dtdz
395 sva =wavg*dsdz
396 uva =wavg*dudz
397 vva =wavg*dvdz
398 3323 continue
401 !calculate pressure gradient
402 pgx = (pp(ip,k,j) - pp(im,k,j))*mrdx/2.0
403 pgy = (pp(i,k,jp) - pp(i,k,jm))*mrdy/2.0
405 if (( i .eq. its .and. its .eq. ids) .or. (i .eq. ite .and. ite .eq. ide-1) .or. &
406 ( j .eq. jts .and. jts .eq. jds) .or. (j .eq. jte .and. jte .eq. jde-1)) then
409 tt(i,k,j) = tt(i,k,j)
410 st(i,k,j) = st(i,k,j)
411 ut(i,k,j) = ut(i,k,j)
412 vt(i,k,j) = vt(i,k,j)
413 else
416 tt(i,k,j) = tt(i,k,j) + tha*hascale + tva*vascale
417 st(i,k,j) = st(i,k,j) + sha*hascale + sva*vascale
418 ut(i,k,j) = ut(i,k,j) + uha*hascale + uva*vascale - pgx*pscale/1023.0
419 vt(i,k,j) = vt(i,k,j) + vha*hascale + vva*vascale - pgy*pscale/1023.0
420 endif
421 420 continue
423 410 continue
427 do 700 k = 1, okme
428 if (jts .eq. jds)then
429 ja = 1
430 do 794 i = its, ite
431 crbc = -cp*rdy*msfv(i,j)
432 j = jts
433 tt(i,k,j) = tt(i,k,j) + crbc*(om_tmp(i,k,j) - om_tmp(i, k, j+ja))
434 st(i,k,j) = st(i,k,j) + crbc*(om_s(i,k,j) - om_s(i, k, j+ja))
435 ut(i,k,j) = ut(i,k,j) + crbc*(om_u(i,k,j) - om_u(i, k, j+ja))
436 vt(i,k,j) = vt(i,k,j) + crbc*(om_v(i,k,j) - om_v(i, k, j+ja))
437 794 continue
438 endif
439 ja = -1
440 if (jte .eq. jde-1)then
441 do 795 i = its, ite
442 crbc = -cp*rdy*msfv(i,j)
443 j = jte
444 tt(i,k,j) = tt(i,k,j) + crbc*(om_tmp(i,k,j) - om_tmp(i, k, j+ja))
445 st(i,k,j) = st(i,k,j) + crbc*(om_s(i,k,j) - om_s(i, k, j+ja))
446 ut(i,k,j) = ut(i,k,j) + crbc*(om_u(i,k,j) - om_u(i, k, j+ja))
447 vt(i,k,j) = vt(i,k,j) + crbc*(om_v(i,k,j) - om_v(i, k, j+ja))
448 795 continue
449 endif
450 if (its .eq. ids)then
451 do 796 j = jts , jte
452 ia = 1
453 crbc = -cp*(rdx*msfu(i,j))
454 i = its
455 tt(i,k,j) = tt(i,k,j) + crbc*(om_tmp(i,k,j) - om_tmp(i+ia, k, j))
456 st(i,k,j) = st(i,k,j) + crbc*(om_s(i,k,j) - om_s(i+ia, k, j))
457 ut(i,k,j) = ut(i,k,j) + crbc*(om_u(i,k,j) - om_u(i+ia, k, j))
458 vt(i,k,j) = vt(i,k,j) + crbc*(om_v(i,k,j) - om_v(i+ia, k, j))
459 796 continue
460 endif
461 if (ite .eq. ide-1)then
462 do 797 j = jts , jte
463 ia = -1
464 crbc = -cp*(rdx*msfu(i,j))
465 i = ite
466 tt(i,k,j) = tt(i,k,j) + crbc*(om_tmp(i,k,j) - om_tmp(i+ia, k, j))
467 st(i,k,j) = st(i,k,j) + crbc*(om_s(i,k,j) - om_s(i+ia, k, j))
468 ut(i,k,j) = ut(i,k,j) + crbc*(om_u(i,k,j) - om_u(i+ia, k, j))
469 vt(i,k,j) = vt(i,k,j) + crbc*(om_v(i,k,j) - om_v(i+ia, k, j))
470 797 continue
471 endif
473 if (its .eq. ids .and. jts .eq. jds) then
474 tt(its,k,jts) = tt(its+1,k,jts+1)
475 endif
476 if (its .eq. ids .and. jte .eq. jde-1) then
477 tt(its,k,jte) = tt(its+1,k,jte-1)
478 endif
479 if (ite .eq. ide-1 .and. jte .eq. jde-1) then
480 tt(ite,k,jte) = tt(ite-1,k,jte-1)
481 endif
482 if (ite .eq. ide-1 .and. jts .eq. jds) then
483 tt(ite,k,jts) = tt(ite-1,k,jts+1)
484 endif
486 700 enddo
490 !!end radiation bc
491 ! apply the wind stress and air-sea heat hluxes, and deepen the mixed layer by 1Dpwp
493 ! do 510 i = i_start, i_end
494 ! do 510 j = j_start, j_end
495 do 510 i = its, ite
496 do 510 j = jts, jte
497 if (ii .eq. 1) then
498 do k = 1, okme
499 tb5(k) = om_tmp(i,k,j)
500 sb5(k) = om_s(i,k,j)
501 ub5(k) = om_u(i,k,j)
502 vb5(k) = om_v(i,k,j)
503 tb4(k) = tb5(k)
504 sb4(k) = sb5(k)
505 ub4(k) = ub5(k)
506 vb4(k) = vb5(k)
507 enddo
508 endif
509 if (ii .eq. 2) then
510 do k = 1, okme
511 tb5(k) =t0(i,k,j)
512 sb5(k) =s0(i,k,j)
513 ub5(k) =u0(i,k,j)
514 vb5(k) =v0(i,k,j)
515 tb4(k) = tb5(k)
516 sb4(k) = sb5(k)
517 ub4(k) = ub5(k)
518 vb4(k) = vb5(k)
519 enddo
520 endif
522 ! vertical mixing
523 call PWP(i,j,k,ims, ime, jms, jme, okms,okme,tb4,&
524 sb4,om_density(i,:,j),om_depth(i,:,j),ub4,&
525 vb4,om_ml(i,j),om_lat(i,j),&
526 om_lon(i,j),HFX(i,j), QFX(i,j),LH(i,j), GSW(i,j),&
527 GLW(i,j), UST(i,j),U_PHY(i,kts,j),V_PHY(i,kts,j),&
528 STBOLT,DELTSM,aalpha(i,:,j),bbeta(i,:,j),om_tini(i,:,j),&
529 om_sini(i,:,j),trr(i,:,j),srr(i,:,j),omdt)
530 ! add mixing tendencies
531 do k = 1 , okme
532 tt(i,k,j) = tt(i,k,j) + (tb4(k)-tb5(k))/dt
533 st(i,k,j) = st(i,k,j) + (sb4(k)-sb5(k))/dt
534 ut(i,k,j) = ut(i,k,j) + (ub4(k)-ub5(k))/dt
535 vt(i,k,j) = vt(i,k,j) + (vb4(k)-vb5(k))/dt
538 enddo
540 510 continue
541 do 24 i = its, ite
542 do 24 j = jts, jte
543 do 24 k = 1, okme
544 om_tmp(i,k,j) = t0(i,k,j) + dti*tt(i,k,j)
545 om_s(i,k,j) = s0(i,k,j) + dti*st(i,k,j)
546 om_u(i,k,j) = u0(i,k,j) + dti*ut(i,k,j)
547 om_v(i,k,j) = v0(i,k,j) + dti*vt(i,k,j)
548 24 continue
550 do i = its, ite
551 do j = jts, jte
554 if (XLAND(i,j).GE.1.5)then
557 TSK(i,j) = om_tmp(i, 1, j)
558 endif
561 enddo
562 enddo
563 1111 enddo
564 return
566 END subroutine DPWP
567 ! cyl: end 3D PWP
568 !------------------------------------------------------------------
569 ! cyl: add 1D PWP
570 SUBROUTINE PWP(i,j,k,ims, ime, jms, jme, okms,okme,om_tmp,om_s, &
571 om_density, om_depth,om_u, om_v,om_ml,om_lat, om_lon, &
572 HFX, QFX,LH, GSW, GLW, UST,UAIR,VAIR,STBOLT,DELTSM,aalpha,bbeta,&
573 om_tini, om_sini,trr,srr,omdt)
574 implicit none
575 !----------------------------------------------------------------------------
576 ! Subroutine PWP is a one-dimensional vertical mixing scheme used in DPWP
577 ! based on on Price et al. (1994, JPO).
578 !----------------------------------------------------------------------------
581 integer i, j, k
582 integer ,intent(in) ::ims, ime, jms, jme, okms,okme !dimensions
583 integer ::ml,kml,mml
584 real, dimension(okms:okme), intent(inout):: om_tmp, om_s, &
585 om_density,om_depth
586 real, dimension(okms:okme),intent(inout):: om_u, om_v
588 real,intent(inout) :: om_ml
589 real,intent(inout) :: om_lat, om_lon
590 real TSK
591 real,intent(in) :: HFX, QFX,LH, GSW, GLW, ust,STBOLT,DELTSM
593 real, dimension(okms:okme)::absrb,u,v
594 real dt, dz
595 real :: tauxair, tauyair, taux, tauy, &
596 wspd
598 real d2r, g, hcon, f, tr, sr, dr, alpha, beta, energy, sg, &
599 ang,dml
601 real ro, rhoair, cpw
602 real beta1, beta2, ql, qi
603 real sipe
604 real rb,rv, rvc, delr, duv,rg,du,dv
605 real sa,ca,uair,vair,rc
606 real udrag, ucon
607 real, intent(in) :: omdt
608 real, dimension(okms:okme),intent(in)::aalpha,bbeta, om_tini,om_sini,trr,srr
611 dt = omdt*60.0
612 d2r = 2*3.14159/360.
613 g = 9.8
614 ro = 1.024e3 ! kg/m3
615 rhoair = 1
616 cpw = 4183.3
617 f = 2*7.29e-5*sin(om_lat*d2r)
618 beta1 = .6
619 beta2 = 20.0
620 rb = 0.65
621 rg = 0.25
622 dml = om_ml
623 tr = trr(1)
624 sr = srr(1)
625 dr = sgt(tr,sr,sg) ! unit of density is sigma
626 alpha = sgt(tr+0.5, sr, sg )-sgt(tr-0.5,sr, sg)
627 beta = sgt(tr, sr+0.5, sg )-sgt(tr,sr-0.5, sg)
628 udrag = 9999.
629 ucon = 0.
630 if (udrag .lt. 100) ucon = 1./(86400*udrag)
631 call absorb(i,j, ims, ime,jms,jme,okms,okme, beta1,beta2, &
632 absrb, om_depth,GSW,GLW)
633 ql = -1.*(LH+HFX+STBOLT*om_tmp(1)**4) ! corrected V3.8
634 qi = GLW+GSW
635 ql = ql+absrb(1)*qi
636 om_tmp(1) = om_tmp(1) + dt*ql/(om_depth(1)*ro*cpw)
637 om_s(1)=om_s(1)+om_s(1)*QFX*dt/ro/om_depth(1) ! corrected V3.8
638 do k = 2, okme
639 dz = om_depth(k)-om_depth(k-1)
640 om_tmp(k) = om_tmp(k) + dt*qi*absrb(k)/(dz*ro*cpw)
641 enddo
642 do k = 1, okme
643 om_density(k) = (dr+1000.0) +alpha*(om_tmp(k)-tr)+ &
644 beta*(om_s(k)-sr)
645 enddo
646 ml = 0
647 call mldep(i,j,k,ims,ime,jms,jme,okms,okme, om_density, &
648 om_depth, ml)
649 if (dml .gt. 100) print*,'mldep1',dml
650 call si (i, j, ims,ime,jms,jme,okms,okme,om_density, &
651 om_depth,sipe, ml)
653 if (ml .gt. 1) then
654 call mix1(om_tmp, ml, i, j,ims,ime,jms,jme,okms,okme,om_depth)
655 call mix1(om_s, ml, i, j,ims,ime,jms,jme,okms,okme,om_depth)
656 call mix1(om_u, ml, i, j,ims,ime,jms,jme,okms,okme,om_depth)
657 call mix1(om_v, ml, i, j,ims,ime,jms,jme,okms,okme,om_depth)
658 call mix1(om_density, ml, i, j,ims,ime,jms,jme,okms,okme,om_depth)
659 endif
660 ! above this, should be statically stable
662 ang = f*dt/2.
663 sa = sin(ang)
664 ca = cos(ang)
665 do k = 1, okme
666 call rot(i, j,ims,ime,jms,jme,okms,okme, om_u(k), om_v(k), sa, ca)
667 enddo
668 call mldep(i,j,k, ims,ime,jms,jme,okms,okme, om_density, &
669 om_depth, ml)
670 if (ml .eq. 1 ) then
671 om_ml = om_depth(ml)
672 else
673 om_ml = om_depth(ml)+(om_depth(ml) - om_depth(ml-1))/2.
674 endif
675 dml = om_ml
676 wspd = max(sqrt(uair*uair+vair*vair), 0.1)
677 tauxair = ust*ust*uair/wspd
678 taux = rhoair/ro*tauxair
679 tauyair = ust*ust*vair/wspd
680 tauy = rhoair/ro*tauyair
681 du = taux*dt/dml
682 dv = tauy*dt/dml
683 do k = 1,ml
684 om_u(k) = om_u(k) + du
685 om_v(k) = om_v(k) + dv
686 enddo
687 if(ucon .gt. 1.e-10) then
688 do k = 1,okme
689 om_u(k) = om_u(k)*(1. - dt*ucon)
690 om_v(k) = om_v(k)*(1. - dt*ucon)
691 enddo
692 endif
693 ! In pwp origional code, Dr. Price turns off the diffusion
694 ! To use it, check the diffusion parameter in subroutine first
695 ! call diffusion(okms,okme,dt,om_tmp,om_depth)
696 ! call diffusion(okms,okme,dt,om_s,om_depth)
697 ! call diffusion(okms,okme,dt,om_u,om_depth)
698 ! call diffusion(okms,okme,dt,om_v,om_depth)
699 ! call diffusion(okms,okme,dt,om_density,om_depth)
700 ! print*,'ca,sa',ca,sa
701 if(rb .le. 0.00001) go to 50
702 rvc = rb
703 do 40 k = ml, okme-1
704 if (k .eq. 1) then
705 dz = om_depth(k)
706 else
707 dz = om_depth(k) + (om_depth(k) - om_depth(k-1))/2.
708 endif
709 delr = (om_density(k+1)-om_density(k))/ro
710 duv =(om_u(k+1)-om_u(k))**2+(om_v(k+1) &
711 -om_v(k))**2 + 1.e-8
712 rv = g*delr*dz/duv
713 if (rv .gt. rvc) go to 41
714 mml=k+1
715 call mix1(om_tmp, mml, i, j,ims,ime,jms,jme,okms,okme,om_depth)
716 call mix1(om_s, mml, i, j,ims,ime,jms,jme,okms,okme,om_depth)
717 call mix1(om_u, mml, i, j,ims,ime,jms,jme,okms,okme,om_depth)
718 call mix1(om_v, mml, i, j,ims,ime,jms,jme,okms,okme,om_depth)
719 call mix1(om_density, mml, i, j,ims,ime,jms,jme,okms,okme,om_depth)
720 40 continue
721 41 continue
722 50 continue
723 if (rg .gt. 1.e-4) call rimix(i,j,k,ims,ime,jms,jme,okms, &
724 okme,rc, &
725 rg,om_density,om_depth, om_u, om_v, om_s, om_tmp)
726 do k = 1, okme
727 call rot(i,j,ims,ime,jms,jme,okms,okme, om_u(k), om_v(k), &
728 sa, ca)
729 enddo
730 call mldep(i,j,k, ims,ime,jms,jme,okms,okme, om_density, &
731 om_depth, ml)
732 if (ml .eq. 1) then
733 om_ml = om_depth(ml)
734 else
735 om_ml = om_depth(ml)+ (om_depth(ml) - om_depth(ml-1))/2.
736 endif
737 ! TSK = om_tmp(1)
738 return
739 END SUBROUTINE PWP
740 ! cyl : end PWP
741 ! cyl : PWP's subrouting
742 ! cyl--------------------
743 SUBROUTINE absorb(i,j, ims, ime, jms, jme, okms,okme,beta1, &
744 beta2, absrb, om_depth,GSW,GLW)
745 implicit none
746 integer i,j,k
747 integer ims, ime, jms, jme, okms, okme
748 real, dimension(okms:okme),intent(inout)::absrb, om_depth
749 real abss, rs1, rs2, beta1, beta2, z1b1, z2b1, z2b2, z1b2,&
750 z1,z2,GSW,GLW
751 abss = 0.
752 rs1 = GSW/(GSW+GLW)
753 rs2 = 1.0 - rs1
754 do k = 1, okme
755 absrb(k) = 0.
756 enddo
757 do k =1, okme
758 if (k .eq. 1) then
759 z1 = 0
760 else
761 z1 = om_depth(k-1)
762 endif
763 z2 = om_depth(k)
764 z1b1 = z1/beta1
765 z2b1 = z2/beta1
766 z1b2 = z1/beta2
767 z2b2 = z2/beta2
768 if (z2b1 .lt. 70) absrb(k) = absrb(k) + &
769 rs1*(exp(-z1b1)-exp(-z2b1))
770 if (z2b2 .lt. 70) absrb(k) = absrb(k) + &
771 rs2*(exp(-z1b2)-exp(-z2b2))
772 abss = abss + absrb(k)
773 enddo
774 return
775 END SUBROUTINE absorb
776 !cyl--------------------
777 SUBROUTINE mldep (i,j,k, ims,ime,jms,jme,okms,okme, om_density, &
778 om_depth, ml)
779 integer i,j,k,kk, ims,ime,jms,jme,okms,okme
780 real deps, dml, dd
781 integer ml
782 real, dimension(okms:okme),intent(inout)::om_density,om_depth
783 deps = 1.e-5
784 do 70 k = 1, okme-1
785 dd = abs(om_density(k+1) - om_density(k))
786 if (dd .gt. deps) go to 71
787 70 continue
788 71 continue
789 ml = k
790 return
791 END SUBROUTINE mldep
792 !cyl--------------------
793 SUBROUTINE si(i, j, ims,ime,jms,jme,okms,okme,om_density, &
794 om_depth, &
795 sipe, ml)
796 integer i, j, k, ims,ime,jms,jme,okms,okme, mml,ml
797 real, dimension(okms:okme),intent(inout)::om_density, om_depth
798 real, dimension(okms:okme) :: temp_density
799 real pesum, sipe
800 do k = 1, okme
801 temp_density(k) = om_density(k)
802 enddo
803 mml = 1
804 do 80 k = 2, okme
805 if (om_density(k) .gt. om_density(k-1)) go to 81
806 if (om_density(k) .lt. om_density(k-1))then
807 mml = k
808 call mix1(om_density, mml, i, j,ims,ime,jms,jme,okms,okme,om_depth)
809 endif
810 80 continue
811 81 continue
812 ml = mml
814 return
815 END SUBROUTINE si
816 !cyl--------------------
817 SUBROUTINE mix1(b,mml, i, j,ims,ime,jms,jme,okms,okme,om_depth)
818 integer i, j, k, ims,ime,jms,jme,okms,okme, mml
819 real,dimension(okms:okme),intent(inout) :: b
820 real,dimension(okms:okme),intent(in) :: om_depth
821 real,dimension(okms:okme) :: dz
822 real bs,zs
823 do k = 1, okme
824 if (k .eq. 1) then
825 dz(k) = om_depth(k)
826 else
827 dz(k) = om_depth(k)-om_depth(k-1)
828 endif
829 enddo
830 bs = 0.
831 zs = 0.
832 do k = 1, mml
833 bs = bs + b(k)*dz(k)
834 zs = zs + dz(k)
835 enddo
836 bs = bs/zs
837 do k = 1, mml
838 b(k) = bs
839 enddo
840 return
841 END SUBROUTINE mix1
842 !cyl--------------------
843 SUBROUTINE rot(i, j,ims,ime,jms,jme,okms,okme, om_u, om_v, sa, ca)
844 integer i, j, k, ims,ime,jms,jme,okms,okme
845 real,intent(inout)::om_u, om_v
846 real:: u0, v0
847 real sa, ca
848 u0= om_u
849 v0= om_v
850 om_u = u0*ca+v0*sa
851 om_v = v0*ca-u0*sa
852 return
853 END SUBROUTINE rot
854 !cyl--------------------
855 SUBROUTINE rimix(i,j,k,ims,ime,jms,jme,okms,okme,rc, rg, &
856 om_density,om_depth, om_u, om_v, om_s, om_tmp)
857 implicit none
858 integer i,j,k,ims,ime,jms,jme,okms,okme,kcd,nmix,k1,k2,ks
859 real rc, rg,dd,g,dz,dv,ro,rs
860 real,dimension(okms:okme) :: r
861 real, dimension(okms:okme),intent(inout)::om_density, om_u, &
862 om_depth,om_v, om_s, om_tmp
863 ro = 1000.
864 rc = rg
865 g = 9.8
866 kcd = 0
867 nmix = 0
868 k1 = 1
869 k2 = okme - 1
872 10 continue
873 do 1 k = k1, k2
874 dd = om_density(k+1)-om_density(k)
875 if (dd .lt. 1.e-3) dd = 1.e-3
876 dv = (om_u(k+1)-om_u(k))**2 + (om_v(k+1)-om_v(k))**2
877 if(dv .lt. 1.e-6) dv = 1.e-6
878 if (k .eq. 1) then
879 dz = (om_depth(1) + om_depth(2) - om_depth(1))/2.
880 else
881 dz = (om_depth(k)-om_depth(k-1) + om_depth(k+1) - om_depth(k))/2.
882 endif
883 r(k) = g*dz*dd/(dv*ro)
884 1 continue
885 rs = r(1)
886 ks = 1
887 do 2 k = 2, okme-1
888 if (r(k) .lt. rs) go to 3
889 go to 2
890 3 continue
891 rs =r(k)
892 ks = k
893 2 continue
894 if (rs .ge. rc) go to 99
895 if (ks .ge. kcd ) kcd = ks + 1
896 call stir(i,j,ims,ime,jms,jme,okms,okme,rc,rs, om_tmp, ks,om_depth)
897 call stir(i,j,ims,ime,jms,jme,okms,okme,rc,rs, om_s, ks,om_depth)
898 call stir(i,j,ims,ime,jms,jme,okms,okme,rc,rs, om_density, ks,om_depth)
899 call stir(i,j,ims,ime,jms,jme,okms,okme,rc,rs, om_u, ks,om_depth)
900 call stir(i,j,ims,ime,jms,jme,okms,okme,rc,rs, om_v, ks,om_depth)
901 nmix = nmix +1
902 k1 = ks - 2
903 if (k1 .lt. 1) k1 = 1
904 k2 = ks +2
905 if (k2 .gt. okme-1) k2 = okme -1
906 go to 10
907 99 continue
910 return
911 END SUBROUTINE rimix
912 !cyl--------------------
913 SUBROUTINE diffusion(okms,okme,dt,a,om_depth)
914 implicit none
915 integer okms,okme,k
916 real dt,rkz,dz
917 real, dimension(okms:okme):: a, work,dconst,om_depth
918 rkz =0.3
919 do k = okms, okme
920 if (k .eq. 1) then
921 dz = om_depth(1)
922 else
923 dz = om_depth(k)-om_depth(k-1)
924 endif
925 dconst(k) = dt*rkz/dz**2
926 enddo
927 do k = 2, okme-1
928 work(k) = dconst(k)*(a(k-1)+a(k+1)-2.*a(k))
929 enddo
930 do k = 2, okme-1
931 a(k) = a(k) + work(k)
932 enddo
933 return
934 END SUBROUTINE diffusion
935 !cyl--------------------
936 SUBROUTINE stir(i,j,ims,ime,jms,jme,okms,okme,rc,rs,a,ks,depth)
937 implicit none
938 integer i,j,k,ims,ime,jms,jme,okms,okme,ks,b
939 real, dimension(okms:okme),intent(inout)::a
940 real, dimension(okms:okme),intent(in)::depth
941 real, dimension(okms:okme)::dz
942 real rcon,rc,rs,ff,rnew,da
943 rnew = 0.3
944 ff = 1. - rs/rnew
945 do k = 1,okme
946 if (k.eq. 1) then
947 dz(k) = depth(1)
948 else
949 dz(k) = depth(k) -depth(k-1)
950 endif
951 enddo
952 da = (a(ks+1)- a(ks))*ff/2.
953 b=2./(1.+dz(ks)/dz(ks+1))
954 a(ks+1)= a(ks+1)-da*b*dz(ks)/dz(ks+1)
955 a(ks) = a(ks)+da*b
956 ! print*,'da',da
957 ! print*,'a',a(ks+1),a(ks)
958 return
959 END SUBROUTINE stir
960 !cyl--------------------
961 FUNCTION sg0(s)
962 implicit none
963 real s,sg0
964 sg0 = ((6.76786136e-6*s-4.8249614e-4)*s+0.814876577)*s &
965 -0.0934458632
966 return
967 END FUNCTION sg0
968 FUNCTION sgt(tt,s,sg)
969 implicit none
970 real tt,t, s, sg,sgt
971 t=tt-273.0
972 sg = sg0(s)
973 20 sgt = ((((-1.43803061e-7*t-1.98248399e-3)*t-0.545939111)*t &
974 +4.53168426)*t)/(t+67.26)+((((1.667e-8*t-8.164e-7)*t &
975 +1.803e-5)*t)*sg+((-1.0843e-6*t+9.8185e-5)*t-4.7867e-3)*t &
976 +1.0)*sg
977 return
978 END FUNCTION sgt
980 !--------------------
981 END MODULE module_sf_3dpwp