ROMS
Loading...
Searching...
No Matches
tl_rho_eos.F
Go to the documentation of this file.
1#include "cppdefs.h"
3#if defined TANGENT && defined SOLVE3D
4!
5!git $Id$
6!================================================== Hernan G. Arango ===
7! Copyright (c) 2002-2025 The ROMS Group Andrew M. Moore !
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 :: tl_rho_eos
43!
44 CONTAINS
45!
46!***********************************************************************
47 SUBROUTINE tl_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 tl_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_NOT_YET
79 & grid(ng) % Hz, &
80 & grid(ng) % tl_Hz, &
81# endif
82 & grid(ng) % z_r, &
83 & grid(ng) % tl_z_r, &
84 & grid(ng) % z_w, &
85 & grid(ng) % tl_z_w, &
86 & ocean(ng) % t, &
87 & ocean(ng) % tl_t, &
88# ifdef VAR_RHO_2D_NOT_YET
89 & coupling(ng) % rhoA, &
90 & coupling(ng) % tl_rhoA, &
91 & coupling(ng) % rhoS, &
92 & coupling(ng) % tl_rhoS, &
93# endif
94# ifdef BV_FREQUENCY_NOT_YET
95 & mixing(ng) % tl_bvf, &
96# endif
97# if defined LMD_SKPP_NOT_YET || defined LMD_BKPP_NOT_YET || \
98 defined bulk_fluxes
99 & mixing(ng) % alpha, &
100 & mixing(ng) % tl_alpha, &
101 & mixing(ng) % beta, &
102 & mixing(ng) % tl_beta, &
103# ifdef LMD_DDMIX_NOT_YET
104 & mixing(ng) % tl_alfaobeta, &
105# endif
106# endif
107# if defined LMD_SKPP_NOT_YET || defined LMD_BKPP_NOT_YET
108 & ocean(ng) % tl_pden, &
109# endif
110 & ocean(ng) % rho, &
111 & ocean(ng) % tl_rho)
112# ifdef PROFILE
113 CALL wclock_off (ng, model, 14, __line__, myfile)
114# endif
115!
116 RETURN
117 END SUBROUTINE tl_rho_eos
118
119# ifdef NONLIN_EOS
120!
121!***********************************************************************
122 SUBROUTINE tl_rho_eos_tile (ng, tile, model, &
123 & LBi, UBi, LBj, UBj, &
124 & IminS, ImaxS, JminS, JmaxS, &
125 & nrhs, &
126# ifdef MASKING
127 & rmask, &
128# endif
129# ifdef VAR_RHO_2D_NOT_YET
130 & Hz, tl_Hz, &
131# endif
132 & z_r, tl_z_r, &
133 & z_w, tl_z_w, &
134 & t, tl_t, &
135# ifdef VAR_RHO_2D_NOT_YET
136 & rhoA, tl_rhoA, &
137 & rhoS, tl_rhoS, &
138# endif
139# ifdef BV_FREQUENCY_NOT_YET
140 & tl_bvf, &
141# endif
142# if defined LMD_SKPP_NOT_YET || defined LMD_BKPP_NOT_YET || \
143 defined BULK_FLUXES
144 & alpha, tl_alpha, &
145 & beta, tl_beta, &
146# ifdef LMD_DDMIX_NOT_YET
147 & tl_alfaobeta, &
148# endif
149# endif
150# if defined LMD_SKPP_NOT_YET || defined LMD_BKPP_NOT_YET
151 & tl_pden, &
152# endif
153 & rho, tl_rho)
154!***********************************************************************
155!
156 USE mod_param
157 USE mod_eoscoef
158 USE mod_scalars
159# ifdef SEDIMENT_NOT_YET
160 USE mod_sediment
161# endif
162!
165# ifdef DISTRIBUTE
167# endif
168!
169! Imported variable declarations.
170!
171 integer, intent(in) :: ng, tile, model
172 integer, intent(in) :: LBi, UBi, LBj, UBj
173 integer, intent(in) :: IminS, ImaxS, JminS, JmaxS
174 integer, intent(in) :: nrhs
175!
176# ifdef ASSUMED_SHAPE
177# ifdef MASKING
178 real(r8), intent(in) :: rmask(LBi:,LBj:)
179# endif
180# ifdef VAR_RHO_2D_NOT_YET
181 real(r8), intent(in) :: Hz(LBi:,LBj:,:)
182# endif
183 real(r8), intent(in) :: z_r(LBi:,LBj:,:)
184 real(r8), intent(in) :: z_w(LBi:,LBj:,0:)
185 real(r8), intent(in) :: t(LBi:,LBj:,:,:,:)
186# ifdef VAR_RHO_2D_NOT_YET
187 real(r8), intent(in) :: tl_Hz(LBi:,LBj:,:)
188# endif
189 real(r8), intent(in) :: tl_z_r(LBi:,LBj:,:)
190 real(r8), intent(in) :: tl_z_w(LBi:,LBj:,0:)
191 real(r8), intent(in) :: tl_t(LBi:,LBj:,:,:,:)
192# if defined LMD_SKPP_NOT_YET || defined LMD_BKPP_NOT_YET || \
193 defined bulk_fluxes
194 real(r8), intent(inout) :: alpha(LBi:,LBj:)
195 real(r8), intent(inout) :: beta(LBi:,LBj:)
196# endif
197# ifdef VAR_RHO_2D_NOT_YET
198 real(r8), intent(out) :: rhoA(LBi:,LBj:)
199 real(r8), intent(out) :: rhoS(LBi:,LBj:)
200 real(r8), intent(out) :: tl_rhoA(LBi:,LBj:)
201 real(r8), intent(out) :: tl_rhoS(LBi:,LBj:)
202# endif
203# ifdef BV_FREQUENCY_NOT_YET
204 real(r8), intent(out) :: tl_bvf(LBi:,LBj:,0:)
205# endif
206# if defined LMD_SKPP_NOT_YET || defined LMD_BKPP_NOT_YET || \
207 defined bulk_fluxes
208 real(r8), intent(out) :: tl_alpha(LBi:,LBj:)
209 real(r8), intent(out) :: tl_beta(LBi:,LBj:)
210# ifdef LMD_DDMIX_NOT_YET
211 real(r8), intent(out) :: tl_alfaobeta(LBi:,LBj:,0:)
212# endif
213# endif
214# if defined LMD_SKPP_NOT_YET || defined LMD_BKPP_NOT_YET
215 real(r8), intent(out) :: tl_pden(LBi:,LBj:,:)
216# endif
217 real(r8), intent(out) :: rho(LBi:,LBj:,:)
218 real(r8), intent(out) :: tl_rho(LBi:,LBj:,:)
219# else
220# ifdef MASKING
221 real(r8), intent(in) :: rmask(LBi:UBi,LBj:UBj)
222# endif
223# ifdef VAR_RHO_2D_NOT_YET
224 real(r8), intent(in) :: Hz(LBi:UBi,LBj:UBj,N(ng))
225# endif
226 real(r8), intent(in) :: z_r(LBi:UBi,LBj:UBj,N(ng))
227 real(r8), intent(in) :: z_w(LBi:UBi,LBj:UBj,0:N(ng))
228 real(r8), intent(in) :: t(LBi:UBi,LBj:UBj,N(ng),3,NT(ng))
229# ifdef VAR_RHO_2D_NOT_YET
230 real(r8), intent(in) :: tl_Hz(LBi:UBi,LBj:UBj,N(ng))
231# endif
232 real(r8), intent(in) :: tl_z_r(LBi:UBi,LBj:UBj,N(ng))
233 real(r8), intent(in) :: tl_z_w(LBi:UBi,LBj:UBj,0:N(ng))
234 real(r8), intent(in) :: tl_t(LBi:UBi,LBj:UBj,N(ng),3,NT(ng))
235# if defined LMD_SKPP_NOT_YET || defined LMD_BKPP_NOT_YET || \
236 defined bulk_fluxes
237 real(r8), intent(inout) :: alpha(LBi:UBi,LBj:UBj)
238 real(r8), intent(inout) :: beta(LBi:UBi,LBj:UBj)
239# endif
240# ifdef VAR_RHO_2D_NOT_YET
241 real(r8), intent(out) :: rhoA(LBi:UBi,LBj:UBj)
242 real(r8), intent(out) :: rhoS(LBi:UBi,LBj:UBj)
243 real(r8), intent(out) :: tl_rhoA(LBi:UBi,LBj:UBj)
244 real(r8), intent(out) :: tl_rhoS(LBi:UBi,LBj:UBj)
245# endif
246# ifdef BV_FREQUENCY_NOT_YET
247 real(r8), intent(out) :: tl_bvf(LBi:UBi,LBj:UBj,0:N(ng))
248# endif
249# if defined LMD_SKPP_NOT_YET || defined LMD_BKPP_NOT_YET || \
250 defined bulk_fluxes
251 real(r8), intent(out) :: tl_alpha(LBi:UBi,LBj:UBj)
252 real(r8), intent(out) :: tl_beta(LBi:UBi,LBj:UBj)
253# ifdef LMD_DDMIX_NOT_YET
254 real(r8), intent(out) :: tl_alfaobeta(LBi:UBi,LBj:UBj,0:N(ng))
255# endif
256# endif
257# if defined LMD_SKPP_NOT_YET || defined LMD_BKPP_NOT_YET
258 real(r8), intent(out) :: tl_pden(LBi:UBi,LBj:UBj,N(ng))
259# endif
260 real(r8), intent(out) :: rho(LBi:UBi,LBj:UBj,N(ng))
261 real(r8), intent(out) :: tl_rho(LBi:UBi,LBj:UBj,N(ng))
262# endif
263!
264! Local variable declarations.
265!
266 integer :: i, ised, itrc, j, k
267
268 real(r8) :: SedDen, Tp, Tpr10, Ts, Tt, sqrtTs
269 real(r8) :: tl_SedDen, tl_Tp, tl_Tpr10, tl_Ts, tl_Tt, tl_sqrtTs
270# ifdef BV_FREQUENCY_NOT_YET
271 real(r8) :: bulk_dn, bulk_up, den_dn, den_up
272 real(r8) :: tl_bulk_dn, tl_bulk_up, tl_den_dn, tl_den_up
273# endif
274 real(r8) :: cff, cff1, cff2, cff3
275 real(r8) :: tl_cff, tl_cff1, tl_cff2, tl_cff3
276
277 real(r8), dimension(0:9) :: C
278 real(r8), dimension(0:9) :: tl_C
279# ifdef EOS_TDERIVATIVE
280 real(r8), dimension(0:9) :: dCdT(0:9)
281 real(r8), dimension(0:9) :: tl_dCdT(0:9)
282 real(r8), dimension(0:9) :: d2Cd2T(0:9)
283
284 real(r8), dimension(IminS:ImaxS,N(ng)) :: DbulkDS
285 real(r8), dimension(IminS:ImaxS,N(ng)) :: DbulkDT
286 real(r8), dimension(IminS:ImaxS,N(ng)) :: Dden1DS
287 real(r8), dimension(IminS:ImaxS,N(ng)) :: Dden1DT
288 real(r8), dimension(IminS:ImaxS,N(ng)) :: Scof
289 real(r8), dimension(IminS:ImaxS,N(ng)) :: Tcof
290 real(r8), dimension(IminS:ImaxS,N(ng)) :: wrk
291
292 real(r8), dimension(IminS:ImaxS,N(ng)) :: tl_DbulkDS
293 real(r8), dimension(IminS:ImaxS,N(ng)) :: tl_DbulkDT
294 real(r8), dimension(IminS:ImaxS,N(ng)) :: tl_Dden1DS
295 real(r8), dimension(IminS:ImaxS,N(ng)) :: tl_Dden1DT
296 real(r8), dimension(IminS:ImaxS,N(ng)) :: tl_Scof
297 real(r8), dimension(IminS:ImaxS,N(ng)) :: tl_Tcof
298 real(r8), dimension(IminS:ImaxS,N(ng)) :: tl_wrk
299# endif
300 real(r8), dimension(IminS:ImaxS,N(ng)) :: bulk
301 real(r8), dimension(IminS:ImaxS,N(ng)) :: bulk0
302 real(r8), dimension(IminS:ImaxS,N(ng)) :: bulk1
303 real(r8), dimension(IminS:ImaxS,N(ng)) :: bulk2
304 real(r8), dimension(IminS:ImaxS,N(ng)) :: den
305 real(r8), dimension(IminS:ImaxS,N(ng)) :: den1
306
307 real(r8), dimension(IminS:ImaxS,N(ng)) :: tl_bulk
308 real(r8), dimension(IminS:ImaxS,N(ng)) :: tl_bulk0
309 real(r8), dimension(IminS:ImaxS,N(ng)) :: tl_bulk1
310 real(r8), dimension(IminS:ImaxS,N(ng)) :: tl_bulk2
311 real(r8), dimension(IminS:ImaxS,N(ng)) :: tl_den
312 real(r8), dimension(IminS:ImaxS,N(ng)) :: tl_den1
313
314# include "set_bounds.h"
315!
316!=======================================================================
317! Nonlinear equation of state. Notice that this equation of state
318! is only valid for potential temperature range of -2C to 40C and
319! a salinity range of 0 PSU to 42 PSU.
320!=======================================================================
321!
322 DO j=jstrt,jendt
323 DO k=1,n(ng)
324 DO i=istrt,iendt
325!
326! Check temperature and salinity lower values. Assign depth to the
327! pressure.
328!
329 tt=max(-2.0_r8,t(i,j,k,nrhs,itemp))
330 tl_tt=(0.5_r8-sign(0.5_r8,-2.0_r8- &
331 & t(i,j,k,nrhs,itemp)))* &
332 & tl_t(i,j,k,nrhs,itemp)
333# ifdef SALINITY
334 ts=max(0.0_r8,t(i,j,k,nrhs,isalt))
335 tl_ts=(0.5_r8-sign(0.5_r8,-t(i,j,k,nrhs,isalt)))* &
336 & tl_t(i,j,k,nrhs,isalt)
337 sqrtts=sqrt(ts)
338 IF (ts.ne.0.0_r8) THEN
339 tl_sqrtts=0.5_r8*tl_ts/sqrt(ts)
340 ELSE
341 tl_sqrtts=0.0_r8
342 END IF
343# else
344 ts=0.0_r8
345 tl_ts=0.0_r8
346 sqrtts=0.0_r8
347 tl_sqrtts=0.0_r8
348# endif
349 tp=z_r(i,j,k)
350 tl_tp=tl_z_r(i,j,k)
351 tpr10=0.1_r8*tp
352 tl_tpr10=0.1_r8*tl_tp
353!
354!-----------------------------------------------------------------------
355! Compute BASIC STATE and tangent linear density (kg/m3) at standard
356! one atmosphere pressure.
357!-----------------------------------------------------------------------
358!
359 c(0)=q00+tt*(q01+tt*(q02+tt*(q03+tt*(q04+tt*q05))))
360 c(1)=u00+tt*(u01+tt*(u02+tt*(u03+tt*u04)))
361 c(2)=v00+tt*(v01+tt*v02)
362# ifdef EOS_TDERIVATIVE
363!
364 dcdt(0)=q01+tt*(2.0_r8*q02+tt*(3.0_r8*q03+tt*(4.0_r8*q04+ &
365 & tt*5.0_r8*q05)))
366 dcdt(1)=u01+tt*(2.0_r8*u02+tt*(3.0_r8*u03+tt*4.0_r8*u04))
367 dcdt(2)=v01+tt*2.0_r8*v02
368# endif
369 tl_c(0)=tl_tt*dcdt(0)
370 tl_c(1)=tl_tt*dcdt(1)
371 tl_c(2)=tl_tt*dcdt(2)
372!
373 den1(i,k)=c(0)+ts*(c(1)+sqrtts*c(2)+ts*w00)
374 tl_den1(i,k)=tl_c(0)+ &
375 & tl_ts*(c(1)+sqrtts*c(2)+ts*w00)+ &
376 & ts*(tl_c(1)+tl_sqrtts*c(2)+ &
377 & sqrtts*tl_c(2)+tl_ts*w00)
378
379# ifdef EOS_TDERIVATIVE
380!
381! Compute d(den1)/d(S) and d(den1)/d(T) derivatives used in the
382! computation of thermal expansion and saline contraction
383! coefficients.
384!
385 d2cd2t(0)=2.0_r8*q02+tt*(6.0_r8*q03+tt*(12.0_r8*q04+ &
386 & tt*20.0_r8*q05))
387 d2cd2t(1)=2.0_r8*u02+tt*(6.0_r8*u03+tt*12.0_r8*u04)
388 d2cd2t(2)=2.0_r8*v02
389!
390 tl_dcdt(0)=tl_tt*d2cd2t(0)
391 tl_dcdt(1)=tl_tt*d2cd2t(1)
392 tl_dcdt(2)=tl_tt*d2cd2t(2)
393!
394 dden1ds(i,k)=c(1)+1.5_r8*c(2)*sqrtts+2.0_r8*w00*ts
395 dden1dt(i,k)=dcdt(0)+ts*(dcdt(1)+sqrtts*dcdt(2))
396!
397 tl_dden1ds(i,k)=tl_c(1)+ &
398 & 1.5_r8*(tl_c(2)*sqrtts+ &
399 & c(2)*tl_sqrtts)+ &
400 & 2.0_r8*w00*tl_ts
401 tl_dden1dt(i,k)=tl_dcdt(0)+ &
402 & tl_ts*(dcdt(1)+sqrtts*dcdt(2))+ &
403 & ts*(tl_dcdt(1)+tl_sqrtts*dcdt(2)+ &
404 & sqrtts*tl_dcdt(2))
405# endif
406!
407!-----------------------------------------------------------------------
408! Compute BASIC STATE and tangent linear secant bulk modulus.
409!-----------------------------------------------------------------------
410!
411 c(3)=a00+tt*(a01+tt*(a02+tt*(a03+tt*a04)))
412 c(4)=b00+tt*(b01+tt*(b02+tt*b03))
413 c(5)=d00+tt*(d01+tt*d02)
414 c(6)=e00+tt*(e01+tt*(e02+tt*e03))
415 c(7)=f00+tt*(f01+tt*f02)
416 c(8)=g01+tt*(g02+tt*g03)
417 c(9)=h00+tt*(h01+tt*h02)
418# ifdef EOS_TDERIVATIVE
419!
420 dcdt(3)=a01+tt*(2.0_r8*a02+tt*(3.0_r8*a03+tt*4.0_r8*a04))
421 dcdt(4)=b01+tt*(2.0_r8*b02+tt*3.0_r8*b03)
422 dcdt(5)=d01+tt*2.0_r8*d02
423 dcdt(6)=e01+tt*(2.0_r8*e02+tt*3.0_r8*e03)
424 dcdt(7)=f01+tt*2.0_r8*f02
425 dcdt(8)=g02+tt*2.0_r8*g03
426 dcdt(9)=h01+tt*2.0_r8*h02
427# endif
428!
429 tl_c(3)=tl_tt*dcdt(3)
430 tl_c(4)=tl_tt*dcdt(4)
431 tl_c(5)=tl_tt*dcdt(5)
432 tl_c(6)=tl_tt*dcdt(6)
433 tl_c(7)=tl_tt*dcdt(7)
434 tl_c(8)=tl_tt*dcdt(8)
435 tl_c(9)=tl_tt*dcdt(9)
436!
437 bulk0(i,k)=c(3)+ts*(c(4)+sqrtts*c(5))
438 bulk1(i,k)=c(6)+ts*(c(7)+sqrtts*g00)
439 bulk2(i,k)=c(8)+ts*c(9)
440 bulk(i,k)=bulk0(i,k)-tp*(bulk1(i,k)-tp*bulk2(i,k))
441!
442 tl_bulk0(i,k)=tl_c(3)+ &
443 & tl_ts*(c(4)+sqrtts*c(5))+ &
444 & ts*(tl_c(4)+tl_sqrtts*c(5)+ &
445 & sqrtts*tl_c(5))
446 tl_bulk1(i,k)=tl_c(6)+ &
447 & tl_ts*(c(7)+sqrtts*g00)+ &
448 & ts*(tl_c(7)+tl_sqrtts*g00)
449 tl_bulk2(i,k)=tl_c(8)+tl_ts*c(9)+ts*tl_c(9)
450 tl_bulk(i,k)=tl_bulk0(i,k)- &
451 & tl_tp*(bulk1(i,k)-tp*bulk2(i,k))- &
452 & tp*(tl_bulk1(i,k)- &
453 & tl_tp*bulk2(i,k)- &
454 & tp*tl_bulk2(i,k))
455
456# if defined LMD_SKPP_NOT_YET || defined LMD_BKPP_NOT_YET || \
457 defined bulk_fluxes
458!
459! Compute d(bulk)/d(S) and d(bulk)/d(T) derivatives used
460! in the computation of thermal expansion and saline contraction
461! coefficients.
462!
463 d2cd2t(3)=2.0_r8*a02+tt*(6.0_r8*a03+tt*12.0_r8*a04)
464 d2cd2t(4)=2.0_r8*b02+tt*6.0_r8*b03
465 d2cd2t(5)=2.0_r8*d02
466 d2cd2t(6)=2.0_r8*e02+tt*6.0_r8*e03
467 d2cd2t(7)=2.0_r8*f02
468 d2cd2t(8)=2.0_r8*g03
469 d2cd2t(9)=2.0_r8*h02
470!
471 tl_dcdt(3)=tl_tt*d2cd2t(3)
472 tl_dcdt(4)=tl_tt*d2cd2t(4)
473 tl_dcdt(5)=tl_tt*d2cd2t(5)
474 tl_dcdt(6)=tl_tt*d2cd2t(6)
475 tl_dcdt(7)=tl_tt*d2cd2t(7)
476 tl_dcdt(8)=tl_tt*d2cd2t(8)
477 tl_dcdt(9)=tl_tt*d2cd2t(9)
478!
479 dbulkds(i,k)=c(4)+sqrtts*1.5_r8*c(5)- &
480 & tp*(c(7)+sqrtts*1.5_r8*g00-tp*c(9))
481 dbulkdt(i,k)=dcdt(3)+ts*(dcdt(4)+sqrtts*dcdt(5))- &
482 & tp*(dcdt(6)+ts*dcdt(7)- &
483 & tp*(dcdt(8)+ts*dcdt(9)))
484!
485 tl_dbulkds(i,k)=tl_c(4)+ &
486 & 1.5_r8*(tl_sqrtts*c(5)+ &
487 & sqrtts*tl_c(5))- &
488 & tl_tp*(c(7)+sqrtts*1.5_r8*g00- &
489 & tp*c(9))- &
490 & tp*(tl_c(7)+tl_sqrtts*1.5_r8*g00- &
491 & tl_tp*c(9)-tp*tl_c(9))
492 tl_dbulkdt(i,k)=tl_dcdt(3)+ &
493 & tl_ts*(dcdt(4)+sqrtts*dcdt(5))+ &
494 & ts*(tl_dcdt(4)+tl_sqrtts*dcdt(5)+ &
495 & sqrtts*tl_dcdt(5))- &
496 & tl_tp*(dcdt(6)+ts*dcdt(7)- &
497 & tp*(dcdt(8)+ts*dcdt(9)))- &
498 & tp*(tl_dcdt(6)+tl_ts*dcdt(7)+ts*tl_dcdt(7)- &
499 & tl_tp*(dcdt(8)+ts*dcdt(9))- &
500 & tp*(tl_dcdt(8)+tl_ts*dcdt(9)+ &
501 & ts*tl_dcdt(9)))
502# endif
503!
504!-----------------------------------------------------------------------
505! Compute local "in situ" density anomaly (kg/m3 - 1000).
506!-----------------------------------------------------------------------
507!
508 cff=1.0_r8/(bulk(i,k)+tpr10)
509 tl_cff=-cff*cff*(tl_bulk(i,k)+tl_tpr10)
510 den(i,k)=den1(i,k)*bulk(i,k)*cff
511 tl_den(i,k)=tl_den1(i,k)*bulk(i,k)*cff+ &
512 & den1(i,k)*(tl_bulk(i,k)*cff+ &
513 & bulk(i,k)*tl_cff)
514# if defined SEDIMENT_NOT_YET && defined SED_DENS_NOT_YET
515 sedden=0.0_r8
516 tl_sedden=0.0_r8
517 DO ised=1,nst
518 itrc=idsed(ised)
519 cff1=1.0_r8/srho(ised,ng)
520 sedden=sedden+ &
521 & t(i,j,k,nrhs,itrc)* &
522 & (srho(ised,ng)-den(i,k))*cff1
523 tl_sedden=tl_sedden+ &
524 & (tl_t(i,j,k,nrhs,itrc)* &
525 & (srho(ised,ng)-den(i,k))- &
526 & t(i,j,k,nrhs,itrc)* &
527 & tl_den(i,k))*cff1
528 END DO
529 den(i,k)=den(i,k)+sedden
530 tl_den(i,k)=tl_den(i,k)+tl_sedden
531# endif
532 den(i,k)=den(i,k)-1000.0_r8
533# ifdef MASKING
534 den(i,k)=den(i,k)*rmask(i,j)
535 tl_den(i,k)=tl_den(i,k)*rmask(i,j)
536# endif
537 END DO
538 END DO
539
540# ifdef VAR_RHO_2D_NOT_YET
541!
542!-----------------------------------------------------------------------
543! Compute vertical averaged density (rhoA) and density perturbation
544! (rhoS) used in barotropic pressure gradient.
545!-----------------------------------------------------------------------
546!
547 DO i=istrt,iendt
548 cff1=den(i,n(ng))*hz(i,j,n(ng))
549 tl_cff1=tl_den(i,n(ng))*hz(i,j,n(ng))+ &
550 & den(i,n(ng))*tl_hz(i,j,n(ng))
551 rhos(i,j)=0.5_r8*cff1*hz(i,j,n(ng))
552 tl_rhos(i,j)=0.5_r8*(tl_cff1*hz(i,j,n(ng))+ &
553 & cff1*tl_hz(i,j,n(ng)))
554 rhoa(i,j)=cff1
555 tl_rhoa(i,j)=tl_cff1
556 END DO
557 DO k=n(ng)-1,1,-1
558 DO i=istrt,iendt
559 cff1=den(i,k)*hz(i,j,k)
560 tl_cff1=tl_den(i,k)*hz(i,j,k)+ &
561 & den(i,k)*tl_hz(i,j,k)
562 rhos(i,j)=rhos(i,j)+hz(i,j,k)*(rhoa(i,j)+0.5_r8*cff1)
563 tl_rhos(i,j)=tl_rhos(i,j)+ &
564 & tl_hz(i,j,k)*(rhoa(i,j)+0.5_r8*cff1)+ &
565 & hz(i,j,k)*(tl_rhoa(i,j)+0.5_r8*tl_cff1)
566 rhoa(i,j)=rhoa(i,j)+cff1
567 tl_rhoa(i,j)=tl_rhoa(i,j)+tl_cff1
568 END DO
569 END DO
570 cff2=1.0_r8/rho0
571 DO i=istrt,iendt
572 cff1=1.0_r8/(z_w(i,j,n(ng))-z_w(i,j,0))
573 tl_cff1=-cff1*cff1*(tl_z_w(i,j,n(ng))-tl_z_w(i,j,0))
574!
575! Here we reverse the order of the NL and TL operations since an
576! intermeridiate value of rhoA and rhoS is needed because they are
577! recursive.
578!
579 tl_rhoa(i,j)=cff2*(tl_cff1*rhoa(i,j)+cff1*tl_rhoa(i,j))
580 rhoa(i,j)=cff2*cff1*rhoa(i,j)
581 tl_rhos(i,j)=2.0_r8*cff2* &
582 & cff1*(2.0_r8*tl_cff1*rhos(i,j)+ &
583 & cff1*tl_rhos(i,j))
584 rhos(i,j)=2.0_r8*cff1*cff1*cff2*rhos(i,j)
585 END DO
586# endif
587
588# if defined BV_FREQUENCY_NOT_YET
589!
590!-----------------------------------------------------------------------
591! Compute Brunt-Vaisala frequency (1/s2) at horizontal RHO-points
592! and vertical W-points:
593!
594! bvf = - g/rho d(rho)/d(z).
595!
596! The density anomaly difference is computed by lowering/rising the
597! water parcel above/below adiabatically at W-point depth "z_w".
598!-----------------------------------------------------------------------
599!
600 DO k=1,n(ng)-1
601 DO i=istrt,iendt
602 bulk_up=bulk0(i,k+1)- &
603 & z_w(i,j,k)*(bulk1(i,k+1)- &
604 & bulk2(i,k+1)*z_w(i,j,k))
605 tl_bulk_up=tl_bulk0(i,k+1)- &
606 & tl_z_w(i,j,k)*(bulk1(i,k+1)- &
607 & bulk2(i,k+1)*z_w(i,j,k))- &
608 & z_w(i,j,k)*(tl_bulk1(i,k+1)- &
609 & tl_bulk2(i,k+1)*z_w(i,j,k)- &
610 & bulk2(i,k+1)*tl_z_w(i,j,k))
611 bulk_dn=bulk0(i,k )- &
612 & z_w(i,j,k)*(bulk1(i,k )- &
613 & bulk2(i,k )*z_w(i,j,k))
614 tl_bulk_dn=tl_bulk0(i,k )- &
615 & tl_z_w(i,j,k)*(bulk1(i,k )- &
616 & bulk2(i,k )*z_w(i,j,k))- &
617 & z_w(i,j,k)*(tl_bulk1(i,k )- &
618 & tl_bulk2(i,k )*z_w(i,j,k)- &
619 & bulk2(i,k )*tl_z_w(i,j,k))
620 cff1=1.0_r8/(bulk_up+0.1_r8*z_w(i,j,k))
621 cff2=1.0_r8/(bulk_dn+0.1_r8*z_w(i,j,k))
622 tl_cff1=-cff1*cff1*(tl_bulk_up+0.1_r8*tl_z_w(i,j,k))
623 tl_cff2=-cff2*cff2*(tl_bulk_dn+0.1_r8*tl_z_w(i,j,k))
624 den_up=cff1*(den1(i,k+1)*bulk_up)
625 den_dn=cff2*(den1(i,k )*bulk_dn)
626 tl_den_up=tl_cff1*(den1(i,k+1)*bulk_up)+ &
627 & cff1*(tl_den1(i,k+1)*bulk_up+ &
628 & den1(i,k+1)*tl_bulk_up)
629 tl_den_dn=tl_cff2*(den1(i,k )*bulk_dn)+ &
630 & cff2*(tl_den1(i,k )*bulk_dn+ &
631 & den1(i,k )*tl_bulk_dn)
632!^ bvf(i,j,k)=-g*(den_up-den_dn)/ &
633!^ & (0.5_r8*(den_up+den_dn)* &
634!^ & (z_r(i,j,k+1)-z_r(i,j,k)))
635!^
636 cff3=1.0_r8/(0.5_r8*(den_up+den_dn)* &
637 & (z_r(i,j,k+1)-z_r(i,j,k)))
638 tl_cff3=-cff3*cff3* &
639 & 0.5_r8*((tl_den_up+tl_den_dn)* &
640 & (z_r(i,j,k+1)-z_r(i,j,k))+ &
641 & (den_up+den_dn)* &
642 & (tl_z_r(i,j,k+1)-tl_z_r(i,j,k)))
643 tl_bvf(i,j,k)=-g*((tl_den_up-tl_den_dn)*cff3+ &
644 & (den_up-den_dn)*tl_cff3)
645 END DO
646 END DO
647 DO i=istrt,iendt
648!^ bvf(i,j,0)=0.0_r8
649!^
650 tl_bvf(i,j,0)=0.0_r8
651!^ bvf(i,j,N(ng))=0.0_r8
652!^
653 tl_bvf(i,j,n(ng))=0.0_r8
654 END DO
655# endif
656
657# if defined LMD_SKPP_NOT_YET || defined LMD_BKPP_NOT_YET || \
658 defined bulk_fluxes
659!
660!-----------------------------------------------------------------------
661! Compute thermal expansion (1/Celsius) and saline contraction
662! (1/PSU) coefficients.
663!-----------------------------------------------------------------------
664!
665# ifdef LMD_DDMIX_NOT_YET
666 DO k=1,n(ng)
667# else
668 DO k=n(ng),n(ng)
669# endif
670 DO i=istrt,iendt
671 tpr10=0.1_r8*z_r(i,j,k)
672 tl_tpr10=0.1_r8*tl_z_r(i,j,k)
673!
674! Compute thermal expansion and saline contraction coefficients.
675!
676 cff=bulk(i,k)+tpr10
677 tl_cff=tl_bulk(i,k)+tl_tpr10
678 cff1=tpr10*den1(i,k)
679 tl_cff1=tl_tpr10*den1(i,k)+tpr10*tl_den1(i,k)
680 cff2=bulk(i,k)*cff
681 tl_cff2=tl_bulk(i,k)*cff+bulk(i,k)*tl_cff
682 wrk(i,k)=(den(i,k)+1000.0_r8)*cff*cff
683 tl_wrk(i,k)=cff*(cff*tl_den(i,k)+ &
684 & 2.0_r8*tl_cff*(den(i,k)+1000.0_r8))
685 tcof(i,k)=-(dbulkdt(i,k)*cff1+ &
686 & dden1dt(i,k)*cff2)
687 tl_tcof(i,k)=-(tl_dbulkdt(i,k)*cff1+ &
688 & dbulkdt(i,k)*tl_cff1+ &
689 & tl_dden1dt(i,k)*cff2+ &
690 & dden1dt(i,k)*tl_cff2)
691 scof(i,k)= (dbulkds(i,k)*cff1+ &
692 & dden1ds(i,k)*cff2)
693 tl_scof(i,k)= (tl_dbulkds(i,k)*cff1+ &
694 & dbulkds(i,k)*tl_cff1+ &
695 & tl_dden1ds(i,k)*cff2+ &
696 & dden1ds(i,k)*tl_cff2)
697# ifdef LMD_DDMIX_NOT_YET
698!^ alfaobeta(i,j,k)=Tcof(i,k)/Scof(i,k)
699!^
700 tl_alfaobeta(i,j,k)=(tl_tcof(i,k)*scof(i,k)- &
701 & tcof(i,k)*tl_scof(i,k))/ &
702 & (scof(i,k)*scof(i,k))
703# endif
704 END DO
705 IF (k.eq.n(ng)) THEN
706 DO i=istrt,iendt
707 cff=1.0_r8/wrk(i,n(ng))
708 tl_cff=-cff*cff*tl_wrk(i,n(ng))
709 alpha(i,j)=cff*tcof(i,n(ng))
710 tl_alpha(i,j)=tl_cff*tcof(i,n(ng))+cff*tl_tcof(i,n(ng))
711 beta(i,j)=cff*scof(i,n(ng))
712 tl_beta(i,j)=tl_cff*scof(i,n(ng))+cff*tl_scof(i,n(ng))
713 END DO
714 END IF
715 END DO
716# endif
717!
718!-----------------------------------------------------------------------
719! Load "in situ" density anomaly (kg/m3 - 1000) and potential
720! density anomaly (kg/m3 - 1000) referenced to the surface into global
721! arrays. Notice that this is done in a separate (i,k) DO-loops to
722! facilitate the adjoint.
723!-----------------------------------------------------------------------
724!
725 DO k=1,n(ng)
726 DO i=istrt,iendt
727 rho(i,j,k)=den(i,k)
728 tl_rho(i,j,k)=tl_den(i,k)
729# if defined LMD_SKPP_NOT_YET || defined LMD_BKPP_NOT_YET
730!^ pden(i,j,k)=(den1(i,k)-1000.0_r8) ! This gives a fatal
731!^ ! result in 4D-Var
732 tl_pden(i,j,k)=tl_den1(i,k) ! posterior error...
733# ifdef MASKING
734!^ pden(i,j,k)=pden(i,k)*rmask(i,j)
735!^
736 tl_pden(i,j,k)=tl_pden(i,k)*rmask(i,j)
737# endif
738# endif
739 END DO
740 END DO
741 END DO
742!
743!-----------------------------------------------------------------------
744! Exchange boundary data.
745!-----------------------------------------------------------------------
746!
747 IF (ewperiodic(ng).or.nsperiodic(ng)) THEN
748 CALL exchange_r3d_tile (ng, tile, &
749 & lbi, ubi, lbj, ubj, 1, n(ng), &
750 & rho)
751 CALL exchange_r3d_tile (ng, tile, &
752 & lbi, ubi, lbj, ubj, 1, n(ng), &
753 & tl_rho)
754
755# if defined LMD_SKPP_NOT_YET || defined LMD_BKPP_NOT_YET
756!^ CALL exchange_r3d_tile (ng, tile, &
757!^ & LBi, UBi, LBj, UBj, 1, N(ng), &
758!^ & pden)
759!^
760 CALL exchange_r3d_tile (ng, tile, &
761 & lbi, ubi, lbj, ubj, 1, n(ng), &
762 & tl_pden)
763# endif
764
765# if defined LMD_SKPP_NOT_YET || defined LMD_BKPP_NOT_YET || \
766 defined bulk_fluxes_not_yet
767# ifdef LMD_DDMIX_NOT_YET
768!^ CALL exchange_w3d_tile (ng, tile, &
769!^ & LBi, UBi, LBj, UBj, 0, N(ng), &
770!^ & alfaobeta)
771!^
772 CALL exchange_w3d_tile (ng, tile, &
773 & lbi, ubi, lbj, ubj, 0, n(ng), &
774 & tl_alfaobeta)
775# endif
776 CALL exchange_r2d_tile (ng, tile, &
777 & lbi, ubi, lbj, ubj, &
778 & alpha)
779 CALL exchange_r2d_tile (ng, tile, &
780 & lbi, ubi, lbj, ubj, &
781 & tl_alpha)
782 CALL exchange_r2d_tile (ng, tile, &
783 & lbi, ubi, lbj, ubj, &
784 & beta)
785 CALL exchange_r2d_tile (ng, tile, &
786 & lbi, ubi, lbj, ubj, &
787 & tl_beta)
788# endif
789
790# ifdef VAR_RHO_2D_NOT_YET
791 CALL exchange_r2d_tile (ng, tile, &
792 & lbi, ubi, lbj, ubj, &
793 & rhoa)
794 CALL exchange_r2d_tile (ng, tile, &
795 & lbi, ubi, lbj, ubj, &
796 & tl_rhoa)
797 CALL exchange_r2d_tile (ng, tile, &
798 & lbi, ubi, lbj, ubj, &
799 & rhos)
800 CALL exchange_r2d_tile (ng, tile, &
801 & lbi, ubi, lbj, ubj, &
802 & tl_rhos)
803# endif
804
805# ifdef BV_FREQUENCY_NOT_YET
806!^ CALL exchange_w3d_tile (ng, tile, &
807!^ & LBi, UBi, LBj, UBj, 0, N(ng), &
808!^ & bvf)
809!^
810 CALL exchange_w3d_tile (ng, tile, &
811 & lbi, ubi, lbj, ubj, 0, n(ng), &
812 & tl_bvf)
813# endif
814 END IF
815
816# ifdef DISTRIBUTE
817!
818 CALL mp_exchange3d (ng, tile, model, 1, &
819 & lbi, ubi, lbj, ubj, 1, n(ng), &
820 & nghostpoints, &
821 & ewperiodic(ng), nsperiodic(ng), &
822 & rho)
823 CALL mp_exchange3d (ng, tile, model, 1, &
824 & lbi, ubi, lbj, ubj, 1, n(ng), &
825 & nghostpoints, &
826 & ewperiodic(ng), nsperiodic(ng), &
827 & tl_rho)
828
829# if defined LMD_SKPP_NOT_YET || defined LMD_BKPP_NOT_YET
830!^ CALL mp_exchange3d (ng, tile, model, 1, &
831!^ & LBi, UBi, LBj, UBj, 1, N(ng), &
832!^ & NghostPoints, &
833!^ & EWperiodic(ng), NSperiodic(ng), &
834!^ & pden)
835!^
836 CALL mp_exchange3d (ng, tile, model, 1, &
837 & lbi, ubi, lbj, ubj, 1, n(ng), &
838 & nghostpoints, &
839 & ewperiodic(ng), nsperiodic(ng), &
840 & tl_pden)
841# endif
842
843# if defined LMD_SKPP_NOT_YET || defined LMD_BKPP_NOT_YET || \
844 defined bulk_fluxes
845# ifdef LMD_DDMIX_NOT_YET
846!^ CALL mp_exchange3d (ng, tile, model, 1, &
847!^ & LBi, UBi, LBj, UBj, 0, N(ng), &
848!^ & NghostPoints, &
849!^ & EWperiodic(ng), NSperiodic(ng), &
850!^ & alfaobeta)
851!^
852 CALL mp_exchange3d (ng, tile, model, 1, &
853 & lbi, ubi, lbj, ubj, 0, n(ng), &
854 & nghostpoints, &
855 & ewperiodic(ng), nsperiodic(ng), &
856 & tl_alfaobeta)
857# endif
858 CALL mp_exchange2d (ng, tile, model, 2, &
859 & lbi, ubi, lbj, ubj, &
860 & nghostpoints, &
861 & ewperiodic(ng), nsperiodic(ng), &
862 & alpha, beta)
863 CALL mp_exchange2d (ng, tile, model, 2, &
864 & lbi, ubi, lbj, ubj, &
865 & nghostpoints, &
866 & ewperiodic(ng), nsperiodic(ng), &
867 & tl_alpha, tl_beta)
868# endif
869
870# ifdef VAR_RHO_2D_NOT_YET
871 CALL mp_exchange2d (ng, tile, model, 2, &
872 & lbi, ubi, lbj, ubj, &
873 & nghostpoints, &
874 & ewperiodic(ng), nsperiodic(ng), &
875 & rhoa, rhos)
876 CALL mp_exchange2d (ng, tile, model, 2, &
877 & lbi, ubi, lbj, ubj, &
878 & nghostpoints, &
879 & ewperiodic(ng), nsperiodic(ng), &
880 & tl_rhoa, tl_rhos)
881# endif
882
883# ifdef BV_FREQUENCY_NOT_YET
884!^ CALL mp_exchange3d (ng, tile, model, 1, &
885!^ & LBi, UBi, LBj, UBj, 0, N(ng), &
886!^ & NghostPoints, &
887!^ & EWperiodic(ng), NSperiodic(ng), &
888!^ & bvf)
889!^
890 CALL mp_exchange3d (ng, tile, model, 1, &
891 & lbi, ubi, lbj, ubj, 0, n(ng), &
892 & nghostpoints, &
893 & ewperiodic(ng), nsperiodic(ng), &
894 & tl_bvf)
895# endif
896# endif
897!
898 RETURN
899 END SUBROUTINE tl_rho_eos_tile
900# endif
901
902# ifndef NONLIN_EOS
903!
904!***********************************************************************
905 SUBROUTINE tl_rho_eos_tile (ng, tile, model, &
906 & LBi, UBi, LBj, UBj, &
907 & IminS, ImaxS, JminS, JmaxS, &
908 & nrhs, &
909# ifdef MASKING
910 & rmask, &
911# endif
912# ifdef VAR_RHO_2D_NOT_YET
913 & Hz, tl_Hz, &
914# endif
915 & z_r, tl_z_r, &
916 & z_w, tl_z_w, &
917 & t, tl_t, &
918# ifdef VAR_RHO_2D_NOT_YET
919 & rhoA, tl_rhoA, &
920 & rhoS, tl_rhoS, &
921# endif
922# ifdef BV_FREQUENCY_NOT_YET
923 & tl_bvf, &
924# endif
925# if defined LMD_SKPP_NOT_YET || defined LMD_BKPP_NOT_YET || \
926 defined BULK_FLUXES
927 & alpha, tl_alpha, &
928 & beta, tl_beta, &
929# ifdef LMD_DDMIX_NOT_YET
930 & tl_alfaobeta, &
931# endif
932# endif
933# if defined LMD_SKPP_NOT_YET || defined LMD_BKPP_NOT_YET
934 & tl_pden, &
935# endif
936 & rho, tl_rho)
937!***********************************************************************
938!
939 USE mod_param
940 USE mod_scalars
941# ifdef SEDIMENT_NOT_YET
942 USE mod_sediment
943# endif
944!
947# ifdef DISTRIBUTE
949# endif
950!
951! Imported variable declarations.
952!
953 integer, intent(in) :: ng, tile, model
954 integer, intent(in) :: LBi, UBi, LBj, UBj
955 integer, intent(in) :: IminS, ImaxS, JminS, JmaxS
956 integer, intent(in) :: nrhs
957!
958# ifdef ASSUMED_SHAPE
959# ifdef MASKING
960 real(r8), intent(in) :: rmask(LBi:,LBj:)
961# endif
962# ifdef VAR_RHO_2D_NOT_YET
963 real(r8), intent(in) :: Hz(LBi:,LBj:,:)
964# endif
965 real(r8), intent(in) :: z_r(LBi:,LBj:,:)
966 real(r8), intent(in) :: z_w(LBi:,LBj:,0:)
967 real(r8), intent(in) :: t(LBi:,LBj:,:,:,:)
968# ifdef VAR_RHO_2D_NOT_YET
969 real(r8), intent(in) :: tl_Hz(LBi:,LBj:,:)
970# endif
971 real(r8), intent(in) :: tl_z_r(LBi:,LBj:,:)
972 real(r8), intent(in) :: tl_z_w(LBi:,LBj:,0:)
973 real(r8), intent(in) :: tl_t(LBi:,LBj:,:,:,:)
974# if defined LMD_SKPP_NOT_YET || defined LMD_BKPP_NOT_YET || \
975 defined bulk_fluxes
976 real(r8), intent(inout) :: alpha(LBi:,LBj:)
977 real(r8), intent(inout) :: beta(LBi:,LBj:)
978# endif
979# ifdef VAR_RHO_2D_NOT_YET
980 real(r8), intent(out) :: rhoA(LBi:,LBj:)
981 real(r8), intent(out) :: rhoS(LBi:,LBj:)
982 real(r8), intent(out) :: tl_rhoA(LBi:,LBj:)
983 real(r8), intent(out) :: tl_rhoS(LBi:,LBj:)
984# endif
985# ifdef BV_FREQUENCY_NOT_YET
986 real(r8), intent(out) :: tl_bvf(LBi:,LBj:,0:)
987# endif
988# if defined LMD_SKPP_NOT_YET || defined LMD_BKPP_NOT_YET || \
989 defined bulk_fluxes
990 real(r8), intent(out) :: tl_alpha(LBi:,LBj:)
991 real(r8), intent(out) :: tl_beta(LBi:,LBj:)
992# ifdef LMD_DDMIX_NOT_YET
993 real(r8), intent(out) :: tl_alfaobeta(LBi:,LBj:,0:)
994# endif
995# endif
996# if defined LMD_SKPP_NOT_YET || defined LMD_BKPP_NOT_YET
997 real(r8), intent(out) :: tl_pden(LBi:,LBj:,:)
998# endif
999 real(r8), intent(out) :: rho(LBi:,LBj:,:)
1000 real(r8), intent(out) :: tl_rho(LBi:,LBj:,:)
1001# else
1002# ifdef MASKING
1003 real(r8), intent(in) :: rmask(LBi:UBi,LBj:UBj)
1004# endif
1005# ifdef VAR_RHO_2D_NOT_YET
1006 real(r8), intent(in) :: Hz(LBi:UBi,LBj:UBj,N(ng))
1007# endif
1008 real(r8), intent(in) :: z_r(LBi:UBi,LBj:UBj,N(ng))
1009 real(r8), intent(in) :: z_w(LBi:UBi,LBj:UBj,0:N(ng))
1010 real(r8), intent(in) :: t(LBi:UBi,LBj:UBj,N(ng),3,NT(ng))
1011# ifdef VAR_RHO_2D_NOT_YET
1012 real(r8), intent(in) :: tl_Hz(LBi:UBi,LBj:UBj,N(ng))
1013# endif
1014 real(r8), intent(in) :: tl_z_r(LBi:UBi,LBj:UBj,N(ng))
1015 real(r8), intent(in) :: tl_z_w(LBi:UBi,LBj:UBj,0:N(ng))
1016 real(r8), intent(in) :: tl_t(LBi:UBi,LBj:UBj,N(ng),3,NT(ng))
1017# if defined LMD_SKPP_NOT_YET || defined LMD_BKPP_NOT_YET || \
1018 defined bulk_fluxes
1019 real(r8), intent(inout) :: alpha(LBi:UBi,LBj:UBj)
1020 real(r8), intent(inout) :: beta(LBi:UBi,LBj:UBj)
1021# endif
1022# ifdef VAR_RHO_2D_NOT_YET
1023 real(r8), intent(out) :: rhoA(LBi:UBi,LBj:UBj)
1024 real(r8), intent(out) :: rhoS(LBi:UBi,LBj:UBj)
1025 real(r8), intent(out) :: tl_rhoA(LBi:UBi,LBj:UBj)
1026 real(r8), intent(out) :: tl_rhoS(LBi:UBi,LBj:UBj)
1027# endif
1028# ifdef BV_FREQUENCY_NOT_YET
1029 real(r8), intent(out) :: tl_bvf(LBi:UBi,LBj:UBj,0:N(ng))
1030# endif
1031# if defined LMD_SKPP_NOT_YET || defined LMD_BKPP_NOT_YET || \
1032 defined bulk_fluxes
1033 real(r8), intent(out) :: tl_alpha(LBi:UBi,LBj:UBj)
1034 real(r8), intent(out) :: tl_beta(LBi:UBi,LBj:UBj)
1035# ifdef LMD_DDMIX_NOT_YET
1036 real(r8), intent(out) :: tl_alfaobeta(LBi:UBi,LBj:UBj,0:N(ng))
1037# endif
1038# endif
1039# if defined LMD_SKPP_NOT_YET || defined LMD_BKPP_NOT_YET
1040 real(r8), intent(out) :: tl_pden(LBi:UBi,LBj:UBj,N(ng))
1041# endif
1042 real(r8), intent(out) :: rho(LBi:UBi,LBj:UBj,N(ng))
1043 real(r8), intent(out) :: tl_rho(LBi:UBi,LBj:UBj,N(ng))
1044# endif
1045!
1046! Local variable declarations.
1047!
1048 integer :: i, ised, itrc, j, k
1049
1050 real(r8) :: SedDen, cff, cff1, cff2
1051 real(r8) :: tl_SedDen, tl_cff, tl_cff1
1052
1053# include "set_bounds.h"
1054!
1055!=======================================================================
1056! Tangent linear, linear equation of state.
1057!=======================================================================
1058!
1059!-----------------------------------------------------------------------
1060! Compute "in situ" density anomaly (kg/m3 - 1000) using the linear
1061! equation of state.
1062!-----------------------------------------------------------------------
1063!
1064 DO j=jstrt,jendt
1065 DO k=1,n(ng)
1066 DO i=istrt,iendt
1067 rho(i,j,k)=r0(ng)- &
1068 & r0(ng)*tcoef(ng)*(t(i,j,k,nrhs,itemp)-t0(ng))
1069 tl_rho(i,j,k)=-r0(ng)*tcoef(ng)*tl_t(i,j,k,nrhs,itemp)
1070# ifdef SALINITY
1071 rho(i,j,k)=rho(i,j,k)+ &
1072 & r0(ng)*scoef(ng)*(t(i,j,k,nrhs,isalt)-s0(ng))
1073 tl_rho(i,j,k)=tl_rho(i,j,k)+ &
1074 & r0(ng)*scoef(ng)*tl_t(i,j,k,nrhs,isalt)
1075# endif
1076# if defined SEDIMENT_NOT_YET && defined SED_DENS_NOT_YET
1077 sedden=0.0_r8
1078 tl_sedden=0.0_r8
1079 DO ised=1,nst
1080 itrc=idsed(ised)
1081 cff1=1.0_r8/srho(ised,ng)
1082 sedden=sedden+ &
1083 & t(i,j,k,nrhs,itrc)* &
1084 & (srho(ised,ng)-rho(i,j,k))*cff1
1085 tl_sedden=tl_sedden+ &
1086 & (tl_t(i,j,k,nrhs,itrc)* &
1087 & (srho(ised,ng)-rho(i,j,k))- &
1088 & t(i,j,k,nrhs,itrc)* &
1089 & tl_rho(i,j,k))*cff1
1090 END DO
1091 rho(i,j,k)=rho(i,j,k)+sedden
1092 tl_rho(i,j,k)=tl_rho(i,j,k)+tl_sedden
1093# endif
1094 rho(i,j,k)=rho(i,j,k)-1000.0_r8
1095# ifdef MASKING
1096 rho(i,j,k)=rho(i,j,k)*rmask(i,j)
1097 tl_rho(i,j,k)=tl_rho(i,j,k)*rmask(i,j)
1098# endif
1099# if defined LMD_SKPP_NOT_YET || defined LMD_BKPP_NOT_YET
1100!^ pden(i,j,k)=rho(i,j,k)
1101!^
1102 tl_pden(i,j,k)=tl_rho(i,j,k)
1103# endif
1104 END DO
1105 END DO
1106
1107# ifdef VAR_RHO_2D_NOT_YET
1108!
1109!-----------------------------------------------------------------------
1110! Compute vertical averaged density (rhoA) and density perturbation
1111! used (rhoS) in barotropic pressure gradient.
1112!-----------------------------------------------------------------------
1113!
1114 DO i=istrt,iendt
1115 cff1=rho(i,j,n(ng))*hz(i,j,n(ng))
1116 tl_cff1=tl_rho(i,j,n(ng))*hz(i,j,n(ng))+ &
1117 & rho(i,j,n(ng))*tl_hz(i,j,n(ng))
1118 rhos(i,j)=0.5_r8*cff1*hz(i,j,n(ng))
1119 tl_rhos(i,j)=0.5_r8*(tl_cff1*hz(i,j,n(ng))+ &
1120 & cff1*tl_hz(i,j,n(ng)))
1121 rhoa(i,j)=cff1
1122 tl_rhoa(i,j)=tl_cff1
1123 END DO
1124 DO k=n(ng)-1,1,-1
1125 DO i=istrt,iendt
1126 cff1=rho(i,j,k)*hz(i,j,k)
1127 tl_cff1=tl_rho(i,j,k)*hz(i,j,k)+ &
1128 & rho(i,j,k)*tl_hz(i,j,k)
1129 rhos(i,j)=rhos(i,j)+hz(i,j,k)*(rhoa(i,j)+0.5_r8*cff1)
1130 tl_rhos(i,j)=tl_rhos(i,j)+ &
1131 & tl_hz(i,j,k)*(rhoa(i,j)+0.5_r8*cff1)+ &
1132 & hz(i,j,k)*(tl_rhoa(i,j)+0.5_r8*tl_cff1)
1133 rhoa(i,j)=rhoa(i,j)+cff1
1134 tl_rhoa(i,j)=tl_rhoa(i,j)+tl_cff1
1135 END DO
1136 END DO
1137 cff2=1.0_r8/rho0
1138 DO i=istrt,iendt
1139 cff1=1.0_r8/(z_w(i,j,n(ng))-z_w(i,j,0))
1140 tl_cff1=-cff1*cff1*(tl_z_w(i,j,n(ng))-tl_z_w(i,j,0))
1141!
1142! Here we reverse the order of the NL and TL operations since an
1143! intermeridiate value of rhoA and rhoS is needed because they are
1144! recursive.
1145!
1146 tl_rhoa(i,j)=cff2*(tl_cff1*rhoa(i,j)+cff1*tl_rhoa(i,j))
1147 rhoa(i,j)=cff2*cff1*rhoa(i,j)
1148 tl_rhos(i,j)=2.0_r8*cff2* &
1149 & cff1*(2.0_r8*tl_cff1*rhos(i,j)+ &
1150 & cff1*tl_rhos(i,j))
1151 rhos(i,j)=2.0_r8*cff1*cff1*cff2*rhos(i,j)
1152 END DO
1153# endif
1154
1155# ifdef BV_FREQUENCY_NOT_YET
1156!
1157!-----------------------------------------------------------------------
1158! Compute Brunt-Vaisala frequency (1/s2) at horizontal RHO-points
1159! and vertical W-points.
1160!-----------------------------------------------------------------------
1161!
1162 DO k=1,n(ng)-1
1163 DO i=istrt,iendt
1164!^ bvf(i,j,k)=-gorho0*(rho(i,j,k+1)-rho(i,j,k))/ &
1165!^ & (z_r(i,j,k+1)-z_r(i,j,k))
1166!^
1167 cff=1.0_r8/(z_r(i,j,k+1)-z_r(i,j,k))
1168 tl_cff=-cff*cff*(tl_z_r(i,j,k+1)-tl_z_r(i,j,k))
1169 tl_bvf(i,j,k)=-gorho0* &
1170 & (cff*(tl_rho(i,j,k+1)-tl_rho(i,j,k))+ &
1171 & tl_cff*(rho(i,j,k+1)-rho(i,j,k)))
1172 END DO
1173 END DO
1174# endif
1175
1176# if defined LMD_SKPP_NOT_YET || defined LMD_BKPP_NOT_YET || \
1177 defined bulk_fluxes
1178!
1179!-----------------------------------------------------------------------
1180! Compute thermal expansion (1/Celsius) and saline contraction
1181! (1/PSU) coefficients.
1182!-----------------------------------------------------------------------
1183!
1184 DO i=istrt,iendt
1185 alpha(i,j)=abs(tcoef(ng))
1186 tl_alpha(i,j)=0.0_r8
1187# ifdef SALINITY
1188 beta(i,j)=abs(scoef(ng))
1189 tl_beta(i,j)=0.0_r8
1190# else
1191 beta(i,j)=0.0_r8
1192 tl_beta(i,j)=0.0_r8
1193# endif
1194 END DO
1195# ifdef LMD_DDMIX_NOT_YET
1196!
1197! Compute ratio of thermal expansion and saline contraction
1198! coefficients.
1199!
1200 IF (scoef(ng).eq.0.0_r8) THEN
1201 cff=1.0_r8
1202 ELSE
1203 cff=1.0_r8/scoef(ng)
1204 END IF
1205 DO k=1,n(ng)
1206 DO i=istrt,iendt
1207!^ alfaobeta(i,j,k)=cff*Tcoef(ng)
1208!^
1209 tl_alfaobeta(i,j,k)=0.0_r8
1210 END DO
1211 END DO
1212# endif
1213# endif
1214 END DO
1215!
1216!-----------------------------------------------------------------------
1217! Exchange boundary data.
1218!-----------------------------------------------------------------------
1219!
1220 IF (ewperiodic(ng).or.nsperiodic(ng)) THEN
1221 CALL exchange_r3d_tile (ng, tile, &
1222 & lbi, ubi, lbj, ubj, 1, n(ng), &
1223 & rho)
1224 CALL exchange_r3d_tile (ng, tile, &
1225 & lbi, ubi, lbj, ubj, 1, n(ng), &
1226 & tl_rho)
1227
1228# if defined LMD_SKPP_NOT_YET || defined LMD_BKPP_NOT_YET
1229!^ CALL exchange_r3d_tile (ng, tile, &
1230!^ & LBi, UBi, LBj, UBj, 1, N(ng), &
1231!^ & pden)
1232!^
1233 CALL exchange_r3d_tile (ng, tile, &
1234 & lbi, ubi, lbj, ubj, 1, n(ng), &
1235 & tl_pden)
1236# endif
1237
1238# if defined LMD_SKPP_NOT_YET || defined LMD_BKPP_NOT_YET || \
1239 defined bulk_fluxes_not_yet
1240# ifdef LMD_DDMIX_NOT_YET
1241!^ CALL exchange_w3d_tile (ng, tile, &
1242!^ & LBi, UBi, LBj, UBj, 0, N(ng), &
1243!^ & alfaobeta)
1244!^
1245 CALL exchange_w3d_tile (ng, tile, &
1246 & lbi, ubi, lbj, ubj, 0, n(ng), &
1247 & tl_alfaobeta)
1248# endif
1249 CALL exchange_r2d_tile (ng, tile, &
1250 & lbi, ubi, lbj, ubj, &
1251 & alpha)
1252 CALL exchange_r2d_tile (ng, tile, &
1253 & lbi, ubi, lbj, ubj, &
1254 & tl_alpha)
1255 CALL exchange_r2d_tile (ng, tile, &
1256 & lbi, ubi, lbj, ubj, &
1257 & beta)
1258 CALL exchange_r2d_tile (ng, tile, &
1259 & lbi, ubi, lbj, ubj, &
1260 & tl_beta)
1261# endif
1262
1263# ifdef VAR_RHO_2D_NOT_YET
1264 CALL exchange_r2d_tile (ng, tile, &
1265 & lbi, ubi, lbj, ubj, &
1266 & rhoa)
1267 CALL exchange_r2d_tile (ng, tile, &
1268 & lbi, ubi, lbj, ubj, &
1269 & tl_rhoa)
1270 CALL exchange_r2d_tile (ng, tile, &
1271 & lbi, ubi, lbj, ubj, &
1272 & rhos)
1273 CALL exchange_r2d_tile (ng, tile, &
1274 & lbi, ubi, lbj, ubj, &
1275 & tl_rhos)
1276# endif
1277
1278# ifdef BV_FREQUENCY_NOT_YET
1279!^ CALL exchange_w3d_tile (ng, tile, &
1280!^ & LBi, UBi, LBj, UBj, 0, N(ng), &
1281!^ & bvf)
1282!^
1283 CALL exchange_w3d_tile (ng, tile, &
1284 & lbi, ubi, lbj, ubj, 0, n(ng), &
1285 & tl_bvf)
1286# endif
1287 END IF
1288
1289# ifdef DISTRIBUTE
1290!
1291 CALL mp_exchange3d (ng, tile, model, 1, &
1292 & lbi, ubi, lbj, ubj, 1, n(ng), &
1293 & nghostpoints, &
1294 & ewperiodic(ng), nsperiodic(ng), &
1295 & rho)
1296 CALL mp_exchange3d (ng, tile, model, 1, &
1297 & lbi, ubi, lbj, ubj, 1, n(ng), &
1298 & nghostpoints, &
1299 & ewperiodic(ng), nsperiodic(ng), &
1300 & tl_rho)
1301
1302# if defined LMD_SKPP_NOT_YET || defined LMD_BKPP_NOT_YET
1303!^ CALL mp_exchange3d (ng, tile, model, 1, &
1304!^ & LBi, UBi, LBj, UBj, 1, N(ng), &
1305!^ & NghostPoints, &
1306!^ & EWperiodic(ng), NSperiodic(ng), &
1307!^ & pden)
1308!^
1309 CALL mp_exchange3d (ng, tile, model, 1, &
1310 & lbi, ubi, lbj, ubj, 1, n(ng), &
1311 & nghostpoints, &
1312 & ewperiodic(ng), nsperiodic(ng), &
1313 & tl_pden)
1314# endif
1315
1316# if defined LMD_SKPP_NOT_YET || defined LMD_BKPP_NOT_YET || \
1317 defined bulk_fluxes
1318# ifdef LMD_DDMIX_NOT_YET
1319!^ CALL mp_exchange3d (ng, tile, model, 1, &
1320!^ & LBi, UBi, LBj, UBj, 0, N(ng), &
1321!^ & NghostPoints, &
1322!^ & EWperiodic(ng), NSperiodic(ng), &
1323!^ & alfaobeta)
1324!^
1325 CALL mp_exchange3d (ng, tile, model, 1, &
1326 & lbi, ubi, lbj, ubj, 0, n(ng), &
1327 & nghostpoints, &
1328 & ewperiodic(ng), nsperiodic(ng), &
1329 & tl_alfaobeta)
1330# endif
1331 CALL mp_exchange2d (ng, tile, model, 2, &
1332 & lbi, ubi, lbj, ubj, &
1333 & nghostpoints, &
1334 & ewperiodic(ng), nsperiodic(ng), &
1335 & alpha, beta)
1336 CALL mp_exchange2d (ng, tile, model, 2, &
1337 & lbi, ubi, lbj, ubj, &
1338 & nghostpoints, &
1339 & ewperiodic(ng), nsperiodic(ng), &
1340 & tl_alpha, tl_beta)
1341# endif
1342
1343# ifdef VAR_RHO_2D_NOT_YET
1344 CALL mp_exchange2d (ng, tile, model, 2, &
1345 & lbi, ubi, lbj, ubj, &
1346 & nghostpoints, &
1347 & ewperiodic(ng), nsperiodic(ng), &
1348 & rhoa, rhos)
1349 CALL mp_exchange2d (ng, tile, model, 2, &
1350 & lbi, ubi, lbj, ubj, &
1351 & nghostpoints, &
1352 & ewperiodic(ng), nsperiodic(ng), &
1353 & tl_rhoa, tl_rhos)
1354# endif
1355
1356# ifdef BV_FREQUENCY_NOT_YET
1357!^ CALL mp_exchange3d (ng, tile, model, 1, &
1358!^ & LBi, UBi, LBj, UBj, 0, N(ng), &
1359!^ & NghostPoints, &
1360!^ & EWperiodic(ng), NSperiodic(ng), &
1361!^ & bvf)
1362!^
1363 CALL mp_exchange3d (ng, tile, model, 1, &
1364 & lbi, ubi, lbj, ubj, 0, n(ng), &
1365 & nghostpoints, &
1366 & ewperiodic(ng), nsperiodic(ng), &
1367 & tl_bvf)
1368# endif
1369# endif
1370!
1371 RETURN
1372 END SUBROUTINE tl_rho_eos_tile
1373# endif
1374#endif
1375 END MODULE tl_rho_eos_mod
subroutine exchange_r2d_tile(ng, tile, lbi, ubi, lbj, ubj, a)
subroutine exchange_r3d_tile(ng, tile, lbi, ubi, lbj, ubj, lbk, ubk, a)
subroutine exchange_w3d_tile(ng, tile, lbi, ubi, lbj, ubj, lbk, ubk, a)
type(t_coupling), dimension(:), allocatable coupling
real(r8), parameter b03
Definition mod_eoscoef.F:32
real(r8), parameter e02
Definition mod_eoscoef.F:38
real(r8), parameter h02
Definition mod_eoscoef.F:49
real(r8), parameter q03
Definition mod_eoscoef.F:53
real(r8), parameter a03
Definition mod_eoscoef.F:27
real(r8), parameter q04
Definition mod_eoscoef.F:54
real(r8), parameter q05
Definition mod_eoscoef.F:55
real(r8), parameter u01
Definition mod_eoscoef.F:57
real(r8), parameter f00
Definition mod_eoscoef.F:40
real(r8), parameter v01
Definition mod_eoscoef.F:62
real(r8), parameter g01
Definition mod_eoscoef.F:44
real(r8), parameter a04
Definition mod_eoscoef.F:28
real(r8), parameter f01
Definition mod_eoscoef.F:41
real(r8), parameter g00
Definition mod_eoscoef.F:43
real(r8), parameter u02
Definition mod_eoscoef.F:58
real(r8), parameter d01
Definition mod_eoscoef.F:34
real(r8), parameter w00
Definition mod_eoscoef.F:64
real(r8), parameter b00
Definition mod_eoscoef.F:29
real(r8), parameter u03
Definition mod_eoscoef.F:59
real(r8), parameter d00
Definition mod_eoscoef.F:33
real(r8), parameter d02
Definition mod_eoscoef.F:35
real(r8), parameter g03
Definition mod_eoscoef.F:46
real(r8), parameter b01
Definition mod_eoscoef.F:30
real(r8), parameter u04
Definition mod_eoscoef.F:60
real(r8), parameter u00
Definition mod_eoscoef.F:56
real(r8), parameter v02
Definition mod_eoscoef.F:63
real(r8), parameter q01
Definition mod_eoscoef.F:51
real(r8), parameter e00
Definition mod_eoscoef.F:36
real(r8), parameter h00
Definition mod_eoscoef.F:47
real(r8), parameter g02
Definition mod_eoscoef.F:45
real(r8), parameter e03
Definition mod_eoscoef.F:39
real(r8), parameter a00
Definition mod_eoscoef.F:24
real(r8), parameter a01
Definition mod_eoscoef.F:25
real(r8), parameter a02
Definition mod_eoscoef.F:26
real(r8), parameter b02
Definition mod_eoscoef.F:31
real(r8), parameter h01
Definition mod_eoscoef.F:48
real(r8), parameter e01
Definition mod_eoscoef.F:37
real(r8), parameter f02
Definition mod_eoscoef.F:42
real(r8), parameter v00
Definition mod_eoscoef.F:61
real(r8), parameter q02
Definition mod_eoscoef.F:52
real(r8), parameter q00
Definition mod_eoscoef.F:50
type(t_grid), dimension(:), allocatable grid
Definition mod_grid.F:365
type(t_mixing), dimension(:), allocatable mixing
Definition mod_mixing.F:399
type(t_ocean), dimension(:), allocatable ocean
Definition mod_ocean.F:351
integer nghostpoints
Definition mod_param.F:710
integer nst
Definition mod_param.F:521
real(r8), dimension(:), allocatable t0
logical, dimension(:), allocatable ewperiodic
logical, dimension(:), allocatable nsperiodic
real(r8), dimension(:), allocatable s0
real(r8), dimension(:), allocatable tcoef
real(dp) gorho0
integer isalt
integer itemp
real(dp) g
real(dp) rho0
real(r8), dimension(:), allocatable r0
real(r8), dimension(:), allocatable scoef
real(r8), dimension(:,:), allocatable srho
integer, dimension(:), allocatable idsed
integer, dimension(:), allocatable nrhs
subroutine mp_exchange2d(ng, tile, model, nvar, lbi, ubi, lbj, ubj, nghost, ew_periodic, ns_periodic, a, b, c, d)
subroutine mp_exchange3d(ng, tile, model, nvar, lbi, ubi, lbj, ubj, lbk, ubk, nghost, ew_periodic, ns_periodic, a, b, c, d)
subroutine tl_rho_eos_tile(ng, tile, model, lbi, ubi, lbj, ubj, imins, imaxs, jmins, jmaxs, nrhs, rmask, z_r, tl_z_r, z_w, tl_z_w, t, tl_t, alpha, tl_alpha, beta, tl_beta, rho, tl_rho)
Definition tl_rho_eos.F:154
subroutine, public tl_rho_eos(ng, tile, model)
Definition tl_rho_eos.F:48
recursive subroutine wclock_off(ng, model, region, line, routine)
Definition timers.F:148
recursive subroutine wclock_on(ng, model, region, line, routine)
Definition timers.F:3