New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Changeset 14779 – NEMO

Changeset 14779


Ignore:
Timestamp:
2021-05-04T12:21:30+02:00 (3 years ago)
Author:
smueller
Message:

Upgrade of internal subroutines zdf_osm_external_gradients, zdf_osm_calculate_dhdt, and zdf_osm_timestep_hbl to module subroutines of module zdfosm (ticket #2353)

File:
1 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2021/dev_r14122_HPC-08_Mueller_OSMOSIS_streamlining/src/OCE/ZDF/zdfosm.F90

    r14775 r14779  
    4545   !!         zdf_osm_velocity_rotation_3d :    rotation of 3d fields 
    4646   !!      zdf_osm_osbl_state              : determine the state of the OSBL 
     47   !!      zdf_osm_external_gradients      : calculate gradients below the OSBL 
     48   !!      zdf_osm_calculate_dhdt          : calculate rate of change of hbl 
     49   !!      zdf_osm_timestep_hbl            : hbl timestep 
    4750   !!      zdf_osm_diffusivity_viscosity   : compute eddy diffusivity and viscosity profiles 
    4851   !!      zdf_osm_fgr_terms               : compute flux-gradient relationship terms 
     
    318321 
    319322      REAL(wp), DIMENSION(jpi,jpj) :: zdh   ! pycnocline depth - grid 
    320       REAL(wp), DIMENSION(jpi,jpj) :: zdhdt ! BL depth tendency 
    321       REAL(wp), DIMENSION(jpi,jpj) :: zdtdz_bl_ext,zdsdz_bl_ext,zdbdz_bl_ext              ! external temperature/salinity and buoyancy gradients 
    322       REAL(wp), DIMENSION(jpi,jpj) :: zdtdz_mle_ext,zdsdz_mle_ext,zdbdz_mle_ext              ! external temperature/salinity and buoyancy gradients 
     323      REAL(wp), DIMENSION(A2D(0)) ::   zdhdt                                         ! BL depth tendency 
     324      REAL(wp), DIMENSION(A2D(0)) ::   zdtdz_bl_ext,  zdsdz_bl_ext,  zdbdz_bl_ext    ! external temperature/salinity and buoyancy gradients 
    323325      REAL(wp), DIMENSION(jpi,jpj) :: zdtdx, zdtdy, zdsdx, zdsdy      ! horizontal gradients for Fox-Kemper parametrization. 
    324326 
     
    376378      ! mixed layer 
    377379      ! no initialization of zhbl or zhml (or zdh?) 
    378       zhbl(:,:)   = zlarge ; zhml(:,:)   = zlarge  ; zdh(:,:)        = zlarge ; zdhdt(:,:)     = zlarge 
     380      zhbl(:,:)   = zlarge ; zhml(:,:)   = zlarge  ; zdh(:,:)        = zlarge 
    379381      zt_bl(:,:)  = zlarge ; zs_bl(:,:)  = zlarge  ; zu_bl(:,:)      = zlarge 
    380382      zv_bl(:,:)  = zlarge ; zb_bl(:,:)  = zlarge 
     
    390392      zdbdz_pyc(A2D(0),:) = 0.0_wp 
    391393      ! 
    392       zdtdz_bl_ext(:,:) = zlarge ; zdsdz_bl_ext(:,:) = zlarge ; zdbdz_bl_ext(:,:) = zlarge 
    393394 
    394395      IF ( ln_osm_mle ) THEN  ! only initialise arrays if needed 
     
    398399         zhmle(:,:)  = zlarge ; zmld(:,:)     = zlarge 
    399400      ENDIF 
    400       zwb_fk_b(:,:) = zlarge   ! must be initialised even with ln_osm_mle=F as used in zdf_osm_calculate_dhdt 
    401  
    402401      zhbl_t(:,:)   = zlarge 
    403402 
     
    746745         !! Calculate fairly-well-mixed depth zmld & its index mld_prof + lateral zmld-averaged gradients 
    747746         CALL zdf_osm_zmld_horizontal_gradients( zmld, zdtdx, zdtdy, zdsdx, zdsdy, dbdx_mle, dbdy_mle, zdbds_mle ) 
    748          !! Calculate vertical gradients immediately below zmld 
    749          CALL zdf_osm_external_gradients( mld_prof, zdtdz_mle_ext, zdsdz_mle_ext, zdbdz_mle_ext ) 
    750747         !! Calculate max vertical FK flux zwb_fk & set logical descriptors 
    751748         CALL zdf_osm_osbl_state_fk( lpyc, lflux, lmle, zwb_fk ) 
     
    761758               & 'dbdx_mle+=', dbdx_mle(ji,jj),' dbdx_mle-=', dbdx_mle(ji-1,jj),& 
    762759               & 'dbdy_mle+=', dbdy_mle(ji,jj),' dbdy_mle-=',dbdy_mle(ji,jj-1),' zdbds_mle=',zdbds_mle(ji,jj), & 
    763                & 'zdtdz_mle_ext=', zdtdz_mle_ext(ji,jj), ' zdsdz_mle_ext=', zdsdz_mle_ext(ji,jj), & 
    764                & ' zdbdz_mle_ext=', zdbdz_mle_ext(ji,jj), & 
    765760               & 'After updating hmle: mld_prof=', mld_prof(ji,jj),' hmle=', hmle(ji,jj), ' zhmle=', zhmle(ji,jj),& 
    766761               & 'zvel_mle =', zvel_mle(ji,jj), ' zdiff_mle=', zdiff_mle(ji,jj), ' zwb_fk=', zwb_fk(ji,jj) 
     
    779774 
    780775      !! External gradient below BL needed both with and w/o FK 
    781       CALL zdf_osm_external_gradients( ibld+1, zdtdz_bl_ext, zdsdz_bl_ext, zdbdz_bl_ext )   ! ag 19/03 
     776      CALL zdf_osm_external_gradients( Kmm, ibld+1, zdtdz_bl_ext, zdsdz_bl_ext, zdbdz_bl_ext )   ! ag 19/03 
    782777 
    783778      ! Test if pycnocline well resolved 
     
    834829 
    835830! Rate of change of hbl 
    836       CALL zdf_osm_calculate_dhdt( zdhdt ) 
     831      CALL zdf_osm_calculate_dhdt( lconv, lshear, j_ddh, zdhdt, zhbl, zdh, zwb_ent, zwb_min,   & 
     832         &                         zdb_bl, zdbdz_bl_ext, zdb_ml, zwb_fk_b, zwb_fk, zvel_mle ) 
    837833      ! Test if surface boundary layer coupled to bottom 
    838834      lcoup(:,:) = .FALSE.   ! ag 19/03 
     
    869865      ! Step through model levels taking account of buoyancy change to determine the effect on dhdt 
    870866      ! 
    871       CALL zdf_osm_timestep_hbl( zdhdt ) 
     867      CALL zdf_osm_timestep_hbl( Kmm, ibld, imld, lconv, lpyc, zdhdt, zhbl, zhbl_t, zt_bl, zs_bl, zwb_ent, zwb_fk_b ) 
    872868      ! is external level in bounds? 
    873869 
     
    914910         &                           jp_ext, zdt_ml, zds_ml, zdb_ml, zdu_ml, zdv_ml ) 
    915911 
    916       CALL zdf_osm_external_gradients( ibld+1, zdtdz_bl_ext, zdsdz_bl_ext, zdbdz_bl_ext ) 
     912      CALL zdf_osm_external_gradients( Kmm, ibld+1, zdtdz_bl_ext, zdsdz_bl_ext, zdbdz_bl_ext ) 
    917913#ifdef key_osm_debug 
    918914      IF(narea==nn_narea_db) THEN 
     
    989985      CALL zdf_osm_fgr_terms( Kmm, ibld, imld, jp_ext, lconv, lpyc, j_ddh, zhbl, zhml, zdh, zdhdt, zshear,          & 
    990986         &                    zdt_bl, zds_bl, zdb_bl, zdu_bl, zdv_bl, zdt_ml, zds_ml, zdb_ml, zdu_ml, zdv_ml,       & 
    991          &                    zdtdz_bl_ext, zdsdz_bl_ext, zdbdz_bl_ext, zdbdz_pyc, zalpha_pyc, zdiffut, zviscos ) 
     987         &                    zdtdz_bl_ext, zdsdz_bl_ext, zdbdz_pyc, zalpha_pyc, zdiffut, zviscos ) 
    992988 
    993989      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
     
    13401336      END SUBROUTINE zdf_osm_osbl_state_fk 
    13411337 
    1342       SUBROUTINE zdf_osm_external_gradients(jbase, zdtdz, zdsdz, zdbdz ) 
    1343          !!--------------------------------------------------------------------- 
    1344          !!                   ***  ROUTINE zdf_osm_external_gradients  *** 
    1345          !! 
    1346          !! ** Purpose : Calculates the gradients below the OSBL 
    1347          !! 
    1348          !! ** Method  : Uses ibld and ibld_ext to determine levels to calculate the gradient. 
    1349          !! 
    1350          !!---------------------------------------------------------------------- 
    1351           
    1352          INTEGER, DIMENSION(jpi,jpj)  :: jbase 
    1353          REAL(wp), DIMENSION(jpi,jpj) :: zdtdz, zdsdz, zdbdz   ! External gradients of temperature, salinity and buoyancy. 
    1354  
    1355          INTEGER :: jj, ji, jkb, jkb1 
    1356          REAL(wp) :: zthermal, zbeta 
    1357  
    1358  
    1359          IF( ln_timing ) CALL timing_start('zdf_osm_eg') 
    1360          DO_2D( 0, 0, 0, 0 ) 
    1361             IF ( jbase(ji,jj)+1 < mbkt(ji,jj) ) THEN 
    1362                zthermal = rab_n(ji,jj,1,jp_tem) !ideally use ibld not 1?? 
    1363                zbeta    = rab_n(ji,jj,1,jp_sal) 
    1364                jkb = jbase(ji,jj) 
    1365                jkb1 = MIN(jkb + 1, mbkt(ji,jj)) 
    1366                zdtdz(ji,jj) = - ( ts(ji,jj,jkb1,jp_tem,Kmm) - ts(ji,jj,jkb,jp_tem,Kmm ) ) & 
    1367                   &   / e3w(ji,jj,jkb1,Kmm) 
    1368                zdsdz(ji,jj) = - ( ts(ji,jj,jkb1,jp_sal,Kmm) - ts(ji,jj,jkb,jp_sal,Kmm ) ) & 
    1369                   &   / e3w(ji,jj,jkb1,Kmm) 
    1370                zdbdz(ji,jj) = grav * zthermal * zdtdz(ji,jj) - grav * zbeta * zdsdz(ji,jj) 
    1371             ELSE 
    1372                zdtdz(ji,jj) = 0._wp 
    1373                zdsdz(ji,jj) = 0._wp 
    1374                zdbdz(ji,jj) = 0._wp 
    1375             END IF 
    1376          END_2D 
    1377          IF( ln_timing ) CALL timing_stop('zdf_osm_eg') 
    1378       END SUBROUTINE zdf_osm_external_gradients 
    1379  
    13801338      SUBROUTINE zdf_osm_pycnocline_buoyancy_profiles( pdbdz, palpha ) 
    13811339         REAL(wp), DIMENSION(:,:,:), INTENT( inout ) ::   pdbdz   ! Gradients in the pycnocline 
     
    14711429         ! 
    14721430      END SUBROUTINE zdf_osm_pycnocline_buoyancy_profiles 
    1473  
    1474       SUBROUTINE zdf_osm_calculate_dhdt( zdhdt ) 
    1475          !!--------------------------------------------------------------------- 
    1476          !!                   ***  ROUTINE zdf_osm_calculate_dhdt  *** 
    1477          !! 
    1478          !! ** Purpose : Calculates the rate at which hbl changes. 
    1479          !! 
    1480          !! ** Method  : 
    1481          !! 
    1482          !!---------------------------------------------------------------------- 
    1483  
    1484          REAL(wp), DIMENSION(jpi,jpj) :: zdhdt        ! Rate of change of hbl 
    1485  
    1486          INTEGER :: jj, ji 
    1487          REAL(wp) :: zgamma_b_nd, zgamma_dh_nd, zpert, zpsi 
    1488          REAL(wp) :: zvel_max, zddhdt 
    1489          REAL(wp), PARAMETER :: zzeta_m  = 0.3_wp 
    1490          REAL(wp), PARAMETER :: zgamma_c = 2.0_wp 
    1491          REAL(wp), PARAMETER :: zdhoh    = 0.1_wp 
    1492          REAL(wp), PARAMETER :: zalpha_b = 0.3_wp 
    1493          REAL(wp), PARAMETER :: a_ddh    = 2.5_wp, a_ddh_2 = 3.5 ! also in pycnocline_depth 
    1494  
    1495          IF( ln_timing ) CALL timing_start('zdf_osm_cd') 
    1496          DO_2D( 0, 0, 0, 0 ) 
    1497  
    1498             IF ( lshear(ji,jj) ) THEN 
    1499                IF ( lconv(ji,jj) ) THEN    ! Convective 
    1500  
    1501                   IF ( ln_osm_mle ) THEN 
    1502  
    1503                      IF ( hmle(ji,jj) > hbl(ji,jj) ) THEN 
    1504                         ! Fox-Kemper buoyancy flux average over OSBL 
    1505                         zwb_fk_b(ji,jj) = zwb_fk(ji,jj) *  & 
    1506                            (1.0 + hmle(ji,jj) / ( 6.0 * hbl(ji,jj) ) * (-1.0 + ( 1.0 - 2.0 * hbl(ji,jj) / hmle(ji,jj))**3) ) 
    1507                      ELSE 
    1508                         zwb_fk_b(ji,jj) = 0.5 * zwb_fk(ji,jj) * hmle(ji,jj) / hbl(ji,jj) 
    1509                      ENDIF 
    1510                      zvel_max = ( svstr(ji,jj)**3 + 0.5 * swstrc(ji,jj)**3 )**p2third / hbl(ji,jj) 
    1511                      IF ( ( zwb_ent(ji,jj) + 2.0 * zwb_fk_b(ji,jj) ) < 0.0 ) THEN 
    1512                         ! OSBL is deepening, entrainment > restratification 
    1513                         IF ( zdb_bl(ji,jj) > 1e-15 ) THEN 
    1514                            zgamma_b_nd = MAX( zdbdz_bl_ext(ji,jj), 0.0_wp ) * zdh(ji,jj) / ( zvel_max + MAX( zdb_bl(ji,jj), 1e-15 ) ) 
    1515                            zpsi = ( 1.0_wp - 0.5_wp * zdh(ji,jj) / zhbl(ji,jj) ) *   & 
    1516                               &   ( swb0(ji,jj) - MIN( ( zwb_min(ji,jj) + 2.0_wp * zwb_fk_b(ji,jj) ), 0.0_wp ) ) * zdh(ji,jj) / zhbl(ji,jj) 
    1517                            zpsi = zpsi + 1.75_wp * ( 1.0_wp - 0.5_wp * zdh(ji,jj) / zhbl(ji,jj) ) *   & 
    1518                               &   ( zdh(ji,jj) / zhbl(ji,jj) + zgamma_b_nd ) * MIN( ( zwb_min(ji,jj) + 2.0_wp * zwb_fk_b(ji,jj) ), 0.0_wp ) 
    1519                            zpsi = zalpha_b * MAX( zpsi, 0.0_wp ) 
    1520                            zdhdt(ji,jj) = -1.0_wp * ( zwb_ent(ji,jj) + 2.0_wp * zwb_fk_b(ji,jj) ) /   & 
    1521                               &                     ( zvel_max + MAX( zdb_bl(ji,jj), 1e-15_wp ) ) +   & 
    1522                               &           zpsi / ( zvel_max + MAX( zdb_bl(ji,jj), 1e-15 ) ) 
    1523 #ifdef key_osm_debug 
    1524                            IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN 
    1525                               WRITE(narea+100,'(a,g11.3)')'Inside 1st major loop of zdf_osm_calculate_dhdt, OSBL is deepening, entrainment > restratification:  zdhdt=',zdhdt(ji,jj) 
    1526                               WRITE(narea+100,'(3(a,g11.3))') '  zpsi=',zpsi, '  zgamma_b_nd=', zgamma_b_nd, '  zdh=', zdh(ji,jj) 
    1527                               FLUSH(narea+100) 
    1528                            END IF 
    1529 #endif 
    1530                            IF ( j_ddh(ji,jj) == 1 ) THEN 
    1531                               IF ( ( swstrc(ji,jj) / svstr(ji,jj) )**3 <= 0.5 ) THEN 
    1532                                  zari = MIN( 1.5 * zdb_bl(ji,jj) / ( zhbl(ji,jj) * ( MAX(zdbdz_bl_ext(ji,jj),0._wp) + zdb_bl(ji,jj)**2 / MAX(4.5 * svstr(ji,jj)**2 , 1.e-12 )) ), 0.2d0 ) 
    1533                               ELSE 
    1534                                  zari = MIN( 1.5 * zdb_bl(ji,jj) / ( zhbl(ji,jj) * ( MAX(zdbdz_bl_ext(ji,jj),0._wp) + zdb_bl(ji,jj)**2 / MAX(4.5 * swstrc(ji,jj)**2 , 1.e-12 )) ), 0.2d0 ) 
    1535                               ENDIF 
    1536                               ! Relaxation to dh_ref = zari * hbl 
    1537                               zddhdt = -1.0_wp * a_ddh_2 * ( 1.0 - zdh(ji,jj) / ( zari * zhbl(ji,jj) ) ) * zwb_ent(ji,jj) /   & 
    1538                                  &     ( zvel_max + MAX( zdb_bl(ji,jj), 1e-15 ) ) 
    1539 #ifdef key_osm_debug 
    1540                               IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN 
    1541                                  WRITE(narea+100,'(a,g11.3)')'Inside 1st major loop of zdf_osm_calculate_dhdt,j_ddh(ji,jj) == 1:  zari=',zari 
    1542                                  FLUSH(narea+100) 
    1543                               END IF 
    1544 #endif 
    1545  
    1546                            ELSE IF ( j_ddh(ji,jj) == 0 ) THEN 
    1547                               ! Growing shear layer 
    1548                               zddhdt = -1.0_wp * a_ddh * ( 1.0 - 1.6_wp * zdh(ji,jj) / zhbl(ji,jj) ) * zwb_ent(ji,jj) /   & 
    1549                                  &     ( zvel_max + MAX( zdb_bl(ji,jj), 1e-15 ) ) 
    1550                               zddhdt = EXP( -4.0_wp * ABS( ff_t(ji,jj) ) * zhbl(ji,jj) / MAX(sustar(ji,jj), 1e-8_wp ) ) * zddhdt 
    1551                            ELSE 
    1552                               zddhdt = 0.0_wp 
    1553                            ENDIF ! j_ddh 
    1554                            zdhdt(ji,jj) = zdhdt(ji,jj) + zalpha_b * ( 1.0_wp - 0.5_wp * zdh(ji,jj) / zhbl(ji,jj) ) *   & 
    1555                               &                          zdb_ml(ji,jj) * MAX( zddhdt, 0.0_wp ) / ( zvel_max + MAX( zdb_bl(ji,jj), 1e-15 ) ) 
    1556                         ELSE    ! zdb_bl >0 
    1557                            zdhdt(ji,jj) = -( zwb_ent(ji,jj) + 2.0 * zwb_fk_b(ji,jj) ) /  MAX( zvel_max, 1.0e-15) 
    1558                         ENDIF 
    1559                      ELSE   ! zwb_min + 2*zwb_fk_b < 0 
    1560                         ! OSBL shoaling due to restratification flux. This is the velocity defined in Fox-Kemper et al (2008) 
    1561                         zdhdt(ji,jj) = -1.0_wp * MIN( zvel_mle(ji,jj), hbl(ji,jj) / 10800.0_wp ) 
    1562  
    1563  
    1564                      ENDIF 
    1565  
    1566                   ELSE 
    1567                      ! Fox-Kemper not used. 
    1568  
    1569                      zvel_max = - ( 1.0 + 1.0 * ( svstr(ji,jj)**3 + 0.5 * swstrc(ji,jj)**3 )**pthird * rn_Dt / hbl(ji,jj) ) * zwb_ent(ji,jj) / & 
    1570                         &        MAX((svstr(ji,jj)**3 + 0.5 * swstrc(ji,jj)**3 )**pthird, epsln) 
    1571                      zdhdt(ji,jj) = -zwb_ent(ji,jj) / ( zvel_max + MAX(zdb_bl(ji,jj), 1.0e-15) ) 
    1572                      ! added ajgn 23 July as temporay fix 
    1573  
    1574                   ENDIF  ! ln_osm_mle 
    1575  
    1576                ELSE    ! lconv - Stable 
    1577                   zdhdt(ji,jj) = ( 0.06 + 0.52 * shol(ji,jj) / 2.0 ) * svstr(ji,jj)**3 / hbl(ji,jj) + swbav(ji,jj) 
    1578                   IF ( zdhdt(ji,jj) < 0._wp ) THEN 
    1579                      ! For long timsteps factor in brackets slows the rapid collapse of the OSBL 
    1580                      zpert = 2.0 * ( 1.0 + 0.0 * 2.0 * svstr(ji,jj) * rn_Dt / hbl(ji,jj) ) * svstr(ji,jj)**2 / hbl(ji,jj) 
    1581                   ELSE 
    1582                      zpert = MAX( svstr(ji,jj)**2 / hbl(ji,jj), zdb_bl(ji,jj) ) 
    1583                   ENDIF 
    1584                   zdhdt(ji,jj) = 2.0 * zdhdt(ji,jj) / MAX(zpert, epsln) 
    1585                   zdhdt(ji,jj) = MAX( zdhdt(ji,jj), -1.0_wp * hbl(ji,jj) / 5400.0_wp ) 
    1586                ENDIF  ! lconv 
    1587             ELSE ! lshear 
    1588                IF ( lconv(ji,jj) ) THEN    ! Convective 
    1589  
    1590                   IF ( ln_osm_mle ) THEN 
    1591  
    1592                      IF ( hmle(ji,jj) > hbl(ji,jj) ) THEN 
    1593                         ! Fox-Kemper buoyancy flux average over OSBL 
    1594                         zwb_fk_b(ji,jj) = zwb_fk(ji,jj) *  & 
    1595                            (1.0 + hmle(ji,jj) / ( 6.0 * hbl(ji,jj) ) * (-1.0 + ( 1.0 - 2.0 * hbl(ji,jj) / hmle(ji,jj))**3) ) 
    1596                      ELSE 
    1597                         zwb_fk_b(ji,jj) = 0.5 * zwb_fk(ji,jj) * hmle(ji,jj) / hbl(ji,jj) 
    1598                      ENDIF 
    1599                      zvel_max = ( swstrl(ji,jj)**3 + 0.5 * swstrc(ji,jj)**3 )**p2third / hbl(ji,jj) 
    1600                      IF ( ( zwb_ent(ji,jj) + 2.0 * zwb_fk_b(ji,jj) ) < 0.0 ) THEN 
    1601                         ! OSBL is deepening, entrainment > restratification 
    1602                         IF ( zdb_bl(ji,jj) > 0.0 .and. zdbdz_bl_ext(ji,jj) > 0.0 ) THEN 
    1603                            zdhdt(ji,jj) = -( zwb_ent(ji,jj) + 2.0 * zwb_fk_b(ji,jj) ) / ( zvel_max + MAX(zdb_bl(ji,jj), 1.0e-15) ) 
    1604                         ELSE 
    1605                            zdhdt(ji,jj) = -( zwb_ent(ji,jj) + 2.0 * zwb_fk_b(ji,jj) ) /  MAX( zvel_max, 1.0e-15) 
    1606                         ENDIF 
    1607                      ELSE 
    1608                         ! OSBL shoaling due to restratification flux. This is the velocity defined in Fox-Kemper et al (2008) 
    1609                         zdhdt(ji,jj) = -1.0_wp * MIN( zvel_mle(ji,jj), hbl(ji,jj) / 10800.0_wp ) 
    1610  
    1611  
    1612                      ENDIF 
    1613  
    1614                   ELSE 
    1615                      ! Fox-Kemper not used. 
    1616  
    1617                      zvel_max = -zwb_ent(ji,jj) / & 
    1618                         &        MAX((svstr(ji,jj)**3 + 0.5 * swstrc(ji,jj)**3 )**pthird, epsln) 
    1619                      zdhdt(ji,jj) = -zwb_ent(ji,jj) / ( zvel_max + MAX(zdb_bl(ji,jj), 1.0e-15) ) 
    1620                      ! added ajgn 23 July as temporay fix 
    1621  
    1622                   ENDIF  ! ln_osm_mle 
    1623  
    1624                ELSE                        ! Stable 
    1625                   zdhdt(ji,jj) = ( 0.06 + 0.52 * shol(ji,jj) / 2.0 ) * svstr(ji,jj)**3 / hbl(ji,jj) + swbav(ji,jj) 
    1626                   IF ( zdhdt(ji,jj) < 0._wp ) THEN 
    1627                      ! For long timsteps factor in brackets slows the rapid collapse of the OSBL 
    1628                      zpert = 2.0 * svstr(ji,jj)**2 / hbl(ji,jj) 
    1629                   ELSE 
    1630                      zpert = MAX( svstr(ji,jj)**2 / hbl(ji,jj), zdb_bl(ji,jj) ) 
    1631                   ENDIF 
    1632                   zdhdt(ji,jj) = 2.0 * zdhdt(ji,jj) / MAX(zpert, epsln) 
    1633                   zdhdt(ji,jj) = MAX( zdhdt(ji,jj), -1.0_wp * hbl(ji,jj) / 5400.0_wp ) 
    1634                ENDIF  ! lconv 
    1635             ENDIF ! lshear 
    1636 #ifdef key_osm_debug 
    1637             IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN 
    1638                WRITE(narea+100,'(4(a,g11.3))')'end of 1st major loop of zdf_osm_calculate_dhdt:  zdhdt=',zdhdt(ji,jj), & 
    1639                   &  '  zpert=', zpert, '  zddhdt=', zddhdt, '  zvel_max=', zvel_max 
    1640  
    1641                IF ( ln_osm_mle ) THEN 
    1642                   WRITE(narea+100,'(3(a,g11.3),/)') 'zvel_mle=',zvel_mle(ji,jj), ' zwb_fk_b=', zwb_fk_b(ji,jj), & 
    1643                      & '  zwb_ent + 2*zwb_fk_b =', zwb_ent(ji,jj) + 2.0 * zwb_fk_b(ji,jj) 
    1644                   FLUSH(narea+100) 
    1645                END IF 
    1646             END IF 
    1647 #endif 
    1648          END_2D 
    1649          IF( ln_timing ) CALL timing_stop('zdf_osm_cd') 
    1650       END SUBROUTINE zdf_osm_calculate_dhdt 
    1651  
    1652       SUBROUTINE zdf_osm_timestep_hbl( zdhdt ) 
    1653          !!--------------------------------------------------------------------- 
    1654          !!                   ***  ROUTINE zdf_osm_timestep_hbl  *** 
    1655          !! 
    1656          !! ** Purpose : Increments hbl. 
    1657          !! 
    1658          !! ** Method  : If thechange in hbl exceeds one model level the change is 
    1659          !!              is calculated by moving down the grid, changing the buoyancy 
    1660          !!              jump. This is to ensure that the change in hbl does not 
    1661          !!              overshoot a stable layer. 
    1662          !! 
    1663          !!---------------------------------------------------------------------- 
    1664  
    1665  
    1666          REAL(wp), DIMENSION(jpi,jpj) :: zdhdt   ! rates of change of hbl. 
    1667  
    1668          INTEGER :: jk, jj, ji, jm 
    1669          REAL(wp) :: zhbl_s, zvel_max, zdb 
    1670          REAL(wp) :: zthermal, zbeta 
    1671  
    1672          IF( ln_timing ) CALL timing_start('zdf_osm_th') 
    1673          DO_2D( 0, 0, 0, 0 ) 
    1674 #ifdef key_osm_debug 
    1675             IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN 
    1676                WRITE(narea+100,'(2(a,i7))')'start of zdf_osm_timestep_hbl: old ibld=',imld(ji,jj),' trial ibld=', ibld(ji,jj) 
    1677                FLUSH(narea+100) 
    1678             END IF 
    1679 #endif 
    1680             IF ( ibld(ji,jj) - imld(ji,jj) > 1 ) THEN 
    1681                ! 
    1682                ! If boundary layer changes by more than one level, need to check for stable layers between initial and final depths. 
    1683                ! 
    1684                zhbl_s = hbl(ji,jj) 
    1685                jm = imld(ji,jj) 
    1686                zthermal = rab_n(ji,jj,1,jp_tem) 
    1687                zbeta = rab_n(ji,jj,1,jp_sal) 
    1688  
    1689  
    1690                IF ( lconv(ji,jj) ) THEN 
    1691                   !unstable 
    1692  
    1693                   IF( ln_osm_mle ) THEN 
    1694                      zvel_max = ( swstrl(ji,jj)**3 + swstrc(ji,jj)**3 )**p2third / hbl(ji,jj) 
    1695                   ELSE 
    1696  
    1697                      zvel_max = -( 1.0 + 1.0 * ( svstr(ji,jj)**3 + 0.5 * swstrc(ji,jj)**3 )**pthird * rn_Dt / hbl(ji,jj) ) * zwb_ent(ji,jj) / & 
    1698                         &      ( svstr(ji,jj)**3 + 0.5 * swstrc(ji,jj)**3 )**pthird 
    1699  
    1700                   ENDIF 
    1701 #ifdef key_osm_debug 
    1702                   IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN 
    1703                      WRITE(narea+100,'(a,g11.3)')'In zdf_osm_timestep_hbl, ibld - imld > 1, lconv=T: zvel_max=',zvel_max 
    1704                      FLUSH(narea+100) 
    1705                   END IF 
    1706 #endif 
    1707  
    1708                   DO jk = imld(ji,jj), ibld(ji,jj) 
    1709                      zdb = MAX( grav * ( zthermal * ( zt_bl(ji,jj) - ts(ji,jj,jm,jp_tem,Kmm) ) & 
    1710                         & - zbeta * ( zs_bl(ji,jj) - ts(ji,jj,jm,jp_sal,Kmm) ) ), & 
    1711                         &  0.0 ) + zvel_max 
    1712  
    1713  
    1714                      IF ( ln_osm_mle ) THEN 
    1715                         zhbl_s = zhbl_s + MIN( & 
    1716                            & rn_Dt * ( ( -zwb_ent(ji,jj) - 2.0 * zwb_fk_b(ji,jj) )/ zdb ) / FLOAT(ibld(ji,jj) - imld(ji,jj) ), & 
    1717                            & e3w(ji,jj,jm,Kmm) ) 
    1718                      ELSE 
    1719                         zhbl_s = zhbl_s + MIN( & 
    1720                            & rn_Dt * (  -zwb_ent(ji,jj) / zdb ) / FLOAT(ibld(ji,jj) - imld(ji,jj) ), & 
    1721                            & e3w(ji,jj,jm,Kmm) ) 
    1722                      ENDIF 
    1723  
    1724                      !                    zhbl_s = MIN(zhbl_s,  gdepw(ji,jj, mbkt(ji,jj) + 1,Kmm) - depth_tol) 
    1725                      IF ( zhbl_s >= gdepw(ji,jj,mbkt(ji,jj) + 1,Kmm) ) THEN 
    1726                         zhbl_s = MIN(zhbl_s,  gdepw(ji,jj, mbkt(ji,jj) + 1,Kmm) - depth_tol) 
    1727                         lpyc(ji,jj) = .FALSE. 
    1728                      ENDIF 
    1729                      IF ( zhbl_s >= gdepw(ji,jj,jm+1,Kmm) ) jm = jm + 1 
    1730 #ifdef key_osm_debug 
    1731                      IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN 
    1732                         WRITE(narea+100,'(2(a,i7))')' jk=',jk,' jm=', jm 
    1733                         WRITE(narea+100,'(2(a,g11.3),a,l7)')'zdb=',zdb,' zhbl_s=', zhbl_s,' lpyc=',lpyc(ji,jj) 
    1734                         FLUSH(narea+100) 
    1735                      END IF 
    1736 #endif 
    1737                   END DO 
    1738                   hbl(ji,jj) = zhbl_s 
    1739                   ibld(ji,jj) = jm 
    1740                ELSE 
    1741                   ! stable 
    1742 #ifdef key_osm_debug 
    1743                   IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN 
    1744                      WRITE(narea+100,'(a)')'In zdf_osm_timestep_hbl, ibld - imld > 1, lconv=F' 
    1745                      FLUSH(narea+100) 
    1746                   END IF 
    1747 #endif 
    1748                   DO jk = imld(ji,jj), ibld(ji,jj) 
    1749                      zdb = MAX( & 
    1750                         & grav * ( zthermal * ( zt_bl(ji,jj) - ts(ji,jj,jm,jp_tem,Kmm) )& 
    1751                         &           - zbeta * ( zs_bl(ji,jj) - ts(ji,jj,jm,jp_sal,Kmm) ) ),& 
    1752                         & 0.0 ) + & 
    1753                         &       2.0 * svstr(ji,jj)**2 / zhbl_s 
    1754  
    1755                      ! Alan is thuis right? I have simply changed hbli to hbl 
    1756                      shol(ji,jj) = -zhbl_s / ( ( svstr(ji,jj)**3 + epsln )/ swbav(ji,jj) ) 
    1757                      zdhdt(ji,jj) = -( swbav(ji,jj) - 0.04 / 2.0 * swstrl(ji,jj)**3 / zhbl_s - 0.15 / 2.0 * ( 1.0 - EXP( -1.5 * sla(ji,jj) ) ) * & 
    1758                         &                  sustar(ji,jj)**3 / zhbl_s ) * ( 0.725 + 0.225 * EXP( -7.5 * shol(ji,jj) ) ) 
    1759                      zdhdt(ji,jj) = zdhdt(ji,jj) + swbav(ji,jj) 
    1760                      zhbl_s = zhbl_s + MIN( zdhdt(ji,jj) / zdb * rn_Dt / FLOAT( ibld(ji,jj) - imld(ji,jj) ), e3w(ji,jj,jm,Kmm) ) 
    1761  
    1762                      !                    zhbl_s = MIN(zhbl_s, gdepw(ji,jj, mbkt(ji,jj) + 1,Kmm) - depth_tol) 
    1763                      IF ( zhbl_s >= mbkt(ji,jj) + 1 ) THEN 
    1764                         zhbl_s = MIN(zhbl_s,  gdepw(ji,jj, mbkt(ji,jj) + 1,Kmm) - depth_tol) 
    1765                         lpyc(ji,jj) = .FALSE. 
    1766                      ENDIF 
    1767                      IF ( zhbl_s >= gdepw(ji,jj,jm,Kmm) ) jm = jm + 1 
    1768 #ifdef key_osm_debug 
    1769                      IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN 
    1770                         WRITE(narea+100,'(2(a,i7))')' jk=',jk,' jm=', jm 
    1771                         WRITE(narea+100,'(4(a,g11.3),a,l7)')'zdb=',zdb,' shol',shol(ji,jj),' zdhdt',zdhdt(ji,jj),' zhbl_s=', zhbl_s,' lpyc=',lpyc(ji,jj) 
    1772                         FLUSH(narea+100) 
    1773                      END IF 
    1774 #endif 
    1775                   END DO 
    1776                ENDIF   ! IF ( lconv ) 
    1777                hbl(ji,jj) = MAX(zhbl_s, gdepw(ji,jj,4,Kmm) ) 
    1778                ibld(ji,jj) = MAX(jm, 4 ) 
    1779             ELSE 
    1780                ! change zero or one model level. 
    1781                hbl(ji,jj) = MAX(zhbl_t(ji,jj), gdepw(ji,jj,4,Kmm) ) 
    1782             ENDIF 
    1783             zhbl(ji,jj) = gdepw(ji,jj,ibld(ji,jj),Kmm) 
    1784 #ifdef key_osm_debug 
    1785             IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN 
    1786                WRITE(narea+100,'(2(a,g11.3),a,i7,/)')'end of zdf_osm_timestep_hbl: hbl=', hbl(ji,jj),' zhbl=', zhbl(ji,jj),' ibld=', ibld(ji,jj) 
    1787                FLUSH(narea+100) 
    1788             END IF 
    1789 #endif 
    1790          END_2D 
    1791          IF( ln_timing ) CALL timing_stop('zdf_osm_th') 
    1792  
    1793       END SUBROUTINE zdf_osm_timestep_hbl 
    17941431 
    17951432      SUBROUTINE zdf_osm_pycnocline_thickness( dh, zdh ) 
     
    25222159   END SUBROUTINE zdf_osm_osbl_state 
    25232160 
     2161   SUBROUTINE zdf_osm_external_gradients( Kmm, kbase, pdtdz, pdsdz, pdbdz ) 
     2162      !!--------------------------------------------------------------------- 
     2163      !!                   ***  ROUTINE zdf_osm_external_gradients  *** 
     2164      !! 
     2165      !! ** Purpose : Calculates the gradients below the OSBL 
     2166      !! 
     2167      !! ** Method  : Uses ibld and ibld_ext to determine levels to calculate the gradient. 
     2168      !! 
     2169      !!----------------------------------------------------------------------    
     2170      INTEGER,                      INTENT(in   ) ::   Kmm                   ! Ocean time-level index 
     2171      INTEGER,  DIMENSION(jpi,jpj), INTENT(in   ) ::   kbase                 ! OSBL base layer index 
     2172      REAL(wp), DIMENSION(A2D(0)),  INTENT(  out) ::   pdtdz, pdsdz, pdbdz   ! External gradients of temperature, salinity and buoyancy 
     2173      ! 
     2174      ! Local variables 
     2175      INTEGER  ::   ji, jj, jkb, jkb1 
     2176      REAL(wp) ::   zthermal, zbeta 
     2177      ! 
     2178      REAL(wp), PARAMETER ::   zlarge = -1e10_wp 
     2179      ! 
     2180      IF( ln_timing ) CALL timing_start('zdf_osm_eg') 
     2181      ! 
     2182      pdtdz(:,:) = zlarge 
     2183      pdsdz(:,:) = zlarge 
     2184      pdbdz(:,:) = zlarge 
     2185      ! 
     2186      DO_2D( 0, 0, 0, 0 ) 
     2187         IF ( kbase(ji,jj)+1 < mbkt(ji,jj) ) THEN 
     2188            zthermal = rab_n(ji,jj,1,jp_tem)   ! Ideally use ibld not 1?? 
     2189            zbeta    = rab_n(ji,jj,1,jp_sal) 
     2190            jkb = kbase(ji,jj) 
     2191            jkb1 = MIN( jkb + 1, mbkt(ji,jj) ) 
     2192            pdtdz(ji,jj) = -1.0_wp * ( ts(ji,jj,jkb1,jp_tem,Kmm) - ts(ji,jj,jkb,jp_tem,Kmm ) ) / e3w(ji,jj,jkb1,Kmm) 
     2193            pdsdz(ji,jj) = -1.0_wp * ( ts(ji,jj,jkb1,jp_sal,Kmm) - ts(ji,jj,jkb,jp_sal,Kmm ) ) / e3w(ji,jj,jkb1,Kmm) 
     2194            pdbdz(ji,jj) = grav * zthermal * pdtdz(ji,jj) - grav * zbeta * pdsdz(ji,jj) 
     2195         ELSE 
     2196            pdtdz(ji,jj) = 0.0_wp 
     2197            pdsdz(ji,jj) = 0.0_wp 
     2198            pdbdz(ji,jj) = 0.0_wp 
     2199         END IF 
     2200      END_2D 
     2201      ! 
     2202      IF( ln_timing ) CALL timing_stop('zdf_osm_eg') 
     2203      ! 
     2204   END SUBROUTINE zdf_osm_external_gradients 
     2205 
     2206   SUBROUTINE zdf_osm_calculate_dhdt( ldconv, ldshear, k_ddh, pdhdt, phbl, pdh, pwb_ent, pwb_min,   & 
     2207      &                               pdb_bl, pdbdz_bl_ext, pdb_ml, pwb_fk_b, pwb_fk, pvel_mle ) 
     2208      !!--------------------------------------------------------------------- 
     2209      !!                   ***  ROUTINE zdf_osm_calculate_dhdt  *** 
     2210      !! 
     2211      !! ** Purpose : Calculates the rate at which hbl changes. 
     2212      !! 
     2213      !! ** Method  : 
     2214      !! 
     2215      !!---------------------------------------------------------------------- 
     2216      LOGICAL,  DIMENSION(:,:),    INTENT(in   ) ::   ldconv         ! BL stability flags 
     2217      LOGICAL,  DIMENSION(:,:),    INTENT(in   ) ::   ldshear        ! Shear layers 
     2218      INTEGER,  DIMENSION(:,:),    INTENT(in   ) ::   k_ddh          ! k_ddh=0: active shear layer 
     2219      !                                                              ! k_ddh=1: shear layer not active 
     2220      !                                                              ! k_ddh=2: shear production low 
     2221      REAL(wp), DIMENSION(A2D(0)), INTENT(  out) ::   pdhdt          ! Rate of change of hbl 
     2222      REAL(wp), DIMENSION(:,:),    INTENT(in   ) ::   phbl           ! BL depth 
     2223      REAL(wp), DIMENSION(:,:),    INTENT(in   ) ::   pdh            ! Pycnocline depth 
     2224      REAL(wp), DIMENSION(:,:),    INTENT(in   ) ::   pwb_ent        ! Buoyancy entrainment flux 
     2225      REAL(wp), DIMENSION(:,:),    INTENT(in   ) ::   pwb_min 
     2226      REAL(wp), DIMENSION(:,:),    INTENT(in   ) ::   pdb_bl         ! Buoyancy diff. between BL average and basal value 
     2227      REAL(wp), DIMENSION(A2D(0)), INTENT(in   ) ::   pdbdz_bl_ext   ! external buoyancy gradients 
     2228      REAL(wp), DIMENSION(:,:),    INTENT(in   ) ::   pdb_ml         ! Buoyancy diff. between mixed-layer average and basal value 
     2229      REAL(wp), DIMENSION(:,:),    INTENT(  out) ::   pwb_fk_b       ! MLE buoyancy flux averaged over OSBL 
     2230      REAL(wp), DIMENSION(:,:),    INTENT(in   ) ::   pwb_fk         ! Max MLE buoyancy flux 
     2231      REAL(wp), DIMENSION(:,:),    INTENT(in   ) ::   pvel_mle       ! Vvelocity scale for dhdt with stable ML and FK 
     2232      ! 
     2233      ! Local variables 
     2234      INTEGER  ::   jj, ji 
     2235      REAL(wp) ::   zgamma_b_nd, zgamma_dh_nd, zpert, zpsi, zari 
     2236      REAL(wp) ::   zvel_max, zddhdt 
     2237      ! 
     2238      REAL(wp), PARAMETER ::   zzeta_m  = 0.3_wp 
     2239      REAL(wp), PARAMETER ::   zgamma_c = 2.0_wp 
     2240      REAL(wp), PARAMETER ::   zdhoh    = 0.1_wp 
     2241      REAL(wp), PARAMETER ::   zalpha_b = 0.3_wp 
     2242      REAL(wp), PARAMETER ::   a_ddh    = 2.5_wp, a_ddh_2 = 3.5   ! Also in pycnocline_depth 
     2243      REAL(wp), PARAMETER ::   zlarge   = -1e10_wp 
     2244      ! 
     2245      IF( ln_timing ) CALL timing_start('zdf_osm_cd') 
     2246      ! 
     2247      pdhdt(:,:)    = zlarge 
     2248      pwb_fk_b(:,:) = zlarge 
     2249      ! 
     2250      DO_2D( 0, 0, 0, 0 ) 
     2251         ! 
     2252         IF ( ldshear(ji,jj) ) THEN 
     2253            ! 
     2254            IF ( ldconv(ji,jj) ) THEN   ! Convective 
     2255               ! 
     2256               IF ( ln_osm_mle ) THEN 
     2257                  IF ( hmle(ji,jj) > hbl(ji,jj) ) THEN   ! Fox-Kemper buoyancy flux average over OSBL 
     2258                     pwb_fk_b(ji,jj) = pwb_fk(ji,jj) * ( 1.0_wp + hmle(ji,jj) / ( 6.0_wp * hbl(ji,jj) ) *   & 
     2259                        &                                         ( -1.0_wp + ( 1.0_wp - 2.0_wp * hbl(ji,jj) / hmle(ji,jj) )**3 ) ) 
     2260                  ELSE 
     2261                     pwb_fk_b(ji,jj) = 0.5_wp * pwb_fk(ji,jj) * hmle(ji,jj) / hbl(ji,jj) 
     2262                  ENDIF 
     2263                  zvel_max = ( svstr(ji,jj)**3 + 0.5_wp * swstrc(ji,jj)**3 )**p2third / hbl(ji,jj) 
     2264                  IF ( ( pwb_ent(ji,jj) + 2.0_wp * pwb_fk_b(ji,jj) ) < 0.0_wp ) THEN   ! OSBL is deepening, 
     2265                     !                                                                 !    entrainment > restratification 
     2266                     IF ( pdb_bl(ji,jj) > 1e-15_wp ) THEN 
     2267                        zgamma_b_nd = MAX( pdbdz_bl_ext(ji,jj), 0.0_wp ) * pdh(ji,jj) / ( zvel_max + MAX( pdb_bl(ji,jj), 1e-15_wp ) ) 
     2268                        zpsi = ( 1.0_wp - 0.5_wp * pdh(ji,jj) / phbl(ji,jj) ) *   & 
     2269                           &   ( swb0(ji,jj) - MIN( ( pwb_min(ji,jj) + 2.0_wp * pwb_fk_b(ji,jj) ), 0.0_wp ) ) * pdh(ji,jj) / phbl(ji,jj) 
     2270                        zpsi = zpsi + 1.75_wp * ( 1.0_wp - 0.5_wp * pdh(ji,jj) / phbl(ji,jj) ) *   & 
     2271                           &   ( pdh(ji,jj) / phbl(ji,jj) + zgamma_b_nd ) * MIN( ( pwb_min(ji,jj) + 2.0_wp * pwb_fk_b(ji,jj) ), 0.0_wp ) 
     2272                        zpsi = zalpha_b * MAX( zpsi, 0.0_wp ) 
     2273                        pdhdt(ji,jj) = -1.0_wp * ( pwb_ent(ji,jj) + 2.0_wp * pwb_fk_b(ji,jj) ) /   & 
     2274                           &                      ( zvel_max + MAX( pdb_bl(ji,jj), 1e-15_wp ) ) +   & 
     2275                           &            zpsi / ( zvel_max + MAX( pdb_bl(ji,jj), 1e-15_wp ) ) 
     2276#ifdef key_osm_debug 
     2277                        IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN 
     2278                           WRITE(narea+100,'(a,g11.3)')'Inside 1st major loop of zdf_osm_calculate_dhdt, OSBL is deepening, entrainment > restratification:  zdhdt=',pdhdt(ji,jj) 
     2279                           WRITE(narea+100,'(3(a,g11.3))') '  zpsi=',zpsi, '  zgamma_b_nd=', zgamma_b_nd, '  zdh=', pdh(ji,jj) 
     2280                           FLUSH(narea+100) 
     2281                        END IF 
     2282#endif 
     2283                        IF ( k_ddh(ji,jj) == 1 ) THEN 
     2284                           IF ( ( swstrc(ji,jj) / svstr(ji,jj) )**3 <= 0.5_wp ) THEN 
     2285                              zari = MIN( 1.5_wp * pdb_bl(ji,jj) /                                                   & 
     2286                                 &        ( phbl(ji,jj) * ( MAX( pdbdz_bl_ext(ji,jj), 0.0_wp ) +                     & 
     2287                                 &                               pdb_bl(ji,jj)**2 / MAX( 4.5_wp * svstr(ji,jj)**2,   & 
     2288                                 &                                                       1e-12_wp ) ) ), 0.2_wp ) 
     2289                           ELSE 
     2290                              zari = MIN( 1.5_wp * pdb_bl(ji,jj) /                                                    & 
     2291                                 &        ( phbl(ji,jj) * ( MAX( pdbdz_bl_ext(ji,jj), 0.0_wp ) +                      & 
     2292                                 &                               pdb_bl(ji,jj)**2 / MAX( 4.5_wp * swstrc(ji,jj)**2,   & 
     2293                                 &                                                       1e-12_wp ) ) ), 0.2_wp ) 
     2294                           ENDIF 
     2295                           ! Relaxation to dh_ref = zari * hbl 
     2296                           zddhdt = -1.0_wp * a_ddh_2 * ( 1.0_wp - pdh(ji,jj) / ( zari * phbl(ji,jj) ) ) * pwb_ent(ji,jj) /   & 
     2297                              &     ( zvel_max + MAX( pdb_bl(ji,jj), 1e-15_wp ) ) 
     2298#ifdef key_osm_debug 
     2299                           IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN 
     2300                              WRITE(narea+100,'(a,g11.3)')'Inside 1st major loop of zdf_osm_calculate_dhdt,j_ddh(ji,jj) == 1:  zari=',zari 
     2301                              FLUSH(narea+100) 
     2302                           END IF 
     2303#endif 
     2304                        ELSE IF ( k_ddh(ji,jj) == 0 ) THEN   ! Growing shear layer 
     2305                           zddhdt = -1.0_wp * a_ddh * ( 1.0_wp - 1.6_wp * pdh(ji,jj) / phbl(ji,jj) ) * pwb_ent(ji,jj) /   & 
     2306                              &     ( zvel_max + MAX( pdb_bl(ji,jj), 1e-15_wp ) ) 
     2307                           zddhdt = EXP( -4.0_wp * ABS( ff_t(ji,jj) ) * phbl(ji,jj) / MAX(sustar(ji,jj), 1e-8_wp ) ) * zddhdt 
     2308                        ELSE 
     2309                           zddhdt = 0.0_wp 
     2310                        ENDIF   ! k_ddh 
     2311                        pdhdt(ji,jj) = pdhdt(ji,jj) + zalpha_b * ( 1.0_wp - 0.5_wp * pdh(ji,jj) / phbl(ji,jj) ) *   & 
     2312                           &                            pdb_ml(ji,jj) * MAX( zddhdt, 0.0_wp ) /   & 
     2313                           &                            ( zvel_max + MAX( pdb_bl(ji,jj), 1e-15_wp ) ) 
     2314                     ELSE   ! pdb_bl >0 
     2315                        pdhdt(ji,jj) = -1.0_wp * ( pwb_ent(ji,jj) + 2.0_wp * pwb_fk_b(ji,jj) ) /  MAX( zvel_max, 1e-15_wp ) 
     2316                     ENDIF 
     2317                  ELSE   ! pwb_min + 2*pwb_fk_b < 0 
     2318                     ! OSBL shoaling due to restratification flux. This is the velocity defined in Fox-Kemper et al (2008) 
     2319                     pdhdt(ji,jj) = -1.0_wp * MIN( pvel_mle(ji,jj), hbl(ji,jj) / 10800.0_wp ) 
     2320                  ENDIF 
     2321               ELSE   ! Fox-Kemper not used. 
     2322                  zvel_max = -1.0_wp * ( 1.0_wp + 1.0_wp * ( svstr(ji,jj)**3 + 0.5_wp * swstrc(ji,jj)**3 )**pthird *     & 
     2323                     &                                                         rn_Dt / hbl(ji,jj) ) * pwb_ent(ji,jj) /   & 
     2324                     &       MAX( ( svstr(ji,jj)**3 + 0.5_wp * swstrc(ji,jj)**3 )**pthird, epsln ) 
     2325                  pdhdt(ji,jj) = -1.0_wp * pwb_ent(ji,jj) / ( zvel_max + MAX( pdb_bl(ji,jj), 1e-15_wp ) ) 
     2326                  ! added ajgn 23 July as temporay fix 
     2327               ENDIF   ! ln_osm_mle 
     2328               ! 
     2329            ELSE   ! ldconv - Stable 
     2330               ! 
     2331               pdhdt(ji,jj) = ( 0.06_wp + 0.52_wp * shol(ji,jj) / 2.0_wp ) * svstr(ji,jj)**3 / hbl(ji,jj) + swbav(ji,jj) 
     2332               IF ( pdhdt(ji,jj) < 0.0_wp ) THEN   ! For long timsteps factor in brackets slows the rapid collapse of the OSBL 
     2333                  zpert = 2.0_wp * ( 1.0_wp + 0.0_wp * 2.0_wp * svstr(ji,jj) * rn_Dt / hbl(ji,jj) ) * svstr(ji,jj)**2 / hbl(ji,jj) 
     2334               ELSE 
     2335                  zpert = MAX( svstr(ji,jj)**2 / hbl(ji,jj), pdb_bl(ji,jj) ) 
     2336               ENDIF 
     2337               pdhdt(ji,jj) = 2.0_wp * pdhdt(ji,jj) / MAX( zpert, epsln ) 
     2338               pdhdt(ji,jj) = MAX( pdhdt(ji,jj), -1.0_wp * hbl(ji,jj) / 5400.0_wp ) 
     2339               ! 
     2340            ENDIF   ! ldconv 
     2341            ! 
     2342         ELSE   ! ldshear 
     2343            ! 
     2344            IF ( ldconv(ji,jj) ) THEN   ! Convective 
     2345               ! 
     2346               IF ( ln_osm_mle ) THEN 
     2347                  IF ( hmle(ji,jj) > hbl(ji,jj) ) THEN   ! Fox-Kemper buoyancy flux average over OSBL 
     2348                     pwb_fk_b(ji,jj) = pwb_fk(ji,jj) *                       & 
     2349                        ( 1.0_wp + hmle(ji,jj) / ( 6.0_wp * hbl(ji,jj) ) *   & 
     2350                        &          ( -1.0_wp + ( 1.0_wp - 2.0_wp * hbl(ji,jj) / hmle(ji,jj))**3) ) 
     2351                  ELSE 
     2352                     pwb_fk_b(ji,jj) = 0.5_wp * pwb_fk(ji,jj) * hmle(ji,jj) / hbl(ji,jj) 
     2353                  ENDIF 
     2354                  zvel_max = ( swstrl(ji,jj)**3 + 0.5_wp * swstrc(ji,jj)**3 )**p2third / hbl(ji,jj) 
     2355                  IF ( ( pwb_ent(ji,jj) + 2.0_wp * pwb_fk_b(ji,jj) ) < 0.0_wp ) THEN   ! OSBL is deepening, 
     2356                     !                                                                 !    entrainment > restratification 
     2357                     IF ( pdb_bl(ji,jj) > 0.0_wp .AND. pdbdz_bl_ext(ji,jj) > 0.0_wp ) THEN 
     2358                        pdhdt(ji,jj) = -1.0_wp * ( pwb_ent(ji,jj) + 2.0_wp * pwb_fk_b(ji,jj) ) /   & 
     2359                           &            ( zvel_max + MAX( pdb_bl(ji,jj), 1e-15_wp ) ) 
     2360                     ELSE 
     2361                        pdhdt(ji,jj) = -1.0_wp * ( pwb_ent(ji,jj) + 2.0_wp * pwb_fk_b(ji,jj) ) / MAX( zvel_max, 1e-15_wp ) 
     2362                     ENDIF 
     2363                  ELSE   ! OSBL shoaling due to restratification flux. This is the velocity defined in Fox-Kemper et al (2008) 
     2364                     pdhdt(ji,jj) = -1.0_wp * MIN( pvel_mle(ji,jj), hbl(ji,jj) / 10800.0_wp ) 
     2365                  ENDIF 
     2366               ELSE   ! Fox-Kemper not used 
     2367                  zvel_max = -1.0_wp * pwb_ent(ji,jj) / MAX( ( svstr(ji,jj)**3 + 0.5_wp * swstrc(ji,jj)**3 )**pthird, epsln ) 
     2368                  pdhdt(ji,jj) = -1.0_wp * pwb_ent(ji,jj) / ( zvel_max + MAX( pdb_bl(ji,jj), 1e-15_wp ) ) 
     2369                  ! added ajgn 23 July as temporay fix 
     2370               ENDIF  ! ln_osm_mle 
     2371               ! 
     2372            ELSE                        ! Stable 
     2373               ! 
     2374               pdhdt(ji,jj) = ( 0.06_wp + 0.52_wp * shol(ji,jj) / 2.0_wp ) * svstr(ji,jj)**3 / hbl(ji,jj) + swbav(ji,jj) 
     2375               IF ( pdhdt(ji,jj) < 0.0_wp ) THEN 
     2376                  ! For long timsteps factor in brackets slows the rapid collapse of the OSBL 
     2377                  zpert = 2.0_wp * svstr(ji,jj)**2 / hbl(ji,jj) 
     2378               ELSE 
     2379                  zpert = MAX( svstr(ji,jj)**2 / hbl(ji,jj), pdb_bl(ji,jj) ) 
     2380               ENDIF 
     2381               pdhdt(ji,jj) = 2.0_wp * pdhdt(ji,jj) / MAX(zpert, epsln) 
     2382               pdhdt(ji,jj) = MAX( pdhdt(ji,jj), -1.0_wp * hbl(ji,jj) / 5400.0_wp ) 
     2383               ! 
     2384            ENDIF  ! ldconv 
     2385            ! 
     2386         ENDIF ! ldshear 
     2387#ifdef key_osm_debug 
     2388         IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN 
     2389            WRITE(narea+100,'(4(a,g11.3))')'end of 1st major loop of zdf_osm_calculate_dhdt:  zdhdt=',pdhdt(ji,jj), & 
     2390               &  '  zpert=', zpert, '  zddhdt=', zddhdt, '  zvel_max=', zvel_max 
     2391            IF ( ln_osm_mle ) THEN 
     2392               WRITE(narea+100,'(3(a,g11.3),/)') 'zvel_mle=',pvel_mle(ji,jj), ' zwb_fk_b=', pwb_fk_b(ji,jj), & 
     2393                  & '  zwb_ent + 2*zwb_fk_b =', pwb_ent(ji,jj) + 2.0 * pwb_fk_b(ji,jj) 
     2394               FLUSH(narea+100) 
     2395            END IF 
     2396         END IF 
     2397#endif 
     2398         ! 
     2399      END_2D 
     2400      ! 
     2401      IF( ln_timing ) CALL timing_stop('zdf_osm_cd') 
     2402      ! 
     2403   END SUBROUTINE zdf_osm_calculate_dhdt 
     2404 
     2405   SUBROUTINE zdf_osm_timestep_hbl( Kmm, kbld, kmld, ldconv, ldpyc, pdhdt, phbl, phbl_t, pt_bl, ps_bl, pwb_ent, pwb_fk_b ) 
     2406      !!--------------------------------------------------------------------- 
     2407      !!                ***  ROUTINE zdf_osm_timestep_hbl  *** 
     2408      !! 
     2409      !! ** Purpose : Increments hbl. 
     2410      !! 
     2411      !! ** Method  : If the change in hbl exceeds one model level the change is 
     2412      !!              is calculated by moving down the grid, changing the 
     2413      !!              buoyancy jump. This is to ensure that the change in hbl 
     2414      !!              does not overshoot a stable layer. 
     2415      !! 
     2416      !!---------------------------------------------------------------------- 
     2417      INTEGER,                     INTENT(in   ) ::   Kmm        ! Ocean time-level index 
     2418      INTEGER,  DIMENSION(:,:),    INTENT(inout) ::   kbld       ! BL base layer 
     2419      INTEGER,  DIMENSION(:,:),    INTENT(in   ) ::   kmld       ! ML base layer 
     2420      LOGICAL,  DIMENSION(:,:),    INTENT(in   ) ::   ldconv     ! BL stability flags 
     2421      LOGICAL,  DIMENSION(:,:),    INTENT(inout) ::   ldpyc      ! Pycnocline flags 
     2422      REAL(wp), DIMENSION(A2D(0)), INTENT(inout) ::   pdhdt      ! Rates of change of hbl 
     2423      REAL(wp), DIMENSION(:,:),    INTENT(inout) ::   phbl       ! BL depth 
     2424      REAL(wp), DIMENSION(:,:),    INTENT(in   ) ::   phbl_t     ! BL depth 
     2425      REAL(wp), DIMENSION(:,:),    INTENT(in   ) ::   pt_bl      ! Temperature average over the depth of the blayer 
     2426      REAL(wp), DIMENSION(:,:),    INTENT(in   ) ::   ps_bl      ! Salinity average over the depth of the blayer 
     2427      REAL(wp), DIMENSION(:,:),    INTENT(in   ) ::   pwb_ent    ! Buoyancy entrainment flux 
     2428      REAL(wp), DIMENSION(:,:),    INTENT(in   ) ::   pwb_fk_b   ! MLE buoyancy flux averaged over OSBL 
     2429      ! 
     2430      ! Local variables 
     2431      INTEGER  :: jk, jj, ji, jm 
     2432      REAL(wp) :: zhbl_s, zvel_max, zdb 
     2433      REAL(wp) :: zthermal, zbeta 
     2434      ! 
     2435      IF( ln_timing ) CALL timing_start('zdf_osm_th') 
     2436      ! 
     2437      DO_2D( 0, 0, 0, 0 ) 
     2438#ifdef key_osm_debug 
     2439         IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN 
     2440            WRITE(narea+100,'(2(a,i7))')'start of zdf_osm_timestep_hbl: old ibld=',kmld(ji,jj),' trial ibld=', kbld(ji,jj) 
     2441            FLUSH(narea+100) 
     2442         END IF 
     2443#endif 
     2444         IF ( kbld(ji,jj) - kmld(ji,jj) > 1 ) THEN 
     2445            ! 
     2446            ! If boundary layer changes by more than one level, need to check for stable layers between initial and final depths. 
     2447            ! 
     2448            zhbl_s   = hbl(ji,jj) 
     2449            jm       = kmld(ji,jj) 
     2450            zthermal = rab_n(ji,jj,1,jp_tem) 
     2451            zbeta    = rab_n(ji,jj,1,jp_sal) 
     2452            ! 
     2453            IF ( ldconv(ji,jj) ) THEN   ! Unstable 
     2454               ! 
     2455               IF( ln_osm_mle ) THEN 
     2456                  zvel_max = ( swstrl(ji,jj)**3 + swstrc(ji,jj)**3 )**p2third / hbl(ji,jj) 
     2457               ELSE 
     2458                  zvel_max = -1.0_wp * ( 1.0_wp + 1.0_wp * ( svstr(ji,jj)**3 + 0.5_wp * swstrc(ji,jj)**3 )**pthird * rn_Dt /   & 
     2459                     &                                     hbl(ji,jj) ) * pwb_ent(ji,jj) /                                     & 
     2460                     &       ( svstr(ji,jj)**3 + 0.5_wp * swstrc(ji,jj)**3 )**pthird 
     2461               ENDIF 
     2462#ifdef key_osm_debug 
     2463               IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN 
     2464                  WRITE(narea+100,'(a,g11.3)')'In zdf_osm_timestep_hbl, ibld - imld > 1, lconv=T: zvel_max=',zvel_max 
     2465                  FLUSH(narea+100) 
     2466               END IF 
     2467#endif 
     2468               DO jk = kmld(ji,jj), kbld(ji,jj) 
     2469                  zdb = MAX( grav * ( zthermal * ( pt_bl(ji,jj) - ts(ji,jj,jm,jp_tem,Kmm) ) -   & 
     2470                     &                zbeta    * ( ps_bl(ji,jj) - ts(ji,jj,jm,jp_sal,Kmm) ) ), 0.0_wp ) + zvel_max 
     2471                  ! 
     2472                  IF ( ln_osm_mle ) THEN 
     2473                     zhbl_s = zhbl_s + MIN( rn_Dt * ( ( -1.0_wp * pwb_ent(ji,jj) - 2.0_wp * pwb_fk_b(ji,jj) ) / zdb ) /   & 
     2474                        &                   REAL( kbld(ji,jj) - kmld(ji,jj), KIND=wp ), e3w(ji,jj,jm,Kmm) ) 
     2475                  ELSE 
     2476                     zhbl_s = zhbl_s + MIN( rn_Dt * ( -1.0_wp * pwb_ent(ji,jj) / zdb ) /   & 
     2477                        &                   REAL( kbld(ji,jj) - kmld(ji,jj), KIND=wp ), e3w(ji,jj,jm,Kmm) ) 
     2478                  ENDIF 
     2479                  !                    zhbl_s = MIN(zhbl_s,  gdepw(ji,jj, mbkt(ji,jj) + 1,Kmm) - depth_tol) 
     2480                  IF ( zhbl_s >= gdepw(ji,jj,mbkt(ji,jj) + 1,Kmm) ) THEN 
     2481                     zhbl_s = MIN( zhbl_s, gdepw(ji,jj, mbkt(ji,jj) + 1, Kmm ) - depth_tol ) 
     2482                     ldpyc(ji,jj) = .FALSE. 
     2483                  ENDIF 
     2484                  IF ( zhbl_s >= gdepw(ji,jj,jm+1,Kmm) ) jm = jm + 1 
     2485#ifdef key_osm_debug 
     2486                  IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN 
     2487                     WRITE(narea+100,'(2(a,i7))')' jk=',jk,' jm=', jm 
     2488                     WRITE(narea+100,'(2(a,g11.3),a,l7)')'zdb=',zdb,' zhbl_s=', zhbl_s,' lpyc=',ldpyc(ji,jj) 
     2489                     FLUSH(narea+100) 
     2490                  END IF 
     2491#endif 
     2492               END DO 
     2493               hbl(ji,jj)  = zhbl_s 
     2494               kbld(ji,jj) = jm 
     2495            ELSE   ! Stable 
     2496#ifdef key_osm_debug 
     2497               IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN 
     2498                  WRITE(narea+100,'(a)')'In zdf_osm_timestep_hbl, ibld - imld > 1, lconv=F' 
     2499                  FLUSH(narea+100) 
     2500               END IF 
     2501#endif 
     2502               DO jk = kmld(ji,jj), kbld(ji,jj) 
     2503                  zdb = MAX(  grav * ( zthermal * ( pt_bl(ji,jj) - ts(ji,jj,jm,jp_tem,Kmm) ) -               & 
     2504                     &                 zbeta    * ( ps_bl(ji,jj) - ts(ji,jj,jm,jp_sal,Kmm) ) ), 0.0_wp ) +   & 
     2505                     &  2.0 * svstr(ji,jj)**2 / zhbl_s 
     2506                  ! 
     2507                  ! Alan is thuis right? I have simply changed hbli to hbl 
     2508                  shol(ji,jj)  = -1.0_wp * zhbl_s / ( ( svstr(ji,jj)**3 + epsln ) / swbav(ji,jj) ) 
     2509                  pdhdt(ji,jj) = -1.0_wp * ( swbav(ji,jj) - 0.04_wp / 2.0_wp * swstrl(ji,jj)**3 / zhbl_s -   & 
     2510                     &                       0.15_wp / 2.0_wp * ( 1.0_wp - EXP( -1.5_wp * sla(ji,jj) ) ) *   & 
     2511                     &                                 sustar(ji,jj)**3 / zhbl_s ) *                         & 
     2512                     &           ( 0.725_wp + 0.225_wp * EXP( -7.5_wp * shol(ji,jj) ) ) 
     2513                  pdhdt(ji,jj) = pdhdt(ji,jj) + swbav(ji,jj) 
     2514                  zhbl_s = zhbl_s + MIN( pdhdt(ji,jj) / zdb * rn_Dt / REAL( kbld(ji,jj) - kmld(ji,jj), KIND=wp ),   & 
     2515                     &                   e3w(ji,jj,jm,Kmm) ) 
     2516                   
     2517                  !                    zhbl_s = MIN(zhbl_s, gdepw(ji,jj, mbkt(ji,jj) + 1,Kmm) - depth_tol) 
     2518                  IF ( zhbl_s >= mbkt(ji,jj) + 1 ) THEN 
     2519                     zhbl_s      = MIN( zhbl_s,  gdepw(ji,jj,mbkt(ji,jj)+1,Kmm) - depth_tol ) 
     2520                     ldpyc(ji,jj) = .FALSE. 
     2521                  ENDIF 
     2522                  IF ( zhbl_s >= gdepw(ji,jj,jm,Kmm) ) jm = jm + 1 
     2523#ifdef key_osm_debug 
     2524                  IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN 
     2525                     WRITE(narea+100,'(2(a,i7))')' jk=',jk,' jm=', jm 
     2526                     WRITE(narea+100,'(4(a,g11.3),a,l7)')'zdb=',zdb,' shol',shol(ji,jj),' zdhdt',pdhdt(ji,jj),' zhbl_s=', zhbl_s,' lpyc=',ldpyc(ji,jj) 
     2527                     FLUSH(narea+100) 
     2528                  END IF 
     2529#endif 
     2530               END DO 
     2531            ENDIF   ! IF ( ldconv ) 
     2532            hbl(ji,jj)  = MAX( zhbl_s, gdepw(ji,jj,4,Kmm) ) 
     2533            kbld(ji,jj) = MAX( jm, 4 ) 
     2534         ELSE 
     2535            ! change zero or one model level. 
     2536            hbl(ji,jj) = MAX( phbl_t(ji,jj), gdepw(ji,jj,4,Kmm) ) 
     2537         ENDIF 
     2538         phbl(ji,jj) = gdepw(ji,jj,kbld(ji,jj),Kmm) 
     2539#ifdef key_osm_debug 
     2540         IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN 
     2541            WRITE(narea+100,'(2(a,g11.3),a,i7,/)')'end of zdf_osm_timestep_hbl: hbl=', hbl(ji,jj),' zhbl=', phbl(ji,jj),' ibld=', kbld(ji,jj) 
     2542            FLUSH(narea+100) 
     2543         END IF 
     2544#endif 
     2545      END_2D 
     2546      ! 
     2547      IF( ln_timing ) CALL timing_stop('zdf_osm_th') 
     2548      ! 
     2549   END SUBROUTINE zdf_osm_timestep_hbl 
     2550 
    25242551   SUBROUTINE zdf_osm_diffusivity_viscosity( Kbb, Kmm, kbld, kmld, ldconv, ldshear, ldpyc, ldcoup, k_ddh, pdiffut, pviscos,       & 
    25252552      &                                      phbl, phml, pdh, pdhdt, pshear, pb_bl, pu_bl, pv_bl, pb_ml, pu_ml, pv_ml, pwb_ent,   & 
     
    25342561      !! 
    25352562      !!---------------------------------------------------------------------- 
    2536       INTEGER,                    INTENT(in   ) ::   Kbb, Kmm       ! Ocean time-level indices 
    2537       INTEGER,  DIMENSION(:,:),   INTENT(in   ) ::   kbld           ! BL base layer 
    2538       INTEGER,  DIMENSION(:,:),   INTENT(in   ) ::   kmld           ! ML base layer 
    2539       LOGICAL,  DIMENSION(:,:),   INTENT(in   ) ::   ldconv         ! BL stability flags 
    2540       LOGICAL,  DIMENSION(:,:),   INTENT(in   ) ::   ldshear        ! Shear layers 
    2541       LOGICAL,  DIMENSION(:,:),   INTENT(in   ) ::   ldpyc          ! Pycnocline flags 
    2542       LOGICAL,  DIMENSION(:,:),   INTENT(in   ) ::   ldcoup         ! Coupling to bottom 
    2543       INTEGER,  DIMENSION(:,:),   INTENT(in   ) ::   k_ddh          ! Type of shear layer 
    2544       REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   pdiffut        ! t-diffusivity 
    2545       REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   pviscos        ! Viscosity 
    2546       REAL(wp), DIMENSION(:,:),   INTENT(in   ) ::   phbl           ! BL depth 
    2547       REAL(wp), DIMENSION(:,:),   INTENT(in   ) ::   phml           ! ML depth 
    2548       REAL(wp), DIMENSION(:,:),   INTENT(in   ) ::   pdh            ! Pycnocline depth 
    2549       REAL(wp), DIMENSION(:,:),  INTENT(in   ) ::   pdhdt          ! BL depth tendency 
    2550       REAL(wp), DIMENSION(:,:),   INTENT(in   ) ::   pshear         ! Shear production 
    2551       REAL(wp), DIMENSION(:,:),   INTENT(in   ) ::   pb_bl          ! Buoyancy average over the depth of the boundary layer 
    2552       REAL(wp), DIMENSION(:,:),   INTENT(in   ) ::   pu_bl          ! Velocity average over the depth of the boundary layer 
    2553       REAL(wp), DIMENSION(:,:),   INTENT(in   ) ::   pv_bl          ! Velocity average over the depth of the boundary layer 
    2554       REAL(wp), DIMENSION(:,:),   INTENT(in   ) ::   pb_ml          ! Buoyancy average over the depth of the mixed layer 
    2555       REAL(wp), DIMENSION(:,:),   INTENT(in   ) ::   pu_ml          ! Velocity average over the depth of the mixed layer 
    2556       REAL(wp), DIMENSION(:,:),   INTENT(in   ) ::   pv_ml          ! Velocity average over the depth of the mixed layer 
    2557       REAL(wp), DIMENSION(:,:),   INTENT(in   ) ::   pwb_ent        ! Buoyancy entrainment flux 
    2558       REAL(wp), DIMENSION(:,:),   INTENT(in   ) ::   pwb_min 
     2563      INTEGER,                     INTENT(in   ) ::   Kbb, Kmm       ! Ocean time-level indices 
     2564      INTEGER,  DIMENSION(:,:),    INTENT(in   ) ::   kbld           ! BL base layer 
     2565      INTEGER,  DIMENSION(:,:),    INTENT(in   ) ::   kmld           ! ML base layer 
     2566      LOGICAL,  DIMENSION(:,:),    INTENT(in   ) ::   ldconv         ! BL stability flags 
     2567      LOGICAL,  DIMENSION(:,:),    INTENT(in   ) ::   ldshear        ! Shear layers 
     2568      LOGICAL,  DIMENSION(:,:),    INTENT(in   ) ::   ldpyc          ! Pycnocline flags 
     2569      LOGICAL,  DIMENSION(:,:),    INTENT(in   ) ::   ldcoup         ! Coupling to bottom 
     2570      INTEGER,  DIMENSION(:,:),    INTENT(in   ) ::   k_ddh          ! Type of shear layer 
     2571      REAL(wp), DIMENSION(:,:,:),  INTENT(inout) ::   pdiffut        ! t-diffusivity 
     2572      REAL(wp), DIMENSION(:,:,:),  INTENT(inout) ::   pviscos        ! Viscosity 
     2573      REAL(wp), DIMENSION(:,:),    INTENT(in   ) ::   phbl           ! BL depth 
     2574      REAL(wp), DIMENSION(:,:),    INTENT(in   ) ::   phml           ! ML depth 
     2575      REAL(wp), DIMENSION(:,:),    INTENT(in   ) ::   pdh            ! Pycnocline depth 
     2576      REAL(wp), DIMENSION(A2D(0)), INTENT(in   ) ::   pdhdt          ! BL depth tendency 
     2577      REAL(wp), DIMENSION(:,:),    INTENT(in   ) ::   pshear         ! Shear production 
     2578      REAL(wp), DIMENSION(:,:),    INTENT(in   ) ::   pb_bl          ! Buoyancy average over the depth of the boundary layer 
     2579      REAL(wp), DIMENSION(:,:),    INTENT(in   ) ::   pu_bl          ! Velocity average over the depth of the boundary layer 
     2580      REAL(wp), DIMENSION(:,:),    INTENT(in   ) ::   pv_bl          ! Velocity average over the depth of the boundary layer 
     2581      REAL(wp), DIMENSION(:,:),    INTENT(in   ) ::   pb_ml          ! Buoyancy average over the depth of the mixed layer 
     2582      REAL(wp), DIMENSION(:,:),    INTENT(in   ) ::   pu_ml          ! Velocity average over the depth of the mixed layer 
     2583      REAL(wp), DIMENSION(:,:),    INTENT(in   ) ::   pv_ml          ! Velocity average over the depth of the mixed layer 
     2584      REAL(wp), DIMENSION(:,:),    INTENT(in   ) ::   pwb_ent        ! Buoyancy entrainment flux 
     2585      REAL(wp), DIMENSION(:,:),    INTENT(in   ) ::   pwb_min 
    25592586      ! 
    25602587      ! Local variables 
     
    27952822   SUBROUTINE zdf_osm_fgr_terms( Kmm, kbld, kmld, kp_ext, ldconv, ldpyc, k_ddh, phbl, phml, pdh, pdhdt, pshear,             & 
    27962823      &                          pdt_bl, pds_bl, pdb_bl, pdu_bl, pdv_bl, pdt_ml, pds_ml, pdb_ml, pdu_ml, pdv_ml,                  & 
    2797       &                          pdtdz_bl_ext, pdsdz_bl_ext, pdbdz_bl_ext, pdbdz_pyc, palpha_pyc, pdiffut, pviscos ) 
     2824      &                          pdtdz_bl_ext, pdsdz_bl_ext, pdbdz_pyc, palpha_pyc, pdiffut, pviscos ) 
    27982825      !!--------------------------------------------------------------------- 
    27992826      !!                 ***  ROUTINE zdf_osm_fgr_terms *** 
     
    28042831      !! 
    28052832      !!---------------------------------------------------------------------- 
    2806       INTEGER,                    INTENT(in   ) ::   Kmm            ! Time-level index 
    2807       INTEGER,  DIMENSION(:,:),   INTENT(in   ) ::   kbld           ! BL base layer 
    2808       INTEGER,  DIMENSION(:,:),   INTENT(in   ) ::   kmld           ! ML base layer 
    2809       INTEGER,  DIMENSION(:,:),   INTENT(in   ) ::   kp_ext         ! Offset for external level 
    2810       LOGICAL,  DIMENSION(:,:),   INTENT(in   ) ::   ldconv         ! BL stability flags 
    2811       LOGICAL,  DIMENSION(:,:),   INTENT(in   ) ::   ldpyc          ! Pycnocline flags 
    2812       INTEGER,  DIMENSION(:,:),   INTENT(in   ) ::   k_ddh          ! Type of shear layer 
    2813       REAL(wp), DIMENSION(:,:),   INTENT(in   ) ::   phbl           ! BL depth 
    2814       REAL(wp), DIMENSION(:,:),   INTENT(in   ) ::   phml           ! ML depth 
    2815       REAL(wp), DIMENSION(:,:),   INTENT(in   ) ::   pdh            ! Pycnocline depth 
    2816       REAL(wp), DIMENSION(:,:),   INTENT(in   ) ::   pdhdt          ! BL depth tendency 
    2817       REAL(wp), DIMENSION(:,:),   INTENT(in   ) ::   pshear         ! Shear production 
    2818       REAL(wp), DIMENSION(:,:),   INTENT(in   ) ::   pdt_bl         ! Temperature diff. between BL average and basal value 
    2819       REAL(wp), DIMENSION(:,:),   INTENT(in   ) ::   pds_bl         ! Salinity diff. between BL average and basal value 
    2820       REAL(wp), DIMENSION(:,:),   INTENT(in   ) ::   pdb_bl         ! Buoyancy diff. between BL average and basal value 
    2821       REAL(wp), DIMENSION(:,:),   INTENT(in   ) ::   pdu_bl         ! Velocity diff. (u) between BL average and basal value 
    2822       REAL(wp), DIMENSION(:,:),   INTENT(in   ) ::   pdv_bl         ! Velocity diff. (u) between BL average and basal value 
    2823       REAL(wp), DIMENSION(:,:),   INTENT(in   ) ::   pdt_ml         ! Temperature diff. between mixed-layer average and basal value 
    2824       REAL(wp), DIMENSION(:,:),   INTENT(in   ) ::   pds_ml         ! Salinity diff. between mixed-layer average and basal value 
    2825       REAL(wp), DIMENSION(:,:),   INTENT(in   ) ::   pdb_ml         ! Buoyancy diff. between mixed-layer average and basal value 
    2826       REAL(wp), DIMENSION(:,:),   INTENT(in   ) ::   pdu_ml         ! Velocity diff. (u) between mixed-layer average and basal value 
    2827       REAL(wp), DIMENSION(:,:),   INTENT(in   ) ::   pdv_ml         ! Velocity diff. (v) between mixed-layer average and basal value 
    2828       REAL(wp), DIMENSION(:,:),   INTENT(in   ) ::   pdtdz_bl_ext   ! External temperature gradients 
    2829       REAL(wp), DIMENSION(:,:),   INTENT(in   ) ::   pdsdz_bl_ext   ! External salinity gradients 
    2830       REAL(wp), DIMENSION(:,:),   INTENT(in   ) ::   pdbdz_bl_ext   ! External buoyancy gradients 
    2831       REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   pdbdz_pyc      ! Pycnocline buoyancy gradients 
    2832       REAL(wp), DIMENSION(:,:),   INTENT(in   ) ::   palpha_pyc     ! 
    2833       REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   pdiffut        ! t-diffusivity 
    2834       REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   pviscos        ! Viscosity 
     2833      INTEGER,                     INTENT(in   ) ::   Kmm            ! Time-level index 
     2834      INTEGER,  DIMENSION(:,:),    INTENT(in   ) ::   kbld           ! BL base layer 
     2835      INTEGER,  DIMENSION(:,:),    INTENT(in   ) ::   kmld           ! ML base layer 
     2836      INTEGER,  DIMENSION(:,:),    INTENT(in   ) ::   kp_ext         ! Offset for external level 
     2837      LOGICAL,  DIMENSION(:,:),    INTENT(in   ) ::   ldconv         ! BL stability flags 
     2838      LOGICAL,  DIMENSION(:,:),    INTENT(in   ) ::   ldpyc          ! Pycnocline flags 
     2839      INTEGER,  DIMENSION(:,:),    INTENT(in   ) ::   k_ddh          ! Type of shear layer 
     2840      REAL(wp), DIMENSION(:,:),    INTENT(in   ) ::   phbl           ! BL depth 
     2841      REAL(wp), DIMENSION(:,:),    INTENT(in   ) ::   phml           ! ML depth 
     2842      REAL(wp), DIMENSION(:,:),    INTENT(in   ) ::   pdh            ! Pycnocline depth 
     2843      REAL(wp), DIMENSION(A2D(0)), INTENT(in   ) ::   pdhdt          ! BL depth tendency 
     2844      REAL(wp), DIMENSION(:,:),    INTENT(in   ) ::   pshear         ! Shear production 
     2845      REAL(wp), DIMENSION(:,:),    INTENT(in   ) ::   pdt_bl         ! Temperature diff. between BL average and basal value 
     2846      REAL(wp), DIMENSION(:,:),    INTENT(in   ) ::   pds_bl         ! Salinity diff. between BL average and basal value 
     2847      REAL(wp), DIMENSION(:,:),    INTENT(in   ) ::   pdb_bl         ! Buoyancy diff. between BL average and basal value 
     2848      REAL(wp), DIMENSION(:,:),    INTENT(in   ) ::   pdu_bl         ! Velocity diff. (u) between BL average and basal value 
     2849      REAL(wp), DIMENSION(:,:),    INTENT(in   ) ::   pdv_bl         ! Velocity diff. (v) between BL average and basal value 
     2850      REAL(wp), DIMENSION(:,:),    INTENT(in   ) ::   pdt_ml         ! Temperature diff. between mixed-layer average and basal value 
     2851      REAL(wp), DIMENSION(:,:),    INTENT(in   ) ::   pds_ml         ! Salinity diff. between mixed-layer average and basal value 
     2852      REAL(wp), DIMENSION(:,:),    INTENT(in   ) ::   pdb_ml         ! Buoyancy diff. between mixed-layer average and basal value 
     2853      REAL(wp), DIMENSION(:,:),    INTENT(in   ) ::   pdu_ml         ! Velocity diff. (u) between mixed-layer average and basal value 
     2854      REAL(wp), DIMENSION(:,:),    INTENT(in   ) ::   pdv_ml         ! Velocity diff. (v) between mixed-layer average and basal value 
     2855      REAL(wp), DIMENSION(A2D(0)), INTENT(in   ) ::   pdtdz_bl_ext   ! External temperature gradients 
     2856      REAL(wp), DIMENSION(A2D(0)), INTENT(in   ) ::   pdsdz_bl_ext   ! External salinity gradients 
     2857      REAL(wp), DIMENSION(:,:,:),  INTENT(in   ) ::   pdbdz_pyc      ! Pycnocline buoyancy gradients 
     2858      REAL(wp), DIMENSION(:,:),    INTENT(in   ) ::   palpha_pyc     ! 
     2859      REAL(wp), DIMENSION(:,:,:),  INTENT(in   ) ::   pdiffut        ! t-diffusivity 
     2860      REAL(wp), DIMENSION(:,:,:),  INTENT(in   ) ::   pviscos        ! Viscosity 
    28352861      ! 
    28362862      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   z3ddz_pyc_1, z3ddz_pyc_2   ! Pycnocline gradient/shear profiles 
Note: See TracChangeset for help on using the changeset viewer.