SUBROUTINE DENSCOEF( 1,6@DYALLOC.4429
*CALL ARGSIZE 
SF011193.46
* ZDZZ,XLAT) ! ############################################# @DYALLOC.4430
DENSCOEF.25
DENSCOEF.26
C RH141293.39
IMPLICIT NONE RH141293.40
C RH141293.41
INTEGER DENSCOEF.27
& LM ! NUMBER OF SALINITY POINTS IN RANGE DENSCOEF.28
& ,LMM ! NUMBER OF TEMPERATURE POINTS IN RANGE DENSCOEF.29
DENSCOEF.30
PARAMETER (LM=5,LMM=2*LM) DENSCOEF.31
DENSCOEF.32
DENSCOEF.33
*CALL TYPSIZE @DYALLOC.4431
*CALL COCSTATE DENSCOEF.35
DENSCOEF.36
INTEGER DENSCOEF.37
& I ! LOOP INDEX DENSCOEF.38
& ,K ! LOOP INDEX DENSCOEF.39
& ,MP ! LOOP INDEX DENSCOEF.40
& ,KK ! TOTAL NUMBER OF TS POINTS (LN*LMM) DENSCOEF.41
& ,J ! LOOP INDEX DENSCOEF.42
& ,N ! NUMBER OF COLUMNS OF A FOR LEAST SQUARES DENSCOEF.43
& ,IN ! OPTION NUMBER FOR LEAST SQUARES DENSCOEF.44
& ,ITMAX ! MAXIMUM NUMBER OF ITERATIONS DENSCOEF.45
& ,IT ! VARIABLE FOR LEAST SQUARES DENSCOEF.46
& ,IEQ ! VARIABLE FOR LEAST SQUARES DENSCOEF.47
& ,IRANK ! VARIABLE FOR LEAST SQUARES DENSCOEF.48
& ,NHDIM ! VARIABLE FOR LEAST SQUARES DENSCOEF.49
& ,L ! LOOP INDEX DENSCOEF.50
DENSCOEF.51
REAL DENSCOEF.52
& ZDZZ(KM+1) ! DEPTH OF LAYER MIDPOINTS IN CM DENSCOEF.53
& ,TS(33,4) ! RANGE OF T/S OVER WHICH POLYNOMIAL DENSCOEF.54
& ,DT(KM) ! MAXIMUM TEMPERATURE DENSCOEF.55
& ,DS(KM) ! MAXIMUM SALINITY FOR MODEL LEVEL DENSCOEF.56
& ,DD(KM) ! INCREMENT IN T DENSCOEF.57
& ,SS(KM) ! INCREMENT IN S DENSCOEF.58
& ,TA(LMM) ! POINTS IN T RANGE DENSCOEF.59
& ,SA(LMM) ! POINTS IN S RANGE @DYALLOC.4432
& ,TP(LM*LMM) ! POINTS OF T IN T*S RANGE DENSCOEF.61
& ,SP(LM*LMM) ! POINTS OF S IN T*S RANGE DENSCOEF.62
& ,T1 ! TEMPORARY ACCUMULATION OF TEMPERATURE DENSCOEF.63
& ,S1 ! TEMPORARY ACCUMULATION OF SALINITY DENSCOEF.64
& ,TOT ! ACCUMULATION OF DENSITY DENSCOEF.65
& ,TH1 ! ACCUMULATION OF POTENTIAL TEMPERATURE DENSCOEF.66
& ,FK ! REAL (KK) DENSCOEF.67
& ,D ! DEPTH IN METRES DENSCOEF.68
& ,S ! SALINITY DENSCOEF.69
& ,T ! TEMPERATURE DENSCOEF.70
& ,PIN ! PRESSURE AT DEPTH Z DENSCOEF.71
& ,XLAT ! REFERENCE LATITUDE DENSCOEF.72
& ,B(LM*LMM) ! DENSITY AT T*S POINTS DENSCOEF.73
& ,Delta_T(LM*LMM)! temperature (difference from reference) OMB1F404.9
! at T*S points OMB1F404.10
& ,PEND ! SURFACE PRESSURE DENSCOEF.74
& ,DPP ! PRESSURE INCREMENT DENSCOEF.75
& ,TH(LM*LMM) ! POTENTIAL TEMPERATURES AT T*S POINTS DENSCOEF.76
& ,TT ! T DIFFERENCE DENSCOEF.78
& ,SL ! S DIFFERENCE DENSCOEF.79
& ,A(LM*LMM,9) ! THIRD ORDER POLYNOMIAL DENSCOEF.80
& ,EPS ! PERMISSIBLE ERROR DENSCOEF.81
& ,X(9) ! KNUDSEN COEFICIENTS FOR LEVEL DENSCOEF.82
& ,Y(9) ! temperature coeffs for a given level OMB1F404.11
& ,ENORM ! VARIABLE FOR LEAST SQUARES DENSCOEF.83
& ,H ! VARIABLE FOR LEAST SQUARES DENSCOEF.84
& ,CARRAY(LM*LMM,LM*LMM) !WORK SPACE FOR LEAST SQUARES DENSCOEF.85
& ,R(4*LM*LMM) ! WORK SPACE FOR LEAST SQUARES DENSCOEF.86
& ,SB(4*LM*LMM) ! WORK SPACE FOR LEAST SQUARES DENSCOEF.87
DENSCOEF.88
REAL DN ! DENSITY DENSCOEF.89
DENSCOEF.90
REAL DENSCOEF.91
& fnztop ! function to convert depth to pressure DENSCOEF.92
&,POTTEM ! function to calculate potential temperature DENSCOEF.93
DENSCOEF.94
EXTERNAL DENSCOEF.95
& fnztop DENSCOEF.96
&,POTTEM DENSCOEF.97
&,UNESCO DENSCOEF.98
DENSCOEF.99
C DENSCOEF.100
C DENSCOEF.101
C BOUNDS FOR POLYNOMIAL FIT: DENSCOEF.102
C TS(K,1)=LOWER BND OF T AT Z=(K-1)*250 METERS DENSCOEF.103
C TS(K,2)=UPPER BND OF T " DENSCOEF.104
C TS(K,3)=LOWER BND OF S " DENSCOEF.105
C TS(K,4)=UPPER BND OF S " DENSCOEF.106
C DENSCOEF.107
DATA TS/4*-2.,15*-1.,14*0.,29.,19.,14.,11.,9.,28*7. DENSCOEF.108
1,28.5,33.7,34.0,34.1,34.2,34.4,2*34.5,15*34.6,10*34.7,37.0,36.6, DENSCOEF.109
235.8,35.7,35.3,2*35.1,26*35.0/ DENSCOEF.110
C DENSCOEF.111
C NOTE: THE MAXIMUM ALLOWABLE DEPTH IS 8000 METERS DENSCOEF.112
C DENSCOEF.113
C DENSCOEF.114
C SET UPPER & LOWER BOUNDS ON T & S RANGES FOR EACH MODEL LEVEL DENSCOEF.115
C DENSCOEF.116
DO 200 I=1,KM DENSCOEF.117
C-- ALSO CONVERT TO METRES BY /100 DENSCOEF.118
K=INT(ZDZZ(I)/250./100.0)+1 DENSCOEF.119
TO(I)=TS(K,1) DENSCOEF.120
DT(I)=TS(K,2) DENSCOEF.121
SO(I)=TS(K,3) DENSCOEF.122
DS(I)=TS(K,4) DENSCOEF.123
200 CONTINUE DENSCOEF.124
C DENSCOEF.125
C CALCULATE INCREMENTS IN T & S BETWEEN LOWER & UPPER BOUNDS DENSCOEF.126
C DENSCOEF.127
DO 701 MP=1,KM DENSCOEF.128
DD(MP)=(DT(MP)-TO(MP))/(2.*LM-1.0) DENSCOEF.129
SS(MP)=(DS(MP)-SO(MP))/(LM-1.0) DENSCOEF.130
701 CONTINUE DENSCOEF.131
KK=LM*LMM DENSCOEF.132
C DENSCOEF.133
C LOOP OVER LEVELS DENSCOEF.134
C DENSCOEF.135
DO 1000 MP=1,KM DENSCOEF.136
C DENSCOEF.137
C CALCULATE EVENLY-SPACED T & S POINTS DENSCOEF.138
C (SA INCLUDES REDUNDANT POINTS GREATER THAN SMAX) DENSCOEF.139
C DENSCOEF.140
DO 900 I=1,LMM DENSCOEF.141
TA(I)=TO(MP)+(FLOAT(I)-1.0)*DD(MP) DENSCOEF.142
SA(I)=SO(MP)+(FLOAT(I)-1.0)*SS(MP) DENSCOEF.143
900 CONTINUE DENSCOEF.144
C DENSCOEF.145
C CREATE 1-D ARRAY WITH ALL COMBINATIONS OF T & S POINTS DENSCOEF.146
C DENSCOEF.147
DO 800 I=1,LMM DENSCOEF.148
DO 800 J=1,LM DENSCOEF.149
K=LM*I+J-LM DENSCOEF.150
TP(K)=TA(I) DENSCOEF.151
SP(K)=SA(J) DENSCOEF.152
800 CONTINUE DENSCOEF.153
C DENSCOEF.154
C INITIALISE MEAN QUANTITIES DENSCOEF.155
C DENSCOEF.156
T1=0.0 DENSCOEF.157
S1=0.0 DENSCOEF.158
TOT=0.0 DENSCOEF.159
TH1=0.0 DENSCOEF.160
FK=KK DENSCOEF.161
C DENSCOEF.162
C CALCULATE DENSITY AND POTENTIAL TEMPERATURE DENSCOEF.163
C DENSCOEF.164
DO 700 K=1,KK DENSCOEF.165
D=ZDZZ(MP)/100.0 DENSCOEF.166
S=SP(K) DENSCOEF.167
T=TP(K) DENSCOEF.168
DENSCOEF.169
C-- CALCULATE THE PRESSURE AT GIVEN DEPTH DENSCOEF.170
PIN=fnztop(D,XLAT) DENSCOEF.171
DENSCOEF.172
C-- CALCULATE IN-SITU DENSITY DENSCOEF.173
CALL UNESCO(T,S,PIN,DN) DENSCOEF.174
B(K)=DN DENSCOEF.175
DENSCOEF.176
C-- CALCULATE POTENTIAL TEMPERATURE DENSCOEF.177
PEND=0.0 DENSCOEF.178
DPP=1.0 DENSCOEF.179
TH(K)=POTTEM(T,S,PIN,PEND,DPP) DENSCOEF.180
C DENSCOEF.181
C SUM T,S,DENSITY AND POTENTIAL TEMPERATURE DENSCOEF.182
C DENSCOEF.183
T1=T1+TP(K) DENSCOEF.184
S1=S1+SP(K) DENSCOEF.185
TOT=TOT+B(K) DENSCOEF.186
TH1=TH1+TH(K) DENSCOEF.187
700 CONTINUE DENSCOEF.188
C DENSCOEF.189
C CALCULATE "MID-POINT" VALUES T,S,DENSITY AND POTENTIAL TEMPERATURE DENSCOEF.190
C AND DENSITY FOR MID-POINT T,S DENSCOEF.191
C DENSCOEF.192
T1=T1/FK DENSCOEF.193
S1=S1/FK DENSCOEF.194
TOT=TOT/FK DENSCOEF.195
DENSCOEF.196
C-- CALCULATE IN-SITU DENSITY DENSCOEF.197
CALL UNESCO(T1,S1,PIN,DN) DENSCOEF.198
TOT=DN DENSCOEF.199
TH1=TH1/FK DENSCOEF.200
C DENSCOEF.201
C PUT MID-POINT VALUES IN OUTPUT ARRAY DENSCOEF.202
C DENSCOEF.203
SIGO(MP)=TOT DENSCOEF.204
TO(MP)=TH1 DENSCOEF.205
SO(MP)=S1 DENSCOEF.206
TempO(MP) = T1 OMB1F404.12
C DENSCOEF.207
C CALCULATE DIFFERENCE FROM MID-POINT VALUES DENSCOEF.208
C DENSCOEF.209
DO 600 K=1,KK DENSCOEF.210
TT=TH(K)-TH1 DENSCOEF.211
SL=SP(K)-S1 DENSCOEF.212
B(K)=B(K)-TOT DENSCOEF.213
Delta_T(K)=TP(K)-T1 OMB1F404.13
C DENSCOEF.214
C CALCULATE VALUES OF THE 9 TERMS OF THE 3RD-ORDER POLYNOMIAL DENSCOEF.215
C DENSCOEF.216
A(K,1)=TT DENSCOEF.217
A(K,2)=SL DENSCOEF.218
A(K,3)=TT*TT DENSCOEF.219
A(K,4)=TT*SL DENSCOEF.220
A(K,5)=SL*SL DENSCOEF.221
A(K,6)=A(K,3)*TT DENSCOEF.222
A(K,7)=A(K,5)*TT DENSCOEF.223
A(K,8)=A(K,3)*SL DENSCOEF.224
A(K,9)=A(K,5)*SL DENSCOEF.225
600 CONTINUE DENSCOEF.226
C DENSCOEF.227
C SET THE ARGUMENTS IN CALL TO LSQSL2 DENSCOEF.228
C NUMBER OF COLUMNS OF A DENSCOEF.229
N=9 DENSCOEF.230
C OPTION NUMBER OF LSQSL2 DENSCOEF.231
IN=1 DENSCOEF.232
C DENSCOEF.233
C ITMAX=NUMBER OF ITERATIONS DENSCOEF.234
ITMAX=4 DENSCOEF.235
C DENSCOEF.236
IT=0 DENSCOEF.237
IEQ=2 DENSCOEF.238
IRANK=0 DENSCOEF.239
EPS=1.0E-7 DENSCOEF.240
NHDIM=9 DENSCOEF.241
C DENSCOEF.242
C CALL LEAST SQUARES ROUTINE TO MINIMISE (AX-B) DENSCOEF.243
C (finding X which holds polynomial coefficients to calculate density) OMB1F404.14
C DENSCOEF.244
CALL LSQSL2(KK,A,KK,N,B,X,IRANK,IN,ITMAX,IT,IEQ,ENORM,EPS, DENSCOEF.245
* NHDIM,H,CARRAY,R,SB) DENSCOEF.246
C OMB1F404.15
C CALL LEAST SQUARES ROUTINE TO MINIMISE (AY-Delta_T) OMB1F404.16
C (finding Y which holds polynomial coeffs to calculate temperature) OMB1F404.17
C Note that IN=2 option no longer works properly so IN=1 used again OMB1F404.18
C CARRAY holds copy of A; so CARRAY and A are interchanged in this call OMB1F404.19
C OMB1F404.20
OMB1F404.21
IN = 1 OMB1F404.22
CALL LSQSL2(KK,CARRAY,KK,N,Delta_T,Y,IRANK,IN,ITMAX,IT,IEQ, OMB1F404.23
* ENORM,EPS,NHDIM,H,A,R,SB) OMB1F404.24
OMB1F404.25
C DENSCOEF.247
C PUT COEFFICIENTS INTO OUTPUT ARRAY DENSCOEF.248
C DENSCOEF.249
DO 324 I=1,N DENSCOEF.250
C(MP,I)=X(I) DENSCOEF.251
coeff_T(MP,I)=Y(I) OMB1F404.26
324 CONTINUE DENSCOEF.252
C DENSCOEF.253
1000 CONTINUE DENSCOEF.254
C DENSCOEF.255
C DENSCOEF.256
C OMIT NEXT SECTION IF OUTPUT NOT REQUIRED DENSCOEF.257
C DENSCOEF.258
C DENSCOEF.259
C SCALE COEFFICIENTS DENSCOEF.260
C DENSCOEF.261
DO 1250 L=1,KM DENSCOEF.262
SIGO(L)=SIGO(L)*1.E-3 DENSCOEF.263
SO(L)=1.E-3*SO(L)-0.035 DENSCOEF.264
DENSCOEF.265
C(L,1)=1.E-3*C(L,1) DENSCOEF.266
C(L,3)=1.E-3*C(L,3) DENSCOEF.267
C(L,6)=1.E-3*C(L,6) DENSCOEF.268
C(L,5)=1.E+3*C(L,5) DENSCOEF.269
C(L,7)=1.E+3*C(L,7) DENSCOEF.270
C(L,9)=1.E+6*C(L,9) DENSCOEF.271
coeff_T(L,2)=1.E+3*coeff_T(L,2) OMB1F404.27
coeff_T(L,4)=1.E+3*coeff_T(L,4) OMB1F404.28
coeff_T(L,8)=1.E+3*coeff_T(L,8) OMB1F404.29
coeff_T(L,5)=1.E+6*coeff_T(L,5) OMB1F404.30
coeff_T(L,7)=1.E+6*coeff_T(L,7) OMB1F404.31
coeff_T(L,9)=1.E+9*coeff_T(L,9) OMB1F404.32
DENSCOEF.272
1250 CONTINUE DENSCOEF.273
C DENSCOEF.274
1200 CONTINUE DENSCOEF.275
1288 CONTINUE DENSCOEF.276
C DENSCOEF.277
C-------------- PRINT OUTS -------------- DENSCOEF.278
WRITE (6,*)'----- PRINTOUT FOR SUBROUTINE DENSCOEF -----' DENSCOEF.279
WRITE (6,*) DENSCOEF.280
write (6,*)'Reference Latitude =',xlat DENSCOEF.281
write (6,*) DENSCOEF.282
DENSCOEF.283
DO K=1,KM DENSCOEF.284
WRITE (6,*)'LEVEL=',K,' TO=',TO(K) DENSCOEF.285
& ,' SO=',SO(K),' SIGO=',SIGO(K),' TempO=',TempO(K) OMB1F404.33
WRITE (6,*)'Density polynomial coeffs are:',(C(K,I),I=1,9) OMB1F404.34
WRITE (6,*)'Temperature poly coeffs are:',(coeff_T(K,I),I=1,9) OMB1F404.35
WRITE (6,*) DENSCOEF.288
END DO DENSCOEF.289
DENSCOEF.290
RETURN DENSCOEF.291
END DENSCOEF.292
*ENDIF @DYALLOC.4433