SUBROUTINE VERTCOFC 2,14VERTCOFC.15
& ( J,IMT,KM,KMM1,NT, VERTCOFC.16
& JMT, ORH7F404.45
& TB,TBP, UB,VB, VERTCOFC.17
& DZZ2RQ,DZ2RQ, OLA3F403.48
& NERGY,GRAV_SI,GM, VERTCOFC.19
& RHOSRN,RHOSRNA,RHOSRNB, OOM1F405.900
& WSX, WSY, OLA3F403.49
& ZDZZ,ZDZ,DZ,DZZ, OOM1F405.901
& Rim,hm,max_qLarge_depth,crit_Ri, OLA0F404.202
& L_M,MAX_LARGE_DEPTH,MAX_LARGE_LEVELS,RHO_WATER_SI, OOM1F405.902
& MLD_LARGE,MLD_LARGEP,HTN,HTNP,PME,PMEP,SOL,SOLP, OOM1F405.903
& WATERFLUX_ICE,WATERFLUX_ICEP,LAMBDA_LARGE,SPECIFIC_HEAT_SI, OOM1F405.904
& WME,WMEP,L_OWINDMIX,L_OBULKMAXMLD, OOM1F405.905
& PHI, OOM1F405.906
& GNU,FNU0_SI,FNUB_SI,STABLM_SI,GNUMINC_SI OOM1F405.907
&,OCEANHEATFLUX,OCEANHEATFLUXP OOM1F405.908
&,CARYHEAT,CARYHEATP OOM1F405.909
&,FLXTOICE,FLXTOICEP) OOM1F405.910
C VERTCOFC.22
C VERTCOFC.23
IMPLICIT NONE VERTCOFC.24
*CALL CNTLOCN 
ORH1F305.5358
*CALL OARRYSIZ ORH1F305.5359
C VERTCOFC.25
*CALL C_VKMAN OLA0F404.203
*CALL C_PI OOM1F405.913
*CALL C_OMEGA OOM1F405.914
INTEGER VERTCOFC.26
& I,J,K,IMT,KM,KMM1,NT,NERGY VERTCOFC.27
&,JMT ! IN No of rows ORH1F305.5360
C VERTCOFC.28
REAL VERTCOFC.29
& TB(IMT,KM,NT),TBP(IMT,KM,NT),UB(IMT,KM),VB(IMT,KM), VERTCOFC.30
& GM(IMT,KM),DZZ2RQ(IMT,KM),DZ2RQ(IMT,KM),GRAV_SI, OLA3F403.52
& rhosrn(IMT,KM), ! OUT Density on TS row to S of vel row J VERTCOFC.32
& RHOSRNA(IMT,KM+1), ! OUT DENSITY ON TS ROW TO S OF VEL ROW J OOM1F405.911
& RHOSRNB(IMT,KM+1), ! OUT DENSITY ON TS ROW TO S OF VEL ROW J OOM1F405.912
& gnu(IMT,KM), ! OUT Vertical viscosity (cm2/s) at TOP of box VERTCOFC.33
& FNU0_SI, ! IN Max value of stability-dependent term VERTCOFC.34
C of vert. momentum diffusivity (m2/s) VERTCOFC.35
& FNUB_SI, ! IN Background value of momentum VERTCOFC.36
C diffusivity (m2/s) VERTCOFC.37
& STABLM_SI ! IN Max value of diffusivity (m2/s) VERTCOFC.38
& ,GNUMINC_SI ! IN Min value of momentum diffusivity VERTCOFC.39
C between top two levels (m2/s) VERTCOFC.40
REAL OOM1F405.920
& MLD_LARGE(IMT) ! IN MIXED LAYER DEPTH ON T GRID, ROW J (CM) OOM1F405.921
&, MLD_LARGEP(IMT) ! IN MIXED LAYER DEPTH ON T GRID, ROW J+1 (CM) OOM1F405.922
&, WME(IMT),WMEP(IMT) ! IN WIND MIXING ENERGY FOR ROW J,J+1 OOM1F405.923
&, DZZ(KM+1) ! IN DISTANCE BETWEEN MIDPTS OF LEVELS (CM) OOM1F405.924
&, MAX_LARGE_DEPTH ! IN MAX DEPTH LARGE CAN BE APPLIED TO (M) OOM1F405.925
&, RHO_WATER_SI ! IN DENSITY OF SEA WATER = 1026. KG/M^3 OOM1F405.926
&, L_M(IMT) ! OUT MONIN OBUKHOV LENGTH LARGE SCHEME (MOMENTUM) OOM1F405.927
&, ALPHAI(IMT,1),BETAI(IMT,1) ! EXPANSION COEFF-TRACER GRID (J) OOM1F405.928
&, ALPHAIP(IMT,1),BETAIP(IMT,1) ! EXPANSION COEFF-TRACER GRID (J+1) OOM1F405.929
REAL OOM1F405.930
& HTN(IMT) ! IN NON-PENETRATING HEAT FLUX (W/M2) ON ROW J OOM1F405.931
&, HTNP(IMT) ! IN NON-PENETRATING HEAT FLUX (W/M2) ON ROW J+1 OOM1F405.932
&, PME(IMT) ! IN PRECIP MINUS EVAP (KG/M2/S) ON ROW J OOM1F405.933
&, PMEP(IMT) ! IN PRECIP MINUS EVAP (KG/M2/S) ON ROW J+1 OOM1F405.934
&, SOL(IMT) ! IN SOLAR IRRADIANCE (W/M2) AT SURFACE ON ROW J OOM1F405.935
&, SOLP(IMT) ! IN SOLAR IRRADIANCE (W/M2) AT SURFACE ON ROW J+1 OOM1F405.936
&, WATERFLUX_ICE(IMT) ! IN WATER FLUX FROM ICE (KG/M2/S) ,ROW J OOM1F405.937
&, WATERFLUX_ICEP(IMT) ! IN WATER FLUX FROM ICE (KG/M2/S) ,ROW J+1 OOM1F405.938
&, LAMBDA_LARGE,SPECIFIC_HEAT_SI ! IN FOR CALCULATING MINIMUM MLD OOM1F405.939
&, PHI ! IN LATITUDE OF ROW J ON UV GRID (RADIANS) OOM1F405.940
INTEGER OOM1F405.941
& MAX_LARGE_LEVELS ! IN MAX NO OF LEVS FOR LARGE SCHEME OOM1F405.942
LOGICAL L_OWINDMIX,L_OBULKMAXMLD OOM1F405.943
REAL WSX(IMT),WSY(IMT) ! surface stress (in SI units, N m^-2) OLA3F403.53
REAL ZDZZ(KM+1) ! IN DEPTH OF MIDDLE OF LAYER (CM) OOM1F405.947
REAL ZDZ ( KM) ! depth of bottom of layer (cm) OLA3F403.55
REAL DZ(KM) ! IN distance between half-levels (cm) OLA0F404.210
REAL OLA0F404.204
& hm(IMT_QLARGE) ! OUT depth quadratic Large scheme applied to OLA0F404.205
&, Rim(IMT,KMM1) ! OUT Richardson no at BOTTOM of box k OLA0F404.206
&, max_qLarge_depth ! IN max depth allowed for quadratic Large OLA0F404.207
&, crit_Ri ! IN critical Richardson no used for depth of OLA0F404.208
! quadratic Large scheme OLA0F404.209
REAL OOM1F405.915
& OCEANHEATFLUX(IMT),OCEANHEATFLUXP(IMT) OOM1F405.916
& !HTN:NON-PENETRATIVE SURFACE HEATFLUX INTO OCEAN BUDGET OOM1F405.917
&,CARYHEAT(IMT),CARYHEATP(IMT) !MISCELLANEOUS HEATFLUX FROM ICE OOM1F405.918
&,FLXTOICE(IMT),FLXTOICEP(IMT) !OCEAN TO ICE HEATFLUX OOM1F405.919
C VERTCOFC.41
C VERTCOFC.42
C Declare local variables and arrays VERTCOFC.43
C VERTCOFC.44
REAL VERTCOFC.45
& ripq(IMT,KMM1), ! d(rho)/dz at BOTTOM of box k VERTCOFC.46
& dvel2q(IMT,KMM1),! Uz*Uz + Vz*Vz at BOTTOM of box k VERTCOFC.47
& rhonrn(IMT,KM), ! Density on TS row to N of vel row J VERTCOFC.48
& RHONRNA(IMT,KM+1), ! DENSITY ON TS ROW TO N OF VEL ROW J OOM1F405.944
& RHONRNB(IMT,KM+1), ! DENSITY ON TS ROW TO N OF VEL ROW J OOM1F405.945
& worka(IMT,KM), ! Workspace VERTCOFC.49
& workb(IMT,KM), ! Workspace VERTCOFC.50
& uz, ! Uz at BOTTOM of current box VERTCOFC.51
& vz, ! Vz at BOTTOM of current box VERTCOFC.52
& ustar(IMT_QLARGE),! surface friction velocity on velocity grid OLA3F403.58
& Kath, ! K at z=h OLA3F403.59
& b2, ! value used in calculation of gnu for Large. OLA3F403.60
& sigma, ! depth/h OLA3F403.61
& RHOZ,RHOZA,RHOZB, OOM1F405.946
& fx, VERTCOFC.54
& fxa, VERTCOFC.55
& fxb VERTCOFC.56
C VERTCOFC.57
REAL OOM1F405.948
& GNUATHM ! VERTICAL DIFFUSION COEFF (CM2/S) AT LEVEL HM OOM1F405.949
&, GNUZATHM ! D(GNU)/DZ AT HM OOM1F405.950
&, ZDZ_AT_MLKM1 ! ZDZ AT LEVEL MLK-1 OOM1F405.951
&, GNUOFMLK ! GNU(I,MLK(I)) OOM1F405.952
&, GNUZ(2) ! D(GNU)/DZ AT MIDPOINTS OF LEVELS MLK,MLK+1 OOM1F405.953
&, COEFF ! PROPORTION OF DEPTH BELOW TOP OF BOX FOR MLD OOM1F405.954
&, RHO_FRESHWATER ! DENSITY OF FRESH WATER = 1.0 G/CM^3 OOM1F405.955
&, RHO0 ! DENSITY OF SEA WATER = 1.026 G/CM^3 OOM1F405.956
&, WATERFLUX_MOM,SFC_HEATFLUX_MOM ! WATERFLUX ON UV GRID IN CGS OOM1F405.957
&, SOL_MOM ! SOL ON UV GRID IN CGS UNITS OOM1F405.958
&, WME_MOM ! WIND MIXING ENERGY ON UV GRID IN CGS UNITS OOM1F405.959
&, GRAV ! GRAV_SI IN CGS UNITS OOM1F405.960
&, ALPHA,BETA ! THERMODYNAMIC EXPANSION COEFFS FOR T,S OOM1F405.961
&, CP0 ! SPECIFIC HEAT CAPACITY AT CONST PRESSURE AT SURFACE(CGS) OOM1F405.962
&, BF ! BUOYANCY FREQUENCY OOM1F405.963
REAL OOM1F405.964
& W1M ! TURBULENT VELOCITY SCALE AT MLD OOM1F405.965
&, DW1M ! DERIV OF TURBULENT VELOCITY SCALE WRT SIGMA AT MLD OOM1F405.966
&, STABIL_PARAM ! STABILITY PARAMETER OOM1F405.967
&, G1 ! SHAPE FUNCTION WHEN SIGMA=Z/H=1 OOM1F405.968
&, DG1 ! DG/DSIGMA EVALUATED AT SIGMA=1 OOM1F405.969
&, A(3) ! COEFFICIENTS OF CUBIC SHAPE FUNCTION OOM1F405.970
&, WM(KM) ! TURBULENT VELOCITY SCALE PROFILE AT TOP OF BOX K OOM1F405.971
&, EPSILON ! =0.1 (LARGE P365) OOM1F405.972
&, G ! SHAPE FUNCTION AT TOP OF BOX K. OOM1F405.973
INTEGER OOM1F405.974
& MLK(IMT) ! LEVEL WITHIN WHICH MLD_LARGE IS CONTAINED OOM1F405.975
&, IP1 ! I+1 OOM1F405.976
INTEGER I2 OOM1F405.977
C DEFINE FUNCTIONS OOM1F405.978
REAL OOM1F405.979
& SOL_IRRA_1B OOM1F405.980
&, FLUX_PROFILEM OOM1F405.981
C PARAMETERS FOR PETERS SCHEME OOM1F405.982
C OOM1F405.983
REAL RICRITM,RICRITM100 ! CRITICAL NUMBERS FOR MOMENTUM OOM1F405.984
PARAMETER(RICRITM = 0.39331063,RICRITM100 = 0.2288417) OOM1F405.985
C OOM1F405.986
INTEGER l ! level within which hm is contained OLA0F404.211
! for quadratic Large scheme OLA0F404.212
CL Parameters OLA3F403.65
C VERTCOFC.58
C Zero out gnu VERTCOFC.59
C VERTCOFC.60
fx=0. VERTCOFC.61
DO 110 K=1,KM VERTCOFC.62
DO 110 I=1,IMT VERTCOFC.63
gnu(I,K)=fx VERTCOFC.64
110 CONTINUE VERTCOFC.65
C CALCULATE SURFACE FRICTION VELOCITY IN CGS UNITS OOM1F405.1116
IF (L_OFULARGE) THEN OOM1F405.1117
RHO_FRESHWATER=1.0 ! IN G/CM^3 OOM1F405.1118
RHO0=RHO_WATER_SI*0.001 ! CONVERT TO CGS UNITS OOM1F405.1119
GRAV=GRAV_SI*100. ! CONVERT TO CGS UNITS OOM1F405.1120
EPSILON=0.1 OOM1F405.1121
IF (L_OUSTARWME) THEN OOM1F405.1122
C 1000.0 CONVERTS FROM SI TO CGS OOM1F405.1123
C USTAR^3 = WME/RHO0 OOM1F405.1124
DO I=1,IMT-1 OOM1F405.1125
WME_MOM=1000.0*0.25*(WME(I)+WME(I+1)+WMEP(I)+WMEP(I+1)) OOM1F405.1126
IF (GM(I,1).GT.0.5) THEN OOM1F405.1127
USTAR(I)=(WME_MOM/RHO0)**(1.0/3.0) OOM1F405.1128
ELSE OOM1F405.1129
USTAR(I)=0.0 OOM1F405.1130
ENDIF OOM1F405.1131
END DO OOM1F405.1132
IF (L_OCYCLIC) THEN OOM1F405.1133
USTAR(IMT)=USTAR(2) OOM1F405.1134
USTAR(IMT-1)=USTAR(1) OOM1F405.1135
ELSE OOM1F405.1136
WME_MOM=1000.0*0.5*(WME(IMT)+WMEP(IMT)) OOM1F405.1137
IF (GM(IMT,1).GT.0.5) THEN OOM1F405.1138
USTAR(IMT)=(WME_MOM/RHO0)**(1.0/3.0) OOM1F405.1139
ELSE OOM1F405.1140
USTAR(IMT)=0.0 OOM1F405.1141
ENDIF OOM1F405.1142
ENDIF OOM1F405.1143
ELSE OOM1F405.1144
C 10.0 CONVERTS FROM N M^-2 TO DYNES CM^-2 OOM1F405.1145
C USTAR^2 = TAU/RHO0 OOM1F405.1146
DO I=1,IMT OOM1F405.1147
USTAR(I)= SQRT(WSX(I)*WSX(I)+WSY(I)*WSY(I)) OOM1F405.1148
USTAR(I)= 10.0 * USTAR(I) / RHO0 OOM1F405.1149
USTAR(I)= SQRT(USTAR(I)) OOM1F405.1150
END DO OOM1F405.1151
ENDIF OOM1F405.1152
C DO PRELIMINARY CALCULATION OF HM AS NEEDED TO CALCULATE BUOYANCY OOM1F405.1153
C FREQUENCY OOM1F405.1154
DO I=1,IMT-1 OOM1F405.1155
HM(I)=0.25*(MLD_LARGE(I)+MLD_LARGE(I+1)+ OOM1F405.1156
& MLD_LARGEP(I)+MLD_LARGEP(I+1)) OOM1F405.1157
HM(I)=MIN(HM(I),MAX_LARGE_DEPTH*100.) OOM1F405.1158
HM(I)=MAX(HM(I),1.5*ZDZ(1)) OOM1F405.1159
ENDDO OOM1F405.1160
IF (L_OCYCLIC) THEN OOM1F405.1161
HM(IMT)=HM(2) OOM1F405.1162
HM(IMT-1)=HM(1) OOM1F405.1163
ELSE OOM1F405.1164
HM(IMT)=0.5*(MLD_LARGE(IMT)+MLD_LARGEP(IMT)) OOM1F405.1165
HM(IMT)=MIN(HM(IMT),MAX_LARGE_DEPTH*100.) OOM1F405.1166
HM(IMT)=MAX(HM(IMT),1.5*ZDZ(1)) OOM1F405.1167
ENDIF OOM1F405.1168
C OOM1F405.1169
C CALCULATE SURFACE BUOYANCY FORCING , BF, IN CM^2/S^3 OOM1F405.1170
C ************************************************************** OOM1F405.1171
C SET UP ALPHA AND BETA AT THE SURFACE FOR ROWS J AND J+1 OOM1F405.1172
C OOM1F405.1173
CALL DRODT(TB,TB(1,1,2),ALPHAI,IMT,1) OOM1F405.1174
CALL DRODS(TB,TB(1,1,2),BETAI,IMT,1) OOM1F405.1175
CALL DRODT(TBP,TBP(1,1,2),ALPHAIP,IMT,1) OOM1F405.1176
CALL DRODS(TBP,TBP(1,1,2),BETAIP,IMT,1) OOM1F405.1177
C OOM1F405.1178
CP0 = SPECIFIC_HEAT_SI*1.E4 ! IN CM^2 S^-2 K^-1 OOM1F405.1179
C ************************************************************** OOM1F405.1180
DO I=1,IMT OOM1F405.1181
L_M(I)=0.0 OOM1F405.1182
IF (GM(I,1).GT.0.5) THEN OOM1F405.1183
IP1=I+1 OOM1F405.1184
IF (I.EQ.IMT) THEN OOM1F405.1185
C---- I EQ IMT IS EQUIVALENT TO I=2 ON VELOCITY GRID, HENCE... OOM1F405.1186
IP1=3 OOM1F405.1187
IF (.NOT.L_OCYCLIC) IP1=IMT OOM1F405.1188
ENDIF OOM1F405.1189
ALPHA=-0.25*( ALPHAI(I,1)+ALPHAI(IP1,1)+ OOM1F405.1190
& ALPHAIP(I,1)+ALPHAIP(IP1,1) ) OOM1F405.1191
BETA=0.25*( BETAI(I,1)+BETAI(IP1,1)+ OOM1F405.1192
& BETAIP(I,1)+BETAIP(IP1,1) ) OOM1F405.1193
C FACTORS CONVERT TO CGS UNITS. W/M^2=KG/S^3=1000G/S^3 OOM1F405.1194
C KG/M^2/S=0.1G/CM^2/S OOM1F405.1195
C OOM1F405.1196
IF(L_SEAICE)THEN OOM1F405.1197
SFC_HEATFLUX_MOM=1000.0*0.25*( OOM1F405.1198
& OCEANHEATFLUX(I)+OCEANHEATFLUX(IP1)+ OOM1F405.1199
& OCEANHEATFLUXP(I)+OCEANHEATFLUXP(IP1)- OOM1F405.1200
& FLXTOICE(I)-FLXTOICE(IP1)- OOM1F405.1201
& FLXTOICEP(I)-FLXTOICEP(IP1)+ OOM1F405.1202
& CARYHEAT(I)+CARYHEAT(IP1)+ OOM1F405.1203
& CARYHEATP(I)+CARYHEATP(IP1) ) OOM1F405.1204
ELSE OOM1F405.1205
SFC_HEATFLUX_MOM=1000.0*0.25*( OOM1F405.1206
& HTN(I)+HTN(IP1)+ OOM1F405.1207
& HTNP(I)+HTNP(IP1) ) OOM1F405.1208
ENDIF OOM1F405.1209
C OOM1F405.1210
SOL_MOM=0.25*(SOL(I)+SOL(IP1)+SOLP(I)+SOLP(IP1))*1000. OOM1F405.1211
WATERFLUX_MOM=0.25*(PME(I)+WATERFLUX_ICE(I) OOM1F405.1212
& +PME(IP1)+WATERFLUX_ICE(IP1) OOM1F405.1213
& +PMEP(I)+WATERFLUX_ICEP(I) OOM1F405.1214
& +PMEP(IP1)+WATERFLUX_ICEP(IP1))*0.1 OOM1F405.1215
C ******************************************************************* OOM1F405.1216
BF=GRAV * ( ALPHA*SFC_HEATFLUX_MOM/(RHO0*CP0) OOM1F405.1217
& + BETA*WATERFLUX_MOM*0.035/RHO_FRESHWATER ) OOM1F405.1218
& + GRAV*ALPHA* ( SOL_MOM/(RHO0*CP0) OOM1F405.1219
& - SOL_IRRA_1B(SOL_MOM,HM(I))/(RHO0*CP0) ) OOM1F405.1220
C ******************************************************************* OOM1F405.1221
C CALCULATE MONIN-OBUKHOV LENGTH, L_M, CM = (CM/S)^3 / (CM^2/S^3) OOM1F405.1222
IF (ABS(BF).LT.1.0E-12) THEN OOM1F405.1223
IF (BF.GE.0.0) BF=1.0E-12 OOM1F405.1224
IF (BF.LT.0.0) BF=-1.0E-12 OOM1F405.1225
ENDIF OOM1F405.1226
L_M(I) = USTAR(I)**3/(VKMAN*BF) OOM1F405.1227
IF (ABS(L_M(I)).LT.1.0E-12) THEN OOM1F405.1228
IF (L_M(I).GE.0.0) L_M(I)=1.0E-12 OOM1F405.1229
IF (L_M(I).LT.0.0) L_M(I)=-1.0E-12 OOM1F405.1230
ENDIF OOM1F405.1231
ENDIF OOM1F405.1232
ENDDO OOM1F405.1233
C CALCULATE DEPTH TO APPLY LARGE SCHEME TO (100 CONVERTS TO CM) OOM1F405.1234
DO I=1,IMT-1 OOM1F405.1235
HM(I)=0.25*(MLD_LARGE(I)+MLD_LARGE(I+1)+ OOM1F405.1236
& MLD_LARGEP(I)+MLD_LARGEP(I+1)) OOM1F405.1237
C OOM1F405.1238
IF (L_OWINDMIX.AND.L_M(I).GT.0.0) THEN OOM1F405.1239
HM(I)=MAX(HM(I),2*LAMBDA_LARGE*VKMAN*L_M(I)) OOM1F405.1240
ENDIF OOM1F405.1241
C OOM1F405.1242
IF (L_OBULKMAXMLD.AND.L_M(I).GT.0.0) THEN OOM1F405.1243
HM(I)=MIN(HM(I),L_M(I)) OOM1F405.1244
IF (.NOT.L_OROTATE) THEN OOM1F405.1245
C CALCULATE CORIOLIS PARAMETER F. IF PHI<5 DEG, ASSUME PHI=5 OOM1F405.1246
FX=2*OMEGA*SIN(MAX(PHI,5.0*PI_OVER_180)) OOM1F405.1247
HM(I)=MIN(HM(I),0.7*USTAR(I)/FX) OOM1F405.1248
ENDIF OOM1F405.1249
ENDIF OOM1F405.1250
C OOM1F405.1251
C SOME FINAL CHECKS ON HM OOM1F405.1252
C OOM1F405.1253
HM(I)=MIN(HM(I),MAX_LARGE_DEPTH*100.) OOM1F405.1254
HM(I)=MAX(HM(I),1.5*ZDZ(1)) OOM1F405.1255
C OOM1F405.1256
ENDDO OOM1F405.1257
IF (L_OCYCLIC) THEN OOM1F405.1258
HM(IMT)=HM(2) OOM1F405.1259
HM(IMT-1)=HM(1) OOM1F405.1260
ELSE OOM1F405.1261
HM(IMT)=0.5*(MLD_LARGE(IMT)+MLD_LARGEP(IMT)) OOM1F405.1262
C OOM1F405.1263
IF (L_OWINDMIX.AND.L_M(IMT).GT.0.0) THEN OOM1F405.1264
HM(IMT)=MAX(HM(IMT),2*LAMBDA_LARGE*VKMAN*L_M(IMT)) OOM1F405.1265
ENDIF OOM1F405.1266
C OOM1F405.1267
IF (L_OBULKMAXMLD.AND.L_M(IMT).GT.0.0) THEN OOM1F405.1268
HM(IMT)=MIN(HM(IMT),L_M(IMT)) OOM1F405.1269
IF (.NOT.L_OROTATE) THEN OOM1F405.1270
C CALCULATE CORIOLIS PARAMETER F. IF PHI<5 DEG, ASSUME PHI=5 OOM1F405.1271
FX=2*OMEGA*SIN(MAX(PHI,5.0*PI_OVER_180)) OOM1F405.1272
HM(IMT)=MIN(HM(IMT),0.7*USTAR(IMT)/FX) OOM1F405.1273
ENDIF OOM1F405.1274
ENDIF OOM1F405.1275
C OOM1F405.1276
C SOME FINAL CHECKS ON HM OOM1F405.1277
C OOM1F405.1278
HM(IMT)=MAX(HM(IMT),1.5*ZDZ(1)) OOM1F405.1279
HM(IMT)=MIN(HM(IMT),MAX_LARGE_DEPTH*100.) OOM1F405.1280
C OOM1F405.1281
ENDIF OOM1F405.1282
C CALCULATE MIN NO OF LEVELS WITHIN WHICH MIXED LAYER IS CONTAINED OOM1F405.1283
C I.E. LEVEL WITHIN WHICH THE MIXED LAYER LIES OOM1F405.1284
DO I=1,IMT OOM1F405.1285
MLK(I)=1 OOM1F405.1286
DO K=MAX_LARGE_LEVELS,1,-1 OOM1F405.1287
IF (MLK(I).EQ.1.AND.ZDZ(K).LE.HM(I)) MLK(I)=K+1 OOM1F405.1288
ENDDO OOM1F405.1289
ENDDO OOM1F405.1290
ELSE OOM1F405.1291
C For quadratic Large scheme calculate surface friction velocity OLA3F403.68
C in cgs units. 10.0 converts from N m^-2 to dynes cm^-2 OLA3F403.69
C GM(I,1) is land-sea mask for velocity grid OLA3F403.70
C RHO0 = 1 g cm^-3 used in last line; tau = rho0 u*^2 OLA3F403.71
IF (L_OQLARGE) THEN OLA3F403.72
RHO0=RHO_WATER_SI*0.001 ! CONVERT TO CGS UNITS OOM1F405.1388
IF (L_OUSTARWME) THEN OOM1F405.1389
C 1000.0 CONVERTS FROM SI TO CGS OOM1F405.1390
C USTAR^3 = WME/RHO0 OOM1F405.1391
DO I=1,IMT-1 OOM1F405.1392
WME_MOM=1000.0*0.25*(WME(I)+WME(I+1)+WMEP(I)+WMEP(I+1)) OOM1F405.1393
IF (GM(I,1).GT.0.5) THEN OOM1F405.1394
USTAR(I)=(WME_MOM/RHO0)**(1.0/3.0) OOM1F405.1395
ELSE OOM1F405.1396
USTAR(I)=0.0 OOM1F405.1397
ENDIF OOM1F405.1398
END DO OOM1F405.1399
IF (L_OCYCLIC) THEN OOM1F405.1400
USTAR(IMT)=USTAR(2) OOM1F405.1401
USTAR(IMT-1)=USTAR(1) OOM1F405.1402
ELSE OOM1F405.1403
WME_MOM=1000.0*0.5*(WME(IMT)+WMEP(IMT)) OOM1F405.1404
IF (GM(IMT,1).GT.0.5) THEN OOM1F405.1405
USTAR(IMT)=(WME_MOM/RHO0)**(1.0/3.0) OOM1F405.1406
ELSE OOM1F405.1407
USTAR(IMT)=0.0 OOM1F405.1408
ENDIF OOM1F405.1409
ENDIF OOM1F405.1410
ELSE OOM1F405.1411
do i=1,imt OLA3F403.73
ustar(i)= sqrt(WSX(i)*WSX(i)+WSY(i)*WSY(i)) OLA3F403.74
ustar(i)= 10.0 * GM(I,1) * ustar(i) OLA3F403.75
ustar(i)= sqrt(ustar(i)) OLA3F403.76
end do OLA3F403.77
ENDIF OOM1F405.1412
C OOM1F405.1413
C Initialise hm to zero OLA3F403.79
do i=1,imt OLA0F404.213
hm(i)=0.0 OLA3F403.81
enddo OLA3F403.82
ENDIF ! L_OQLARGE OLA0F404.214
C Initialise Rim to zero OLA3F403.83
do k=1,kmm1 OLA3F403.84
do i=1,imt OLA3F403.85
Rim(i,k)=0.0 OLA3F403.86
enddo OLA3F403.87
enddo OLA3F403.88
ENDIF OOM1F405.1292
C VERTCOFC.66
C -------------------------------------------------------------- VERTCOFC.67
C Calculate density on TS rows to north and south of UV row VERTCOFC.68
C -------------------------------------------------------------- VERTCOFC.69
C VERTCOFC.70
C OOM1F405.987
IF(L_OSTATEC)THEN OOM1F405.988
C OOM1F405.989
C--------------- USE STATEC TO FIND VALUES FOR RHO OOM1F405.990
CALL STATEC(TB(1,1,1),TB(1,1,2),RHOSRNA,WORKA,WORKB,1, OOM1F405.991
& IMT,KM,J,JMT) OOM1F405.992
CALL STATEC(TB(1,1,1),TB(1,1,2),RHOSRNB,WORKA,WORKB,2, OOM1F405.993
& IMT,KM,J,JMT) OOM1F405.994
CALL STATEC(TBP(1,1,1),TBP(1,1,2),RHONRNA,WORKA,WORKB,1, OOM1F405.995
& IMT,KM,J+1,JMT) OOM1F405.996
CALL STATEC(TBP(1,1,1),TBP(1,1,2),RHONRNB,WORKA,WORKB,2, OOM1F405.997
& IMT,KM,J+1,JMT) OOM1F405.998
C OOM1F405.999
DO I=1,IMT OOM1F405.1000
RHOSRNA(I,KM+1)=RHOSRNA(I,KM) OOM1F405.1001
RHOSRNB(I,KM+1)=RHOSRNB(I,KM) OOM1F405.1002
RHONRNA(I,KM+1)=RHONRNA(I,KM) OOM1F405.1003
RHONRNB(I,KM+1)=RHONRNB(I,KM) OOM1F405.1004
ENDDO OOM1F405.1005
C OOM1F405.1006
ELSE OOM1F405.1007
C OOM1F405.1008
C--------------- USE STATED TO FIND VALUES FOR RHO OOM1F405.1009
JA121293.275
CALL STATED(TB(1,1,1),TB(1,1,2),rhosrn,worka,workb,IMT,KM,J JA121293.276
&,KM JA121293.277
& , JMT ORH7F404.46
&) JA121293.280
CALL STATED(TBP(1,1,1),TBP(1,1,2),rhonrn,worka,workb,IMT,KM,J+1 JA121293.281
&,KM JA121293.282
& , JMT ORH7F404.47
&) JA121293.285
JA121293.286
C OOM1F405.1010
ENDIF ! L_OSTATEC OOM1F405.1011
C OOM1F405.1012
C VERTCOFC.73
C -------------------------------------------------------------- VERTCOFC.74
C Evaluate Richardson No. VERTCOFC.75
C and thereby coefficient of viscosity for mid-levels VERTCOFC.76
C VERTCOFC.77
C Note Ri = (g/rho0)*drhodz/(dudz*dudz + dvdz*dvdz). VERTCOFC.78
C This code assumes rho0=1 (OK in cgs units!) VERTCOFC.79
C -------------------------------------------------------------- VERTCOFC.80
C VERTCOFC.81
fxa=0.5 VERTCOFC.82
fxb=2. VERTCOFC.83
C OOM1F405.1013
IF (L_OSTATEC) THEN OOM1F405.1014
C---------- AT PRESENT HARDWIRE CODING FOR KM EVEN AND STATEC OOM1F405.1015
DO K=1,KM-3,2 OOM1F405.1016
DO I=1,IMT-1 OOM1F405.1017
RHOZA=((RHONRNA(I,K)-RHONRNA(I,K+1)) OOM1F405.1018
& +(RHONRNA(I+1,K)-RHONRNA(I+1,K+1)) OOM1F405.1019
& +( RHOSRNA(I,K)-RHOSRNA(I,K+1)) OOM1F405.1020
& +(RHOSRNA(I+1,K)-RHOSRNA(I+1,K+1))) OOM1F405.1021
& *FXA*DZZ2RQ(I,K+1) OOM1F405.1022
RHOZB=((RHONRNB(I,K+1)-RHONRNB(I,K+2)) OOM1F405.1023
& +(RHONRNB(I+1,K+1)-RHONRNB(I+1,K+2)) OOM1F405.1024
& +( RHOSRNB(I,K+1)-RHOSRNB(I,K+2)) OOM1F405.1025
& +(RHOSRNB(I+1,K+1)-RHOSRNB(I+1,K+2))) OOM1F405.1026
& *FXA*DZZ2RQ(I,K+2) OOM1F405.1027
UZ=(UB(I,K)-UB(I,K+1))*FXB*DZZ2RQ(I,K+1) OOM1F405.1028
VZ=(VB(I,K)-VB(I,K+1))*FXB*DZZ2RQ(I,K+1) OOM1F405.1029
DVEL2Q(I,K)=UZ*UZ+VZ*VZ OOM1F405.1030
RIPQ(I,K)=-GRAV_SI*100.*RHOZA ! 100. CONVERTS GRAV_SI TO CGS OOM1F405.1031
UZ=(UB(I,K+1)-UB(I,K+2))*FXB*DZZ2RQ(I,K+2) OOM1F405.1032
VZ=(VB(I,K+1)-VB(I,K+2))*FXB*DZZ2RQ(I,K+2) OOM1F405.1033
DVEL2Q(I,K+1)=UZ*UZ+VZ*VZ OOM1F405.1034
RIPQ(I,K+1)=-GRAV_SI*100.*RHOZB ! 100. CONVERTS GRAV_SI TO CGS OOM1F405.1035
ENDDO OOM1F405.1036
ENDDO OOM1F405.1037
I=IMT OOM1F405.1038
DO K=1,KM-3,2 OOM1F405.1039
RHOZA=((RHONRNA(I,K)-RHONRNA(I,K+1)) OOM1F405.1040
& +( RHOSRNA(I,K)-RHOSRNA(I,K+1))) OOM1F405.1041
& *DZZ2RQ(I,K+1) OOM1F405.1042
RHOZB=((RHONRNB(I,K+1)-RHONRNB(I,K)) OOM1F405.1043
& +( RHOSRNB(I,K+1)-RHOSRNB(I,K))) OOM1F405.1044
& *DZZ2RQ(I,K+2) OOM1F405.1045
UZ=(UB(I,K)-UB(I,K+1))*FXB*DZZ2RQ(I,K+1) OOM1F405.1046
VZ=(VB(I,K)-VB(I,K+1))*FXB*DZZ2RQ(I,K+1) OOM1F405.1047
DVEL2Q(I,K)=UZ*UZ+VZ*VZ OOM1F405.1048
RIPQ(I,K)=-GRAV_SI*100.*RHOZA ! 100. CONVERTS GRAV_SI TO CGS OOM1F405.1049
UZ=(UB(I,K+1)-UB(I,K+2))*FXB*DZZ2RQ(I,K+2) OOM1F405.1050
VZ=(VB(I,K+1)-VB(I,K+2))*FXB*DZZ2RQ(I,K+2) OOM1F405.1051
DVEL2Q(I,K+1)=UZ*UZ+VZ*VZ OOM1F405.1052
RIPQ(I,K+1)=-GRAV_SI*100.*RHOZB ! 100. CONVERTS GRAV_SI TO CGS OOM1F405.1053
ENDDO OOM1F405.1054
C OOM1F405.1055
K=KMM1 OOM1F405.1056
DO I=1,IMT-1 OOM1F405.1057
RHOZA=((RHONRNA(I,K)-RHONRNA(I,K+1)) OOM1F405.1058
& +(RHONRNA(I+1,K)-RHONRNA(I+1,K+1)) OOM1F405.1059
& +( RHOSRNA(I,K)-RHOSRNA(I,K+1)) OOM1F405.1060
& +(RHOSRNA(I+1,K)-RHOSRNA(I+1,K+1))) OOM1F405.1061
& *FXA*DZZ2RQ(I,K+1) OOM1F405.1062
UZ=(UB(I,K)-UB(I,K+1))*FXB*DZZ2RQ(I,K+1) OOM1F405.1063
VZ=(VB(I,K)-VB(I,K+1))*FXB*DZZ2RQ(I,K+1) OOM1F405.1064
DVEL2Q(I,K)=UZ*UZ+VZ*VZ OOM1F405.1065
RIPQ(I,K)=-GRAV_SI*100.*RHOZA ! 100. CONVERTS GRAV_SI TO CGS OOM1F405.1066
ENDDO OOM1F405.1067
I=IMT OOM1F405.1068
RHOZA=((RHONRNA(I,K)-RHONRNA(I,K+1)) OOM1F405.1069
& +( RHOSRNA(I,K)-RHOSRNA(I,K+1))) OOM1F405.1070
& *DZZ2RQ(I,K+1) OOM1F405.1071
UZ=(UB(I,K)-UB(I,K+1))*FXB*DZZ2RQ(I,K+1) OOM1F405.1072
VZ=(VB(I,K)-VB(I,K+1))*FXB*DZZ2RQ(I,K+1) OOM1F405.1073
DVEL2Q(I,K)=UZ*UZ+VZ*VZ OOM1F405.1074
RIPQ(I,K)=-GRAV_SI*100.*RHOZA ! 100. CONVERTS GRAV_SI TO CGS OOM1F405.1075
ELSE OOM1F405.1076
DO 120 K=1,KMM1 VERTCOFC.84
DO 120 I=1,IMT VERTCOFC.85
IF (I.NE.IMT) THEN VERTCOFC.86
rhoz=((rhonrn(I,K)-rhonrn(I,K+1)) VERTCOFC.87
& +(rhonrn(I+1,K)-rhonrn(I+1,K+1)) VERTCOFC.88
& +( rhosrn(I,K)-rhosrn(I,K+1)) VERTCOFC.89
& +(rhosrn(I+1,K)-rhosrn(I+1,K+1))) VERTCOFC.90
& *fxa*DZZ2RQ(I,K+1) VERTCOFC.91
ELSE VERTCOFC.92
rhoz=((rhonrn(I,K)-rhonrn(I,K+1)) VERTCOFC.93
& +( rhosrn(I,K)-rhosrn(I,K+1))) VERTCOFC.94
& *DZZ2RQ(I,K+1) VERTCOFC.95
ENDIF VERTCOFC.96
uz=(UB(I,K)-UB(I,K+1))*fxb*DZZ2RQ(I,K+1) VERTCOFC.97
vz=(VB(I,K)-VB(I,K+1))*fxb*DZZ2RQ(I,K+1) VERTCOFC.98
dvel2q(I,K)=uz*uz+vz*vz VERTCOFC.99
ripq(I,K)=-GRAV_SI*100.*rhoz ! 100. converts GRAV_SI to cgs VERTCOFC.100
120 CONTINUE VERTCOFC.101
ENDIF ! L_OSTATEC OOM1F405.1077
C OOM1F405.1078
C VERTCOFC.102
C Convective adjustment takes care of unstable profiles VERTCOFC.103
C VERTCOFC.104
DO 150 K=1,KMM1 VERTCOFC.105
DO 150 I=1,IMT VERTCOFC.106
IF (ripq(I,K).LT.fx) ripq(I,K)=fx VERTCOFC.107
150 CONTINUE VERTCOFC.108
C VERTCOFC.109
C If currents are zero set dvel2q to a small number. This is done VERTCOFC.110
C to prevent a divide check in the calculation of gnu VERTCOFC.111
C VERTCOFC.112
fx=1.E-12 VERTCOFC.113
DO 160 K=1,KMM1 VERTCOFC.114
DO 160 I=1,IMT VERTCOFC.115
IF (dvel2q(I,K).LT.fx) dvel2q(I,K)=fx VERTCOFC.116
160 CONTINUE VERTCOFC.117
C VERTCOFC.118
C OOM1F405.1079
IF(L_OPANDP)THEN OOM1F405.1080
C OOM1F405.1081
C Now calculate viscosity. 1.E4 converts external consts. to cgs VERTCOFC.119
C Note k's are offset by 1 since gnu(*,k) is the viscostiy at the VERTCOFC.120
C TOP of box k. gnu(*,1) is set at the end of the routine. VERTCOFC.121
C VERTCOFC.122
fx=1. VERTCOFC.123
fxa=5. VERTCOFC.124
DO 170 K=1,KMM1 VERTCOFC.125
DO 170 I=1,IMT VERTCOFC.126
IF (GM(I,K+1).GT.0.9) THEN VERTCOFC.127
gnu(I,K+1)=FNU0_SI*1.E4*dvel2q(I,K)*dvel2q(I,K)/ VERTCOFC.128
& ((dvel2q(I,K)+fxa*ripq(I,K))* VERTCOFC.129
& (dvel2q(I,K)+fxa*ripq(I,K)))+FNUB_SI*1.E4 VERTCOFC.130
ENDIF VERTCOFC.131
170 CONTINUE VERTCOFC.132
C OOM1F405.1082
ENDIF OOM1F405.1083
C OOM1F405.1084
OLA3F403.89
C Calculate Richardson number OLA3F403.90
do K=1,KMM1 OLA3F403.91
do I=1,IMT OLA3F403.92
Rim(I,K)=ripq(I,K)/dvel2q(I,K) OLA3F403.93
end do OLA3F403.94
end do OLA3F403.95
C OOM1F405.1085
IF(.NOT.L_OPANDP)THEN OOM1F405.1086
C OOM1F405.1087
C CALCULATE GNU - MOMENTUM DIFFUSIVITY USING PETERS SCHEME. OOM1F405.1088
C Modify to limit dvel2q setting gnu if below a certain level. OOM1F405.1089
C If dvel2q too small, then regardless of ripq, set gnu to its OOM1F405.1090
C background value there. OOM1F405.1091
C GJR July 1998 OOM1F405.1092
C OOM1F405.1093
DO K=1,KMM1 OOM1F405.1094
DO I=1,IMT OOM1F405.1095
IF (GM(I,K+1).GT.0.9) THEN OOM1F405.1096
IF (DVEL2Q(I,K).GT.1.E-6)THEN OOM1F405.1097
IF (RIM(I,K).LT.RICRITM100) OOM1F405.1098
& GNU(I,K+1)=100. OOM1F405.1099
& +FNUB_SI*1.E4 OOM1F405.1100
IF (RIM(I,K).LT.RICRITM.AND.RIM(I,K).GT.RICRITM100) OOM1F405.1101
& GNU(I,K+1)=5.6E-4/RIM(I,K)**8.2 OOM1F405.1102
& +FNUB_SI*1.E4 OOM1F405.1103
IF (RIM(I,K).GE.RICRITM) OOM1F405.1104
& GNU(I,K+1)=5.0/((1+5*RIM(I,K))**1.5) OOM1F405.1105
& +FNUB_SI*1.E4 OOM1F405.1106
ELSE OOM1F405.1107
GNU(I,K+1)=FNUB_SI*1.E4 OOM1F405.1108
ENDIF OOM1F405.1109
ENDIF OOM1F405.1110
ENDDO OOM1F405.1111
ENDDO OOM1F405.1112
C OOM1F405.1113
ENDIF OOM1F405.1114
C OOM1F405.1115
OLA3F403.96
IF (L_OFULARGE) THEN OOM1F405.1293
C CALCULATE VERTICAL DIFFUSIVITY, GNU, AT HM OOM1F405.1294
DO I=1,IMT OOM1F405.1295
IF (GM(I,1).GT.0.5.AND.HM(I).GE.ZDZ(1)) THEN OOM1F405.1296
ZDZ_AT_MLKM1=0.0 OOM1F405.1297
GNUOFMLK=0.0 OOM1F405.1298
IF (MLK(I).GT.1) THEN OOM1F405.1299
ZDZ_AT_MLKM1=ZDZ(MLK(I)-1) OOM1F405.1300
GNUOFMLK=GNU(I,MLK(I)) OOM1F405.1301
ENDIF OOM1F405.1302
COEFF=(HM(I)-ZDZ_AT_MLKM1)/DZ(MLK(I)) OOM1F405.1303
GNUATHM=(1-COEFF)*GNUOFMLK+COEFF*GNU(I,MLK(I)+1) OOM1F405.1304
C CALCULATE D(GNU)/DZ AT HM AND RESTRICT TO MINIMUM OF 0.0 AS OTHERWISE OOM1F405.1305
C MAY OBTAIN POSSIBLE NEGATIVE VALUES OF GNU FROM LARGE SCHEME DUE OOM1F405.1306
C TO FITTING OF CUBIC FUNCTION. OOM1F405.1307
GNUZ(1)=(GNUOFMLK-GNU(I,MLK(I)+1))/DZ(MLK(I)) OOM1F405.1308
IF (GNUZ(1).GT.0.0) THEN OOM1F405.1309
IF (COEFF.LE.0.5) GNUZATHM=GNUZ(1) OOM1F405.1310
IF (COEFF.GT.0.5) THEN OOM1F405.1311
GNUZ(2)=(GNU(I,MLK(I)+1)-GNU(I,MLK(I)+2))/DZ(MLK(I)+1) OOM1F405.1312
COEFF=(HM(I)-ZDZZ(MLK(I)))/DZZ(MLK(I)+1) OOM1F405.1313
GNUZATHM=(1-COEFF)*GNUZ(1)+COEFF*GNUZ(2) OOM1F405.1314
GNUZATHM=MAX(GNUZATHM,0.0) OOM1F405.1315
ENDIF OOM1F405.1316
ELSE OOM1F405.1317
GNUZATHM=0.0 OOM1F405.1318
ENDIF OOM1F405.1319
C CALCULATE DEPTH DEPENDENT TURBULENT VELOCITY SCALE AT HM, OOM1F405.1320
C W1M, (LARGE EQN 13 WITH SIGMA=1) AND ITS DERIVATIVE WRT OOM1F405.1321
C SIGMA, DW1M. OOM1F405.1322
C OOM1F405.1323
IF (L_M(I).LT.0.0)THEN OOM1F405.1324
C OOM1F405.1325
STABIL_PARAM=EPSILON*HM(I)/L_M(I) OOM1F405.1326
W1M=VKMAN*USTAR(I)/FLUX_PROFILEM(STABIL_PARAM) OOM1F405.1327
DW1M=0.0 OOM1F405.1328
C OOM1F405.1329
ELSE OOM1F405.1330
STABIL_PARAM=HM(I)/L_M(I) OOM1F405.1331
W1M=VKMAN*USTAR(I)/FLUX_PROFILEM(STABIL_PARAM) OOM1F405.1332
IF (0.0.LE.STABIL_PARAM) THEN OOM1F405.1333
DW1M = -5.0*VKMAN*USTAR(I)*STABIL_PARAM OOM1F405.1334
& /((1.0+5.0*STABIL_PARAM)**2) OOM1F405.1335
ELSE OOM1F405.1336
IF (-0.2.LE.STABIL_PARAM) THEN OOM1F405.1337
DW1M = -4.0*VKMAN*USTAR(I)*STABIL_PARAM OOM1F405.1338
& /((1.0-16.0*STABIL_PARAM)**0.75) OOM1F405.1339
ELSE OOM1F405.1340
DW1M = - 8.38/3.0 *VKMAN*USTAR(I)*STABIL_PARAM OOM1F405.1341
& /((1.26-8.38*STABIL_PARAM)**(2.0/3.0)) OOM1F405.1342
ENDIF OOM1F405.1343
ENDIF OOM1F405.1344
C OOM1F405.1345
ENDIF ! L_M(I).LT.0.0 OOM1F405.1346
C OOM1F405.1347
C CALCULATE, WHEN SIGMA=Z/H=1, THE SHAPE FUNCTION AND ITS DERIVATIVE OOM1F405.1348
C WRT SIGMA (LARGE EQN 18) OOM1F405.1349
C N.B. G1=G(1)=SHAPE FUNCTION WHEN SIGMA=1 OOM1F405.1350
C DG1= DG(1)/DSIGMA OOM1F405.1351
IF (ABS(HM(I)).LT.1.0E-12.OR.ABS(W1M).LT.1.0E-12) THEN OOM1F405.1352
IF (W1M.GE.0.) W1M=1.0E-12 OOM1F405.1353
IF (W1M.LT.0.) W1M=-1.0E-12 OOM1F405.1354
ENDIF OOM1F405.1355
G1 = GNUATHM/(HM(I)*W1M) OOM1F405.1356
DG1 = - GNUZATHM/W1M - (GNUATHM*DW1M)/(HM(I)*W1M*W1M) OOM1F405.1357
C CALCULATE COEFFICIENTS OF SHAPE FUNCTION (LARGE EQN 17) OOM1F405.1358
A(1)=1. OOM1F405.1359
A(2)=-2.0+3.0*G1-DG1 OOM1F405.1360
A(3)=1.0-2.0*G1+DG1 OOM1F405.1361
C CALCULATE WM(K) - DEPTH DEPENDENT TURBULENT VELOCITY SCALE PROFILE OOM1F405.1362
C AT TOP OF LEVEL K (LARGE, EQN 13) OOM1F405.1363
DO K=2,MLK(I) OOM1F405.1364
STABIL_PARAM=ZDZ(K-1)/L_M(I) OOM1F405.1365
SIGMA=ZDZ(K-1)/HM(I) OOM1F405.1366
IF (STABIL_PARAM.LT.0.0.AND.SIGMA.GT.EPSILON)THEN OOM1F405.1367
STABIL_PARAM=EPSILON*HM(I)/L_M(I) OOM1F405.1368
ENDIF OOM1F405.1369
WM(K)=VKMAN*USTAR(I)/FLUX_PROFILEM(STABIL_PARAM) OOM1F405.1370
ENDDO OOM1F405.1371
C CALCULATE SHAPE FUNCTION G THEN VERTICAL DIFFUSION COEFF PROFILE, OOM1F405.1372
C GNU(K) AT TOP OF LEVEL K IN THE MIXED LAYER (LARGE EQNS 11 AND 10) OOM1F405.1373
GNU(I,1)=0. OOM1F405.1374
DO K=2,MLK(I) OOM1F405.1375
SIGMA=ZDZ(K-1)/HM(I) OOM1F405.1376
G=SIGMA*(A(1)+SIGMA*(A(2)+A(3)*SIGMA)) OOM1F405.1377
GNU(I,K)=HM(I)*WM(K)*G OOM1F405.1378
ENDDO OOM1F405.1379
ELSE OOM1F405.1380
DO K=1,KM OOM1F405.1381
GNU(I,K)=0.0 OOM1F405.1382
ENDDO OOM1F405.1383
ENDIF OOM1F405.1384
ENDDO OOM1F405.1385
ELSE OOM1F405.1386
IF (L_OQLARGE) THEN OLA3F403.97
C Simple quadratic form of Large scheme - OLA3F403.98
C applies - neutral Large in 'mixed layer' - defined to be where OLA3F403.99
C Richardson number is less than 0.3. OLA3F403.100
C - normal model Packonowski and Philander below mixed layer OLA3F403.101
C - functions specified so that vertical diffusion coeff is OLA3F403.102
C continuous at z=h OLA3F403.103
C - Constants chosen to give a quadratic form for G(z/h) OLA3F403.104
C Reference: OLA3F403.105
C W.G.Large et al 1994, Oceanic Vertical Mixing : A review and a OLA3F403.106
C model with a nonlocal boundary layer parametrisation, Rev Geophys, OLA3F403.107
C 32, 363-403. OLA3F403.108
C This scheme is a simplified neutral case of Large OLA3F403.109
C Vertical diffusion coeff, OLA3F403.110
C K(z)=h*kk*ustar*G(z/h) OLA3F403.111
C where z=depth (positive downwards) OLA3F403.112
C h=mixed layer depth OLA3F403.113
C kk=0.4 , von Karman's constant OLA0F404.215
C ustar= surface friction velocity OLA3F403.115
C = sqrt(wsx**2+wsy**2)/rho_w OLA3F403.116
C G(z/h)=z/h+a2*(z/h)**2 OLA3F403.117
C where a2 = Kath/(kk*ustar*h) - 1 OLA3F403.118
C Kath=K(Ri) at z=h OLA3F403.119
C To avoid problems where ustar is zero this is actually implemented as OLA3F403.120
C K(z)=h*kk*(z/h)*ustar*b2*(z/h)**2 OLA3F403.121
C where b2 = Kath/(kk*h) - ustar OLA3F403.122
C Scheme - for both momentum and tracers OLA3F403.123
C 1. Calculate gnu as usual OLA3F403.124
C 2. Find 'mixed layer depth' for Large scheme, h OLA3F403.125
C 3. For depths shallower than h, recalculate gnu using Large. OLA3F403.126
C OLA3F403.127
do i=1,imt OLA0F404.216
C Find 'mixed layer depth', hm, for quadratic Large scheme defined as OLA0F404.217
C min(depth where the critical Richardson no reached,max_qLarge_depth) OLA0F404.218
IF (GM(I,1).GT.0.9) THEN OLA0F404.219
do k=1,km OLA0F404.220
l=k OLA0F404.221
if (Rim(I,K).ge.crit_Ri.or.zdz(k).gt.max_qLarge_depth*100. OLA0F404.222
& .or.GM(i,min(k+1,km)).lt.0.5) goto 10 OLA0F404.223
enddo OLA0F404.224
10 continue OLA0F404.225
if (l.gt.1) then OLA0F404.226
if (Rim(i,l).ge.crit_Ri) then OLA0F404.227
if (Rim(i,l).ne.Rim(i,l-1)) then OLA0F404.228
hm(i)=zdz(l-1) + dz(l) OLA0F404.229
& *(crit_Ri-Rim(i,l-1))/(Rim(i,l)-Rim(i,l-1)) OLA0F404.230
hm(i)=min(hm(i),max_qLarge_depth*100.) OLA0F404.231
else OLA0F404.232
hm(i)=zdz(l-1) OLA0F404.233
endif OLA0F404.234
else OLA0F404.235
hm(i)=min(zdz(l),max_qLarge_depth*100.) OLA0F404.236
endif OLA0F404.237
C Calculate altered gnu i.e. applying quadratic Large where depth < hm OLA0F404.238
Kath=gnu(i,l-1) + (gnu(i,l)-gnu(i,l-1)) OLA0F404.239
& *(hm(i)-zdz(l-1))/dz(l) OLA0F404.240
b2 = Kath/(vkman*hm(i)) - ustar(i) OLA0F404.241
do k=1,l-1 OLA0F404.242
sigma=zdz(k)/hm(i) OLA0F404.243
gnu(i,k+1)=hm(i)*vkman*sigma*(ustar(i)+b2*sigma) OLA0F404.244
end do OLA0F404.245
endif OLA0F404.246
ENDIF OLA0F404.247
end do OLA3F403.152
ELSE OLA3F403.153
C The following occurs if the quadratic Large scheme is not chosen. OLA3F403.154
C VERTCOFC.133
C Reset gnu if value is greater than the stability limit. VERTCOFC.134
C 1.E4 converts STABLM_SI to cgs VERTCOFC.135
C VERTCOFC.136
DO 180 K=2,KM VERTCOFC.137
DO 180 I=1,IMT VERTCOFC.138
IF (gnu(I,K).GT.STABLM_SI*1.E4) gnu(I,K)=STABLM_SI*1.E4 VERTCOFC.139
180 CONTINUE VERTCOFC.140
ENDIF OLA3F403.155
C VERTCOFC.141
C Reset mixing coeff. between first two levels to a min of GNUMINC_SI VERTCOFC.142
C The surface value gnu(*,1) is also set to GNUMINC_SI here. VERTCOFC.143
C VERTCOFC.144
fx=GNUMINC_SI*1.E4 ! 1.E4 converts to cgs VERTCOFC.145
DO 190 I=1,IMT VERTCOFC.146
IF (gnu(I,2).LT.fx) gnu(I,2)=fx*GM(I,2) VERTCOFC.147
gnu(I,1)=fx*GM(I,1) VERTCOFC.148
190 CONTINUE VERTCOFC.149
C VERTCOFC.150
ENDIF OOM1F405.1387
C Finally, multiply by the land-sea mask for safety. VERTCOFC.151
C VERTCOFC.152
DO 200 K=1,KM VERTCOFC.153
DO 210 I=1,IMT VERTCOFC.154
gnu(I,K)=gnu(I,K)*GM(I,K) VERTCOFC.155
210 CONTINUE VERTCOFC.156
200 CONTINUE VERTCOFC.157
C VERTCOFC.158
RETURN VERTCOFC.159
END VERTCOFC.160
*ENDIF @DYALLOC.4686