/[MITgcm]/MITgcm/model/src/find_rho.F
ViewVC logotype

Annotation of /MITgcm/model/src/find_rho.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph


Revision 1.23 - (hide annotations) (download)
Fri Nov 15 03:01:21 2002 UTC (21 years, 6 months ago) by heimbach
Branch: MAIN
CVS Tags: checkpoint47e_post, checkpoint47c_post, checkpoint48e_post, checkpoint48b_post, checkpoint48c_pre, checkpoint47d_pre, checkpoint47a_post, checkpoint48d_pre, checkpoint47i_post, checkpoint47d_post, checkpoint48d_post, checkpoint48f_post, checkpoint47g_post, checkpoint48a_post, checkpoint47j_post, branch-exfmods-tag, checkpoint48c_post, checkpoint47b_post, checkpoint47f_post, checkpoint47, checkpoint48, checkpoint47h_post
Branch point for: branch-exfmods-curt
Changes since 1.22: +14 -2 lines
differentiable version of checkpoint46n_post
o external_fields_load now part of differentiation list
o pressure needs multiple storing;
  would be nice to have store_pressure at beginning or
  end of forward_step, e.g. by having phiHyd global (5-dim.)
  (NB: pressure is needed for certain cases in find_rho,
  which is also invoked through convective_adjustment).
o recomputations in find_rho for cases
 'JMD95'/'UNESCO' or 'MDJWF' are OK.
o #define ATMOSPHERIC_LOADING should be differentiable
o ini_forcing shifted to begining of initialise_varia

1 heimbach 1.23 C $Header: /u/gcmpack/MITgcm/model/src/find_rho.F,v 1.22 2002/11/01 22:00:33 mlosch Exp $
2 cnh 1.15 C $Name: $
3 cnh 1.1
4 cnh 1.9 #include "CPP_OPTIONS.h"
5 adcroft 1.5 #define USE_FACTORIZED_POLY
6 cnh 1.1
7 cnh 1.15 CBOP
8     C !ROUTINE: FIND_RHO
9     C !INTERFACE:
10 adcroft 1.4 subroutine FIND_RHO(
11 mlosch 1.20 I bi, bj, iMin, iMax, jMin, jMax, k, kRef,
12 adcroft 1.13 I tFld, sFld,
13 adcroft 1.4 O rholoc,
14     I myThid )
15 cnh 1.15
16     C !DESCRIPTION: \bv
17     C *==========================================================*
18     C | o SUBROUTINE FIND_RHO
19 mlosch 1.21 C | Calculates [rho(S,T,z)-rhoConst] of a slice
20 cnh 1.15 C *==========================================================*
21     C |
22     C | k - is the Theta/Salt level
23     C | kRef - determines pressure reference level
24     C | (not used in 'LINEAR' mode)
25     C |
26     C *==========================================================*
27     C \ev
28    
29     C !USES:
30 cnh 1.1 implicit none
31 heimbach 1.12 C == Global variables ==
32 cnh 1.1 #include "SIZE.h"
33 cnh 1.7 #include "EEPARAMS.h"
34 cnh 1.1 #include "PARAMS.h"
35 mlosch 1.16 #include "EOS.h"
36     #include "GRID.h"
37 heimbach 1.12
38 cnh 1.15 C !INPUT/OUTPUT PARAMETERS:
39 heimbach 1.12 C == Routine arguments ==
40 adcroft 1.4 integer bi,bj,iMin,iMax,jMin,jMax
41     integer k ! Level of Theta/Salt slice
42     integer kRef ! Pressure reference level
43 adcroft 1.13 _RL tFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
44     _RL sFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
45 adcroft 1.4 _RL rholoc(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
46     integer myThid
47 heimbach 1.12
48 cnh 1.15 C !LOCAL VARIABLES:
49 heimbach 1.12 C == Local variables ==
50 adcroft 1.4 integer i,j
51 mlosch 1.21 _RL refTemp,refSalt,sigRef,tP,sP,deltaSig,dRho
52 mlosch 1.16 _RL rhoP0(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
53     _RL bulkMod(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
54 mlosch 1.19 _RL rhoNum(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
55     _RL rhoDen(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
56 adcroft 1.10 character*(max_len_mbuf) msgbuf
57 cnh 1.15 CEOP
58 heimbach 1.11
59     #ifdef ALLOW_AUTODIFF_TAMC
60     DO j=1-OLy,sNy+OLy
61     DO i=1-OLx,sNx+OLx
62 mlosch 1.16 rholoc(i,j) = 0. _d 0
63     rhoP0(i,j) = 0. _d 0
64     bulkMod(i,j) = 0. _d 0
65 heimbach 1.11 ENDDO
66     ENDDO
67     #endif
68 cnh 1.1
69 mlosch 1.22 #ifdef CHECK_SALINITY_FOR_NEGATIVE_VALUES
70     CALL LOOK_FOR_NEG_SALINITY( bi, bj, iMin, iMax, jMin, jMax, k,
71     & sFld, myThid )
72     #endif
73    
74 mlosch 1.19 if (equationOfState.eq.'LINEAR') then
75 adcroft 1.4
76     C ***NOTE***
77     C In the linear EOS, to make the static stability calculation meaningful
78     C we alway calculate the perturbation with respect to the surface level.
79     C **********
80 adcroft 1.6 refTemp=tRef(kRef)
81     refSalt=sRef(kRef)
82 adcroft 1.4
83 mlosch 1.21 dRho = rhoNil-rhoConst
84    
85 cnh 1.1 do j=jMin,jMax
86     do i=iMin,iMax
87 mlosch 1.21 rholoc(i,j)=rhoNil*(
88 adcroft 1.13 & sBeta*(sFld(i,j,k,bi,bj)-refSalt)
89     & -tAlpha*(tFld(i,j,k,bi,bj)-refTemp) )
90 mlosch 1.21 & + dRho
91 cnh 1.1 enddo
92     enddo
93 cnh 1.8
94 mlosch 1.19 elseif (equationOfState.eq.'POLY3') then
95 adcroft 1.4
96     refTemp=eosRefT(kRef)
97     refSalt=eosRefS(kRef)
98 mlosch 1.21 sigRef=eosSig0(kRef) + (1000.-rhoConst)
99 adcroft 1.4
100     do j=jMin,jMax
101     do i=iMin,iMax
102 adcroft 1.13 tP=tFld(i,j,k,bi,bj)-refTemp
103     sP=sFld(i,j,k,bi,bj)-refSalt
104 adcroft 1.5 #ifdef USE_FACTORIZED_POLY
105 adcroft 1.4 deltaSig=
106 adcroft 1.5 & (( eosC(9,kRef)*sP + eosC(5,kRef) )*sP + eosC(2,kRef) )*sP
107     & + ( ( eosC(6,kRef)
108     & *tP
109     & +eosC(7,kRef)*sP + eosC(3,kRef)
110     & )*tP
111     & +(eosC(8,kRef)*sP + eosC(4,kRef) )*sP + eosC(1,kRef)
112     & )*tP
113     #else
114     deltaSig=
115     & eosC(1,kRef)*tP
116     & +eosC(2,kRef) *sP
117     & +eosC(3,kRef)*tP*tP
118     & +eosC(4,kRef)*tP *sP
119     & +eosC(5,kRef) *sP*sP
120     & +eosC(6,kRef)*tP*tP*tP
121     & +eosC(7,kRef)*tP*tP *sP
122     & +eosC(8,kRef)*tP *sP*sP
123     & +eosC(9,kRef) *sP*sP*sP
124     #endif
125     rholoc(i,j)=sigRef+deltaSig
126 adcroft 1.4 enddo
127     enddo
128    
129 mlosch 1.19 elseif ( equationOfState(1:5).eq.'JMD95'
130     & .or. equationOfState.eq.'UNESCO' ) then
131 mlosch 1.16 C nonlinear equation of state in pressure coordinates
132    
133     CALL FIND_RHOP0(
134     I bi, bj, iMin, iMax, jMin, jMax, k,
135     I tFld, sFld,
136     O rhoP0,
137     I myThid )
138    
139     CALL FIND_BULKMOD(
140 adcroft 1.18 I bi, bj, iMin, iMax, jMin, jMax, k, kRef,
141 mlosch 1.16 I tFld, sFld,
142     O bulkMod,
143     I myThid )
144    
145 heimbach 1.23 #ifdef ALLOW_AUTODIFF_TAMC
146     cph can't do storing here since find_rho is called multiple times;
147     cph additional recomp. should be acceptable
148     cphCADJ STORE rhoP0(:,:) = comlev1_bibj_k , key=kkey , byte=isbyte
149     cphCADJ STORE bulkMod(:,:) = comlev1_bibj_k , key=kkey , byte=isbyte
150     #endif /* ALLOW_AUTODIFF_TAMC */
151 mlosch 1.16 do j=jMin,jMax
152     do i=iMin,iMax
153    
154     C density of sea water at pressure p
155     rholoc(i,j) = rhoP0(i,j)
156     & /(1. _d 0 -
157 adcroft 1.18 & pressure(i,j,kRef,bi,bj)*SItoBar/bulkMod(i,j))
158 mlosch 1.21 & - rhoConst
159 mlosch 1.16
160     end do
161     end do
162    
163 mlosch 1.19 elseif ( equationOfState.eq.'MDJWF' ) then
164    
165     CALL FIND_RHONUM( bi, bj, iMin, iMax, jMin, jMax, k, kRef,
166     & tFld, sFld, rhoNum, myThid )
167    
168     CALL FIND_RHODEN( bi, bj, iMin, iMax, jMin, jMax, k, kRef,
169     & tFld, sFld, rhoDen, myThid )
170 heimbach 1.23 #ifdef ALLOW_AUTODIFF_TAMC
171     cph can't do storing here since find_rho is called multiple times;
172     cph additional recomp. should be acceptable
173     cphCADJ STORE rhoNum(:,:) = comlev1_bibj_k , key=kkey , byte=isbyte
174     cphCADJ STORE rhoDen(:,:) = comlev1_bibj_k , key=kkey , byte=isbyte
175     #endif /* ALLOW_AUTODIFF_TAMC */
176 mlosch 1.19 do j=jMin,jMax
177     do i=iMin,iMax
178    
179 mlosch 1.21 rholoc(i,j) = rhoNum(i,j)*rhoDen(i,j) - rhoConst
180 mlosch 1.19
181     end do
182     end do
183 heimbach 1.23
184 mlosch 1.19 elseif( equationOfState .eq. 'IDEALG' ) then
185     C
186 adcroft 1.4 else
187 mlosch 1.19 write(msgbuf,'(3a)')
188     & ' FIND_RHO: equationOfState = "',equationOfState,'"'
189 adcroft 1.10 call print_error( msgbuf, mythid )
190     stop 'ABNORMAL END: S/R FIND_RHO'
191 adcroft 1.4 endif
192 cnh 1.1
193     return
194     end
195 mlosch 1.16
196     CBOP
197     C !ROUTINE: FIND_RHOP0
198     C !INTERFACE:
199     subroutine FIND_RHOP0(
200     I bi, bj, iMin, iMax, jMin, jMax, k,
201     I tFld, sFld,
202     O rhoP0,
203     I myThid )
204    
205     C !DESCRIPTION: \bv
206     C *==========================================================*
207     C | o SUBROUTINE FIND_RHOP0
208     C | Calculates rho(S,T,0) of a slice
209     C *==========================================================*
210     C |
211     C | k - is the surface level
212     C |
213     C *==========================================================*
214     C \ev
215    
216     C !USES:
217     implicit none
218     C == Global variables ==
219     #include "SIZE.h"
220     #include "EEPARAMS.h"
221     #include "PARAMS.h"
222     #include "EOS.h"
223    
224     C !INPUT/OUTPUT PARAMETERS:
225     C == Routine arguments ==
226     integer bi,bj,iMin,iMax,jMin,jMax
227 mlosch 1.19 integer k ! Level of surface slice
228 mlosch 1.16 _RL tFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
229     _RL sFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
230     _RL rhoP0(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
231     _RL rfresh, rsalt
232     integer myThid
233    
234     C !LOCAL VARIABLES:
235     C == Local variables ==
236     integer i,j
237     _RL t, t2, t3, t4, s, s3o2
238     CEOP
239    
240     do j=jMin,jMax
241     do i=iMin,iMax
242     C abbreviations
243     t = tFld(i,j,k,bi,bj)
244     t2 = t*t
245     t3 = t2*t
246     t4 = t3*t
247    
248 adcroft 1.17 s = sFld(i,j,k,bi,bj)
249 mlosch 1.16 s3o2 = s*sqrt(s)
250    
251     C density of freshwater at the surface
252     rfresh =
253     & eosJMDCFw(1)
254     & + eosJMDCFw(2)*t
255     & + eosJMDCFw(3)*t2
256     & + eosJMDCFw(4)*t3
257     & + eosJMDCFw(5)*t4
258     & + eosJMDCFw(6)*t4*t
259     C density of sea water at the surface
260     rsalt =
261     & s*(
262     & eosJMDCSw(1)
263     & + eosJMDCSw(2)*t
264     & + eosJMDCSw(3)*t2
265     & + eosJMDCSw(4)*t3
266     & + eosJMDCSw(5)*t4
267     & )
268     & + s3o2*(
269     & eosJMDCSw(6)
270     & + eosJMDCSw(7)*t
271     & + eosJMDCSw(8)*t2
272     & )
273     & + eosJMDCSw(9)*s*s
274    
275     rhoP0(i,j) = rfresh + rsalt
276    
277     end do
278     end do
279    
280     return
281     end
282    
283     C !ROUTINE: FIND_BULKMOD
284     C !INTERFACE:
285     subroutine FIND_BULKMOD(
286 adcroft 1.18 I bi, bj, iMin, iMax, jMin, jMax, k, kRef,
287 mlosch 1.16 I tFld, sFld,
288     O bulkMod,
289     I myThid )
290    
291     C !DESCRIPTION: \bv
292     C *==========================================================*
293     C | o SUBROUTINE FIND_BULKMOD
294     C | Calculates the secant bulk modulus K(S,T,p) of a slice
295     C *==========================================================*
296     C |
297 adcroft 1.18 C | k - is the surface level
298     C | kRef - is the level to use to determine the pressure
299 mlosch 1.16 C |
300     C *==========================================================*
301     C \ev
302    
303     C !USES:
304     implicit none
305     C == Global variables ==
306     #include "SIZE.h"
307     #include "EEPARAMS.h"
308     #include "PARAMS.h"
309     #include "EOS.h"
310    
311     C !INPUT/OUTPUT PARAMETERS:
312     C == Routine arguments ==
313     integer bi,bj,iMin,iMax,jMin,jMax
314     integer k ! Level of surface slice
315 adcroft 1.18 integer kRef ! Reference pressure level
316 mlosch 1.16 _RL tFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
317     _RL sFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
318     _RL bulkMod(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
319     _RL bMfresh, bMsalt, bMpres
320     integer myThid
321    
322     C !LOCAL VARIABLES:
323     C == Local variables ==
324     integer i,j
325     _RL t, t2, t3, t4, s, s3o2, p, p2
326     CEOP
327    
328     do j=jMin,jMax
329     do i=iMin,iMax
330     C abbreviations
331     t = tFld(i,j,k,bi,bj)
332     t2 = t*t
333     t3 = t2*t
334     t4 = t3*t
335    
336 adcroft 1.17 s = sFld(i,j,k,bi,bj)
337 mlosch 1.16 s3o2 = s*sqrt(s)
338     C
339 adcroft 1.18 p = pressure(i,j,kRef,bi,bj)*SItoBar
340 mlosch 1.16 p2 = p*p
341     C secant bulk modulus of fresh water at the surface
342     bMfresh =
343     & eosJMDCKFw(1)
344     & + eosJMDCKFw(2)*t
345     & + eosJMDCKFw(3)*t2
346     & + eosJMDCKFw(4)*t3
347     & + eosJMDCKFw(5)*t4
348     C secant bulk modulus of sea water at the surface
349     bMsalt =
350     & s*( eosJMDCKSw(1)
351     & + eosJMDCKSw(2)*t
352     & + eosJMDCKSw(3)*t2
353     & + eosJMDCKSw(4)*t3
354     & )
355     & + s3o2*( eosJMDCKSw(5)
356     & + eosJMDCKSw(6)*t
357     & + eosJMDCKSw(7)*t2
358     & )
359     C secant bulk modulus of sea water at pressure p
360     bMpres =
361     & p*( eosJMDCKP(1)
362     & + eosJMDCKP(2)*t
363     & + eosJMDCKP(3)*t2
364     & + eosJMDCKP(4)*t3
365     & )
366     & + p*s*( eosJMDCKP(5)
367     & + eosJMDCKP(6)*t
368     & + eosJMDCKP(7)*t2
369     & )
370     & + p*s3o2*eosJMDCKP(8)
371     & + p2*( eosJMDCKP(9)
372     & + eosJMDCKP(10)*t
373     & + eosJMDCKP(11)*t2
374     & )
375     & + p2*s*( eosJMDCKP(12)
376     & + eosJMDCKP(13)*t
377     & + eosJMDCKP(14)*t2
378     & )
379    
380     bulkMod(i,j) = bMfresh + bMsalt + bMpres
381    
382     end do
383     end do
384    
385     return
386     end
387    
388 mlosch 1.19 CBOP
389     C !ROUTINE: FIND_RHONUM
390     C !INTERFACE:
391     subroutine FIND_RHONUM(
392     I bi, bj, iMin, iMax, jMin, jMax, k, kRef,
393     I tFld, sFld,
394     O rhoNum,
395     I myThid )
396    
397     C !DESCRIPTION: \bv
398     C *==========================================================*
399     C | o SUBROUTINE FIND_RHONUM
400     C | Calculates the numerator of the McDougall et al.
401     C | equation of state
402     C | - the code is more or less a copy of MOM4
403     C *==========================================================*
404     C |
405     C | k - is the surface level
406     C | kRef - is the level to use to determine the pressure
407     C |
408     C *==========================================================*
409     C \ev
410    
411     C !USES:
412     implicit none
413     C == Global variables ==
414     #include "SIZE.h"
415     #include "EEPARAMS.h"
416     #include "PARAMS.h"
417     #include "EOS.h"
418    
419     C !INPUT/OUTPUT PARAMETERS:
420     C == Routine arguments ==
421     integer bi,bj,iMin,iMax,jMin,jMax
422     integer k ! Level of surface slice
423     integer kRef ! Reference pressure level
424     _RL tFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
425     _RL sFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
426     _RL rhoNum(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
427     integer myThid
428    
429     C !LOCAL VARIABLES:
430     C == Local variables ==
431     integer i,j
432     _RL t1, t2, s1, p1
433     CEOP
434     do j=jMin,jMax
435     do i=iMin,iMax
436     C abbreviations
437     t1 = tFld(i,j,k,bi,bj)
438     t2 = t1*t1
439     s1 = sFld(i,j,k,bi,bj)
440    
441     p1 = pressure(i,j,kRef,bi,bj)*SItodBar
442    
443     rhoNum(i,j) = eosMDJWFnum(0)
444     & + t1*(eosMDJWFnum(1)
445     & + t1*(eosMDJWFnum(2) + eosMDJWFnum(3)*t1) )
446     & + s1*(eosMDJWFnum(4)
447     & + eosMDJWFnum(5)*t1 + eosMDJWFnum(6)*s1)
448     & + p1*(eosMDJWFnum(7) + eosMDJWFnum(8)*t2
449     & + eosMDJWFnum(9)*s1
450     & + p1*(eosMDJWFnum(10) + eosMDJWFnum(11)*t2) )
451    
452     end do
453     end do
454    
455     return
456     end
457    
458     CBOP
459     C !ROUTINE: FIND_RHODEN
460     C !INTERFACE:
461     subroutine FIND_RHODEN(
462     I bi, bj, iMin, iMax, jMin, jMax, k, kRef,
463     I tFld, sFld,
464     O rhoDen,
465     I myThid )
466    
467     C !DESCRIPTION: \bv
468     C *==========================================================*
469     C | o SUBROUTINE FIND_RHODEN
470     C | Calculates the denominator of the McDougall et al.
471     C | equation of state
472     C | - the code is more or less a copy of MOM4
473     C *==========================================================*
474     C |
475     C | k - is the surface level
476     C | kRef - is the level to use to determine the pressure
477     C |
478     C *==========================================================*
479     C \ev
480    
481     C !USES:
482     implicit none
483     C == Global variables ==
484     #include "SIZE.h"
485     #include "EEPARAMS.h"
486     #include "PARAMS.h"
487     #include "EOS.h"
488    
489     C !INPUT/OUTPUT PARAMETERS:
490     C == Routine arguments ==
491     integer bi,bj,iMin,iMax,jMin,jMax
492     integer k ! Level of surface slice
493     integer kRef ! Reference pressure level
494     _RL tFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
495     _RL sFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
496     _RL rhoDen(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
497     integer myThid
498    
499     C !LOCAL VARIABLES:
500     C == Local variables ==
501     integer i,j
502     _RL t1, t2, s1, sp5, p1, p1t1
503     _RL den, epsln
504     parameter ( epsln = 0. _d 0 )
505     CEOP
506     do j=jMin,jMax
507     do i=iMin,iMax
508     C abbreviations
509     t1 = tFld(i,j,k,bi,bj)
510     t2 = t1*t1
511     s1 = sFld(i,j,k,bi,bj)
512     sp5 = sqrt(s1)
513    
514     p1 = pressure(i,j,kRef,bi,bj)*SItodBar
515     p1t1 = p1*t1
516    
517     den = eosMDJWFden(0)
518     & + t1*(eosMDJWFden(1)
519     & + t1*(eosMDJWFden(2)
520     & + t1*(eosMDJWFden(3) + t1*eosMDJWFden(4) ) ) )
521     & + s1*(eosMDJWFden(5)
522     & + t1*(eosMDJWFden(6)
523     & + eosMDJWFden(7)*t2)
524     & + sp5*(eosMDJWFden(8) + eosMDJWFden(9)*t2) )
525     & + p1*(eosMDJWFden(10)
526     & + p1t1*(eosMDJWFden(11)*t2 + eosMDJWFden(12)*p1) )
527    
528     rhoDen(i,j) = 1.0/(epsln+den)
529    
530     end do
531     end do
532    
533     return
534     end
535 mlosch 1.20
536     subroutine find_rho_scalar(
537     I tLoc, sLoc, pLoc,
538     O rhoLoc,
539     I myThid )
540    
541     C !DESCRIPTION: \bv
542     C *==========================================================*
543     C | o SUBROUTINE FIND_RHO_SCALAR
544 mlosch 1.21 C | Calculates [rho(S,T,p)-rhoConst]
545 mlosch 1.20 C *==========================================================*
546     C \ev
547    
548     C !USES:
549     implicit none
550     C == Global variables ==
551     #include "SIZE.h"
552     #include "EEPARAMS.h"
553     #include "PARAMS.h"
554     #include "EOS.h"
555    
556     C !INPUT/OUTPUT PARAMETERS:
557     C == Routine arguments ==
558     _RL sLoc, tLoc, pLoc
559     _RL rhoLoc
560     integer myThid
561    
562     C !LOCAL VARIABLES:
563     C == Local variables ==
564    
565     _RL t1, t2, t3, t4, s1, s3o2, p1, p2, sp5, p1t1
566     _RL rfresh, rsalt, rhoP0
567     _RL bMfresh, bMsalt, bMpres, BulkMod
568     _RL rhoNum, rhoDen, den, epsln
569     parameter ( epsln = 0. _d 0 )
570    
571     character*(max_len_mbuf) msgbuf
572     CEOP
573    
574     rhoLoc = 0. _d 0
575     rhoP0 = 0. _d 0
576     bulkMod = 0. _d 0
577     rfresh = 0. _d 0
578     rsalt = 0. _d 0
579     bMfresh = 0. _d 0
580     bMsalt = 0. _d 0
581     bMpres = 0. _d 0
582     rhoNum = 0. _d 0
583     rhoDen = 0. _d 0
584     den = 0. _d 0
585    
586     t1 = tLoc
587     t2 = t1*t1
588     t3 = t2*t1
589     t4 = t3*t1
590    
591     s1 = sLoc
592     if ( s1 .lt. 0. _d 0 ) then
593     C issue a warning
594     write(*,'(a,i3,a,i3,a,i3,a,e13.5)')
595     & ' FIND_RHO_SCALAR: WARNING, salinity = ', s1
596     s1 = 0. _d 0
597     end if
598    
599     if (equationOfState.eq.'LINEAR') then
600    
601     rhoLoc = 0. _d 0
602    
603     elseif (equationOfState.eq.'POLY3') then
604    
605     C this is not correct, there is a field eosSig0 which should be use here
606     C but I do not intent to include the reference level in this routine
607     write(*,'(a)')
608     & ' FIND_RHO_SCALAR: for POLY3, the density is not'
609     write(*,'(a)')
610     & ' computed correctly in this routine'
611     rhoLoc = 0. _d 0
612    
613     elseif ( equationOfState(1:5).eq.'JMD95'
614     & .or. equationOfState.eq.'UNESCO' ) then
615     C nonlinear equation of state in pressure coordinates
616    
617     s3o2 = s1*sqrt(s1)
618    
619     p1 = pLoc*SItoBar
620     p2 = p1*p1
621    
622     C density of freshwater at the surface
623     rfresh =
624     & eosJMDCFw(1)
625     & + eosJMDCFw(2)*t1
626     & + eosJMDCFw(3)*t2
627     & + eosJMDCFw(4)*t3
628     & + eosJMDCFw(5)*t4
629     & + eosJMDCFw(6)*t4*t1
630     C density of sea water at the surface
631     rsalt =
632     & s1*(
633     & eosJMDCSw(1)
634     & + eosJMDCSw(2)*t1
635     & + eosJMDCSw(3)*t2
636     & + eosJMDCSw(4)*t3
637     & + eosJMDCSw(5)*t4
638     & )
639     & + s3o2*(
640     & eosJMDCSw(6)
641     & + eosJMDCSw(7)*t1
642     & + eosJMDCSw(8)*t2
643     & )
644     & + eosJMDCSw(9)*s1*s1
645    
646     rhoP0 = rfresh + rsalt
647    
648     C secant bulk modulus of fresh water at the surface
649     bMfresh =
650     & eosJMDCKFw(1)
651     & + eosJMDCKFw(2)*t1
652     & + eosJMDCKFw(3)*t2
653     & + eosJMDCKFw(4)*t3
654     & + eosJMDCKFw(5)*t4
655     C secant bulk modulus of sea water at the surface
656     bMsalt =
657     & s1*( eosJMDCKSw(1)
658     & + eosJMDCKSw(2)*t1
659     & + eosJMDCKSw(3)*t2
660     & + eosJMDCKSw(4)*t3
661     & )
662     & + s3o2*( eosJMDCKSw(5)
663     & + eosJMDCKSw(6)*t1
664     & + eosJMDCKSw(7)*t2
665     & )
666     C secant bulk modulus of sea water at pressure p
667     bMpres =
668     & p1*( eosJMDCKP(1)
669     & + eosJMDCKP(2)*t1
670     & + eosJMDCKP(3)*t2
671     & + eosJMDCKP(4)*t3
672     & )
673     & + p1*s1*( eosJMDCKP(5)
674     & + eosJMDCKP(6)*t1
675     & + eosJMDCKP(7)*t2
676     & )
677     & + p1*s3o2*eosJMDCKP(8)
678     & + p2*( eosJMDCKP(9)
679     & + eosJMDCKP(10)*t1
680     & + eosJMDCKP(11)*t2
681     & )
682     & + p2*s1*( eosJMDCKP(12)
683     & + eosJMDCKP(13)*t1
684     & + eosJMDCKP(14)*t2
685     & )
686    
687     bulkMod = bMfresh + bMsalt + bMpres
688    
689     C density of sea water at pressure p
690 mlosch 1.21 rhoLoc = rhoP0/(1. _d 0 - p1/bulkMod) - rhoConst
691 mlosch 1.20
692     elseif ( equationOfState.eq.'MDJWF' ) then
693    
694     sp5 = sqrt(s1)
695    
696     p1 = pLoc*SItodBar
697     p1t1 = p1*t1
698    
699     rhoNum = eosMDJWFnum(0)
700     & + t1*(eosMDJWFnum(1)
701     & + t1*(eosMDJWFnum(2) + eosMDJWFnum(3)*t1) )
702     & + s1*(eosMDJWFnum(4)
703     & + eosMDJWFnum(5)*t1 + eosMDJWFnum(6)*s1)
704     & + p1*(eosMDJWFnum(7) + eosMDJWFnum(8)*t2
705     & + eosMDJWFnum(9)*s1
706     & + p1*(eosMDJWFnum(10) + eosMDJWFnum(11)*t2) )
707    
708    
709     den = eosMDJWFden(0)
710     & + t1*(eosMDJWFden(1)
711     & + t1*(eosMDJWFden(2)
712     & + t1*(eosMDJWFden(3) + t1*eosMDJWFden(4) ) ) )
713     & + s1*(eosMDJWFden(5)
714     & + t1*(eosMDJWFden(6)
715     & + eosMDJWFden(7)*t2)
716     & + sp5*(eosMDJWFden(8) + eosMDJWFden(9)*t2) )
717     & + p1*(eosMDJWFden(10)
718     & + p1t1*(eosMDJWFden(11)*t2 + eosMDJWFden(12)*p1) )
719    
720     rhoDen = 1.0/(epsln+den)
721    
722 mlosch 1.21 rhoLoc = rhoNum*rhoDen - rhoConst
723 mlosch 1.20
724     elseif( equationOfState .eq. 'IDEALG' ) then
725     C
726     else
727     write(msgbuf,'(3a)')
728     & ' FIND_RHO_SCALAR : equationOfState = "',
729     & equationOfState,'"'
730     call print_error( msgbuf, mythid )
731     stop 'ABNORMAL END: S/R FIND_RHO_SCALAR'
732     endif
733 mlosch 1.22
734     return
735     end
736    
737     CBOP
738     C !ROUTINE: LOOK_FOR_NEG_SALINITY
739     C !INTERFACE:
740     subroutine LOOK_FOR_NEG_SALINITY(
741     I bi, bj, iMin, iMax, jMin, jMax, k,
742     I sFld,
743     I myThid )
744    
745     C !DESCRIPTION: \bv
746     C *==========================================================*
747     C | o SUBROUTINE LOOK_FOR_NEG_SALINITY
748     C | looks for and fixes negative salinity values
749     C | this is necessary if the equation of state uses
750     C | the square root of salinity
751     C *==========================================================*
752     C |
753     C | k - is the Salt level
754     C |
755     C *==========================================================*
756     C \ev
757    
758     C !USES:
759     implicit none
760     C == Global variables ==
761     #include "SIZE.h"
762     #include "EEPARAMS.h"
763     #include "PARAMS.h"
764     #include "GRID.h"
765    
766     C !INPUT/OUTPUT PARAMETERS:
767     C == Routine arguments ==
768     integer bi,bj,iMin,iMax,jMin,jMax
769     integer k ! Level of Salt slice
770     _RL sFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
771     integer myThid
772    
773     C !LOCAL VARIABLES:
774     C == Local variables ==
775     integer i,j, localWarning
776     character*(max_len_mbuf) msgbuf
777     CEOP
778    
779     localWarning = 0
780     do j=jMin,jMax
781     do i=iMin,iMax
782     C abbreviations
783     if ( sFld(i,j,k,bi,bj) .lt. 0. _d 0 ) then
784     localWarning = localWarning + 1
785     sFld(i,j,k,bi,bj) = 0. _d 0
786     end if
787     end do
788     end do
789     C issue a warning
790     if ( localWarning .gt. 0 ) then
791     write(*,'(a,a)')
792     & 'S/R LOOK_FOR_NEG_SALINITY: found negative salinity',
793     & 'values and reset them to zero.'
794     write(*,'(a,I3)')
795     & 'S/R LOOK_FOR_NEG_SALINITY: current level is k = ', k
796     end if
797 mlosch 1.20
798     return
799     end

  ViewVC Help
Powered by ViewVC 1.1.22