/[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.5 - (hide annotations) (download)
Thu May 28 16:19:50 1998 UTC (26 years ago) by adcroft
Branch: MAIN
CVS Tags: checkpoint4
Changes since 1.4: +24 -12 lines
Changes to:
 o read in coefficients for POLY3 EOS.
 o find_rho() polynomial evaluation has been factorized.
 o additional density field needed in calc_iso_slopes() with non-linear EOS.

This EOS must use the appropriate version of KNUDSEN to generate the
coefficients file (POLY3.COEFFS). C7 and C8 were back to front in
all previous versions of the model (compare01).

1 adcroft 1.5 C $Header: /u/gcmpack/models/MITgcmUV/model/src/find_rho.F,v 1.4 1998/04/29 21:30:18 adcroft Exp $
2 cnh 1.1
3     #include "CPP_EEOPTIONS.h"
4 adcroft 1.5 #define USE_FACTORIZED_POLY
5 cnh 1.1
6     ! ==============================================================================
7 adcroft 1.4 subroutine FIND_RHO(
8     I bi, bj, iMin, iMax, jMin, jMax, k, kRef, eqn,
9     O rholoc,
10     I myThid )
11     C /==========================================================\
12     C | o SUBROUTINE FIND_RHO |
13     C | Calculates [rho(S,T,z)-Rhonil] of a slice |
14     C |==========================================================|
15     C | |
16     C | k - is the Theta/Salt level |
17     C | kRef - determines pressure reference level |
18     C | (not used in 'LINEAR' mode) |
19     C | eqn - determines the eqn. of state: 'LINEAR' or 'POLY3' |
20     C | |
21     C \==========================================================/
22 cnh 1.1 implicit none
23     ! Common
24     #include "SIZE.h"
25     #include "DYNVARS.h"
26     #include "PARAMS.h"
27 adcroft 1.4 ! Arguments
28     integer bi,bj,iMin,iMax,jMin,jMax
29     integer k ! Level of Theta/Salt slice
30     integer kRef ! Pressure reference level
31     character*(*) eqn
32     _RL rholoc(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
33     integer myThid
34 cnh 1.1 ! Local
35 adcroft 1.4 integer i,j
36     _RL refTemp,refSalt,sigRef,tP,sP,deltaSig
37 cnh 1.1 ! ------------------------------------------------------------------------------
38    
39 adcroft 1.4 if (eqn.eq.'LINEAR') then
40    
41     C ***NOTE***
42     C In the linear EOS, to make the static stability calculation meaningful
43     C we alway calculate the perturbation with respect to the surface level.
44     C **********
45     refTemp=tRef(1)
46     refSalt=sRef(1)
47    
48 cnh 1.1 do j=jMin,jMax
49     do i=iMin,iMax
50 adcroft 1.4 rholoc(i,j)=rhonil*(
51     & sBeta*( salt(i,j,k,bi,bj)-refSalt)
52     & -tAlpha*(theta(i,j,k,bi,bj)-refTemp) )
53 cnh 1.1 enddo
54     enddo
55 adcroft 1.4
56     elseif (eqn.eq.'POLY3') then
57    
58     refTemp=eosRefT(kRef)
59     refSalt=eosRefS(kRef)
60 adcroft 1.5 sigRef=eosSig0(kRef) + (1000.-Rhonil)
61 adcroft 1.4
62     do j=jMin,jMax
63     do i=iMin,iMax
64     tP=theta(i,j,k,bi,bj)-refTemp
65     sP=salt(i,j,k,bi,bj)-refSalt
66 adcroft 1.5 #ifdef USE_FACTORIZED_POLY
67 adcroft 1.4 deltaSig=
68 adcroft 1.5 & (( eosC(9,kRef)*sP + eosC(5,kRef) )*sP + eosC(2,kRef) )*sP
69     & + ( ( eosC(6,kRef)
70     & *tP
71     & +eosC(7,kRef)*sP + eosC(3,kRef)
72     & )*tP
73     & +(eosC(8,kRef)*sP + eosC(4,kRef) )*sP + eosC(1,kRef)
74     & )*tP
75     #else
76     deltaSig=
77     & eosC(1,kRef)*tP
78     & +eosC(2,kRef) *sP
79     & +eosC(3,kRef)*tP*tP
80     & +eosC(4,kRef)*tP *sP
81     & +eosC(5,kRef) *sP*sP
82     & +eosC(6,kRef)*tP*tP*tP
83     & +eosC(7,kRef)*tP*tP *sP
84     & +eosC(8,kRef)*tP *sP*sP
85     & +eosC(9,kRef) *sP*sP*sP
86     #endif
87     rholoc(i,j)=sigRef+deltaSig
88 adcroft 1.4 enddo
89     enddo
90    
91     else
92     write(0,*) 'FIND_RHO: eqn = ',eqn
93     stop 'FIND_RHO: argument "eqn" has illegal value'
94     endif
95 cnh 1.1
96     ! ------------------------------------------------------------------------------
97     return
98     end
99     ! ==============================================================================

  ViewVC Help
Powered by ViewVC 1.1.22