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

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

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


Revision 1.4 - (hide annotations) (download)
Tue Apr 10 22:35:25 2001 UTC (23 years, 2 months ago) by heimbach
Branch: MAIN
CVS Tags: checkpoint40pre3, checkpoint40pre1, checkpoint40pre7, checkpoint40pre6, checkpoint40pre9, checkpoint40pre8, checkpoint38, checkpoint40pre2, checkpoint40pre4, checkpoint39, checkpoint40pre5, checkpoint40
Changes since 1.3: +4 -8 lines
See doc/tag-index and doc/notes_c37_adj.txt
Preparation for stand-alone autodifferentiability.

1 heimbach 1.4 C $Header: /u/gcmpack/models/MITgcmUV/model/src/find_alpha.F,v 1.3 2001/02/04 14:38:47 cnh Exp $
2     C $Name: checkpoint37 $
3 adcroft 1.1
4     #include "CPP_OPTIONS.h"
5     #define USE_FACTORIZED_POLY
6    
7     subroutine FIND_ALPHA (
8     I bi, bj, iMin, iMax, jMin, jMax, k, kRef, eqn,
9     O alphaloc )
10     C /==========================================================\
11     C | o SUBROUTINE FIND_ALPHA |
12     C | Calculates [drho(S,T,z) / dT] of a horizontal slice |
13     C |==========================================================|
14     C | |
15     C | k - is the Theta/Salt level |
16     C | kRef - determines pressure reference level |
17     C | (not used in 'LINEAR' mode) |
18     C | eqn - determines the eqn. of state: 'LINEAR' or 'POLY3' |
19     C | |
20     C | alphaloc - drho / dT (kg/m^3/C) |
21     C | |
22     C \==========================================================/
23     implicit none
24    
25     c Common
26     #include "SIZE.h"
27     #include "DYNVARS.h"
28     #include "EEPARAMS.h"
29     #include "PARAMS.h"
30    
31     c Arguments
32     integer bi,bj,iMin,iMax,jMin,jMax
33     integer k ! Level of Theta/Salt slice
34     integer kRef ! Pressure reference level
35     character*(*) eqn
36     _RL alphaloc(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
37    
38     c Local
39     integer i,j
40     _RL refTemp,refSalt,tP,sP
41    
42     if (eqn.eq.'LINEAR') then
43    
44     do j=jMin,jMax
45     do i=iMin,iMax
46     alphaloc(i,j) = -rhonil * tAlpha
47     enddo
48     enddo
49    
50     elseif (eqn.eq.'POLY3') then
51    
52     refTemp=eosRefT(kRef)
53     refSalt=eosRefS(kRef)
54    
55     do j=jMin,jMax
56     do i=iMin,iMax
57     tP=theta(i,j,k,bi,bj)-refTemp
58     sP=salt(i,j,k,bi,bj)-refSalt
59     #ifdef USE_FACTORIZED_POLY
60     alphaloc(i,j) =
61     & ( eosC(6,kRef)
62     & *tP*3.
63     & +(eosC(7,kRef)*sP + eosC(3,kRef))*2.
64     & )*tP
65     & +(eosC(8,kRef)*sP + eosC(4,kRef) )*sP + eosC(1,kRef)
66     &
67     #else
68     alphaloc(i,j) =
69     & eosC(1,kRef) +
70     & eosC(3,kRef)*tP*2. +
71     & eosC(4,kRef) *sP +
72     & eosC(6,kRef)*tP*tP*3. +
73     & eosC(7,kRef)*tP*2. *sP +
74     & eosC(8,kRef) *sP*sP
75     #endif
76     enddo
77     enddo
78    
79     else
80 heimbach 1.4 write(*,*) 'FIND_ALPHA: eqn = ',eqn
81 adcroft 1.1 stop 'FIND_ALPHA: argument "eqn" has illegal value'
82     endif
83    
84     return
85     end
86    
87     subroutine FIND_BETA (
88     I bi, bj, iMin, iMax, jMin, jMax, k, kRef, eqn,
89     O betaloc )
90     C /==========================================================\
91     C | o SUBROUTINE FIND_BETA |
92     C | Calculates [drho(S,T,z) / dS] of a horizontal slice |
93     C |==========================================================|
94     C | |
95     C | k - is the Theta/Salt level |
96     C | kRef - determines pressure reference level |
97     C | (not used in 'LINEAR' mode) |
98     C | eqn - determines the eqn. of state: 'LINEAR' or 'POLY3' |
99     C | |
100     C | betaloc - drho / dS (kg/m^3/PSU) |
101     C | |
102     C \==========================================================/
103     implicit none
104    
105     c Common
106     #include "SIZE.h"
107     #include "DYNVARS.h"
108     #include "EEPARAMS.h"
109     #include "PARAMS.h"
110    
111     c Arguments
112     integer bi,bj,iMin,iMax,jMin,jMax
113     integer k ! Level of Theta/Salt slice
114     integer kRef ! Pressure reference level
115     character*(*) eqn
116     _RL betaloc(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
117    
118     c Local
119     integer i,j
120     _RL refTemp,refSalt,tP,sP
121    
122     if (eqn.eq.'LINEAR') then
123    
124     do j=jMin,jMax
125     do i=iMin,iMax
126     betaloc(i,j) = rhonil * sBeta
127     enddo
128     enddo
129    
130     elseif (eqn.eq.'POLY3') then
131    
132     refTemp=eosRefT(kRef)
133     refSalt=eosRefS(kRef)
134    
135     do j=jMin,jMax
136     do i=iMin,iMax
137     tP=theta(i,j,k,bi,bj)-refTemp
138     sP=salt(i,j,k,bi,bj)-refSalt
139     #ifdef USE_FACTORIZED_POLY
140     betaloc(i,j) =
141     & ( eosC(9,kRef)*sP*3. + eosC(5,kRef)*2. )*sP + eosC(2,kRef)
142     & + ( eosC(7,kRef)*tP
143     & +eosC(8,kRef)*sP*2. + eosC(4,kRef)
144     & )*tP
145     #else
146     betaloc(i,j) =
147     & eosC(2,kRef) +
148     & eosC(4,kRef)*tP +
149     & eosC(5,kRef) *sP*2. +
150     & eosC(7,kRef)*tP*tP +
151     & eosC(8,kRef)*tP *sP*2. +
152     & eosC(9,kRef) *sP*sP*3.
153     #endif
154     enddo
155     enddo
156    
157     else
158 heimbach 1.4 write(*,*) 'FIND_BETA: eqn = ',eqn
159 adcroft 1.1 stop 'FIND_BETA: argument "eqn" has illegal value'
160     endif
161    
162     return
163     end

  ViewVC Help
Powered by ViewVC 1.1.22