SUBROUTINE VERTCOFT 6,7VERTCOFT.15
& ( J,IMT,KM,KMM1,NT, VERTCOFT.16
& IMT_QLARGE,L_OQLARGE,L_OCYCLIC, OLA3F403.295
& IMT_FULARGE,L_OFULARGE,L_OUSTARWME,L_OWINDMIX,L_OBULKMAXMLD, OOM1F405.1423
& L_OROTATE,L_OSTATEC,L_SEAICE,L_OPANDP,PHIT,TB, OOM1F405.1424
& UB,VB,UBM,VBM, VERTCOFT.17
& WSX,WSY,WSXM,WSYM, ! surface stress fields OLA3F403.296
& ZDZZ,ZDZ,DZ,DZZ, OOM1F405.1425
& DZZ2RQ,DZ2RQ, OLA3F403.298
& NERGY,GRAV_SI,FM, VERTCOFT.19
& RHOSRN,RHOSRNA,RHOSRNB, OOM1F405.1419
& hT,RiT,max_qLarge_depth,crit_Ri, OLA0F404.153
& MLD_LARGE,MAX_LARGE_DEPTH,MAX_LARGE_LEVELS,RHO_WATER_SI, OOM1F405.1426
& HTN,PME,WATERFLUX_ICE,SOL,WME,L_T,LAMBDA_LARGE, OOM1F405.1427
& SPECIFIC_HEAT_SI, OOM1F405.1428
& gnu,FNU0_SI,FNUB_SI,STABLM_SI,FKAPB_SI,GNUMINT_SI RW071293.17
& ,KAPPA_B_SI,OCEANHEATFLUX,CARYHEAT,FLXTOICE,CGFLUX ) OOM1F405.1429
C VERTCOFT.22
C VERTCOFT.23
IMPLICIT NONE VERTCOFT.24
*CALL C_VKMAN 
OLA0F404.154
*CALL C_OMEGA OOM1F405.1431
*CALL C_PI OOM1F405.1432
C VERTCOFT.25
INTEGER VERTCOFT.26
& I,J,K,IMT,KM,KMM1,NT,ITT,NERGY,NDUM OOM1F405.1418
C VERTCOFT.28
REAL VERTCOFT.29
& UB(IMT,KM),VB(IMT,KM),UBM(IMT,KM),VBM(IMT,KM), VERTCOFT.30
& FM(IMT,KM),DZZ2RQ(IMT,KM),DZ2RQ(IMT,KM),GRAV_SI, OLA3F403.300
& rhosrn(IMT,KM), ! IN Density on row J. Calculated using VERTCOFT.32
& RHOSRNA(IMT,KM+1),RHOSRNB(IMT,KM+1), OOM1F405.1420
C ! IN DENSITY ON ROW J. CALCULATED USING STATEC IN VERTCOFC OOM1F405.1421
C ! IN CLINIC AND PASSED THROUGH ARGUMENT LISTS. OOM1F405.1422
C STATED in VERTCOFC in CLINIC and VERTCOFT.33
C passed through argument lists. VERTCOFT.34
& gnu(IMT,KM), ! OUT Vertical diffusion coeff. (cm2/s) VERTCOFT.35
C at TOP of box k VERTCOFT.36
& FNU0_SI, ! IN Max value of stability-dependent term VERTCOFT.37
C of vert. momentum diffusivity (m2/s) VERTCOFT.38
& FNUB_SI, ! IN Background value of momentum VERTCOFT.39
C diffusivity (m2/s) VERTCOFT.40
& STABLM_SI, ! IN Max value of diffusivity (m2/s) VERTCOFT.41
& FKAPB_SI ! IN Background value of tracer VERTCOFT.42
C diffusivity (m2/s) VERTCOFT.43
& ,GNUMINT_SI ! IN Min. value of tracer diffusivity VERTCOFT.44
C between top two levels (m2/s) VERTCOFT.45
& ,KAPPA_B_SI(KM) ! IN Vertically varying background tracer RW071293.19
C diffusivity at TOP of box K (m2/s) RW071293.20
REAL WSX(IMT),WSY(IMT) ! surface wind stress fields on row j OLA3F403.301
REAL WSXM(IMT),WSYM(IMT) ! surface wind stress on row j-1 OLA3F403.302
REAL ZDZZ(KM+1) ! IN DEPTH OF MIDDLE OF LAYER (CM) OOM1F405.1433
REAL ZDZ ( KM) ! depth of bottom of layer (cm) OLA3F403.304
REAL DZ(KM) ! IN distance between half-levels (cm) OLA0F404.155
INTEGER IMT_QLARGE ! imt for large scheme PXORDER.58
REAL OLA0F404.156
& hT(IMT_QLARGE) ! OUT depth quadratic Large scheme applied to OLA0F404.157
&, RiT(IMT,KMM1) ! OUT Richardson no at half levels (+ve) OLA0F404.158
&, max_qLarge_depth ! IN max depth allowed for quadratic Large OLA0F404.159
&, crit_Ri ! IN critical Richardson no used for depth of OLA0F404.160
! quadratic Large scheme OLA0F404.161
LOGICAL L_OQLARGE ! Logical for large scheme OLA3F403.306
LOGICAL L_OCYCLIC ! Logical for cyclic model EW OLA3F403.307
REAL OOM1F405.1434
& OCEANHEATFLUX(IMT) !HTN:NON-PEN SURF HEATFLUX INTO OCEAN BUDGET OOM1F405.1435
&,CARYHEAT(IMT) !MISCELLANEOUS HEATFLUX FROM ICE OOM1F405.1436
&,FLXTOICE(IMT) !OCEAN TO ICE HEATFLUX OOM1F405.1437
C VERTCOFT.46
REAL OOM1F405.1438
& MLD_LARGE(IMT) ! IN MIXED LAYER DEPTH (CM) OOM1F405.1439
&, WME(IMT) ! IN WIND MIXING ENERGY FOR ROW J OOM1F405.1440
&, DZZ(KM+1) ! IN DISTANCE BETWEEN MIDPTS OF LEVELS (CM) OOM1F405.1441
&, MAX_LARGE_DEPTH ! IN MAX DEPTH LARGE SCHEME APPLIED TO (M) OOM1F405.1442
&, RHO_WATER_SI ! IN DENSITY OF SEA WATER = 1026. KG/M^3 OOM1F405.1443
&, L_T(IMT) ! OUT MONIN OBUKHOV LENGTH LARGE SCHEME (MOMENTUM) OOM1F405.1444
&, ALPHAI(IMT,1),BETAI(IMT,1) ! EXPANSION COEFFICIENTS ON ROW J OOM1F405.1445
REAL OOM1F405.1446
& HTN(IMT) ! IN NON-PENETRATING HEAT FLUX (W/M2) ON ROW J OOM1F405.1447
&, PME(IMT) ! IN PRECIP MINUS EVAP (KG/M2/S) ON ROW J OOM1F405.1448
&, SOL(IMT) ! IN SOLAR IRRADIANCE (W/M2) AT SURFACE ON ROW J OOM1F405.1449
&, WATERFLUX_ICE(IMT) ! IN WATER FLUX FROM ICE (KG/M2/S) ON ROW J OOM1F405.1450
&, LAMBDA_LARGE,SPECIFIC_HEAT_SI ! IN FOR CALCULATING MINIMUM MLD OOM1F405.1451
&, PHIT ! IN LATITUDE OF ROW J ON TRACER GRID (RADIANS) OOM1F405.1452
&, TB(IMT,KM,NT) ! IN TEMPERATURE ARRAY FROM PREVIOUS TIMESTEP OOM1F405.1453
&, CGFLUX(IMT,KM,NT) !OUT COUNTER-GRADIENT FLUXES FOR TRACERS OOM1F405.1454
&, WT0(IMT),WS0(IMT),WTR0(IMT),WTRZ !LOCAL VARIABLES OOM1F405.1455
INTEGER OOM1F405.1456
& MAX_LARGE_LEVELS ! IN MAX NO OF LEVS FOR LARGE SCHEME OOM1F405.1457
&, IMT_FULARGE ! IN IMT FOR FULL LARGE SCHEME OOM1F405.1458
LOGICAL OOM1F405.1459
& L_OFULARGE ! IN LOGICAL FOR FULL LARGE SCHEME OOM1F405.1460
&, L_OUSTARWME,L_OWINDMIX,L_OBULKMAXMLD,L_OROTATE,L_OSTATEC OOM1F405.1461
LOGICAL L_OPANDP,L_SEAICE OOM1F405.1462
C VERTCOFT.47
C Decalre local variables and arrays VERTCOFT.48
C VERTCOFT.49
REAL VERTCOFT.50
& ripq(IMT,KMM1), ! d(rho)/dz at BOTTOM of box k VERTCOFT.51
& dvel2q(IMT,KMM1), ! Uz*Uz + Vz*Vz at BOTTOM of box k VERTCOFT.52
& uz, ! Uz at BOTTOM of current box VERTCOFT.53
& vz, ! Vz at BOTTOM of current box VERTCOFT.54
& wsxT,wsyT, ! wind stress on tracer grid OLA3F403.310
& ustar(IMT_QLARGE),! surface friction velocity on tracer grid OLA3F403.311
& Kath, ! K at z=h OLA3F403.312
& b2, ! used in calculation of gnu for Large scheme OLA3F403.313
& sigma, ! depth/h OLA3F403.314
& RHOZ,RHOZA,RHOZB, OOM1F405.1430
& fx, VERTCOFT.56
& fxa, VERTCOFT.57
& fxb VERTCOFT.58
C VERTCOFT.59
REAL OOM1F405.1463
& GNUATHT ! VERTICAL VISCOSITY (CM2/S) AT LEVEL HT OOM1F405.1464
&, GNUZATHT ! D(GNU)/DZV AT LEVEL HT OOM1F405.1465
&, ZDZ_AT_MLKM1 ! ZDZ AT LEVEL MLK-1 OOM1F405.1466
&, GNUOFMLK ! GNU(I,MLK(I)) OOM1F405.1467
&, GNUZ(2) ! D(GNU)/DZ AT MIDPOINTS OF LEVELS MLK,MLK+1 OOM1F405.1468
&, COEFF ! PROPORTION OF DEPTH BELOW TOP OF BOX FOR MLD OOM1F405.1469
&, RHO_FRESHWATER ! DENSITY OF FRESH WATER = 1.0 G/CM^3 OOM1F405.1470
&, RHO0 ! DENSITY OF SEA WATER = 1.026 G/CM^3 OOM1F405.1471
&, GRAV ! GRAV_SI IN CGS UNITS OOM1F405.1472
&, ALPHA,BETA ! THERMODYNAMIC EXPANSION COEFFS FOR T,S OOM1F405.1473
&, CP0 ! SPECIFIC HEAT CAPACITY AT CONST PRESSURE AT SURFACE(CGS) OOM1F405.1474
&, BF ! BUOYANCY FREQUENCY OOM1F405.1475
&, WATERFLUX,SFC_HEATFLUX ! TOTAL WATER FLUX (PME + WATERFLUX_ICE) OOM1F405.1476
REAL OOM1F405.1477
& W1T ! TURBULENT VELOCITY SCALE AT MLD OOM1F405.1478
&, DW1T ! DERIV OF TURBULENT VELOCITY SCALE WRT SIGMA AT MLD OOM1F405.1479
&, STABIL_PARAM ! STABILITY PARAMETER OOM1F405.1480
&, G1 ! SHAPE FUNCTION WHEN SIGMA=Z/H=1 OOM1F405.1481
&, DG1 ! DG/DSIGMA EVALUATED AT SIGMA=1 OOM1F405.1482
&, A(3) ! COEFFICIENTS OF CUBIC SHAPE FUNCTION OOM1F405.1483
&, WT(KM) ! TURBULENT VELOCITY SCALE PROFILE AT TOP OF BOX K OOM1F405.1484
&, EPSILON ! =0.1 (LARGE P365) OOM1F405.1485
&, G ! SHAPE FUNCTION AT TOP OF BOX K. OOM1F405.1486
INTEGER OOM1F405.1487
& MLK(IMT) ! LEVEL WITHIN WHICH MLD_LARGE IS CONTAINED OOM1F405.1488
INTEGER I2 OOM1F405.1489
C DEFINE FUNCTION OOM1F405.1490
REAL OOM1F405.1491
& SOL_IRRA_1B OOM1F405.1492
&, FLUX_PROFILET OOM1F405.1493
C OOM1F405.1494
C PARAMETERS FOR PETERS SCHEME OOM1F405.1495
C OOM1F405.1496
REAL RICRITT,RICRITT100 ! CRITICAL NUMBERS FOR TRACERS OOM1F405.1497
PARAMETER(RICRITT = 0.3752198666334,RICRITT100=0.2091865) OOM1F405.1498
C OOM1F405.1499
INTEGER l ! level within which hT is contained OLA0F404.162
! for quadratic Large scheme OLA0F404.163
CL Parameters OLA3F403.318
C VERTCOFT.60
C Zero out gnu VERTCOFT.61
C VERTCOFT.62
fx=0. VERTCOFT.63
DO 110 K=1,KM VERTCOFT.64
DO 110 I=1,IMT VERTCOFT.65
gnu(I,K)=fx VERTCOFT.66
110 CONTINUE VERTCOFT.67
C OOM1F405.1586
C Zero cgflux to be sure OOM1F405.1587
C OOM1F405.1588
DO NDUM=1,NT OOM1F405.1589
DO K=1,KM OOM1F405.1590
DO I=1,IMT OOM1F405.1591
CGFLUX(I,K,NDUM)=0.0 OOM1F405.1592
ENDDO OOM1F405.1593
ENDDO OOM1F405.1594
ENDDO OOM1F405.1595
C OOM1F405.1596
IF (L_OFULARGE) THEN OOM1F405.1597
C SET CONSTANTS OOM1F405.1598
RHO_FRESHWATER=1.0 ! IN G/CM^3 OOM1F405.1599
RHO0=RHO_WATER_SI*0.001 ! CONVERT TO CGS UNITS OOM1F405.1600
GRAV=GRAV_SI*100. ! CONVERT TO CGS UNITS OOM1F405.1601
EPSILON=0.1 OOM1F405.1602
C CALCULATE SURFACE FRICTION VELOCITY IN CGS UNITS OOM1F405.1603
IF (L_OUSTARWME) THEN OOM1F405.1604
C 1000.0 CONVERTS FROM SI UNITS TO CGS UNITS OOM1F405.1605
C USTAR^3 = WME/RHO0 OOM1F405.1606
DO I=1,IMT OOM1F405.1607
IF (FM(I,1).GT.0.5) THEN OOM1F405.1608
USTAR(I)=(1000.0*WME(I)/RHO0)**(1.0/3.0) OOM1F405.1609
ELSE OOM1F405.1610
USTAR(I)=0.0 OOM1F405.1611
ENDIF OOM1F405.1612
ENDDO OOM1F405.1613
ELSE OOM1F405.1614
C 10.0 CONVERTS FROM N M^-2 TO DYNES CM^-2 OOM1F405.1615
C USTAR^2 = TAU/RHO0 OOM1F405.1616
C OOM1F405.1617
C IF TAKE THIS CHOICE THEN CAREFUL ON OCEAN OOM1F405.1618
C POINTS ADJOINING LAND SINCE WSX/Y ARE SET TO MISSING DATA OOM1F405.1619
C VALUES ON ALL LAND POINTS OOM1F405.1620
C OOM1F405.1621
DO I=2,IMT OOM1F405.1622
WSXT=0.25*(WSX(I-1)+WSX(I)+WSXM(I-1)+WSXM(I)) OOM1F405.1623
WSYT=0.25*(WSY(I-1)+WSY(I)+WSYM(I-1)+WSYM(I)) OOM1F405.1624
USTAR(I)= SQRT(WSXT*WSXT+WSYT*WSYT) OOM1F405.1625
USTAR(I)= 10.0 * USTAR(I) / RHO0 OOM1F405.1626
USTAR(I)= SQRT(USTAR(I)) OOM1F405.1627
END DO OOM1F405.1628
IF (L_OCYCLIC) THEN OOM1F405.1629
USTAR(1)=USTAR(IMT-1) OOM1F405.1630
ELSE OOM1F405.1631
WSXT=0.5*(WSX(1)+WSXM(1)) OOM1F405.1632
WSYT=0.5*(WSY(1)+WSYM(1)) OOM1F405.1633
USTAR(1)= SQRT(WSXT*WSXT+WSYT*WSYT) OOM1F405.1634
USTAR(1)= 10.0 * FM(1,1) * USTAR(1) / RHO0 OOM1F405.1635
USTAR(1)= SQRT(USTAR(1)) OOM1F405.1636
ENDIF OOM1F405.1637
ENDIF OOM1F405.1638
C DO PRELIMINARY CALCULATION OF HT AS NEEDED TO CALC BUOYANCY FREQUENCY OOM1F405.1639
DO I=1,IMT OOM1F405.1640
HT(I)=MIN(MLD_LARGE(I),MAX_LARGE_DEPTH*100.) OOM1F405.1641
HT(I)=MAX(HT(I),1.5*ZDZ(1)) OOM1F405.1642
ENDDO OOM1F405.1643
C OOM1F405.1644
C CALCULATE SURFACE BUOYANCY FORCING , BF, IN CM^2/S^3 OOM1F405.1645
C ************************************************************** OOM1F405.1646
C OOM1F405.1647
C CALCULATE ALPHA AND BETA FOR ROW J OOM1F405.1648
CALL DRODT(TB,TB(1,1,2),ALPHAI,IMT,1) OOM1F405.1649
CALL DRODS(TB,TB(1,1,2),BETAI,IMT,1) OOM1F405.1650
C OOM1F405.1651
CP0 = SPECIFIC_HEAT_SI*1.E4 ! IN CM^2 S^-2 K^-1 OOM1F405.1652
C OOM1F405.1653
DO I=1,IMT OOM1F405.1654
ALPHA=-1.*ALPHAI(I,1) OOM1F405.1655
BETA=BETAI(I,1) OOM1F405.1656
C OOM1F405.1657
L_T(I)=0.0 OOM1F405.1658
C OOM1F405.1659
WS0(I)=0.0 OOM1F405.1660
WT0(I)=0.0 OOM1F405.1661
WTR0(I)=0.0 OOM1F405.1662
C OOM1F405.1663
IF (FM(I,1).GT.0.5) THEN OOM1F405.1664
C FACTORS CONVERT TO CGS UNITS. OOM1F405.1665
C HTN,SOL W/M^2=KG/S^3=1000G/S^3 OOM1F405.1666
C PME,WATERFLUX_ICE KG/M^2/S=0.1G/CM^2/S OOM1F405.1667
C ******************************************************************* OOM1F405.1668
WATERFLUX=PME(I)+WATERFLUX_ICE(I) OOM1F405.1669
C OOM1F405.1670
IF(L_SEAICE)THEN OOM1F405.1671
SFC_HEATFLUX=OCEANHEATFLUX(I)-FLXTOICE(I)+CARYHEAT(I) OOM1F405.1672
ELSE OOM1F405.1673
SFC_HEATFLUX=HTN(I) OOM1F405.1674
ENDIF OOM1F405.1675
C OOM1F405.1676
C ************************************************************** OOM1F405.1677
BF=GRAV * ( (ALPHA*SFC_HEATFLUX*1000.)/(RHO0*CP0) OOM1F405.1678
& + (BETA*WATERFLUX*0.1*0.035)/RHO_FRESHWATER ) OOM1F405.1679
& + GRAV*ALPHA* ( (SOL(I)*1000.)/(RHO0*CP0) OOM1F405.1680
& - (SOL_IRRA_1B(SOL(I)*1000.,HT(I)))/(RHO0*CP0) ) OOM1F405.1681
C ************************************************************** OOM1F405.1682
C OOM1F405.1683
C local surface variables for cgflux calculation OOM1F405.1684
C OOM1F405.1685
WS0(I)=(WATERFLUX*0.1*0.035)/RHO_FRESHWATER OOM1F405.1686
WT0(I)=-(SFC_HEATFLUX*1000.)/(RHO0*CP0) OOM1F405.1687
C note that define WTR=WTRZ-WTR0 OOM1F405.1688
WTR0(I)=(SOL(I)*1000.)/(RHO0*CP0) OOM1F405.1689
C OOM1F405.1690
IF (ABS(BF).LT.1.0E-12) THEN OOM1F405.1691
IF (BF.GE.0.0) BF=1.0E-12 OOM1F405.1692
IF (BF.LT.0.0) BF=-1.0E-12 OOM1F405.1693
ENDIF OOM1F405.1694
C CALCULATE MONIN-OBUKHOV LENGTH, L_T, CM = (CM/S)^3 / (CM^2/S^3) OOM1F405.1695
L_T(I) = USTAR(I)**3/(VKMAN*BF) OOM1F405.1696
IF (ABS(L_T(I)).LT.1.0E-12) THEN OOM1F405.1697
IF (L_T(I).GE.0.0) L_T(I)=1.0E-12 OOM1F405.1698
IF (L_T(I).LT.0.0) L_T(I)=-1.0E-12 OOM1F405.1699
ENDIF OOM1F405.1700
ENDIF OOM1F405.1701
ENDDO OOM1F405.1702
C CALCULATE DEPTH TO APPLY LARGE SCHEME TO (100 CONVERTS TO CM) OOM1F405.1703
DO I=1,IMT OOM1F405.1704
C OOM1F405.1705
HT(I)=MLD_LARGE(I) OOM1F405.1706
C OOM1F405.1707
IF (L_OWINDMIX.AND.L_T(I).GT.0.0) THEN OOM1F405.1708
HT(I)=MAX(HT(I),2*LAMBDA_LARGE*VKMAN*L_T(I)) OOM1F405.1709
ENDIF OOM1F405.1710
C OOM1F405.1711
IF (L_OBULKMAXMLD.AND.L_T(I).GT.0.0) THEN OOM1F405.1712
HT(I)=MIN(HT(I),L_T(I)) OOM1F405.1713
IF (.NOT.L_OROTATE) THEN OOM1F405.1714
C CALCULATE CORIOLIS PARAMETER F. IF PHIT<5 DEG, ASSUME PHIT=5 OOM1F405.1715
FX=2*OMEGA*SIN(MAX(PHIT,5.0*PI_OVER_180)) OOM1F405.1716
HT(I)=MIN(HT(I),0.7*USTAR(I)/FX) OOM1F405.1717
ENDIF OOM1F405.1718
ENDIF OOM1F405.1719
C OOM1F405.1720
C SOME FINAL CHECKS ON HT OOM1F405.1721
C OOM1F405.1722
HT(I)=MAX(HT(I),1.5*ZDZ(1)) OOM1F405.1723
HT(I)=MIN(HT(I),MAX_LARGE_DEPTH*100.) OOM1F405.1724
C OOM1F405.1725
ENDDO OOM1F405.1726
C CALCULATE MIN NO OF LEVELS WITHIN WHICH MIXED LAYER IS CONTAINED OOM1F405.1727
C I.E. LEVEL WITHIN WHICH THE MIXED LAYER LIES OOM1F405.1728
DO I=1,IMT OOM1F405.1729
MLK(I)=1 OOM1F405.1730
DO K=MAX_LARGE_LEVELS,1,-1 OOM1F405.1731
IF ((MLK(I).EQ.1).AND.(ZDZ(K).LE.HT(I))) MLK(I)=K+1 OOM1F405.1732
ENDDO OOM1F405.1733
ENDDO OOM1F405.1734
ELSE OOM1F405.1735
C calculate surface friction velocity in cgs units OLA3F403.321
C 10.0 converts from N m^-2 to dynes cm^-2 OLA3F403.322
C FM(I,1) is land-sea mask for tracer grid OLA3F403.323
C RHO0 = 1 g cm^-3 used in last line; tau = rho0 u*^2 OLA3F403.324
IF (L_OQLARGE) THEN OLA3F403.325
C SET CONSTANTS OOM1F405.1875
RHO0=RHO_WATER_SI*0.001 ! CONVERT TO CGS UNITS OOM1F405.1876
C CALCULATE SURFACE FRICTION VELOCITY IN CGS UNITS OOM1F405.1877
IF (L_OUSTARWME) THEN OOM1F405.1878
C 1000.0 CONVERTS FROM SI UNITS TO CGS UNITS OOM1F405.1879
C USTAR^3 = WME/RHO0 OOM1F405.1880
DO I=1,IMT OOM1F405.1881
IF (FM(I,1).GT.0.5) THEN OOM1F405.1882
USTAR(I)=(1000.0*WME(I)/RHO0)**(1.0/3.0) OOM1F405.1883
ELSE OOM1F405.1884
USTAR(I)=0.0 OOM1F405.1885
ENDIF OOM1F405.1886
ENDDO OOM1F405.1887
ELSE OOM1F405.1888
C OOM1F405.1889
C IF TAKE THIS CHOICE THEN CAREFUL ON OCEAN OOM1F405.1890
C POINTS ADJOINING LAND SINCE WSX/Y ARE SET TO MISSING DATA OOM1F405.1891
C VALUES ON ALL LAND POINTS OOM1F405.1892
C OOM1F405.1893
do i=2,imt OLA3F403.326
wsxT=0.25*(WSX(I-1)+WSX(I)+WSXM(I-1)+WSXM(I)) OLA3F403.327
wsyT=0.25*(WSY(I-1)+WSY(I)+WSYM(I-1)+WSYM(I)) OLA3F403.328
ustar(i)= sqrt(wsxT*wsxT+wsyT*wsyT) OLA3F403.329
ustar(i)= 10.0 * FM(I,1) * ustar(i) OLA3F403.330
ustar(i)= sqrt(ustar(i)) OLA3F403.331
end do OLA3F403.332
IF (L_OCYCLIC) THEN OLA3F403.333
ustar(1)=ustar(imt-1) OLA3F403.334
ELSE OLA3F403.335
wsxT=0.5*(WSX(1)+WSXM(1)) OLA3F403.336
wsyT=0.5*(WSY(1)+WSYM(1)) OLA3F403.337
ustar(1)= sqrt(wsxT*wsxT+wsyT*wsyT) OLA3F403.338
ustar(1)= 10.0 * FM(1,1) * ustar(1) OLA3F403.339
ustar(1)= sqrt(ustar(1)) OLA3F403.340
ENDIF OLA3F403.341
ENDIF !L_OUSTARWME OOM1F405.1894
C OOM1F405.1895
C Initialise hT to zero OLA3F403.343
do i=1,imt OLA0F404.164
hT(i)=0.0 OLA3F403.345
enddo OLA3F403.346
ENDIF ! L_OQLARGE OLA0F404.165
C Initialise RiT to zero OLA3F403.347
do k=1,kmm1 OLA3F403.348
do i=1,imt OLA3F403.349
RiT(i,k)=0.0 OLA3F403.350
enddo OLA3F403.351
enddo OLA3F403.352
ENDIF OOM1F405.1736
C VERTCOFT.68
fx=2. VERTCOFT.69
fxa=0.5 VERTCOFT.70
IF (L_OSTATEC) THEN OOM1F405.1500
C-------------- AT PRESENT HARDWIRE CODING FOR KM EVEN OOM1F405.1501
DO 210 K=1,KM-3,2 OOM1F405.1502
DO 211 I=2,IMT OOM1F405.1503
C OOM1F405.1504
UZ=((UB(I-1,K)-UB(I-1,K+1))+(UB(I,K)-UB(I,K+1))+ OOM1F405.1505
& (UBM(I-1,K)-UBM(I-1,K+1))+(UBM(I,K)-UBM(I,K+1))) OOM1F405.1506
& *FXA*DZZ2RQ(I,K+1) OOM1F405.1507
VZ=((VB(I-1,K)-VB(I-1,K+1))+(VB(I,K)-VB(I,K+1))+ OOM1F405.1508
& (VBM(I-1,K)-VBM(I-1,K+1))+(VBM(I,K)-VBM(I,K+1))) OOM1F405.1509
& *FXA*DZZ2RQ(I,K+1) OOM1F405.1510
DVEL2Q(I,K)=UZ*UZ+VZ*VZ OOM1F405.1511
UZ=((UB(I-1,K+1)-UB(I-1,K+2))+(UB(I,K+1)-UB(I,K+2))+ OOM1F405.1512
& (UBM(I-1,K+1)-UBM(I-1,K+2))+(UBM(I,K+1)-UBM(I,K+2))) OOM1F405.1513
& *FXA*DZZ2RQ(I,K+2) OOM1F405.1514
VZ=((VB(I-1,K+1)-VB(I-1,K+1))+(VB(I,K+1)-VB(I,K+2))+ OOM1F405.1515
& (VBM(I-1,K+1)-VBM(I-1,K+1))+(VBM(I,K+1)-VBM(I,K+2))) OOM1F405.1516
& *FXA*DZZ2RQ(I,K+2) OOM1F405.1517
DVEL2Q(I,K+1)=UZ*UZ+VZ*VZ OOM1F405.1518
211 CONTINUE OOM1F405.1519
DVEL2Q(1,K)=0.0 OOM1F405.1520
DVEL2Q(1,K+1)=0.0 OOM1F405.1521
DO 212 I=1,IMT OOM1F405.1522
RHOZA=(RHOSRNA(I,K)-RHOSRNA(I,K+1))*FX*DZZ2RQ(I,K+1) OOM1F405.1523
RHOZB=(RHOSRNB(I,K+1)-RHOSRNB(I,K+2))*FX*DZZ2RQ(I,K+2) OOM1F405.1524
RIPQ(I,K)=-GRAV_SI*100.*RHOZA ! 100. CONVERTS GRAV_SI TO CGS OOM1F405.1525
RIPQ(I,K+1)=-GRAV_SI*100.*RHOZB ! 100. CONVERTS GRAV_SI TO CGS OOM1F405.1526
212 CONTINUE OOM1F405.1527
210 CONTINUE OOM1F405.1528
C OOM1F405.1529
K=KMM1 OOM1F405.1530
DO I=2,IMT OOM1F405.1531
C OOM1F405.1532
UZ=((UB(I-1,K)-UB(I-1,K+1))+(UB(I,K)-UB(I,K+1))+ OOM1F405.1533
& (UBM(I-1,K)-UBM(I-1,K+1))+(UBM(I,K)-UBM(I,K+1))) OOM1F405.1534
& *FXA*DZZ2RQ(I,K+1) OOM1F405.1535
VZ=((VB(I-1,K)-VB(I-1,K+1))+(VB(I,K)-VB(I,K+1))+ OOM1F405.1536
& (VBM(I-1,K)-VBM(I-1,K+1))+(VBM(I,K)-VBM(I,K+1))) OOM1F405.1537
& *FXA*DZZ2RQ(I,K+1) OOM1F405.1538
DVEL2Q(I,K)=UZ*UZ+VZ*VZ OOM1F405.1539
ENDDO OOM1F405.1540
DVEL2Q(1,K)=0.0 OOM1F405.1541
DO I=1,IMT OOM1F405.1542
RHOZA=(RHOSRNA(I,K)-RHOSRNA(I,K+1))*FX*DZZ2RQ(I,K+1) OOM1F405.1543
RIPQ(I,K)=-GRAV_SI*100.*RHOZA ! 100. CONVERTS GRAV_SI TO CGS OOM1F405.1544
ENDDO OOM1F405.1545
ELSE OOM1F405.1546
DO 410 K=1,KMM1 VERTCOFT.71
DO 411 I=2,IMT VERTCOFT.72
C VERTCOFT.73
C -------------------------------------------------------------- VERTCOFT.74
C Evaluate Richardson No. VERTCOFT.75
C and thereby coefficient of diffusivity for mid-levels. VERTCOFT.76
C Note that the density rhosrn is calculated in CLINIC VERTCOFT.77
C VERTCOFT.78
C Note also that Ri = (g/rho0)*drhodz/(dudz*dudz + dvdz*dvdz). VERTCOFT.79
C This code assumes rho0=1 (OK in cgs units!) VERTCOFT.80
C -------------------------------------------------------------- VERTCOFT.81
C VERTCOFT.82
uz=((UB(I-1,K)-UB(I-1,K+1))+(UB(I,K)-UB(I,K+1))+ VERTCOFT.83
& (UBM(I-1,K)-UBM(I-1,K+1))+(UBM(I,K)-UBM(I,K+1))) VERTCOFT.84
& *fxa*DZZ2RQ(I,K+1) VERTCOFT.85
vz=((VB(I-1,K)-VB(I-1,K+1))+(VB(I,K)-VB(I,K+1))+ VERTCOFT.86
& (VBM(I-1,K)-VBM(I-1,K+1))+(VBM(I,K)-VBM(I,K+1))) VERTCOFT.87
& *fxa*DZZ2RQ(I,K+1) VERTCOFT.88
dvel2q(I,K)=uz*uz+vz*vz VERTCOFT.89
411 CONTINUE VERTCOFT.90
dvel2q(1,K)=0.0 VERTCOFT.91
DO 412 I=1,IMT VERTCOFT.92
rhoz=(rhosrn(I,K)-rhosrn(I,K+1))*fx*DZZ2RQ(I,K+1) VERTCOFT.93
ripq(I,K)=-GRAV_SI*100.*rhoz ! 100. converts GRAV_SI to cgs VERTCOFT.94
412 CONTINUE VERTCOFT.95
410 CONTINUE VERTCOFT.96
ENDIF ! L_OSTATEC OOM1F405.1547
C OOM1F405.1548
C VERTCOFT.97
C Convective adjustment takes care of unstable profiles VERTCOFT.98
C VERTCOFT.99
fx=0. VERTCOFT.100
DO 420 K=1,KMM1 VERTCOFT.101
DO 420 I=1,IMT VERTCOFT.102
IF (ripq(I,K).LT.fx) ripq(I,K)=fx VERTCOFT.103
420 CONTINUE VERTCOFT.104
C VERTCOFT.105
C If currents are zero set dvel2q to a small number.This is done VERTCOFT.106
C to prevent a divide check in the calculation of gnu VERTCOFT.107
C VERTCOFT.108
fx=1.E-12 VERTCOFT.109
DO 430 K=1,KMM1 VERTCOFT.110
DO 430 I=1,IMT VERTCOFT.111
IF (dvel2q(I,K).LT.fx) dvel2q(I,K)=fx VERTCOFT.112
430 CONTINUE VERTCOFT.113
C VERTCOFT.114
C OOM1F405.1549
IF(L_OPANDP)THEN OOM1F405.1550
C OOM1F405.1551
C Now calculate the diffusivity. Note that 1.E4 converts external VERTCOFT.115
C constants to cgs. The k's are offset by 1 since gnu(*,k) is the VERTCOFT.116
C diffusivity at the TOP of box k. gnu(*,1) is set at the end of VERTCOFT.117
C the routine. VERTCOFT.118
C VERTCOFT.119
fx=1. VERTCOFT.120
fxa=5. VERTCOFT.121
DO 440 K=1,KMM1 VERTCOFT.122
DO 440 I=1,IMT VERTCOFT.123
IF (FM(I,K+1).GT.0.9) THEN VERTCOFT.124
gnu(I,K+1)=FNU0_SI*1.E4*dvel2q(I,K)*dvel2q(I,K)/ VERTCOFT.125
& ((dvel2q(I,K)+fxa*ripq(I,K))* VERTCOFT.126
& (dvel2q(I,K)+fxa*ripq(I,K)))+FNUB_SI*1.E4 VERTCOFT.127
gnu(I,K+1)=gnu(I,K+1)*dvel2q(I,K)/ VERTCOFT.128
& (dvel2q(I,K)+fxa*ripq(I,K))+KAPPA_B_SI(K+1)*1.E4 RW071293.21
ENDIF VERTCOFT.130
440 CONTINUE VERTCOFT.131
C VERTCOFT.132
C OOM1F405.1552
ENDIF OOM1F405.1553
C OOM1F405.1554
C Calculate Richardson number OLA3F403.353
do K=1,KMM1 OLA3F403.354
do I=1,IMT OLA3F403.355
RiT(I,K)=ripq(I,K)/dvel2q(I,K) OLA3F403.356
end do OLA3F403.357
end do OLA3F403.358
C OOM1F405.1555
IF(.NOT.L_OPANDP)THEN OOM1F405.1556
C OOM1F405.1557
C CALCULATE GNU - TRACER DIFFUSIVITY USING PETERS. OOM1F405.1558
C Modify to limit dvel2q setting gnu if below a certain level. OOM1F405.1559
C If dvel2q too small, then regardless of ripq, set gnu to its OOM1F405.1560
C background value there. OOM1F405.1561
C GJR July 1998 OOM1F405.1562
C OOM1F405.1563
DO K=1,KMM1 OOM1F405.1564
DO I=1,IMT OOM1F405.1565
IF (FM(I,K+1).GT.0.9) THEN OOM1F405.1566
IF (DVEL2Q(I,K).GT.1.E-6)THEN OOM1F405.1567
IF (RIT(I,K).LT.RICRITT100) OOM1F405.1568
& GNU(I,K+1)=100. OOM1F405.1569
& +KAPPA_B_SI(K+1)*1.E4 OOM1F405.1570
IF (RIT(I,K).LT.RICRITT.AND.RIT(I,K).GT.RICRITT100) OOM1F405.1571
& GNU(I,K+1)=3.0E-5/(RIT(I,K)**9.6) OOM1F405.1572
& +KAPPA_B_SI(K+1)*1.E4 OOM1F405.1573
IF (RIT(I,K).GE.RICRITT) OOM1F405.1574
& GNU(I,K+1)=5.0/((1+5*RIT(I,K))**2.5) OOM1F405.1575
& +KAPPA_B_SI(K+1)*1.E4 OOM1F405.1576
ELSE OOM1F405.1577
GNU(I,K+1)=KAPPA_B_SI(K+1)*1.E4 OOM1F405.1578
ENDIF OOM1F405.1579
ENDIF OOM1F405.1580
ENDDO OOM1F405.1581
ENDDO OOM1F405.1582
C OOM1F405.1583
ENDIF OOM1F405.1584
C OOM1F405.1585
OLA3F403.359
IF (L_OFULARGE) THEN OOM1F405.1737
C CALCULATE VERTICAL DIFFUSIVITY, GNU, AT HT OOM1F405.1738
C OOM1F405.1739
DO I=1,IMT OOM1F405.1740
C OOM1F405.1741
IF (FM(I,1).GT.0.5.AND.HT(I).GE.ZDZ(1)) THEN OOM1F405.1742
C OOM1F405.1743
ZDZ_AT_MLKM1=0.0 OOM1F405.1744
GNUOFMLK=0.0 OOM1F405.1745
IF (MLK(I).GT.1) THEN OOM1F405.1746
ZDZ_AT_MLKM1=ZDZ(MLK(I)-1) OOM1F405.1747
GNUOFMLK=GNU(I,MLK(I)) OOM1F405.1748
ENDIF OOM1F405.1749
COEFF=(HT(I)-ZDZ_AT_MLKM1)/DZ(MLK(I)) OOM1F405.1750
GNUATHT=(1-COEFF)*GNUOFMLK+COEFF*GNU(I,MLK(I)+1) OOM1F405.1751
C CALCULATE D(GNU)/DZ AT HT TO RESTRICT TO MINIMUM OF 0.0 AS OTHERWISE OOM1F405.1752
C MAY OBTAIN POSSIBLE NEGATIVE VALUES OF GNU FROM LARGE SCHEME DUE OOM1F405.1753
C TO FITTING OF CUBIC FUNCTION. OOM1F405.1754
GNUZ(1)=(GNUOFMLK-GNU(I,MLK(I)+1))/DZ(MLK(I)) OOM1F405.1755
IF (GNUZ(1).GT.0.0) THEN OOM1F405.1756
IF (COEFF.LE.0.5) GNUZATHT=GNUZ(1) OOM1F405.1757
IF (COEFF.GT.0.5) THEN OOM1F405.1758
GNUZ(2)=(GNU(I,MLK(I)+1)-GNU(I,MLK(I)+2))/DZ(MLK(I)+1) OOM1F405.1759
COEFF=(HT(I)-ZDZZ(MLK(I)))/DZZ(MLK(I)+1) OOM1F405.1760
GNUZATHT=(1-COEFF)*GNUZ(1)+COEFF*GNUZ(2) OOM1F405.1761
GNUZATHT=MAX(GNUZATHT,0.0) OOM1F405.1762
ENDIF OOM1F405.1763
ELSE OOM1F405.1764
GNUZATHT=0.0 OOM1F405.1765
ENDIF OOM1F405.1766
C CALCULATE DEPTH DEPENDENT TURBULENT VELOCITY SCALE AT HT, OOM1F405.1767
C W1T, (LARGE EQN 13 WITH SIGMA=1) AND ITS DERIVATIVE WRT SIGMA, OOM1F405.1768
C DW1T. OOM1F405.1769
C OOM1F405.1770
IF (L_T(I).LT.0.0)THEN OOM1F405.1771
C OOM1F405.1772
STABIL_PARAM=EPSILON*HT(I)/L_T(I) OOM1F405.1773
W1T=VKMAN*USTAR(I)/FLUX_PROFILET(STABIL_PARAM) OOM1F405.1774
DW1T=0.0 OOM1F405.1775
C OOM1F405.1776
ELSE OOM1F405.1777
C OOM1F405.1778
STABIL_PARAM=HT(I)/L_T(I) OOM1F405.1779
W1T=VKMAN*USTAR(I)/FLUX_PROFILET(STABIL_PARAM) OOM1F405.1780
IF (0.0.LE.STABIL_PARAM) THEN OOM1F405.1781
DW1T = -5.0*VKMAN*USTAR(I)*STABIL_PARAM OOM1F405.1782
& /((1.0+5.0*STABIL_PARAM)**2) OOM1F405.1783
ELSE OOM1F405.1784
IF (-1.0.LE.STABIL_PARAM) THEN OOM1F405.1785
DW1T = -8.0*VKMAN*USTAR(I)*STABIL_PARAM OOM1F405.1786
& /((1.0-16.0*STABIL_PARAM)**0.5) OOM1F405.1787
ELSE OOM1F405.1788
DW1T = - 98.96/3.0 *VKMAN*USTAR(I)*STABIL_PARAM OOM1F405.1789
& /((-28.86-98.96*STABIL_PARAM)**(2.0/3.0)) OOM1F405.1790
ENDIF OOM1F405.1791
ENDIF OOM1F405.1792
C OOM1F405.1793
ENDIF ! L_T(I).LT.0.0 OOM1F405.1794
C OOM1F405.1795
C CALCULATE, WHEN SIGMA=Z/H=1, THE SHAPE FUNCTION AND ITS DERIVATIVE OOM1F405.1796
C WRT SIGMA (LARGE EQN 18) OOM1F405.1797
C N.B. G1=G(1)=SHAPE FUNCTION WHEN SIGMA=1 OOM1F405.1798
C DG1= DG(1)/DSIGMA OOM1F405.1799
IF (ABS(HT(I)).LT.1.0E-12.OR.ABS(W1T).LT.1.0E-12) THEN OOM1F405.1800
IF (W1T.GE.0.) W1T=1.0E-12 OOM1F405.1801
IF (W1T.LT.0.) W1T=-1.0E-12 OOM1F405.1802
ENDIF OOM1F405.1803
G1 = GNUATHT/(HT(I)*W1T) OOM1F405.1804
DG1 = - GNUZATHT/W1T - (GNUATHT*DW1T)/(HT(I)*W1T*W1T) OOM1F405.1805
C CALCULATE COEFFICIENTS OF SHAPE FUNCTION (LARGE EQN 17) OOM1F405.1806
A(1)=1. OOM1F405.1807
A(2)=-2.0+3.0*G1-DG1 OOM1F405.1808
A(3)=1.0-2.0*G1+DG1 OOM1F405.1809
C CALCULATE WT(K) - DEPTH DEPENDENT TURBULENT VELOCITY SCALE PROFILE OOM1F405.1810
C AT TOP OF LEVEL K (LARGE, EQN 13). WILL ALSO INCLUDE OOM1F405.1811
C COUNTER-GRADIENT FLUX CALCULATIONS HERE, TERMS OF WHICH ARE OOM1F405.1812
C DIFFERENT FOR POT. TEMP AND SALINITY, AND WHICH ARE ZERO OOM1F405.1813
C IF STABLE FORCING. THESE FLUXES ARE ALWAYS ZERO FOR MOMENTUM! OOM1F405.1814
C OOM1F405.1815
DO K=2,MLK(I) OOM1F405.1816
C OOM1F405.1817
STABIL_PARAM=ZDZ(K-1)/L_T(I) OOM1F405.1818
SIGMA=ZDZ(K-1)/HT(I) OOM1F405.1819
C OOM1F405.1820
IF(STABIL_PARAM.LT.0.0)THEN OOM1F405.1821
IF (SIGMA.GT.EPSILON)THEN OOM1F405.1822
STABIL_PARAM=(EPSILON*HT(I))/L_T(I) OOM1F405.1823
ENDIF OOM1F405.1824
WT(K)=(VKMAN*USTAR(I))/FLUX_PROFILET(STABIL_PARAM) OOM1F405.1825
C OOM1F405.1826
C FIND NON-ZERO COUNTER-GRAIDENT FLUX TERMS HERE. LARGE et al OOM1F405.1827
C RECOMMEND NOT INCLUDING SOLAR PENETRATION TERMS... OOM1F405.1828
C OOM1F405.1829
C GAMMAT=7.5*(WT0+(WTRZ-WTR0))/(WT(K)*HT(I)) OOM1F405.1830
C OOM1F405.1831
C WTRZ=SOL_IRRA_1B(SOL(I)*1000.,ZDZ(K-1))/(RHO0*CP0) OOM1F405.1832
C CGFLUX(I,K,1)=7.5*(WT0(I)+(WTRZ-WTR0(I)))/(WT(K)*HT(I)) OOM1F405.1833
CGFLUX(I,K,1)=(7.5*WT0(I))/(WT(K)*HT(I)) OOM1F405.1834
C OOM1F405.1835
C GAMMAS=7.5*WS0/(WT(K)*HT(I)) OOM1F405.1836
CGFLUX(I,K,2)=(7.5*WS0(I))/(WT(K)*HT(I)) OOM1F405.1837
C OOM1F405.1838
ELSE OOM1F405.1839
C OOM1F405.1840
WT(K)=(VKMAN*USTAR(I))/FLUX_PROFILET(STABIL_PARAM) OOM1F405.1841
CGFLUX(I,K,1)=0.0 OOM1F405.1842
CGFLUX(I,K,2)=0.0 OOM1F405.1843
C OOM1F405.1844
ENDIF OOM1F405.1845
C OOM1F405.1846
ENDDO OOM1F405.1847
C OOM1F405.1848
C CALCULATE SHAPE FUNCTION G THEN VERTICAL VISCOSITY PROFILE, GNU(K) OOM1F405.1849
C AT TOP OF LEVEL K IN THE MIXED LAYER (LARGE EQNS 11 AND 10) OOM1F405.1850
GNU(I,1)=0. OOM1F405.1851
DO K=2,MLK(I) OOM1F405.1852
SIGMA=ZDZ(K-1)/HT(I) OOM1F405.1853
G=SIGMA*(A(1)+SIGMA*(A(2)+A(3)*SIGMA)) OOM1F405.1854
GNU(I,K)= (HT(I)*WT(K))*G OOM1F405.1855
CGFLUX(I,K,1)=CGFLUX(I,K,1)*GNU(I,K) OOM1F405.1856
CGFLUX(I,K,2)=CGFLUX(I,K,2)*GNU(I,K) OOM1F405.1857
ENDDO OOM1F405.1858
C OOM1F405.1859
ELSE OOM1F405.1860
C OOM1F405.1861
DO K=1,KM OOM1F405.1862
GNU(I,K)=0.0 OOM1F405.1863
ENDDO OOM1F405.1864
C OOM1F405.1865
ENDIF !IF (FM(I,1).GT.0.5.AND.HT(I).GE.ZDZ(1)) OOM1F405.1866
C OOM1F405.1867
ENDDO ! DO I=1,IMT OOM1F405.1868
C OOM1F405.1869
ELSE OOM1F405.1870
IF (L_OQLARGE) THEN OLA3F403.360
C Simple quadratic form of Large scheme - OLA3F403.361
C applies - neutral Large in 'mixed layer' - defined to be where OLA3F403.362
C Richardson number is less than 0.3. OLA3F403.363
C - normal model Packonowski and Philander below mixed layer OLA3F403.364
C - functions specified so that vertical diffusion coeff is OLA3F403.365
C continuous at z=h OLA3F403.366
C - Constants chosen to give a quadratic form for G(z/h) OLA3F403.367
C Reference: OLA3F403.368
C W.G.Large et al 1994, Oceanic Vertical Mixing : A review and a OLA3F403.369
C model with a nonlocal boundary layer parametrisation, Rev Geophys, OLA3F403.370
C 32, 363-403. OLA3F403.371
C This scheme is a simplified neutral case of Large OLA3F403.372
C Vertical diffusion coeff, OLA3F403.373
C K(z)=h*kk*ustar*G(z/h) OLA3F403.374
C where z=depth (positive downwards) OLA3F403.375
C h=mixed layer depth OLA3F403.376
C kk=0.4 , von Karman's constant OLA0F404.166
C ustar= surface friction velocity OLA3F403.378
C = sqrt(wsx**2+wsy**2)/rho_w OLA3F403.379
C G(z/h)=z/h+a2*(z/h)**2 OLA3F403.380
C where a2 = Kath/(kk*ustar*h) - 1 OLA3F403.381
C Kath=K(Ri) at z=h OLA3F403.382
C To avoid problems where ustar is zero this is actually implemented as OLA3F403.383
C K(z)=h*kk*(z/h)*ustar*b2*(z/h)**2 OLA3F403.384
C where b2 = Kath/(kk*h) - ustar OLA3F403.385
C Scheme - for both momentum and tracers OLA3F403.386
C 1. Calculate gnu as usual OLA3F403.387
C 2. Find 'mixed layer depth' for Large scheme, h OLA3F403.388
C 3. For depths shallower than h, recalculate gnu using Large. OLA3F403.389
C OLA3F403.390
do i=1,imt OLA0F404.167
C Find 'mixed layer depth', hT, for quadratic Large scheme defined as OLA0F404.168
C min(depth where the critical Richardson no reached,max_qLarge_depth) OLA0F404.169
IF (FM(I,1).GT.0.9) THEN OLA0F404.170
do k=1,km OLA0F404.171
l=k OLA0F404.172
if (RiT(I,K).ge.crit_Ri.or.zdz(k).gt.max_qLarge_depth*100. OLA0F404.173
& .or.FM(i,min(k+1,km)).lt.0.5) goto 10 OLA0F404.174
enddo OLA0F404.175
10 continue OLA0F404.176
if (l.gt.1) then OLA0F404.177
if (RiT(i,l).ge.crit_Ri) then OLA0F404.178
if (RiT(i,l).ne.RiT(i,l-1)) then OLA0F404.179
hT(i)=zdz(l-1) + dz(l) OLA0F404.180
& *(crit_Ri-RiT(i,l-1))/(RiT(i,l)-RiT(i,l-1)) OLA0F404.181
hT(i)=min(hT(i),max_qLarge_depth*100.) OLA0F404.182
else OLA0F404.183
hT(i)=zdz(l-1) OLA0F404.184
endif OLA0F404.185
else OLA0F404.186
hT(i)=min(zdz(l),max_qLarge_depth*100.) OLA0F404.187
endif OLA0F404.188
C Calculate altered gnu i.e. applying quadratic Large where depth < hT OLA0F404.189
Kath=gnu(i,l-1) + (gnu(i,l)-gnu(i,l-1)) OLA0F404.190
& *(hT(i)-zdz(l-1))/dz(l) OLA0F404.191
b2 = Kath/(vkman*hT(i)) - ustar(i) OLA0F404.192
do k=1,l-1 OLA0F404.193
sigma=zdz(k)/hT(i) OLA0F404.194
gnu(i,k+1)=hT(i)*vkman*sigma*(ustar(i)+b2*sigma) OLA0F404.195
end do OLA0F404.196
endif OLA0F404.197
ENDIF OLA0F404.198
end do OLA3F403.415
ELSE OLA3F403.416
C The following occurs if the quadratic Large scheme is not chosen OLA3F403.417
C Reset gnu if value is greater than the stability limit. VERTCOFT.133
C 1.E4 converts STABLM_SI to cgs. VERTCOFT.134
C VERTCOFT.135
DO 450 K=2,KM VERTCOFT.136
DO 450 I=1,IMT VERTCOFT.137
IF (gnu(I,K).GT.STABLM_SI*1.E4) gnu(I,K)=STABLM_SI*1.E4 VERTCOFT.138
450 CONTINUE VERTCOFT.139
C VERTCOFT.140
ENDIF OLA3F403.418
C Reset diffusivity between top 2 levels to a min of GNUMINT_SI VERTCOFT.141
C The surface value gnu(*,1) is also set to GNUMINT_SI here. VERTCOFT.142
C VERTCOFT.143
fx=GNUMINT_SI*1.E4 ! 1.E4 converts to cgs VERTCOFT.144
DO 460 I=1,IMT VERTCOFT.145
IF (gnu(I,2).LT.fx) gnu(I,2)=fx*FM(I,2) VERTCOFT.146
gnu(I,1)=fx*FM(I,1) VERTCOFT.147
460 CONTINUE VERTCOFT.148
470 CONTINUE VERTCOFT.149
C VERTCOFT.150
ENDIF OOM1F405.1871
C Finally, multiply by the land-sea mask for safety. VERTCOFT.151
C VERTCOFT.152
DO 480 K=1,KM VERTCOFT.153
DO 490 I=1,IMT VERTCOFT.154
gnu(I,K)=gnu(I,K)*FM(I,K) VERTCOFT.155
DO NDUM=1,NT OOM1F405.1872
CGFLUX(I,K,NDUM)=CGFLUX(I,K,NDUM)*FM(I,K) OOM1F405.1873
ENDDO OOM1F405.1874
490 CONTINUE VERTCOFT.156
480 CONTINUE VERTCOFT.157
C VERTCOFT.158
RETURN VERTCOFT.159
END VERTCOFT.160
*ENDIF @DYALLOC.4682