ROMS
Loading...
Searching...
No Matches
rho_eos.F
Go to the documentation of this file.
1#include "cppdefs.h"
2 MODULE rho_eos_mod
3#ifdef SOLVE3D
4!
5!git $Id$
6!================================================== Hernan G. Arango ===
7! Copyright (c) 2002-2025 The ROMS Group !
8! Licensed under a MIT/X style license !
9! See License_ROMS.md !
10!=======================================================================
11! !
12! This routine computes "in situ" density and other associated !
13! quantitites as a function of potential temperature, salinity, !
14! and pressure from a polynomial expression (Jackett and McDougall, !
15! 1992). The polynomial expression was found from fitting to 248 !
16! values in the oceanographic ranges of salinity, potential !
17! temperature, and pressure. It assumes no pressure variation !
18! along geopotential surfaces, that is, depth (meters; negative) !
19! and pressure (dbar; assumed negative here) are interchangeable. !
20! !
21! Check Values: (T=3 C, S=35.5 PSU, Z=-5000 m) !
22! !
23! alpha = 2.1014611551470d-04 (1/Celsius) !
24! beta = 7.2575037309946d-04 (1/PSU) !
25! gamma = 3.9684764511766d-06 (1/Pa) !
26! den = 1050.3639165364 (kg/m3) !
27! den1 = 1028.2845117925 (kg/m3) !
28! sound = 1548.8815240223 (m/s) !
29! bulk = 23786.056026320 (Pa) !
30! !
31! Reference: !
32! !
33! Jackett, D. R. and T. J. McDougall, 1995, Minimal Adjustment of !
34! Hydrostatic Profiles to Achieve Static Stability, J. of Atmos. !
35! and Oceanic Techn., vol. 12, pp. 381-389. !
36! !
37!=======================================================================
38!
39 implicit none
40!
41 PRIVATE
42 PUBLIC :: rho_eos
43!
44 CONTAINS
45!
46!***********************************************************************
47 SUBROUTINE rho_eos (ng, tile, model)
48!***********************************************************************
49!
50 USE mod_param
51 USE mod_coupling
52 USE mod_grid
53 USE mod_mixing
54 USE mod_ocean
55 USE mod_stepping
56!
57! Imported variable declarations.
58!
59 integer, intent(in) :: ng, tile, model
60!
61! Local variable declarations.
62!
63 character (len=*), parameter :: myfile = &
64 & __FILE__
65!
66# include "tile.h"
67!
68# ifdef PROFILE
69 CALL wclock_on (ng, model, 14, __line__, myfile)
70# endif
71 CALL rho_eos_tile (ng, tile, model, &
72 & lbi, ubi, lbj, ubj, &
73 & imins, imaxs, jmins, jmaxs, &
74 & nrhs(ng), &
75# ifdef MASKING
76 & grid(ng) % rmask, &
77# endif
78# ifdef VAR_RHO_2D
79 & grid(ng) % Hz, &
80# endif
81 & grid(ng) % z_r, &
82 & grid(ng) % z_w, &
83 & ocean(ng) % t, &
84# ifdef VAR_RHO_2D
85 & coupling(ng) % rhoA, &
86 & coupling(ng) % rhoS, &
87# endif
88# ifdef BV_FREQUENCY
89 & mixing(ng) % bvf, &
90# endif
91# if defined LMD_SKPP || defined LMD_BKPP || \
92 defined bulk_fluxes || defined balance_operator
93 & mixing(ng) % alpha, &
94 & mixing(ng) % beta, &
95# ifdef LMD_DDMIX
96 & mixing(ng) % alfaobeta, &
97# endif
98# endif
99 & ocean(ng) % pden, &
100 & ocean(ng) % rho)
101# ifdef PROFILE
102 CALL wclock_off (ng, model, 14, __line__, myfile)
103# endif
104!
105 RETURN
106 END SUBROUTINE rho_eos
107
108# ifdef NONLIN_EOS
109!
110!***********************************************************************
111 SUBROUTINE rho_eos_tile (ng, tile, model, &
112 & LBi, UBi, LBj, UBj, &
113 & IminS, ImaxS, JminS, JmaxS, &
114 & nrhs, &
115# ifdef MASKING
116 & rmask, &
117# endif
118# ifdef VAR_RHO_2D
119 & Hz, &
120# endif
121 & z_r, z_w, t, &
122# ifdef VAR_RHO_2D
123 & rhoA, rhoS, &
124# endif
125# ifdef BV_FREQUENCY
126 & bvf, &
127# endif
128# if defined LMD_SKPP || defined LMD_BKPP || \
129 defined BULK_FLUXES || defined BALANCE_OPERATOR
130 & alpha, beta, &
131# ifdef LMD_DDMIX
132 & alfaobeta, &
133# endif
134# endif
135 & pden, &
136 & rho)
137!***********************************************************************
138!
139 USE mod_param
140 USE mod_eoscoef
141 USE mod_scalars
142# ifdef SEDIMENT
143 USE mod_sediment
144# endif
145!
146 USE exchange_2d_mod
147 USE exchange_3d_mod
148# ifdef DISTRIBUTE
149 USE mp_exchange_mod, ONLY : mp_exchange2d, mp_exchange3d
150# endif
151!
152! Imported variable declarations.
153!
154 integer, intent(in) :: ng, tile, model
155 integer, intent(in) :: LBi, UBi, LBj, UBj
156 integer, intent(in) :: IminS, ImaxS, JminS, JmaxS
157 integer, intent(in) :: nrhs
158!
159# ifdef ASSUMED_SHAPE
160# ifdef MASKING
161 real(r8), intent(in) :: rmask(LBi:,LBj:)
162# endif
163# ifdef VAR_RHO_2D
164 real(r8), intent(in) :: Hz(LBi:,LBj:,:)
165# endif
166 real(r8), intent(in) :: z_r(LBi:,LBj:,:)
167 real(r8), intent(in) :: z_w(LBi:,LBj:,0:)
168 real(r8), intent(in) :: t(LBi:,LBj:,:,:,:)
169# ifdef VAR_RHO_2D
170 real(r8), intent(out) :: rhoA(LBi:,LBj:)
171 real(r8), intent(out) :: rhoS(LBi:,LBj:)
172# endif
173# ifdef BV_FREQUENCY
174 real(r8), intent(out) :: bvf(LBi:,LBj:,0:)
175# endif
176# if defined LMD_SKPP || defined LMD_BKPP || \
177 defined bulk_fluxes || defined balance_operator
178 real(r8), intent(out) :: alpha(LBi:,LBj:)
179 real(r8), intent(out) :: beta(LBi:,LBj:)
180# ifdef LMD_DDMIX
181 real(r8), intent(out) :: alfaobeta(LBi:,LBj:,0:)
182# endif
183# endif
184 real(r8), intent(out) :: pden(LBi:,LBj:,:)
185 real(r8), intent(out) :: rho(LBi:,LBj:,:)
186# else
187# ifdef MASKING
188 real(r8), intent(in) :: rmask(LBi:UBi,LBj:UBj)
189# endif
190# ifdef VAR_RHO_2D
191 real(r8), intent(in) :: Hz(LBi:UBi,LBj:UBj,N(ng))
192# endif
193 real(r8), intent(in) :: z_r(LBi:UBi,LBj:UBj,N(ng))
194 real(r8), intent(in) :: z_w(LBi:UBi,LBj:UBj,0:N(ng))
195 real(r8), intent(in) :: t(LBi:UBi,LBj:UBj,N(ng),3,NT(ng))
196# ifdef VAR_RHO_2D
197 real(r8), intent(out) :: rhoA(LBi:UBi,LBj:UBj)
198 real(r8), intent(out) :: rhoS(LBi:UBi,LBj:UBj)
199# endif
200# ifdef BV_FREQUENCY
201 real(r8), intent(out) :: bvf(LBi:UBi,LBj:UBj,0:N(ng))
202# endif
203# if defined LMD_SKPP || defined LMD_BKPP || \
204 defined bulk_fluxes || defined balance_operator
205 real(r8), intent(out) :: alpha(LBi:UBi,LBj:UBj)
206 real(r8), intent(out) :: beta(LBi:UBi,LBj:UBj)
207# ifdef LMD_DDMIX
208 real(r8), intent(out) :: alfaobeta(LBi:UBi,LBj:UBj,0:N(ng))
209# endif
210# endif
211 real(r8), intent(out) :: pden(LBi:UBi,LBj:UBj,N(ng))
212 real(r8), intent(out) :: rho(LBi:UBi,LBj:UBj,N(ng))
213# endif
214!
215! Local variable declarations.
216!
217 integer :: i, ised, itrc, j, k
218
219 real(r8) :: SedDen, Tp, Tpr10, Ts, Tt, sqrtTs
220# ifdef BV_FREQUENCY
221 real(r8) :: bulk_dn, bulk_up, den_dn, den_up
222# endif
223 real(r8) :: cff, cff1, cff2
224
225 real(r8), dimension(0:9) :: C
226# ifdef EOS_TDERIVATIVE
227 real(r8), dimension(0:9) :: dCdT(0:9)
228
229 real(r8), dimension(IminS:ImaxS,N(ng)) :: DbulkDS
230 real(r8), dimension(IminS:ImaxS,N(ng)) :: DbulkDT
231 real(r8), dimension(IminS:ImaxS,N(ng)) :: Dden1DS
232 real(r8), dimension(IminS:ImaxS,N(ng)) :: Dden1DT
233 real(r8), dimension(IminS:ImaxS,N(ng)) :: Scof
234 real(r8), dimension(IminS:ImaxS,N(ng)) :: Tcof
235 real(r8), dimension(IminS:ImaxS,N(ng)) :: wrk
236# endif
237 real(r8), dimension(IminS:ImaxS,N(ng)) :: bulk
238 real(r8), dimension(IminS:ImaxS,N(ng)) :: bulk0
239 real(r8), dimension(IminS:ImaxS,N(ng)) :: bulk1
240 real(r8), dimension(IminS:ImaxS,N(ng)) :: bulk2
241 real(r8), dimension(IminS:ImaxS,N(ng)) :: den
242 real(r8), dimension(IminS:ImaxS,N(ng)) :: den1
243
244# include "set_bounds.h"
245!
246!=======================================================================
247! Nonlinear equation of state. Notice that this equation of state
248! is only valid for potential temperature range of -2C to 40C and
249! a salinity range of 0 PSU to 42 PSU.
250!=======================================================================
251!
252 DO j=jstrt,jendt
253 DO k=1,n(ng)
254 DO i=istrt,iendt
255!
256! Check temperature and salinity lower values. Assign depth to the
257! pressure.
258!
259 tt=max(-2.0_r8,t(i,j,k,nrhs,itemp))
260# ifdef SALINITY
261 ts=max(0.0_r8,t(i,j,k,nrhs,isalt))
262 sqrtts=sqrt(ts)
263# else
264 ts=0.0_r8
265 sqrtts=0.0_r8
266# endif
267 tp=z_r(i,j,k)
268 tpr10=0.1_r8*tp
269!
270!-----------------------------------------------------------------------
271! Compute density (kg/m3) at standard one atmosphere pressure.
272!-----------------------------------------------------------------------
273!
274 c(0)=q00+tt*(q01+tt*(q02+tt*(q03+tt*(q04+tt*q05))))
275 c(1)=u00+tt*(u01+tt*(u02+tt*(u03+tt*u04)))
276 c(2)=v00+tt*(v01+tt*v02)
277# ifdef EOS_TDERIVATIVE
278!
279 dcdt(0)=q01+tt*(2.0_r8*q02+tt*(3.0_r8*q03+tt*(4.0_r8*q04+ &
280 & tt*5.0_r8*q05)))
281 dcdt(1)=u01+tt*(2.0_r8*u02+tt*(3.0_r8*u03+tt*4.0_r8*u04))
282 dcdt(2)=v01+tt*2.0_r8*v02
283# endif
284!
285 den1(i,k)=c(0)+ts*(c(1)+sqrtts*c(2)+ts*w00)
286
287# ifdef EOS_TDERIVATIVE
288!
289! Compute d(den1)/d(S) and d(den1)/d(T) derivatives used in the
290! computation of thermal expansion and saline contraction
291! coefficients.
292!
293 dden1ds(i,k)=c(1)+1.5_r8*c(2)*sqrtts+2.0_r8*w00*ts
294 dden1dt(i,k)=dcdt(0)+ts*(dcdt(1)+sqrtts*dcdt(2))
295# endif
296!
297!-----------------------------------------------------------------------
298! Compute secant bulk modulus.
299!-----------------------------------------------------------------------
300!
301 c(3)=a00+tt*(a01+tt*(a02+tt*(a03+tt*a04)))
302 c(4)=b00+tt*(b01+tt*(b02+tt*b03))
303 c(5)=d00+tt*(d01+tt*d02)
304 c(6)=e00+tt*(e01+tt*(e02+tt*e03))
305 c(7)=f00+tt*(f01+tt*f02)
306 c(8)=g01+tt*(g02+tt*g03)
307 c(9)=h00+tt*(h01+tt*h02)
308# ifdef EOS_TDERIVATIVE
309!
310 dcdt(3)=a01+tt*(2.0_r8*a02+tt*(3.0_r8*a03+tt*4.0_r8*a04))
311 dcdt(4)=b01+tt*(2.0_r8*b02+tt*3.0_r8*b03)
312 dcdt(5)=d01+tt*2.0_r8*d02
313 dcdt(6)=e01+tt*(2.0_r8*e02+tt*3.0_r8*e03)
314 dcdt(7)=f01+tt*2.0_r8*f02
315 dcdt(8)=g02+tt*2.0_r8*g03
316 dcdt(9)=h01+tt*2.0_r8*h02
317# endif
318!
319 bulk0(i,k)=c(3)+ts*(c(4)+sqrtts*c(5))
320 bulk1(i,k)=c(6)+ts*(c(7)+sqrtts*g00)
321 bulk2(i,k)=c(8)+ts*c(9)
322 bulk(i,k)=bulk0(i,k)-tp*(bulk1(i,k)-tp*bulk2(i,k))
323
324# if defined LMD_SKPP || defined LMD_BKPP || \
325 defined bulk_fluxes || defined balance_operator
326!
327! Compute d(bulk)/d(S) and d(bulk)/d(T) derivatives used
328! in the computation of thermal expansion and saline contraction
329! coefficients.
330!
331 dbulkds(i,k)=c(4)+sqrtts*1.5_r8*c(5)- &
332 & tp*(c(7)+sqrtts*1.5_r8*g00-tp*c(9))
333 dbulkdt(i,k)=dcdt(3)+ts*(dcdt(4)+sqrtts*dcdt(5))- &
334 & tp*(dcdt(6)+ts*dcdt(7)- &
335 & tp*(dcdt(8)+ts*dcdt(9)))
336# endif
337!
338!-----------------------------------------------------------------------
339! Compute local "in situ" density anomaly (kg/m3 - 1000).
340!-----------------------------------------------------------------------
341!
342 cff=1.0_r8/(bulk(i,k)+tpr10)
343 den(i,k)=den1(i,k)*bulk(i,k)*cff
344# if defined SEDIMENT && defined SED_DENS
345 sedden=0.0_r8
346 DO ised=1,nst
347 itrc=idsed(ised)
348 cff1=1.0_r8/srho(ised,ng)
349 sedden=sedden+ &
350 & t(i,j,k,nrhs,itrc)* &
351 & (srho(ised,ng)-den(i,k))*cff1
352 END DO
353 den(i,k)=den(i,k)+sedden
354# endif
355 den(i,k)=den(i,k)-1000.0_r8
356# ifdef MASKING
357 den(i,k)=den(i,k)*rmask(i,j)
358# endif
359 END DO
360 END DO
361
362# ifdef VAR_RHO_2D
363!
364!-----------------------------------------------------------------------
365! Compute vertical averaged density (rhoA) and density perturbation
366! (rhoS) used in barotropic pressure gradient. Both values are
367! nondimensional.
368!-----------------------------------------------------------------------
369!
370 DO i=istrt,iendt
371 cff1=den(i,n(ng))*hz(i,j,n(ng))
372 rhos(i,j)=0.5_r8*cff1*hz(i,j,n(ng))
373 rhoa(i,j)=cff1
374 END DO
375 DO k=n(ng)-1,1,-1
376 DO i=istrt,iendt
377 cff1=den(i,k)*hz(i,j,k)
378 rhos(i,j)=rhos(i,j)+hz(i,j,k)*(rhoa(i,j)+0.5_r8*cff1)
379 rhoa(i,j)=rhoa(i,j)+cff1
380 END DO
381 END DO
382 cff2=1.0_r8/rho0
383 DO i=istrt,iendt
384 cff1=1.0_r8/(z_w(i,j,n(ng))-z_w(i,j,0))
385 rhoa(i,j)=cff2*cff1*rhoa(i,j)
386 rhos(i,j)=2.0_r8*cff1*cff1*cff2*rhos(i,j)
387 END DO
388# endif
389
390# if defined BV_FREQUENCY
391!
392!-----------------------------------------------------------------------
393! Compute Brunt-Vaisala frequency (1/s2) at horizontal RHO-points
394! and vertical W-points:
395!
396! bvf = - g/rho d(rho)/d(z).
397!
398! The density anomaly difference is computed by lowering/rising the
399! water parcel above/below adiabatically at W-point depth "z_w".
400!-----------------------------------------------------------------------
401!
402 DO k=1,n(ng)-1
403 DO i=istrt,iendt
404 bulk_up=bulk0(i,k+1)- &
405 & z_w(i,j,k)*(bulk1(i,k+1)- &
406 & bulk2(i,k+1)*z_w(i,j,k))
407 bulk_dn=bulk0(i,k )- &
408 & z_w(i,j,k)*(bulk1(i,k )- &
409 & bulk2(i,k )*z_w(i,j,k))
410 cff1=1.0_r8/(bulk_up+0.1_r8*z_w(i,j,k))
411 cff2=1.0_r8/(bulk_dn+0.1_r8*z_w(i,j,k))
412 den_up=cff1*(den1(i,k+1)*bulk_up)
413 den_dn=cff2*(den1(i,k )*bulk_dn)
414 bvf(i,j,k)=-g*(den_up-den_dn)/ &
415 & (0.5_r8*(den_up+den_dn)* &
416 & (z_r(i,j,k+1)-z_r(i,j,k)))
417 END DO
418 END DO
419 DO i=istrt,iendt
420!! bvf(i,j,0)=bvf(i,j,1)
421!! bvf(i,j,N(ng))=bvf(i,j,N(ng)-1)
422 bvf(i,j,0)=0.0_r8
423 bvf(i,j,n(ng))=0.0_r8
424 END DO
425# endif
426
427# if defined LMD_SKPP || defined LMD_BKPP || \
428 defined bulk_fluxes || defined balance_operator
429!
430!-----------------------------------------------------------------------
431! Compute thermal expansion (1/Celsius) and saline contraction
432! (1/PSU) coefficients.
433!-----------------------------------------------------------------------
434!
435# ifdef LMD_DDMIX
436 DO k=1,n(ng)
437# else
438 DO k=n(ng),n(ng)
439# endif
440 DO i=istrt,iendt
441 tpr10=0.1_r8*z_r(i,j,k)
442!
443! Compute thermal expansion and saline contraction coefficients.
444!
445 cff=bulk(i,k)+tpr10
446 cff1=tpr10*den1(i,k)
447 cff2=bulk(i,k)*cff
448 wrk(i,k)=(den(i,k)+1000.0_r8)*cff*cff
449 tcof(i,k)=-(dbulkdt(i,k)*cff1+ &
450 & dden1dt(i,k)*cff2)
451 scof(i,k)= (dbulkds(i,k)*cff1+ &
452 & dden1ds(i,k)*cff2)
453# ifdef LMD_DDMIX
454 alfaobeta(i,j,k)=tcof(i,k)/scof(i,k)
455# endif
456 END DO
457 IF (k.eq.n(ng)) THEN
458 DO i=istrt,iendt
459 cff=1.0_r8/wrk(i,n(ng))
460 alpha(i,j)=cff*tcof(i,n(ng))
461 beta(i,j)=cff*scof(i,n(ng))
462 END DO
463 END IF
464 END DO
465# endif
466!
467!-----------------------------------------------------------------------
468! Load "in situ" density anomaly (kg/m3 - 1000) and potential
469! density anomaly (kg/m3 - 1000) referenced to the surface into global
470! arrays. Notice that this is done in a separate (i,k) DO-loops to
471! facilitate the adjoint.
472!-----------------------------------------------------------------------
473!
474 DO k=1,n(ng)
475 DO i=istrt,iendt
476 rho(i,j,k)=den(i,k)
477 pden(i,j,k)=(den1(i,k)-1000.0_r8)
478# ifdef MASKING
479 pden(i,j,k)=pden(i,j,k)*rmask(i,j)
480# endif
481 END DO
482 END DO
483 END DO
484!
485!-----------------------------------------------------------------------
486! Exchange boundary data.
487!-----------------------------------------------------------------------
488!
489 IF (ewperiodic(ng).or.nsperiodic(ng)) THEN
490 CALL exchange_r3d_tile (ng, tile, &
491 & lbi, ubi, lbj, ubj, 1, n(ng), &
492 & rho)
493 CALL exchange_r3d_tile (ng, tile, &
494 & lbi, ubi, lbj, ubj, 1, n(ng), &
495 & pden)
496
497# if defined LMD_SKPP || defined LMD_BKPP || \
498 defined bulk_fluxes || defined balance_operator
499# ifdef LMD_DDMIX
500 CALL exchange_w3d_tile (ng, tile, &
501 & lbi, ubi, lbj, ubj, 0, n(ng), &
502 & alfaobeta)
503# endif
504 CALL exchange_r2d_tile (ng, tile, &
505 & lbi, ubi, lbj, ubj, &
506 & alpha)
507 CALL exchange_r2d_tile (ng, tile, &
508 & lbi, ubi, lbj, ubj, &
509 & beta)
510# endif
511
512# ifdef VAR_RHO_2D
513 CALL exchange_r2d_tile (ng, tile, &
514 & lbi, ubi, lbj, ubj, &
515 & rhoa)
516 CALL exchange_r2d_tile (ng, tile, &
517 & lbi, ubi, lbj, ubj, &
518 & rhos)
519# endif
520
521# ifdef BV_FREQUENCY
522 CALL exchange_w3d_tile (ng, tile, &
523 & lbi, ubi, lbj, ubj, 0, n(ng), &
524 & bvf)
525# endif
526 END IF
527
528# ifdef DISTRIBUTE
529!
530 CALL mp_exchange3d (ng, tile, model, 2, &
531 & lbi, ubi, lbj, ubj, 1, n(ng), &
532 & nghostpoints, &
533 & ewperiodic(ng), nsperiodic(ng), &
534 & rho, pden)
535
536# if defined LMD_SKPP || defined LMD_BKPP || \
537 defined bulk_fluxes || defined balance_operator
538# ifdef LMD_DDMIX
539 CALL mp_exchange3d (ng, tile, model, 1, &
540 & lbi, ubi, lbj, ubj, 0, n(ng), &
541 & nghostpoints, &
542 & ewperiodic(ng), nsperiodic(ng), &
543 & alfaobeta)
544# endif
545 CALL mp_exchange2d (ng, tile, model, 2, &
546 & lbi, ubi, lbj, ubj, &
547 & nghostpoints, &
548 & ewperiodic(ng), nsperiodic(ng), &
549 & alpha, beta)
550# endif
551
552# ifdef VAR_RHO_2D
553 CALL mp_exchange2d (ng, tile, model, 2, &
554 & lbi, ubi, lbj, ubj, &
555 & nghostpoints, &
556 & ewperiodic(ng), nsperiodic(ng), &
557 & rhoa, rhos)
558# endif
559
560# ifdef BV_FREQUENCY
561 CALL mp_exchange3d (ng, tile, model, 1, &
562 & lbi, ubi, lbj, ubj, 0, n(ng), &
563 & nghostpoints, &
564 & ewperiodic(ng), nsperiodic(ng), &
565 & bvf)
566# endif
567# endif
568!
569 RETURN
570 END SUBROUTINE rho_eos_tile
571# endif
572
573# ifndef NONLIN_EOS
574!
575!***********************************************************************
576 SUBROUTINE rho_eos_tile (ng, tile, model, &
577 & LBi, UBi, LBj, UBj, &
578 & IminS, ImaxS, JminS, JmaxS, &
579 & nrhs, &
580# ifdef MASKING
581 & rmask, &
582# endif
583# ifdef VAR_RHO_2D
584 & Hz, &
585# endif
586 & z_r, z_w, t, &
587# ifdef VAR_RHO_2D
588 & rhoA, rhoS, &
589# endif
590# ifdef BV_FREQUENCY
591 & bvf, &
592# endif
593# if defined LMD_SKPP || defined LMD_BKPP || \
594 defined BULK_FLUXES || defined BALANCE_OPERATOR
595 & alpha, beta, &
596# ifdef LMD_DDMIX
597 & alfaobeta, &
598# endif
599# endif
600 & pden, &
601 & rho)
602!***********************************************************************
603!
604 USE mod_param
605 USE mod_scalars
606# ifdef SEDIMENT
607 USE mod_sediment
608# endif
609!
610 USE exchange_2d_mod
611 USE exchange_3d_mod
612# ifdef DISTRIBUTE
613 USE mp_exchange_mod, ONLY : mp_exchange2d, mp_exchange3d
614# endif
615!
616! Imported variable declarations.
617!
618 integer, intent(in) :: ng, tile, model
619 integer, intent(in) :: LBi, UBi, LBj, UBj
620 integer, intent(in) :: IminS, ImaxS, JminS, JmaxS
621 integer, intent(in) :: nrhs
622!
623# ifdef ASSUMED_SHAPE
624# ifdef MASKING
625 real(r8), intent(in) :: rmask(LBi:,LBj:)
626# endif
627# ifdef VAR_RHO_2D
628 real(r8), intent(in) :: Hz(LBi:,LBj:,:)
629# endif
630 real(r8), intent(in) :: z_r(LBi:,LBj:,:)
631 real(r8), intent(in) :: z_w(LBi:,LBj:,0:)
632 real(r8), intent(in) :: t(LBi:,LBj:,:,:,:)
633
634# ifdef VAR_RHO_2D
635 real(r8), intent(out) :: rhoA(LBi:,LBj:)
636 real(r8), intent(out) :: rhoS(LBi:,LBj:)
637# endif
638# ifdef BV_FREQUENCY
639 real(r8), intent(out) :: bvf(LBi:,LBj:,0:)
640# endif
641# if defined LMD_SKPP || defined LMD_BKPP || \
642 defined bulk_fluxes || defined balance_operator
643 real(r8), intent(out) :: alpha(LBi:,LBj:)
644 real(r8), intent(out) :: beta(LBi:,LBj:)
645# ifdef LMD_DDMIX
646 real(r8), intent(out) :: alfaobeta(LBi:,LBj:,0:)
647# endif
648# endif
649 real(r8), intent(out) :: pden(LBi:,LBj:,:)
650 real(r8), intent(out) :: rho(LBi:,LBj:,:)
651# else
652# ifdef MASKING
653 real(r8), intent(in) :: rmask(LBi:UBi,LBj:UBj)
654# endif
655# ifdef VAR_RHO_2D
656 real(r8), intent(in) :: Hz(LBi:UBi,LBj:UBj,N(ng))
657# endif
658 real(r8), intent(in) :: z_r(LBi:UBi,LBj:UBj,N(ng))
659 real(r8), intent(in) :: z_w(LBi:UBi,LBj:UBj,0:N(ng))
660 real(r8), intent(in) :: t(LBi:UBi,LBj:UBj,N(ng),3,NT(ng))
661# ifdef VAR_RHO_2D
662 real(r8), intent(out) :: rhoA(LBi:UBi,LBj:UBj)
663 real(r8), intent(out) :: rhoS(LBi:UBi,LBj:UBj)
664# endif
665# ifdef BV_FREQUENCY
666 real(r8), intent(out) :: bvf(LBi:UBi,LBj:UBj,0:N(ng))
667# endif
668# if defined LMD_SKPP || defined LMD_BKPP || \
669 defined bulk_fluxes || defined balance_operator
670 real(r8), intent(out) :: alpha(LBi:UBi,LBj:UBj)
671 real(r8), intent(out) :: beta(LBi:UBi,LBj:UBj)
672# ifdef LMD_DDMIX
673 real(r8), intent(out) :: alfaobeta(LBi:UBi,LBj:UBj,0:N(ng))
674# endif
675# endif
676 real(r8), intent(out) :: pden(LBi:UBi,LBj:UBj,N(ng))
677 real(r8), intent(out) :: rho(LBi:UBi,LBj:UBj,N(ng))
678# endif
679!
680! Local variable declarations.
681!
682 integer :: i, ised, itrc, j, k
683 real(r8) :: SedDen, cff, cff1, cff2
684
685# include "set_bounds.h"
686!
687!=======================================================================
688! Linear equation of state.
689!=======================================================================
690!
691!-----------------------------------------------------------------------
692! Compute "in situ" density anomaly (kg/m3 - 1000) using the linear
693! equation of state.
694!-----------------------------------------------------------------------
695!
696 DO j=jstrt,jendt
697 DO k=1,n(ng)
698 DO i=istrt,iendt
699 rho(i,j,k)=r0(ng)- &
700 & r0(ng)*tcoef(ng)*(t(i,j,k,nrhs,itemp)-t0(ng))
701# ifdef SALINITY
702 rho(i,j,k)=rho(i,j,k)+ &
703 & r0(ng)*scoef(ng)*(t(i,j,k,nrhs,isalt)-s0(ng))
704# endif
705# if defined SEDIMENT && defined SED_DENS
706 sedden=0.0_r8
707 DO ised=1,nst
708 itrc=idsed(ised)
709 cff1=1.0_r8/srho(ised,ng)
710 sedden=sedden+ &
711 & t(i,j,k,nrhs,itrc)* &
712 & (srho(ised,ng)-rho(i,j,k))*cff1
713 END DO
714 rho(i,j,k)=rho(i,j,k)+sedden
715# endif
716 rho(i,j,k)=rho(i,j,k)-1000.0_r8
717# ifdef MASKING
718 rho(i,j,k)=rho(i,j,k)*rmask(i,j)
719# endif
720 pden(i,j,k)=rho(i,j,k)
721 END DO
722 END DO
723
724# ifdef VAR_RHO_2D
725!
726!-----------------------------------------------------------------------
727! Compute vertical averaged density (rhoA) and density perturbation
728! used (rhoS) in barotropic pressure gradient.
729!-----------------------------------------------------------------------
730!
731 DO i=istrt,iendt
732 cff1=rho(i,j,n(ng))*hz(i,j,n(ng))
733 rhos(i,j)=0.5_r8*cff1*hz(i,j,n(ng))
734 rhoa(i,j)=cff1
735 END DO
736 DO k=n(ng)-1,1,-1
737 DO i=istrt,iendt
738 cff1=rho(i,j,k)*hz(i,j,k)
739 rhos(i,j)=rhos(i,j)+hz(i,j,k)*(rhoa(i,j)+0.5_r8*cff1)
740 rhoa(i,j)=rhoa(i,j)+cff1
741 END DO
742 END DO
743 cff2=1.0_r8/rho0
744 DO i=istrt,iendt
745 cff1=1.0_r8/(z_w(i,j,n(ng))-z_w(i,j,0))
746 rhoa(i,j)=cff2*cff1*rhoa(i,j)
747 rhos(i,j)=2.0_r8*cff1*cff1*cff2*rhos(i,j)
748 END DO
749# endif
750
751# ifdef BV_FREQUENCY
752!
753!-----------------------------------------------------------------------
754! Compute Brunt-Vaisala frequency (1/s2) at horizontal RHO-points
755! and vertical W-points.
756!-----------------------------------------------------------------------
757!
758 DO k=1,n(ng)-1
759 DO i=istrt,iendt
760 bvf(i,j,k)=-gorho0*(rho(i,j,k+1)-rho(i,j,k))/ &
761 & (z_r(i,j,k+1)-z_r(i,j,k))
762 END DO
763 END DO
764# endif
765
766# if defined LMD_SKPP || defined LMD_BKPP || \
767 defined bulk_fluxes || defined balance_operator
768!
769!-----------------------------------------------------------------------
770! Compute thermal expansion (1/Celsius) and saline contraction
771! (1/PSU) coefficients.
772!-----------------------------------------------------------------------
773!
774 DO i=istrt,iendt
775 alpha(i,j)=abs(tcoef(ng))
776# ifdef SALINITY
777 beta(i,j)=abs(scoef(ng))
778# else
779 beta(i,j)=0.0_r8
780# endif
781 END DO
782# ifdef LMD_DDMIX
783!
784! Compute ratio of thermal expansion and saline contraction
785! coefficients.
786!
787 IF (scoef(ng).eq.0.0_r8) THEN
788 cff=1.0_r8
789 ELSE
790 cff=1.0_r8/scoef(ng)
791 END IF
792 DO k=1,n(ng)
793 DO i=istrt,iendt
794 alfaobeta(i,j,k)=cff*tcoef(ng)
795 END DO
796 END DO
797# endif
798# endif
799 END DO
800!
801!-----------------------------------------------------------------------
802! Exchange boundary data.
803!-----------------------------------------------------------------------
804!
805 IF (ewperiodic(ng).or.nsperiodic(ng)) THEN
806 CALL exchange_r3d_tile (ng, tile, &
807 & lbi, ubi, lbj, ubj, 1, n(ng), &
808 & rho)
809 CALL exchange_r3d_tile (ng, tile, &
810 & lbi, ubi, lbj, ubj, 1, n(ng), &
811 & pden)
812
813# if defined LMD_SKPP || defined LMD_BKPP || \
814 defined bulk_fluxes || defined balance_operator
815# ifdef LMD_DDMIX
816 CALL exchange_w3d_tile (ng, tile, &
817 & lbi, ubi, lbj, ubj, 0, n(ng), &
818 & alfaobeta)
819# endif
820 CALL exchange_r2d_tile (ng, tile, &
821 & lbi, ubi, lbj, ubj, &
822 & alpha)
823 CALL exchange_r2d_tile (ng, tile, &
824 & lbi, ubi, lbj, ubj, &
825 & beta)
826# endif
827
828# ifdef VAR_RHO_2D
829 CALL exchange_r2d_tile (ng, tile, &
830 & lbi, ubi, lbj, ubj, &
831 & rhoa)
832 CALL exchange_r2d_tile (ng, tile, &
833 & lbi, ubi, lbj, ubj, &
834 & rhos)
835# endif
836
837# ifdef BV_FREQUENCY
838 CALL exchange_w3d_tile (ng, tile, &
839 & lbi, ubi, lbj, ubj, 0, n(ng), &
840 & bvf)
841# endif
842 END IF
843
844# ifdef DISTRIBUTE
845!
846 CALL mp_exchange3d (ng, tile, model, 2, &
847 & lbi, ubi, lbj, ubj, 1, n(ng), &
848 & nghostpoints, &
849 & ewperiodic(ng), nsperiodic(ng), &
850 & rho, pden)
851
852# if defined LMD_SKPP || defined LMD_BKPP || \
853 defined bulk_fluxes || defined balance_operator
854# ifdef LMD_DDMIX
855 CALL mp_exchange3d (ng, tile, model, 1, &
856 & lbi, ubi, lbj, ubj, 0, n(ng), &
857 & nghostpoints, &
858 & ewperiodic(ng), nsperiodic(ng), &
859 & alfaobeta)
860# endif
861 CALL mp_exchange2d (ng, tile, model, 2, &
862 & lbi, ubi, lbj, ubj, &
863 & nghostpoints, &
864 & ewperiodic(ng), nsperiodic(ng), &
865 & alpha, beta)
866# endif
867
868# ifdef VAR_RHO_2D
869 CALL mp_exchange2d (ng, tile, model, 2, &
870 & lbi, ubi, lbj, ubj, &
871 & nghostpoints, &
872 & ewperiodic(ng), nsperiodic(ng), &
873 & rhoa, rhos)
874# endif
875
876# ifdef BV_FREQUENCY
877 CALL mp_exchange3d (ng, tile, model, 1, &
878 & lbi, ubi, lbj, ubj, 0, n(ng), &
879 & nghostpoints, &
880 & ewperiodic(ng), nsperiodic(ng), &
881 & bvf)
882# endif
883# endif
884!
885 RETURN
886 END SUBROUTINE rho_eos_tile
887# endif
888#endif
889 END MODULE rho_eos_mod
exchange_2d_mod
Definition exchange_2d.F:2
exchange_2d_mod::exchange_r2d_tile
subroutine exchange_r2d_tile(ng, tile, lbi, ubi, lbj, ubj, a)
Definition exchange_2d.F:253
exchange_3d_mod
Definition exchange_3d.F:2
exchange_3d_mod::exchange_r3d_tile
subroutine exchange_r3d_tile(ng, tile, lbi, ubi, lbj, ubj, lbk, ubk, a)
Definition exchange_3d.F:283
exchange_3d_mod::exchange_w3d_tile
subroutine exchange_w3d_tile(ng, tile, lbi, ubi, lbj, ubj, lbk, ubk, a)
Definition exchange_3d.F:920
mod_coupling
Definition mod_coupling.F:2
mod_coupling::coupling
type(t_coupling), dimension(:), allocatable coupling
Definition mod_coupling.F:102
mod_eoscoef
Definition mod_eoscoef.F:2
mod_eoscoef::b03
real(r8), parameter b03
Definition mod_eoscoef.F:32
mod_eoscoef::e02
real(r8), parameter e02
Definition mod_eoscoef.F:38
mod_eoscoef::h02
real(r8), parameter h02
Definition mod_eoscoef.F:49
mod_eoscoef::q03
real(r8), parameter q03
Definition mod_eoscoef.F:53
mod_eoscoef::a03
real(r8), parameter a03
Definition mod_eoscoef.F:27
mod_eoscoef::q04
real(r8), parameter q04
Definition mod_eoscoef.F:54
mod_eoscoef::q05
real(r8), parameter q05
Definition mod_eoscoef.F:55
mod_eoscoef::u01
real(r8), parameter u01
Definition mod_eoscoef.F:57
mod_eoscoef::f00
real(r8), parameter f00
Definition mod_eoscoef.F:40
mod_eoscoef::v01
real(r8), parameter v01
Definition mod_eoscoef.F:62
mod_eoscoef::g01
real(r8), parameter g01
Definition mod_eoscoef.F:44
mod_eoscoef::a04
real(r8), parameter a04
Definition mod_eoscoef.F:28
mod_eoscoef::f01
real(r8), parameter f01
Definition mod_eoscoef.F:41
mod_eoscoef::g00
real(r8), parameter g00
Definition mod_eoscoef.F:43
mod_eoscoef::u02
real(r8), parameter u02
Definition mod_eoscoef.F:58
mod_eoscoef::d01
real(r8), parameter d01
Definition mod_eoscoef.F:34
mod_eoscoef::w00
real(r8), parameter w00
Definition mod_eoscoef.F:64
mod_eoscoef::b00
real(r8), parameter b00
Definition mod_eoscoef.F:29
mod_eoscoef::u03
real(r8), parameter u03
Definition mod_eoscoef.F:59
mod_eoscoef::d00
real(r8), parameter d00
Definition mod_eoscoef.F:33
mod_eoscoef::d02
real(r8), parameter d02
Definition mod_eoscoef.F:35
mod_eoscoef::g03
real(r8), parameter g03
Definition mod_eoscoef.F:46
mod_eoscoef::b01
real(r8), parameter b01
Definition mod_eoscoef.F:30
mod_eoscoef::u04
real(r8), parameter u04
Definition mod_eoscoef.F:60
mod_eoscoef::u00
real(r8), parameter u00
Definition mod_eoscoef.F:56
mod_eoscoef::v02
real(r8), parameter v02
Definition mod_eoscoef.F:63
mod_eoscoef::q01
real(r8), parameter q01
Definition mod_eoscoef.F:51
mod_eoscoef::e00
real(r8), parameter e00
Definition mod_eoscoef.F:36
mod_eoscoef::h00
real(r8), parameter h00
Definition mod_eoscoef.F:47
mod_eoscoef::g02
real(r8), parameter g02
Definition mod_eoscoef.F:45
mod_eoscoef::e03
real(r8), parameter e03
Definition mod_eoscoef.F:39
mod_eoscoef::a00
real(r8), parameter a00
Definition mod_eoscoef.F:24
mod_eoscoef::a01
real(r8), parameter a01
Definition mod_eoscoef.F:25
mod_eoscoef::a02
real(r8), parameter a02
Definition mod_eoscoef.F:26
mod_eoscoef::b02
real(r8), parameter b02
Definition mod_eoscoef.F:31
mod_eoscoef::h01
real(r8), parameter h01
Definition mod_eoscoef.F:48
mod_eoscoef::e01
real(r8), parameter e01
Definition mod_eoscoef.F:37
mod_eoscoef::f02
real(r8), parameter f02
Definition mod_eoscoef.F:42
mod_eoscoef::v00
real(r8), parameter v00
Definition mod_eoscoef.F:61
mod_eoscoef::q02
real(r8), parameter q02
Definition mod_eoscoef.F:52
mod_eoscoef::q00
real(r8), parameter q00
Definition mod_eoscoef.F:50
mod_grid
Definition mod_grid.F:2
mod_grid::grid
type(t_grid), dimension(:), allocatable grid
Definition mod_grid.F:365
mod_mixing
Definition mod_mixing.F:2
mod_mixing::mixing
type(t_mixing), dimension(:), allocatable mixing
Definition mod_mixing.F:399
mod_ocean
Definition mod_ocean.F:2
mod_ocean::ocean
type(t_ocean), dimension(:), allocatable ocean
Definition mod_ocean.F:351
mod_param
Definition mod_param.F:2
mod_param::nghostpoints
integer nghostpoints
Definition mod_param.F:710
mod_param::nst
integer nst
Definition mod_param.F:521
mod_scalars
Definition mod_scalars.F:2
mod_scalars::t0
real(r8), dimension(:), allocatable t0
Definition mod_scalars.F:738
mod_scalars::ewperiodic
logical, dimension(:), allocatable ewperiodic
Definition mod_scalars.F:510
mod_scalars::nsperiodic
logical, dimension(:), allocatable nsperiodic
Definition mod_scalars.F:511
mod_scalars::s0
real(r8), dimension(:), allocatable s0
Definition mod_scalars.F:739
mod_scalars::tcoef
real(r8), dimension(:), allocatable tcoef
Definition mod_scalars.F:732
mod_scalars::gorho0
real(dp) gorho0
Definition mod_scalars.F:468
mod_scalars::isalt
integer isalt
Definition mod_scalars.F:127
mod_scalars::itemp
integer itemp
Definition mod_scalars.F:126
mod_scalars::g
real(dp) g
Definition mod_scalars.F:464
mod_scalars::rho0
real(dp) rho0
Definition mod_scalars.F:787
mod_scalars::r0
real(r8), dimension(:), allocatable r0
Definition mod_scalars.F:731
mod_scalars::scoef
real(r8), dimension(:), allocatable scoef
Definition mod_scalars.F:733
mod_sediment
Definition sediment_mod.h:1
mod_sediment::srho
real(r8), dimension(:,:), allocatable srho
Definition sediment_mod.h:174
mod_sediment::idsed
integer, dimension(:), allocatable idsed
Definition sediment_mod.h:88
mod_stepping
Definition mod_stepping.F:2
mod_stepping::nrhs
integer, dimension(:), allocatable nrhs
Definition mod_stepping.F:70
mp_exchange_mod
Definition mp_exchange.F:2
mp_exchange_mod::mp_exchange2d
subroutine mp_exchange2d(ng, tile, model, nvar, lbi, ubi, lbj, ubj, nghost, ew_periodic, ns_periodic, a, b, c, d)
Definition mp_exchange.F:294
mp_exchange_mod::mp_exchange3d
subroutine mp_exchange3d(ng, tile, model, nvar, lbi, ubi, lbj, ubj, lbk, ubk, nghost, ew_periodic, ns_periodic, a, b, c, d)
Definition mp_exchange.F:2037
rho_eos_mod
Definition rho_eos.F:2
rho_eos_mod::rho_eos
subroutine, public rho_eos(ng, tile, model)
Definition rho_eos.F:48
rho_eos_mod::rho_eos_tile
subroutine rho_eos_tile(ng, tile, model, lbi, ubi, lbj, ubj, imins, imaxs, jmins, jmaxs, nrhs, rmask, hz, z_r, z_w, t, rhoa, rhos, alpha, beta, alfaobeta, pden, rho)
Definition rho_eos.F:137
wclock_off
recursive subroutine wclock_off(ng, model, region, line, routine)
Definition timers.F:148
wclock_on
recursive subroutine wclock_on(ng, model, region, line, routine)
Definition timers.F:3