My Project
eqs_of_state.f90
Go to the documentation of this file.
1 
2 
3 
4 
5 
6 
7 
8 
9 
10 
11 
12 !/===========================================================================/
13 ! Copyright (c) 2007, The University of Massachusetts Dartmouth
14 ! Produced at the School of Marine Science & Technology
15 ! Marine Ecosystem Dynamics Modeling group
16 ! All rights reserved.
17 !
18 ! FVCOM has been developed by the joint UMASSD-WHOI research team. For
19 ! details of authorship and attribution of credit please see the FVCOM
20 ! technical manual or contact the MEDM group.
21 !
22 !
23 ! This file is part of FVCOM. For details, see http://fvcom.smast.umassd.edu
24 ! The full copyright notice is contained in the file COPYRIGHT located in the
25 ! root directory of the FVCOM code. This original header must be maintained
26 ! in all distributed versions.
27 !
28 ! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
29 ! AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO,
30 ! THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
31 ! PURPOSE ARE DISCLAIMED.
32 !
33 !/---------------------------------------------------------------------------/
34 ! CVS VERSION INFORMATION
35 ! $Id$
36 ! $Name$
37 ! $Revision$
38 !/===========================================================================/
39 
40 MODULE eqs_of_state
41 
42 
43  ! THIS MODULE CONTAINS THREE EQUATIONS OF STATE WHICH CAN BE ACCESSED
44  ! USING A NUMBER OF DIFFERENT INTERFACES. THE ORIGINAL FVCOM
45  ! INTERFACE (DENS1,DENS2,DENS3) STILL EXISTS, WHERE BOTH RHO AND
46  ! RHO1 ARE UPDATED, OR THE USER CAN ACCESS THE UNDERLIEING
47  ! FUNCTIONS AND PASS AN ARRAY, 1D OR 2D TO THE SUBROUTINE.
48  !
49  ! SEE DETAILS OF METHOD IN THE HEADER FOR EACH SUBROUTINE
50  !
51  ! SUBROUTINES PUBLICLY AVAILABLE IN THIS MODULE:
52  !=================================================================================
53  ! DENS1 - INTERFACE: CALL DENS1
54  ! DENS2 - INTERFACE: CALL DENS2
55  ! DENS3 - INTERFACE: CALL DENS3
56  ! - THESE ROUTINES USE THE VALUES IN S1,T1,ZZ TO UPDATE THE DENSITY IN RHO AND RHO1
57  !==================================================================================
58  ! FOFONOFF_MILLARD - INTERFACE: CALL FOFONOFF_MILLARD_2D(S,T,Z,PREF,RHO)
59  ! REAL(SP), INTENT(IN),DIMENSION(:) :: S,T,P ! SALINITY,TEMPERATURE, PRESSURE
60  ! REAL(SP), INTENT(OUT),DIMENSION(:) :: RHO ! RESULT - DENSITY
61  ! REAL(SP), INTENT(IN) :: PREF ! REFERENCE PRESSURE
62  ! or
63  ! REAL(SP), INTENT(IN),DIMENSION(:,:) :: S,T,P ! SALINITY,TEMPERATURE, PRESSURE
64  ! REAL(SP), INTENT(OUT),DIMENSION(:,:) :: RHO ! RESULT - DENSITY
65  ! REAL(SP), INTENT(IN) :: PREF ! REFERENCE PRESSURE
66  !=================================================================================
67  !
68  !=================================================================================
69  ! DENS2G - INTERFACE: CALL DENS2G(S,T,RHO)
70  ! REAL(SP), INTENT(IN),DIMENSION(:) :: S,T ! SALINITY,TEMPERATURE
71  ! REAL(SP), INTENT(OUT),DIMENSION(:) :: RHO ! DENSITY
72  ! or
73  ! REAL(SP), INTENT(IN),DIMENSION(:,:) :: S,T,P ! SALINITY,TEMPERATURE
74  ! REAL(SP), INTENT(OUT),DIMENSION(:,:) :: RHO ! DENSITY
75  !=================================================================================
76  !
77  !=================================================================================
78  ! JACKET_MCDOUGALL - INTERFACE: CALL JACKET_MCDOUGALL(S,T,Z,RHO)
79  ! REAL(SP), INTENT(IN),DIMENSION(:) :: S,T,P ! SALINITY,TEMPERATURE, PRESSURE
80  ! REAL(SP), INTENT(OUT),DIMENSION(:) :: RHO ! RESULT - DENSITY
81  ! or
82  ! REAL(SP), INTENT(IN),DIMENSION(:,:) :: S,T,P ! SALINITY,TEMPERATURE, PRESSURE
83  ! REAL(SP), INTENT(OUT),DIMENSION(:,:) :: RHO ! RESULT - DENSITY
84  !=================================================================================
85  !=================================================================================
86  !=================================================================================
87 
88  USE all_vars
89  USE mod_utils
90  IMPLICIT NONE
91 
92  PUBLIC
93 
94  INTERFACE fofonoff_millard
95  MODULE PROCEDURE fofonoff_millard_1d
96  MODULE PROCEDURE fofonoff_millard_2d
97  END INTERFACE
98 
99  INTERFACE dens2g
100  MODULE PROCEDURE dens2_1d
101  MODULE PROCEDURE dens2_2d
102  END INTERFACE
103 
104  INTERFACE jacket_mcdougall
105  MODULE PROCEDURE jacket_mcdougall_1d
106  MODULE PROCEDURE jacket_mcdougall_2d
107  END INTERFACE
108 
109  PRIVATE jacket_mcdougall_1d
110  PRIVATE jacket_mcdougall_2d
111 
112  PRIVATE dens2_1d
113  PRIVATE dens2_2d
114 
115  PRIVATE fofonoff_millard_1d
116  PRIVATE fofonoff_millard_2d
117 
118  PRIVATE svan
119  PRIVATE theta
120  PRIVATE atg
121 
122 CONTAINS
123 
124  !==============================================================================|
125  ! Calculate Potential Density Based on Potential Temp and Salinity |
126  ! Pressure effects are incorported (Can Model Fresh Water < 4 Deg C) |
127  ! Ref: algorithms for computation of fundamental properties of |
128  ! seawater , Fofonoff and Millard. |
129  ! |
130  ! calculates: rho1(nnode) density at nodes |
131  ! calculates: rho (ncell) density at elements |
132  !==============================================================================|
133  SUBROUTINE fofonoff_millard_2d(MYS,MYT,MYP,PREF, MYRHO)
134 
135  !------------------------------------------------------------------------------|
136  IMPLICIT NONE
137  REAL(SP), INTENT(IN), DIMENSION(:,:) :: MYS,MYT,MYP
138  REAL(SP), INTENT(IN) :: PREF
139  REAL(SP), INTENT(INOUT), DIMENSION(:,:):: MYRHO
140  INTEGER :: I,K,ub1,lb1,ub2,lb2
141  REAL(SP) ::PT,SVA,SIGMA
142  !==============================================================================|
143 
144  ! SET DIMS USING MY S
145  lb1 = lbound(mys,1)
146  ub1 = ubound(mys,1)
147 
148  lb2 = lbound(mys,2)
149  ub2 = ubound(mys,2)
150 
151  ! CHECK MY T
152  if(lb1 /= lbound(myt,1)) CALL fatal_error("EQS_OF_STATE:: Dimension mismatch!")
153  if(ub1 /= ubound(myt,1)) CALL fatal_error("EQS_OF_STATE:: Dimension mismatch!")
154 
155  if(lb2 /= lbound(myt,2)) CALL fatal_error("EQS_OF_STATE:: Dimension mismatch!")
156  if(ub2 /= ubound(myt,2)) CALL fatal_error("EQS_OF_STATE:: Dimension mismatch!")
157 
158  ! CHECK MY P
159  if(lb1 /= lbound(myp,1)) CALL fatal_error("EQS_OF_STATE:: Dimension mismatch!")
160  if(ub1 /= ubound(myp,1)) CALL fatal_error("EQS_OF_STATE:: Dimension mismatch!")
161 
162  if(lb2 /= lbound(myp,2)) CALL fatal_error("EQS_OF_STATE:: Dimension mismatch!")
163  if(ub2 /= ubound(myp,2)) CALL fatal_error("EQS_OF_STATE:: Dimension mismatch!")
164 
165  ! CHECK MY RHO
166  if(lb1 /= lbound(myrho,1)) CALL fatal_error("EQS_OF_STATE:: Dimension mismatch!")
167  if(ub1 /= ubound(myrho,1)) CALL fatal_error("EQS_OF_STATE:: Dimension mismatch!")
168 
169  if(lb2 /= lbound(myrho,2)) CALL fatal_error("EQS_OF_STATE:: Dimension mismatch!")
170  if(ub2 /= ubound(myrho,2)) CALL fatal_error("EQS_OF_STATE:: Dimension mismatch!")
171 
172 
173 
174  DO i=lb1,ub1
175  DO k=lb2,ub2
176 
177  pt = theta(mys(i,k),myt(i,k),myp(i,k),pref)
178  sva = svan(mys(i,k),pt,myp(i,k),sigma)
179  myrho(i,k) = sigma*1.e-3_sp
180  END DO
181  END DO
182 
183  !----------------TRANSFORM TO FACE CENTER--------------------------------------
184 
185  END SUBROUTINE fofonoff_millard_2d
186  !==============================================================================|
187  SUBROUTINE fofonoff_millard_1d(MYS,MYT,MYP,PREF, MYRHO)
188 
189  !------------------------------------------------------------------------------|
190  IMPLICIT NONE
191  REAL(SP), INTENT(IN), DIMENSION(:) :: MYS,MYT,MYP
192  REAL(SP), INTENT(IN) :: PREF
193  REAL(SP), INTENT(INOUT), DIMENSION(:):: MYRHO
194  INTEGER :: I,K,ub1,lb1
195  REAL(SP) ::PT,SVA,SIGMA
196  !==============================================================================|
197 
198  ! SET DIMS USING MY S
199  lb1 = lbound(mys,1)
200  ub1 = ubound(mys,1)
201 
202  ! CHECK MY T
203  if(lb1 /= lbound(myt,1)) CALL fatal_error("EQS_OF_STATE:: Dimension mismatch!")
204  if(ub1 /= ubound(myt,1)) CALL fatal_error("EQS_OF_STATE:: Dimension mismatch!")
205 
206  ! CHECK MY P
207  if(lb1 /= lbound(myp,1)) CALL fatal_error("EQS_OF_STATE:: Dimension mismatch!")
208  if(ub1 /= ubound(myp,1)) CALL fatal_error("EQS_OF_STATE:: Dimension mismatch!")
209 
210  ! CHECK MY RHO
211  if(lb1 /= lbound(myrho,1)) CALL fatal_error("EQS_OF_STATE:: Dimension mismatch!")
212  if(ub1 /= ubound(myrho,1)) CALL fatal_error("EQS_OF_STATE:: Dimension mismatch!")
213 
214  DO k=lb1,ub1
215 
216  pt = theta(mys(k),myt(k),myp(k),pref)
217  sva = svan(mys(k),pt,myp(k),sigma)
218  myrho(k) = sigma*1.e-3_sp
219  END DO
220 
221  !----------------TRANSFORM TO FACE CENTER--------------------------------------
222 
223  END SUBROUTINE fofonoff_millard_1d
224  !==============================================================================|
225 
226  ! GENERIC CALL FOR FVCOM
227  !==============================================================================|
228  SUBROUTINE dens1
229 
230  !------------------------------------------------------------------------------|
231  IMPLICIT NONE
232  INTEGER :: K
233  REAL(SP), PARAMETER ::PR = 0.0_sp
234  REAL(SP), DIMENSION(0:MT,1:KB) :: RZU
235  !==============================================================================|
236  IF(dbg_set(dbg_sbr)) WRITE(ipt,*) "START: DENS1"
237 
238  ! The thickness of the water column
239  ! Is not the depth realtive to z=0
240  DO k=1,kbm1
241  rzu(:,k) = -grav_n*1.025_sp*(zz(:,k)*d(:))*0.1_sp
242  END DO
243 
244  CALL fofonoff_millard_2d(s1,t1,rzu,pr,rho1)
245  rho1(:,kb)=0.0_sp
246  rho1(0,:)=0.0_sp
247 
248 
249  !----------------TRANSFORM TO FACE CENTER--------------------------------------
250 
251  CALL n2e3d(rho1,rho)
252 
253  IF(dbg_set(dbg_sbr)) WRITE(ipt,*) "END: DENS1"
254  RETURN
255  END SUBROUTINE dens1
256  !==============================================================================|
257 
258  !==============================================================================|
259  ! COMPUTE DENSITY USING SALINITY AND POTENTIAL TEMP |
260  ! APPEARS TO BE BASED ON THE ROMS CODE? |
261  ! CALCULATES: RHO1(M) DENSITY AT NODES |
262  ! CALCULATES: RHO (N) DENSITY AT ELEMENTS |
263  !==============================================================================|
264 
265  SUBROUTINE dens2_2d(MYS,MYT,MYRHO)
266 
267  !==============================================================================|
268  IMPLICIT NONE
269  REAL(SP), INTENT(IN), DIMENSION(:,:) :: MYS,MYT
270  REAL(SP), INTENT(INOUT), DIMENSION(:,:):: MYRHO
271 
272 
273  INTEGER :: I,K,ub1,lb1,ub2,lb2
274  !==============================================================================|
275 
276  ! SET DIMS USING MY S
277  lb1 = lbound(mys,1)
278  ub1 = ubound(mys,1)
279 
280  lb2 = lbound(mys,2)
281  ub2 = ubound(mys,2)
282 
283  ! CHECK MY T
284  if(lb1 /= lbound(myt,1)) CALL fatal_error("EQS_OF_STATE:: Dimension mismatch!")
285  if(ub1 /= ubound(myt,1)) CALL fatal_error("EQS_OF_STATE:: Dimension mismatch!")
286 
287  if(lb2 /= lbound(myt,2)) CALL fatal_error("EQS_OF_STATE:: Dimension mismatch!")
288  if(ub2 /= ubound(myt,2)) CALL fatal_error("EQS_OF_STATE:: Dimension mismatch!")
289 
290  ! CHECK MY RHO
291  if(lb1 /= lbound(myrho,1)) CALL fatal_error("EQS_OF_STATE:: Dimension mismatch!")
292  if(ub1 /= ubound(myrho,1)) CALL fatal_error("EQS_OF_STATE:: Dimension mismatch!")
293 
294  if(lb2 /= lbound(myrho,2)) CALL fatal_error("EQS_OF_STATE:: Dimension mismatch!")
295  if(ub2 /= ubound(myrho,2)) CALL fatal_error("EQS_OF_STATE:: Dimension mismatch!")
296 
297 
298  !
299  ! CALCULATE DENSITY FROM EQUATION OF STATE
300  !
301  myrho =mys*mys*mys*6.76786136e-6_sp &
302  & - mys*mys*4.8249614e-4_sp &
303  & + mys*8.14876577e-1_sp &
304  & - 0.22584586e0_sp
305 
306  myrho = myrho * &
307  & ( myt*myt*myt*1.667e-8_sp &
308  & - myt*myt*8.164e-7_sp &
309  & + myt*1.803e-5_sp &
310  & )
311 
312  myrho = myrho+1.0_sp &
313  & - myt*myt*myt*1.0843e-6_sp &
314  & + myt*myt*9.8185e-5_sp &
315  & - myt*4.786e-3_sp
316 
317  myrho = myrho* &
318  & ( mys*mys*mys*6.76786136e-6_sp &
319  & - mys*mys*4.8249614e-4_sp &
320  & + mys*8.14876577e-1_sp &
321  & +3.895414e-2_sp &
322  & )
323 
324  myrho = myrho &
325  & - (myt-3.98_sp) *(myt-3.98_sp) * (myt+283.0_sp) / (503.57_sp*(myt+67.26_sp))
326 
327  !
328  ! CALCULATE RHO
329  !
330  myrho = myrho*1.e-3_sp
331 
332  END SUBROUTINE dens2_2d
333  !==============================================================================|
334  SUBROUTINE dens2_1d(MYS,MYT,MYRHO)
335 
336  !==============================================================================|
337  IMPLICIT NONE
338  REAL(SP), INTENT(IN), DIMENSION(:) :: MYS,MYT
339  REAL(SP), INTENT(INOUT), DIMENSION(:):: MYRHO
340 
341  INTEGER :: I,K,ub1,lb1,ub2,lb2
342  !==============================================================================|
343 
344  ! SET DIMS USING MY S
345  lb1 = lbound(mys,1)
346  ub1 = ubound(mys,1)
347 
348  ! CHECK MY T
349  if(lb1 /= lbound(myt,1)) CALL fatal_error("EQS_OF_STATE:: Dimension mismatch!")
350  if(ub1 /= ubound(myt,1)) CALL fatal_error("EQS_OF_STATE:: Dimension mismatch!")
351 
352  ! CHECK MY RHO
353  if(lb1 /= lbound(myrho,1)) CALL fatal_error("EQS_OF_STATE:: Dimension mismatch!")
354  if(ub1 /= ubound(myrho,1)) CALL fatal_error("EQS_OF_STATE:: Dimension mismatch!")
355 
356  !
357  ! CALCULATE DENSITY FROM EQUATION OF STATE
358  !
359  myrho =mys*mys*mys*6.76786136e-6_sp &
360  & - mys*mys*4.8249614e-4_sp &
361  & + mys*8.14876577e-1_sp &
362  & - 0.22584586e0_sp
363 
364  myrho = myrho * &
365  & ( myt*myt*myt*1.667e-8_sp &
366  & - myt*myt*8.164e-7_sp &
367  & + myt*1.803e-5_sp &
368  & )
369 
370  myrho = myrho+1.0_sp &
371  & - myt*myt*myt*1.0843e-6_sp &
372  & + myt*myt*9.8185e-5_sp &
373  & - myt*4.786e-3_sp
374 
375  myrho = myrho* &
376  & ( mys*mys*mys*6.76786136e-6_sp &
377  & - mys*mys*4.8249614e-4_sp &
378  & + mys*8.14876577e-1_sp &
379  & +3.895414e-2_sp &
380  & )
381 
382  myrho = myrho &
383  & - (myt-3.98_sp) *(myt-3.98_sp) * (myt+283.0_sp) / (503.57_sp*(myt+67.26_sp))
384 
385  !
386  ! CALCULATE RHO
387  !
388  myrho = myrho*1.e-3_sp
389 
390  END SUBROUTINE dens2_1d
391  !==============================================================================|
392  SUBROUTINE dens2
393 
394  !==============================================================================|
395  IMPLICIT NONE
396  !==============================================================================|
397  IF(dbg_set(dbg_sbr)) WRITE(ipt,*) "START: DENS2"
398 
399  !
400  ! CALCULATE DENSITY FROM EQUATION OF STATE
401  !
402  CALL dens2_2d(s1,t1,rho1)
403  rho1(:,kb)=0.0_sp
404  rho1(0,:)=0.0_sp
405 
406 
407  !
408  ! AVERAGE FROM NODES TO FACE CENTERS
409  !
410  CALL n2e3d(rho1,rho)
411 
412  IF(dbg_set(dbg_sbr)) WRITE(ipt,*) "END: DENS2"
413  RETURN
414  END SUBROUTINE dens2
415  !==============================================================================|
416 
417  !==============================================================================|
418  ! COMPUTE IN SITU DENSITY - 1000 USING SALINITY, POTENTIAL TEMP, |
419  ! AND PRESSURE FROM A POLYNOMIAL EXPRESSION (JACKETT & MCDOUGALL, |
420  ! 1995). IT ASSUMES NO PRESSURE VARIATION ALONG GEOPOTENTIAL |
421  ! SURFACES, THAT IS, DEPTH (METERS; NEGATIVE) AND PRESSURE (DBAR |
422  ! ASSUMED NEGATIVE HERE) ARE INTERCHANGEABLE. |
423  ! |
424  ! check Values: (T=3 C, S=35.5 PSU, Z=-5000 m) |
425  ! RHOF = 1050.3639165364 (kg/m3) |
426  ! DEN1 = 1028.2845117925 (kg/m3) |
427  ! |
428  ! Reference: |
429  ! |
430  ! Jackett, D. R. and T. J. McDougall, 1995, Minimal Adjustment of |
431  ! Hydrostatic Profiles to Achieve Static Stability, J. of Atmos. |
432  ! and Oceanic Techn., vol. 12, pp. 381-389. |
433  ! |
434  ! CALCULATES: RHO1(M) DENSITY AT NODES |
435  ! CALCULATES: RHO (N) DENSITY AT ELEMENTS |
436  !==============================================================================|
437  SUBROUTINE jacket_mcdougall_2d(MYS,MYT,MYP,MYRHO)
438 
439  !==============================================================================|
440  IMPLICIT NONE
441  REAL(SP), INTENT(IN), DIMENSION(:,:) :: MYS,MYT,MYP
442  REAL(SP), INTENT(INOUT), DIMENSION(:,:):: MYRHO
443 
444  REAL(SP) :: TF,SF,sqrtSF,PBAR,TEMP(10),BULK,BULK0,BULK1,BULK2
445  INTEGER :: I,K,ub1,lb1,ub2,lb2
446  !==============================================================================|
447  ! Polynomial expansion coefficients for the computation of in situ |
448  ! density via the nonlinear equation of state for seawater as a |
449  ! function of potential temperature, salinity, and pressure (Jackett |
450  ! and McDougall, 1995). |
451  REAL(SP), PARAMETER :: A00 = +1.965933e+04_sp
452  REAL(SP), PARAMETER :: A01 = +1.444304e+02_sp
453  REAL(SP), PARAMETER :: A02 = -1.706103e+00_sp
454  REAL(SP), PARAMETER :: A03 = +9.648704e-03_sp
455  REAL(SP), PARAMETER :: A04 = -4.190253e-05_sp
456  REAL(SP), PARAMETER :: B00 = +5.284855e+01_sp
457  REAL(SP), PARAMETER :: B01 = -3.101089e-01_sp
458  REAL(SP), PARAMETER :: B02 = +6.283263e-03_sp
459  REAL(SP), PARAMETER :: B03 = -5.084188e-05_sp
460  REAL(SP), PARAMETER :: D00 = +3.886640e-01_sp
461  REAL(SP), PARAMETER :: D01 = +9.085835e-03_sp
462  REAL(SP), PARAMETER :: D02 = -4.619924e-04_sp
463  REAL(SP), PARAMETER :: E00 = +3.186519e+00_sp
464  REAL(SP), PARAMETER :: E01 = +2.212276e-02_sp
465  REAL(SP), PARAMETER :: E02 = -2.984642e-04_sp
466  REAL(SP), PARAMETER :: E03 = +1.956415e-06_sp
467  REAL(SP), PARAMETER :: F00 = +6.704388e-03_sp
468  REAL(SP), PARAMETER :: F01 = -1.847318e-04_sp
469  REAL(SP), PARAMETER :: F02 = +2.059331e-07_sp
470  REAL(SP), PARAMETER :: G00 = +1.480266e-04_sp
471  REAL(SP), PARAMETER :: G01 = +2.102898e-04_sp
472  REAL(SP), PARAMETER :: G02 = -1.202016e-05_sp
473  REAL(SP), PARAMETER :: G03 = +1.394680e-07_sp
474  REAL(SP), PARAMETER :: H00 = -2.040237e-06_sp
475  REAL(SP), PARAMETER :: H01 = +6.128773e-08_sp
476  REAL(SP), PARAMETER :: H02 = +6.207323e-10_sp
477 
478  REAL(SP), PARAMETER :: Q00 = +9.99842594e+02_sp
479  REAL(SP), PARAMETER :: Q01 = +6.793952e-02_sp
480  REAL(SP), PARAMETER :: Q02 = -9.095290e-03_sp
481  REAL(SP), PARAMETER :: Q03 = +1.001685e-04_sp
482  REAL(SP), PARAMETER :: Q04 = -1.120083e-06_sp
483  REAL(SP), PARAMETER :: Q05 = +6.536332e-09_sp
484  REAL(SP), PARAMETER :: U00 = +8.24493e-01_sp
485  REAL(SP), PARAMETER :: U01 = -4.08990e-03_sp
486  REAL(SP), PARAMETER :: U02 = +7.64380e-05_sp
487  REAL(SP), PARAMETER :: U03 = -8.24670e-07_sp
488  REAL(SP), PARAMETER :: U04 = +5.38750e-09_sp
489  REAL(SP), PARAMETER :: V00 = -5.72466e-03_sp
490  REAL(SP), PARAMETER :: V01 = +1.02270e-04_sp
491  REAL(SP), PARAMETER :: V02 = -1.65460e-06_sp
492  REAL(SP), PARAMETER :: W00 = +4.8314e-04_sp
493 
494 
495  ! SET DIMS USING MY S
496  lb1 = lbound(mys,1)
497  ub1 = ubound(mys,1)
498 
499  lb2 = lbound(mys,2)
500  ub2 = ubound(mys,2)
501 
502  ! CHECK MY T
503  if(lb1 /= lbound(myt,1)) CALL fatal_error("EQS_OF_STATE:: Dimension mismatch!")
504  if(ub1 /= ubound(myt,1)) CALL fatal_error("EQS_OF_STATE:: Dimension mismatch!")
505 
506  if(lb2 /= lbound(myt,2)) CALL fatal_error("EQS_OF_STATE:: Dimension mismatch!")
507  if(ub2 /= ubound(myt,2)) CALL fatal_error("EQS_OF_STATE:: Dimension mismatch!")
508 
509  ! CHECK MY P
510  if(lb1 /= lbound(myp,1)) CALL fatal_error("EQS_OF_STATE:: Dimension mismatch!")
511  if(ub1 /= ubound(myp,1)) CALL fatal_error("EQS_OF_STATE:: Dimension mismatch!")
512 
513  if(lb2 /= lbound(myp,2)) CALL fatal_error("EQS_OF_STATE:: Dimension mismatch!")
514  if(ub2 /= ubound(myp,2)) CALL fatal_error("EQS_OF_STATE:: Dimension mismatch!")
515 
516  ! CHECK MY RHO
517  if(lb1 /= lbound(myrho,1)) CALL fatal_error("EQS_OF_STATE:: Dimension mismatch!")
518  if(ub1 /= ubound(myrho,1)) CALL fatal_error("EQS_OF_STATE:: Dimension mismatch!")
519 
520  if(lb2 /= lbound(myrho,2)) CALL fatal_error("EQS_OF_STATE:: Dimension mismatch!")
521  if(ub2 /= ubound(myrho,2)) CALL fatal_error("EQS_OF_STATE:: Dimension mismatch!")
522 
523 
524 
525  !
526  ! CALCULATE DENSITY FROM EQUATION OF STATE
527  !
528  DO i=lb1,ub1
529  DO k=lb2,ub2
530  tf = myt(i,k)
531  sf = mys(i,k)
532  sqrtsf = sqrt(sf)
533 
534  pbar = myp(i,k)
535 
536  ! Compute density (kg/m3) at standard one atmosphere pressure
537  temp(1)=q00+tf*(q01+tf*(q02+tf*(q03+tf*(q04+tf*q05))))
538  temp(2)=u00+tf*(u01+tf*(u02+tf*(u03+tf*u04)))
539  temp(3)=v00+tf*(v01+tf*v02)
540  myrho(i,k)=temp(1)+sf*(temp(2)+sqrtsf*temp(3)+sf*w00)
541 
542  ! Compute secant bulk modulus (BULK = BULK0 + BULK1*PBAR + BULK2*PBAR*PBAR)
543  temp(4)=a00+tf*(a01+tf*(a02+tf*(a03+tf*a04)))
544  temp(5)=b00+tf*(b01+tf*(b02+tf*b03))
545  temp(6)=d00+tf*(d01+tf*d02)
546  temp(7)=e00+tf*(e01+tf*(e02+tf*e03))
547  temp(8)=f00+tf*(f01+tf*f02)
548  temp(9)=g01+tf*(g02+tf*g03)
549  temp(10)=h00+tf*(h01+tf*h02)
550 
551  bulk0=temp(4)+sf*(temp(5)+sqrtsf*temp(6))
552  bulk1=temp(7)+sf*(temp(8)+sqrtsf*g00)
553  bulk2=temp(9)+sf*temp(10)
554  bulk = bulk0 + pbar * (bulk1 + pbar * bulk2)
555 
556  ! Compute "in situ" density anomaly (kg/m3)
557  myrho(i,k)=(myrho(i,k)*bulk)/(bulk-pbar)
558  myrho(i,k)= myrho(i,k)-1000.0_sp
559  END DO
560  END DO
561 
562  !
563  ! CALCULATE RHO1
564  !
565  myrho = myrho*1.e-3_sp
566 
567 
568  END SUBROUTINE jacket_mcdougall_2d
569 !==============================================================================|
570  SUBROUTINE jacket_mcdougall_1d(MYS,MYT,MYP,MYRHO)
571 
572  !==============================================================================|
573  IMPLICIT NONE
574  REAL(SP), INTENT(IN), DIMENSION(:) :: MYS,MYT,MYP
575  REAL(SP), INTENT(INOUT), DIMENSION(:):: MYRHO
576 
577  REAL(SP) :: TF,SF,sqrtSF,PBAR,TEMP(10),BULK,BULK0,BULK1,BULK2
578  INTEGER :: K,ub1,lb1
579  !==============================================================================|
580  ! Polynomial expansion coefficients for the computation of in situ |
581  ! density via the nonlinear equation of state for seawater as a |
582  ! function of potential temperature, salinity, and pressure (Jackett |
583  ! and McDougall, 1995). |
584  REAL(SP), PARAMETER :: A00 = +1.965933e+04_sp
585  REAL(SP), PARAMETER :: A01 = +1.444304e+02_sp
586  REAL(SP), PARAMETER :: A02 = -1.706103e+00_sp
587  REAL(SP), PARAMETER :: A03 = +9.648704e-03_sp
588  REAL(SP), PARAMETER :: A04 = -4.190253e-05_sp
589  REAL(SP), PARAMETER :: B00 = +5.284855e+01_sp
590  REAL(SP), PARAMETER :: B01 = -3.101089e-01_sp
591  REAL(SP), PARAMETER :: B02 = +6.283263e-03_sp
592  REAL(SP), PARAMETER :: B03 = -5.084188e-05_sp
593  REAL(SP), PARAMETER :: D00 = +3.886640e-01_sp
594  REAL(SP), PARAMETER :: D01 = +9.085835e-03_sp
595  REAL(SP), PARAMETER :: D02 = -4.619924e-04_sp
596  REAL(SP), PARAMETER :: E00 = +3.186519e+00_sp
597  REAL(SP), PARAMETER :: E01 = +2.212276e-02_sp
598  REAL(SP), PARAMETER :: E02 = -2.984642e-04_sp
599  REAL(SP), PARAMETER :: E03 = +1.956415e-06_sp
600  REAL(SP), PARAMETER :: F00 = +6.704388e-03_sp
601  REAL(SP), PARAMETER :: F01 = -1.847318e-04_sp
602  REAL(SP), PARAMETER :: F02 = +2.059331e-07_sp
603  REAL(SP), PARAMETER :: G00 = +1.480266e-04_sp
604  REAL(SP), PARAMETER :: G01 = +2.102898e-04_sp
605  REAL(SP), PARAMETER :: G02 = -1.202016e-05_sp
606  REAL(SP), PARAMETER :: G03 = +1.394680e-07_sp
607  REAL(SP), PARAMETER :: H00 = -2.040237e-06_sp
608  REAL(SP), PARAMETER :: H01 = +6.128773e-08_sp
609  REAL(SP), PARAMETER :: H02 = +6.207323e-10_sp
610 
611  REAL(SP), PARAMETER :: Q00 = +9.99842594e+02_sp
612  REAL(SP), PARAMETER :: Q01 = +6.793952e-02_sp
613  REAL(SP), PARAMETER :: Q02 = -9.095290e-03_sp
614  REAL(SP), PARAMETER :: Q03 = +1.001685e-04_sp
615  REAL(SP), PARAMETER :: Q04 = -1.120083e-06_sp
616  REAL(SP), PARAMETER :: Q05 = +6.536332e-09_sp
617  REAL(SP), PARAMETER :: U00 = +8.24493e-01_sp
618  REAL(SP), PARAMETER :: U01 = -4.08990e-03_sp
619  REAL(SP), PARAMETER :: U02 = +7.64380e-05_sp
620  REAL(SP), PARAMETER :: U03 = -8.24670e-07_sp
621  REAL(SP), PARAMETER :: U04 = +5.38750e-09_sp
622  REAL(SP), PARAMETER :: V00 = -5.72466e-03_sp
623  REAL(SP), PARAMETER :: V01 = +1.02270e-04_sp
624  REAL(SP), PARAMETER :: V02 = -1.65460e-06_sp
625  REAL(SP), PARAMETER :: W00 = +4.8314e-04_sp
626 
627 
628  ! SET DIMS USING MY S
629  lb1 = lbound(mys,1)
630  ub1 = ubound(mys,1)
631 
632  ! CHECK MY T
633  if(lb1 /= lbound(myt,1)) CALL fatal_error("EQS_OF_STATE:: Dimension mismatch!")
634  if(ub1 /= ubound(myt,1)) CALL fatal_error("EQS_OF_STATE:: Dimension mismatch!")
635 
636  ! CHECK MY P
637  if(lb1 /= lbound(myp,1)) CALL fatal_error("EQS_OF_STATE:: Dimension mismatch!")
638  if(ub1 /= ubound(myp,1)) CALL fatal_error("EQS_OF_STATE:: Dimension mismatch!")
639 
640  ! CHECK MY RHO
641  if(lb1 /= lbound(myrho,1)) CALL fatal_error("EQS_OF_STATE:: Dimension mismatch!")
642  if(ub1 /= ubound(myrho,1)) CALL fatal_error("EQS_OF_STATE:: Dimension mismatch!")
643 
644  !
645  ! CALCULATE DENSITY FROM EQUATION OF STATE
646  !
647  DO k=lb1,ub1
648  tf = myt(k)
649  sf = mys(k)
650  sqrtsf = sqrt(sf)
651 
652  pbar = myp(k)
653 
654  ! Compute density (kg/m3) at standard one atmosphere pressure
655  temp(1)=q00+tf*(q01+tf*(q02+tf*(q03+tf*(q04+tf*q05))))
656  temp(2)=u00+tf*(u01+tf*(u02+tf*(u03+tf*u04)))
657  temp(3)=v00+tf*(v01+tf*v02)
658  myrho(k)=temp(1)+sf*(temp(2)+sqrtsf*temp(3)+sf*w00)
659 
660  ! Compute secant bulk modulus (BULK = BULK0 + BULK1*PBAR + BULK2*PBAR*PBAR)
661  temp(4)=a00+tf*(a01+tf*(a02+tf*(a03+tf*a04)))
662  temp(5)=b00+tf*(b01+tf*(b02+tf*b03))
663  temp(6)=d00+tf*(d01+tf*d02)
664  temp(7)=e00+tf*(e01+tf*(e02+tf*e03))
665  temp(8)=f00+tf*(f01+tf*f02)
666  temp(9)=g01+tf*(g02+tf*g03)
667  temp(10)=h00+tf*(h01+tf*h02)
668 
669  bulk0=temp(4)+sf*(temp(5)+sqrtsf*temp(6))
670  bulk1=temp(7)+sf*(temp(8)+sqrtsf*g00)
671  bulk2=temp(9)+sf*temp(10)
672  bulk = bulk0 + pbar * (bulk1 + pbar * bulk2)
673 
674  ! Compute "in situ" density anomaly (kg/m3)
675  myrho(k)=(myrho(k)*bulk)/(bulk-pbar)
676  myrho(k)= myrho(k)-1000.0_sp
677  END DO
678 
679  !
680  ! CALCULATE RHO1
681  !
682  myrho = myrho*1.e-3_sp
683 
684 
685  END SUBROUTINE jacket_mcdougall_1d
686 !==============================================================================|
687  SUBROUTINE dens3
688 !==============================================================================|
689  IMPLICIT NONE
690  REAL(SP), DIMENSION(0:MT,KB) :: MYP
691  INTEGER :: K
692 
693  IF(dbg_set(dbg_sbr)) WRITE(ipt,*) "START: DENS3"
694 
695  DO k=1,kbm1
696  myp(:,k) = -grav_n(:)*1.025_sp*zz(:,k)*d(:) *0.01_sp
697  END DO
698 
699  CALL jacket_mcdougall_2d(s1,t1,myp,rho1)
700  rho1(:,kb)=0.0_sp
701  rho1(0,:)=0.0_sp
702 
703  !
704  ! AVERAGE FROM NODES TO FACE CENTERS
705  !
706  CALL n2e3d(rho1,rho)
707 
708  IF(dbg_set(dbg_sbr)) WRITE(ipt,*) "END: DENS3"
709  RETURN
710  END SUBROUTINE dens3
711  !==============================================================================!
712  FUNCTION svan(S4,T4,P04,SIGMA)
713  !==============================================================================!
714  ! specific volume anomaly (steric anomaly) based on 1980 equation |
715  ! of state for seawater and 1978 practerical salinity scale. |
716  ! references: |
717  ! millero, et al (1980) deep-sea res.,27a,255-264 |
718  ! millero and poisson 1981,deep-sea res.,28a pp 625-629. |
719  ! both above references are also found in unesco report 38 (1981) |
720  ! |
721  ! units: |
722  ! pressure p04 decibars |
723  ! temperature t4 deg celsius (ipts-68) |
724  ! salinity s4 (ipss-78) |
725  ! spec. vol. ana. svan m**3/kg *1.0e-8 |
726  ! density ana. sigma kg/m**3 |
727  ! |
728  ! check value: svan=981.3021 e-8 m**3/kg. for s = 40 (ipss-78), |
729  ! t = 40 deg c, p0= 10000 decibars. |
730  ! check value: sigma = 59.82037 kg/m**3. for s = 40 (ipss-78) , |
731  ! t = 40 deg c, p0= 10000 decibars. |
732  !==============================================================================!
733 
734 
735  USE mod_prec
736  IMPLICIT NONE
737  REAL(SP) :: SVAN
738  REAL(SP), INTENT(IN) :: S4,T4,P04
739  REAL(SP), INTENT(OUT) :: SIGMA
740  REAL(SP) P4,SIG,SR,RR1,RR2,RR3,V350P,DK
741  REAL(SP) A4,B4,C4,D4,E4,AA1,BB1,AW,BW,KK,K0,KW,K35,SVA
742  REAL(SP) GAM,PK,DVAN,DR35P
743 
744  REAL(SP), PARAMETER :: R3500 = 1028.1063_sp
745  REAL(SP), PARAMETER :: RR4 = 4.8314e-4_sp
746  REAL(SP), PARAMETER :: DR350 = 28.106331_sp
747 
748 
749  ! rr4 is refered to as c in millero and poisson 1981
750  ! convert pressure to bars and take square root salinity.
751 
752  p4=p04/10.0_sp
753  sr = sqrt(abs(s4))
754 
755  ! pure water density at atmospheric pressure
756  ! bigg p.h.,(1967) br. j. applied physics 8 pp 521-537.
757  !
758 
759  rr1=((((6.536332e-9_sp*t4-1.120083e-6_sp)*t4+1.001685e-4_sp)*t4 &
760  -9.095290e-3_sp)*t4+6.793952e-2_sp)*t4-28.263737_sp
761 
762 
763  ! seawater density atm press.
764  ! coefficients involving salinity
765  ! rr2 = a in notation of millero and poisson 1981
766 
767  rr2=(((5.3875e-9_sp*t4-8.2467e-7_sp)*t4+7.6438e-5_sp)*t4-4.0899e-3_sp)*t4 &
768  +8.24493e-1_sp
769 
770  ! rr3 = b4 in notation of millero and poisson 1981
771 
772  rr3=(-1.6546e-6_sp*t4+1.0227e-4_sp)*t4-5.72466e-3_sp
773 
774  ! international one-atmosphere equation of state of seawater
775 
776  sig=(rr4*s4+rr3*sr+rr2)*s4+rr1
777 
778  ! specific volume at atmospheric pressure
779 
780  v350p = 1.0_sp/r3500
781  sva = -sig*v350p/(r3500+sig)
782  sigma=sig+dr350
783 
784  ! scale specific vol. anamoly to normally reported units
785 
786  svan=sva*1.0e+8_sp
787  IF(p4 == 0.0_sp) RETURN
788 
789  !-------------------------------------------------------------|
790  ! new high pressure equation of sate for seawater |
791  ! |
792  ! millero, el al., 1980 dsr 27a, pp 255-264 |
793  ! constant notation follows article |
794  !-------------------------------------------------------------|
795  ! compute compression terms
796 
797  e4 = (9.1697e-10*t4+2.0816e-8_sp)*t4-9.9348e-7_sp
798  bw = (5.2787e-8_sp*t4-6.12293e-6_sp)*t4+3.47718e-5_sp
799  b4 = bw + e4*s4
800 
801  d4 = 1.91075e-4
802  c4 = (-1.6078e-6_sp*t4-1.0981e-5_sp)*t4+2.2838e-3_sp
803  aw = ((-5.77905e-7_sp*t4+1.16092e-4_sp)*t4+1.43713e-3_sp)*t4 &
804  -0.1194975_sp
805  a4 = (d4*sr + c4)*s4 + aw
806 
807  bb1 = (-5.3009e-4_sp*t4+1.6483e-2_sp)*t4+7.944e-2_sp
808  aa1 = ((-6.1670e-5_sp*t4+1.09987e-2_sp)*t4-0.603459_sp)*t4+54.6746
809  kw = (((-5.155288e-5_sp*t4+1.360477e-2_sp)*t4-2.327105_sp)*t4 &
810  +148.4206_sp)*t4-1930.06_sp
811  k0 = (bb1*sr + aa1)*s4 + kw
812 
813  ! evaluate pressure polynomial
814  !-----------------------------------------------------|
815  ! k equals the secant bulk modulus of seawater |
816  ! dk=k(s,t,p)-k(35,0,p) |
817  ! k35=k(35,0,p) |
818  !-----------------------------------------------------|
819 
820  dk = (b4*p4 + a4)*p4 + k0
821  k35 = (5.03217e-5_sp*p4+3.359406_sp)*p4+21582.27_sp
822  gam=p4/k35
823  pk = 1.0_sp - gam
824  sva = sva*pk + (v350p+sva)*p4*dk/(k35*(k35+dk))
825 
826  ! scale specific vol. anamoly to normally reported units
827 
828  svan=sva*1.0e+8_sp
829  v350p = v350p*pk
830 
831  !----------------------------------------------------------|
832  ! compute density anamoly with respect to 1000.0 kg/m**3 |
833  ! 1) dr350: density anamoly at 35 (ipss-78), |
834  ! 0 deg. c and 0 decibars |
835  ! 2) dr35p: density anamoly at 35 (ipss-78), |
836  ! 0 deg. c, pres. variation |
837  ! 3) dvan : density anamoly variations involving specific |
838  ! volume anamoly |
839  ! |
840  ! check values: sigma = 59.82037 kg/m**3 |
841  ! for s = 40 (ipss-78), t = 40 deg c, p0= 10000 decibars. |
842  !----------------------------------------------------------|
843 
844  dr35p=gam/v350p
845  dvan=sva/(v350p*(v350p+sva))
846  sigma=dr350+dr35p-dvan
847 
848  RETURN
849  END FUNCTION svan
850  !==============================================================================!
851  !==============================================================================!
852 
853  !==============================================================================|
854  FUNCTION theta(S4,T04,P04,PR)
855  !==============================================================================|
856  ! to compute local potential temperature at pr using |
857  ! bryden 1973 polynomial for adiabatic lapse rate and |
858  ! runge-kutta 4th order integration algorithm. |
859  ! ref: bryden,h.,1973,deep-sea res.,20,401-408; |
860  ! fofonoff,n.,1977,deep-sea res.,24,489-491 |
861  ! |
862  ! units: |
863  ! pressure p04 decibars |
864  ! temperature t04 deg celsius (ipts-68) |
865  ! salinity s4 (ipss-78) |
866  ! reference prs pr decibars |
867  ! potential tmp. theta deg celsius |
868  ! checkvalue: |
869  ! theta= 36.89073 c,s=40 (ipss-78), |
870  ! t0=40 deg c,p0=10000 decibars,pr=0 decibars |
871  ! |
872  ! |
873  ! set up intermediate temperature and pressure variables. |
874  !==============================================================================|
875 
876  USE mod_prec
877  IMPLICIT NONE
878  REAL(SP) :: THETA
879  REAL(SP), INTENT(IN) :: S4,T04,P04,PR
880  REAL(SP) :: P4,T4,H4,Q4,XK
881 
882 
883  !==============================================================================|
884 
885  p4 = p04
886  t4 = t04
887  h4 = pr - p4
888  xk = h4*atg(s4,t4,p4)
889  t4 = t4 + .5_sp*xk
890  q4 = xk
891  p4 = p4 + 0.5_sp*h4
892  xk = h4*atg(s4,t4,p4)
893  t4 = t4 + 0.29289322_sp*(xk-q4)
894  q4 = 0.58578644_sp*xk + 0.121320344_sp*q4
895  xk = h4*atg(s4,t4,p4)
896  t4 = t4 + 1.707106781_sp*(xk-q4)
897  q4 = 3.414213562_sp*xk - 4.121320344_sp*q4
898  p4 = p4 + 0.5_sp*h4
899  xk = h4*atg(s4,t4,p4)
900  theta = t4 + (xk-2.0_sp*q4)/6.0_sp
901 
902  RETURN
903  END FUNCTION theta
904  !==============================================================================|
905 
906  !==============================================================================|
907  ! adiabatic temperature gradient deg c per decibar |
908  ! ref: bryden, h., 1973,deep-sea res.,20,401-408 |
909  ! |
910  ! units: |
911  ! pressure P4 decibars |
912  ! temperature T4 deg celsius(ipts-68) |
913  ! salinity s4 (ipss-78) |
914  ! adiabatic atg deg. c/decibar |
915  ! checkvalue: atg=3.255976e-4 c/dbar for s=40 (ipss-78), |
916  ! t=40 deg c,p0=10000 decibars |
917  !==============================================================================|
918 
919  REAL(SP) FUNCTION ATG(S4,T4,P4)
920  USE mod_prec
921 
922  !------------------------------------------------------------------------------|
923 
924  IMPLICIT NONE
925  REAL(SP), INTENT(IN) :: S4,T4,P4
926  REAL(SP) :: DS
927 
928  !==============================================================================|
929 
930  ds = s4 - 35.0_sp
931  atg = (((-2.1687e-16_sp*t4+1.8676e-14_sp)*t4-4.6206e-13_sp)*p4 &
932  +((2.7759e-12_sp*t4-1.1351e-10_sp)*ds+((-5.4481e-14_sp*t4 &
933  +8.733e-12_sp)*t4-6.7795e-10_sp)*t4+1.8741e-8_sp))*p4 &
934  +(-4.2393e-8_sp*t4+1.8932e-6_sp)*ds &
935  +((6.6228e-10_sp*t4-6.836e-8_sp)*t4+8.5258e-6_sp)*t4+3.5803e-5_sp
936 
937  RETURN
938  END FUNCTION atg
939  !==============================================================================|
940 
941 END MODULE eqs_of_state
eqs_of_state::jacket_mcdougall
Definition: eqs_of_state.f90:104
all_vars::d
real(sp), dimension(:), allocatable, target d
Definition: mod_main.f90:1132
mod_utils
Definition: mod_utils.f90:40
mod_utils::dbg_set
logical function dbg_set(vrb)
Definition: mod_utils.f90:182
all_vars::rho1
real(sp), dimension(:,:), allocatable, target rho1
Definition: mod_main.f90:1309
all_vars::t1
real(sp), dimension(:,:), allocatable, target t1
Definition: mod_main.f90:1307
all_vars::rho
real(sp), dimension(:,:), allocatable, target rho
Definition: mod_main.f90:1284
all_vars::s1
real(sp), dimension(:,:), allocatable, target s1
Definition: mod_main.f90:1308
mod_prec
Definition: mod_prec.f90:43
eqs_of_state::dens2
subroutine dens2
Definition: eqs_of_state.f90:393
eqs_of_state::fofonoff_millard
Definition: eqs_of_state.f90:94
all_vars
Definition: mod_main.f90:959
eqs_of_state::dens3
subroutine dens3
Definition: eqs_of_state.f90:688
eqs_of_state
Definition: eqs_of_state.f90:40
mod_utils::fatal_error
subroutine fatal_error(ER1, ER2, ER3, ER4)
Definition: mod_utils.f90:230
eqs_of_state::dens1
subroutine dens1
Definition: eqs_of_state.f90:229
all_vars::grav_n
real(sp), dimension(:), allocatable, target grav_n
Definition: mod_main.f90:1013
all_vars::n2e3d
subroutine n2e3d(NVAR, EVAR)
Definition: mod_main.f90:1370
eqs_of_state::dens2g
Definition: eqs_of_state.f90:99
mod_utils::dbg_sbr
integer, parameter dbg_sbr
Definition: mod_utils.f90:69
all_vars::zz
real(sp), dimension(:,:), allocatable, target zz
Definition: mod_main.f90:1091
Generated by   doxygen 1.8.15