Changeset 12516


Ignore:
Timestamp:
2020-03-05T16:07:35+01:00 (8 months ago)
Author:
cguiavarch
Message:

add modernised diahth w/o key

Location:
NEMO/branches/UKMO/NEMO_4.0_FKOSM/src/OCE
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/UKMO/NEMO_4.0_FKOSM/src/OCE/C1D/step_c1d.F90

    r10888 r12516  
    5858 
    5959                             indic = 0                ! reset to no error condition 
    60       IF( kstp == nit000 )   CALL iom_init( "nemo")   ! iom_put initialization (must be done after nemo_init for AGRIF+XIOS+OASIS) 
     60                             IF( kstp == nit000 )   THEN 
     61                                CALL iom_init( "nemo")   ! iom_put initialization (must be done after nemo_init for AGRIF+XIOS+OASIS) 
     62                                CALL dia_hth_init    ! extra ML depth diagnostics, thermocline depths 
     63                             END IF 
    6164      IF( kstp /= nit000 )   CALL day( kstp )         ! Calendar (day was already called at nit000 in day_init) 
    6265                             CALL iom_setkt( kstp - nit000 + 1, "nemo" )   ! say to iom that we are at time step kstp 
     
    8689      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    8790                         CALL dia_wri( kstp )       ! ocean model: outputs 
    88       IF( lk_diahth  )   CALL dia_hth( kstp )       ! Thermocline depth (20°C) 
     91      IF( ll_diahth  )   CALL dia_hth( kstp )       ! Thermocline depth (20°C) 
    8992 
    9093 
  • NEMO/branches/UKMO/NEMO_4.0_FKOSM/src/OCE/DIA/diahth.F90

    r10888 r12516  
    55   !!====================================================================== 
    66   !! History :  OPA  !  1994-09  (J.-P. Boulanger)  Original code 
    7    !!                 !  1996-11  (E. Guilyardi)  OPA8  
     7   !!                 !  1996-11  (E. Guilyardi)  OPA8 
    88   !!                 !  1997-08  (G. Madec)  optimization 
    9    !!                 !  1999-07  (E. Guilyardi)  hd28 + heat content  
     9   !!                 !  1999-07  (E. Guilyardi)  hd28 + heat content 
    1010   !!   NEMO     1.0  !  2002-06  (G. Madec)  F90: Free form and module 
    1111   !!            3.2  !  2009-07  (S. Masson) hc300 bugfix + cleaning + add new diag 
    12    !!---------------------------------------------------------------------- 
    13 #if   defined key_diahth 
    14    !!---------------------------------------------------------------------- 
    15    !!   'key_diahth' :                              thermocline depth diag. 
     12   !!            4.0.1!  2020-01  (G. Nurser) remove need for key lk_diahth 
    1613   !!---------------------------------------------------------------------- 
    1714   !!   dia_hth      : Compute varius diagnostics associated with the mixed layer 
     
    2421   USE lib_mpp         ! MPP library 
    2522   USE iom             ! I/O library 
    26    USE timing          ! preformance summary 
    2723 
    2824   IMPLICIT NONE 
     
    3026 
    3127   PUBLIC   dia_hth       ! routine called by step.F90 
    32    PUBLIC   dia_hth_alloc ! routine called by nemogcm.F90 
    33  
    34    LOGICAL , PUBLIC, PARAMETER ::   lk_diahth = .TRUE.    !: thermocline-20d depths flag 
    35     
    36    ! note: following variables should move to local variables once iom_put is always used  
    37    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hth    !: depth of the max vertical temperature gradient [m] 
    38    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hd20   !: depth of 20 C isotherm                         [m] 
    39    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hd28   !: depth of 28 C isotherm                         [m] 
    40    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   htc3   !: heat content of first 300 m                    [W] 
     28   PUBLIC   dia_hth_init  ! routine called by nemogcm.F90 
     29 
     30   LOGICAL, PUBLIC ::   ll_diahth   !: Compute further diagnostics of ML and thermocline depth 
    4131 
    4232   !!---------------------------------------------------------------------- 
     
    4737CONTAINS 
    4838 
    49    FUNCTION dia_hth_alloc() 
    50       !!--------------------------------------------------------------------- 
    51       INTEGER :: dia_hth_alloc 
    52       !!--------------------------------------------------------------------- 
     39 
     40   SUBROUTINE dia_hth( kt ) 
     41     !!--------------------------------------------------------------------- 
     42     !!                  ***  ROUTINE dia_hth  *** 
     43     !! 
     44     !! ** Purpose : Computes 
     45     !!      the mixing layer depth (turbocline): avt = 5.e-4 
     46     !!      the depth of strongest vertical temperature gradient 
     47     !!      the mixed layer depth with density     criteria: rho = rho(10m or surf) + 0.03(or 0.01) 
     48     !!      the mixed layer depth with temperature criteria: abs( tn - tn(10m) ) = 0.2 
     49     !!      the top of the thermochine: tn = tn(10m) - ztem2 
     50     !!      the pycnocline depth with density criteria equivalent to a temperature variation 
     51     !!                rho = rho10m + (dr/dT)(T,S,10m)*(-0.2 degC) 
     52     !!      the barrier layer thickness 
     53     !!      the maximal verical inversion of temperature and its depth max( 0, max of tn - tn(10m) ) 
     54     !!      the depth of the 20 degree isotherm (linear interpolation) 
     55     !!      the depth of the 28 degree isotherm (linear interpolation) 
     56     !!      the heat content of first 300 m 
     57     !! 
     58     !! ** Method : 
     59     !!------------------------------------------------------------------- 
     60     INTEGER, INTENT( in ) ::   kt      ! ocean time-step index 
     61     !! 
     62     INTEGER                          ::   ji, jj, jk            ! dummy loop arguments 
     63     INTEGER                          ::   iid, ilevel           ! temporary integers 
     64     INTEGER, DIMENSION(jpi,jpj) ::   ik20, ik28  ! levels 
     65     REAL(wp)                         ::   zavt5 = 5.e-4_wp      ! Kz criterion for the turbocline depth 
     66     REAL(wp)                         ::   zrho3 = 0.03_wp       ! density     criterion for mixed layer depth 
     67     REAL(wp)                         ::   zrho1 = 0.01_wp       ! density     criterion for mixed layer depth 
     68     REAL(wp)                         ::   ztem2 = 0.2_wp        ! temperature criterion for mixed layer depth 
     69     REAL(wp)                         ::   zthick_0, zcoef       ! temporary scalars 
     70     REAL(wp)                         ::   zztmp, zzdep          ! temporary scalars inside do loop 
     71     REAL(wp)                         ::   zu, zv, zw, zut, zvt  ! temporary workspace 
     72     REAL(wp), DIMENSION(jpi,jpj) ::   zabs2      ! MLD: abs( tn - tn(10m) ) = ztem2 
     73     REAL(wp), DIMENSION(jpi,jpj) ::   ztm2       ! Top of thermocline: tn = tn(10m) - ztem2 
     74     REAL(wp), DIMENSION(jpi,jpj) ::   zrho10_3   ! MLD: rho = rho10m + zrho3 
     75     REAL(wp), DIMENSION(jpi,jpj) ::   zpycn      ! pycnocline: rho = rho10m + (dr/dT)(T,S,10m)*(-0.2 degC) 
     76     REAL(wp), DIMENSION(jpi,jpj) ::   ztinv      ! max of temperature inversion 
     77     REAL(wp), DIMENSION(jpi,jpj) ::   zdepinv    ! depth of temperature inversion 
     78     REAL(wp), DIMENSION(jpi,jpj) ::   zrho0_3    ! MLD rho = rho(surf) = 0.03 
     79     REAL(wp), DIMENSION(jpi,jpj) ::   zrho0_1    ! MLD rho = rho(surf) = 0.01 
     80     REAL(wp), DIMENSION(jpi,jpj) ::   zmaxdzT    ! max of dT/dz 
     81     REAL(wp), DIMENSION(jpi,jpj) ::   zthick     ! vertical integration thickness 
     82     REAL(wp), DIMENSION(jpi,jpj) ::   zdelr      ! delta rho equivalent to deltaT = 0.2 
     83   ! note: following variables should move to local variables once iom_put is always used 
     84     REAL(wp), DIMENSION(jpi,jpj) ::   zhth    !: depth of the max vertical temperature gradient [m] 
     85     REAL(wp), DIMENSION(jpi,jpj) ::   zhd20   !: depth of 20 C isotherm                         [m] 
     86     REAL(wp), DIMENSION(jpi,jpj) ::   zhd28   !: depth of 28 C isotherm                         [m] 
     87     REAL(wp), DIMENSION(jpi,jpj) ::   zhtc3   !: heat content of first 300 m                    [W] 
     88 
     89     IF (iom_use("mlddzt") .OR. iom_use("mldr0_3") .OR. iom_use("mldr0_1")) THEN 
     90        ! ------------------------------------------------------------- ! 
     91        ! thermocline depth: strongest vertical gradient of temperature ! 
     92        ! turbocline depth (mixing layer depth): avt = zavt5            ! 
     93        ! MLD: rho = rho(1) + zrho3                                     ! 
     94        ! MLD: rho = rho(1) + zrho1                                     ! 
     95        ! ------------------------------------------------------------- ! 
     96        zmaxdzT(:,:) = 0._wp 
     97        IF( nla10 > 1 ) THEN 
     98           DO jj = 1, jpj 
     99              DO ji = 1, jpi 
     100                 zztmp = gdepw_n(ji,jj,mbkt(ji,jj)+1) 
     101                 zrho0_3(ji,jj) = zztmp 
     102                 zrho0_1(ji,jj) = zztmp 
     103                 zhth(ji,jj) = zztmp 
     104              END DO 
     105           END DO 
     106        ELSE IF (iom_use("mlddzt")) THEN 
     107           DO jj = 1, jpj 
     108              DO ji = 1, jpi 
     109                 zztmp = gdepw_n(ji,jj,mbkt(ji,jj)+1) 
     110                 zhth(ji,jj) = zztmp 
     111              END DO 
     112           END DO 
     113        ELSE 
     114           zhth(:,:) = 0._wp 
     115 
     116        ENDIF 
     117 
     118        DO jk = jpkm1, 2, -1   ! loop from bottom to 2 
     119           DO jj = 1, jpj 
     120              DO ji = 1, jpi 
     121                 ! 
     122                 zzdep = gdepw_n(ji,jj,jk) 
     123                 zztmp = ( tsn(ji,jj,jk-1,jp_tem) - tsn(ji,jj,jk,jp_tem) ) / zzdep * tmask(ji,jj,jk)   ! vertical gradient of temperature (dT/dz) 
     124                 zzdep = zzdep * tmask(ji,jj,1) 
     125 
     126                 IF( zztmp > zmaxdzT(ji,jj) ) THEN 
     127                    zmaxdzT(ji,jj) = zztmp   ;   zhth    (ji,jj) = zzdep                ! max and depth of dT/dz 
     128                 ENDIF 
     129 
     130                 IF( nla10 > 1 ) THEN 
     131                    zztmp = rhop(ji,jj,jk) - rhop(ji,jj,1)                             ! delta rho(1) 
     132                    IF( zztmp > zrho3 )          zrho0_3(ji,jj) = zzdep                ! > 0.03 
     133                    IF( zztmp > zrho1 )          zrho0_1(ji,jj) = zzdep                ! > 0.01 
     134                 ENDIF 
     135 
     136              END DO 
     137           END DO 
     138        END DO 
     139 
     140        IF (iom_use("mlddzt")) CALL iom_put( "mlddzt", zhth*tmask(:,:,1) )            ! depth of the thermocline 
     141        IF( nla10 > 1 ) THEN 
     142           IF (iom_use("mldr0_3")) CALL iom_put( "mldr0_3", zrho0_3*tmask(:,:,1) )   ! MLD delta rho(surf) = 0.03 
     143           IF (iom_use("mldr0_1")) CALL iom_put( "mldr0_1", zrho0_1*tmask(:,:,1) )   ! MLD delta rho(surf) = 0.01 
     144        ENDIF 
     145     ENDIF 
     146 
     147     IF (iom_use("mld_dt02") .OR. iom_use("topthdep") .OR. iom_use("mldr10_3") .OR.  & 
     148          & iom_use("pycndep") .OR. iom_use("tinv") .OR. iom_use("depti")) THEN 
     149        DO jj = 1, jpj 
     150           DO ji = 1, jpi 
     151              zztmp = gdepw_n(ji,jj,mbkt(ji,jj)+1) 
     152              zabs2   (ji,jj) = zztmp 
     153              ztm2    (ji,jj) = zztmp 
     154              zrho10_3(ji,jj) = zztmp 
     155              zpycn   (ji,jj) = zztmp 
     156           END DO 
     157        END DO 
     158        ztinv  (:,:) = 0._wp 
     159        zdepinv(:,:) = 0._wp 
     160 
     161        IF (iom_use("pycndep")) THEN 
     162           ! Preliminary computation 
     163           ! computation of zdelr = (dr/dT)(T,S,10m)*(-0.2 degC) 
     164           DO jj = 1, jpj 
     165              DO ji = 1, jpi 
     166                 IF( tmask(ji,jj,nla10) == 1. ) THEN 
     167                    zu  =  1779.50 + 11.250 * tsn(ji,jj,nla10,jp_tem) - 3.80   * tsn(ji,jj,nla10,jp_sal)                             & 
     168                         &                                              - 0.0745 * tsn(ji,jj,nla10,jp_tem) * tsn(ji,jj,nla10,jp_tem)   & 
     169                         &                                              - 0.0100 * tsn(ji,jj,nla10,jp_tem) * tsn(ji,jj,nla10,jp_sal) 
     170                    zv  =  5891.00 + 38.000 * tsn(ji,jj,nla10,jp_tem) + 3.00   * tsn(ji,jj,nla10,jp_sal)                             & 
     171                         &                                              - 0.3750 * tsn(ji,jj,nla10,jp_tem) * tsn(ji,jj,nla10,jp_tem) 
     172                    zut =    11.25 -  0.149 * tsn(ji,jj,nla10,jp_tem) - 0.01   * tsn(ji,jj,nla10,jp_sal) 
     173                    zvt =    38.00 -  0.750 * tsn(ji,jj,nla10,jp_tem) 
     174                    zw  = (zu + 0.698*zv) * (zu + 0.698*zv) 
     175                    zdelr(ji,jj) = ztem2 * (1000.*(zut*zv - zvt*zu)/zw) 
     176                 ELSE 
     177                    zdelr(ji,jj) = 0._wp 
     178                 ENDIF 
     179              END DO 
     180           END DO 
     181        ELSE 
     182           zdelr(:,:) = 0._wp 
     183        ENDIF 
     184 
     185        ! ------------------------------------------------------------- ! 
     186        ! MLD: abs( tn - tn(10m) ) = ztem2                              ! 
     187        ! Top of thermocline: tn = tn(10m) - ztem2                      ! 
     188        ! MLD: rho = rho10m + zrho3                                     ! 
     189        ! pycnocline: rho = rho10m + (dr/dT)(T,S,10m)*(-0.2 degC)       ! 
     190        ! temperature inversion: max( 0, max of tn - tn(10m) )          ! 
     191        ! depth of temperature inversion                                ! 
     192        ! ------------------------------------------------------------- ! 
     193        DO jk = jpkm1, nlb10, -1   ! loop from bottom to nlb10 
     194           DO jj = 1, jpj 
     195              DO ji = 1, jpi 
     196                 ! 
     197                 zzdep = gdepw_n(ji,jj,jk) * tmask(ji,jj,1) 
     198                 ! 
     199                 zztmp = tsn(ji,jj,nla10,jp_tem) - tsn(ji,jj,jk,jp_tem)  ! - delta T(10m) 
     200                 IF( ABS(zztmp) > ztem2 )      zabs2   (ji,jj) = zzdep   ! abs > 0.2 
     201                 IF(     zztmp  > ztem2 )      ztm2    (ji,jj) = zzdep   ! > 0.2 
     202                 zztmp = -zztmp                                          ! delta T(10m) 
     203                 IF( zztmp >  ztinv(ji,jj) ) THEN                        ! temperature inversion 
     204                    ztinv(ji,jj) = zztmp   ;   zdepinv (ji,jj) = zzdep   ! max value and depth 
     205                 ENDIF 
     206 
     207                 zztmp = rhop(ji,jj,jk) - rhop(ji,jj,nla10)              ! delta rho(10m) 
     208                 IF( zztmp > zrho3        )    zrho10_3(ji,jj) = zzdep   ! > 0.03 
     209                 IF( zztmp > zdelr(ji,jj) )    zpycn   (ji,jj) = zzdep   ! > equi. delta T(10m) - 0.2 
     210                 ! 
     211              END DO 
     212           END DO 
     213        END DO 
     214 
     215        IF (iom_use("mld_dt02")) CALL iom_put( "mld_dt02", zabs2*tmask(:,:,1)        )   ! MLD abs(delta t) - 0.2 
     216        IF (iom_use("topthdep")) CALL iom_put( "topthdep", ztm2*tmask(:,:,1)         )   ! T(10) - 0.2 
     217        IF (iom_use("mldr10_3")) CALL iom_put( "mldr10_3", zrho10_3*tmask(:,:,1)     )   ! MLD delta rho(10m) = 0.03 
     218        IF (iom_use("pycndep")) CALL iom_put( "pycndep" , zpycn*tmask(:,:,1)        )   ! MLD delta rho equi. delta T(10m) = 0.2 
     219        IF (iom_use("tinv")) CALL iom_put( "tinv"    , ztinv*tmask(:,:,1)        )   ! max. temp. inv. (t10 ref) 
     220        IF (iom_use("depti")) CALL iom_put( "depti"   , zdepinv*tmask(:,:,1)      )   ! depth of max. temp. inv. (t10 ref) 
     221     ENDIF 
     222 
     223     IF(iom_use("20d") .OR. iom_use("28d")) THEN 
     224        ! ----------------------------------- ! 
     225        ! search deepest level above 20C/28C  ! 
     226        ! ----------------------------------- ! 
     227        ik20(:,:) = 1 
     228        ik28(:,:) = 1 
     229        DO jk = 1, jpkm1   ! beware temperature is not always decreasing with depth => loop from top to bottom 
     230           DO jj = 1, jpj 
     231              DO ji = 1, jpi 
     232                 zztmp = tsn(ji,jj,jk,jp_tem) 
     233                 IF( zztmp >= 20. )   ik20(ji,jj) = jk 
     234                 IF( zztmp >= 28. )   ik28(ji,jj) = jk 
     235              END DO 
     236           END DO 
     237        END DO 
     238 
     239        ! --------------------------- ! 
     240        !  Depth of 20C/28C isotherm  ! 
     241        ! --------------------------- ! 
     242        DO jj = 1, jpj 
     243           DO ji = 1, jpi 
     244              ! 
     245              zzdep = gdepw_n(ji,jj,mbkt(ji,jj)+1)       ! depth of the oean bottom 
     246              ! 
     247              iid = ik20(ji,jj) 
     248              IF( iid /= 1 ) THEN 
     249                 zztmp =      gdept_n(ji,jj,iid  )   &                     ! linear interpolation 
     250                      &  + (    gdept_n(ji,jj,iid+1) - gdept_n(ji,jj,iid)                       )   & 
     251                      &  * ( 20.*tmask(ji,jj,iid+1) - tsn(ji,jj,iid,jp_tem)                       )   & 
     252                      &  / ( tsn(ji,jj,iid+1,jp_tem) - tsn(ji,jj,iid,jp_tem) + (1.-tmask(ji,jj,1)) ) 
     253                 zhd20(ji,jj) = MIN( zztmp , zzdep) * tmask(ji,jj,1)       ! bound by the ocean depth 
     254              ELSE 
     255                 zhd20(ji,jj) = 0._wp 
     256              ENDIF 
     257              ! 
     258              iid = ik28(ji,jj) 
     259              IF( iid /= 1 ) THEN 
     260                 zztmp =      gdept_n(ji,jj,iid  )   &                     ! linear interpolation 
     261                      &  + (    gdept_n(ji,jj,iid+1) - gdept_n(ji,jj,iid)                       )   & 
     262                      &  * ( 28.*tmask(ji,jj,iid+1) -    tsn(ji,jj,iid,jp_tem)                       )   & 
     263                      &  / (  tsn(ji,jj,iid+1,jp_tem) -    tsn(ji,jj,iid,jp_tem) + (1.-tmask(ji,jj,1)) ) 
     264                 zhd28(ji,jj) = MIN( zztmp , zzdep ) * tmask(ji,jj,1)      ! bound by the ocean depth 
     265              ELSE 
     266                 zhd28(ji,jj) = 0._wp 
     267              ENDIF 
     268 
     269           END DO 
     270        END DO 
     271        CALL iom_put( "20d", zhd20 )   ! depth of the 20 isotherm 
     272        CALL iom_put( "28d", zhd28 )   ! depth of the 28 isotherm 
     273     ENDIF 
     274 
     275     ! ----------------------------- ! 
     276     !  Heat content of first 300 m  ! 
     277     ! ----------------------------- ! 
     278     IF (iom_use("hc300")) THEN 
     279        ! find ilevel with (ilevel+1) the deepest W-level above 300m (we assume we can use e3t_1d to do this search...) 
     280        ilevel   = 0 
     281        zthick_0 = 0._wp 
     282        DO jk = 1, jpkm1 
     283           zthick_0 = zthick_0 + e3t_1d(jk) 
     284           IF( zthick_0 < 300. )   ilevel = jk 
     285        END DO 
     286        ! surface boundary condition 
     287        IF( ln_linssh ) THEN   ;   zthick(:,:) = sshn(:,:)   ;   zhtc3(:,:) = tsn(:,:,1,jp_tem) * sshn(:,:) * tmask(:,:,1) 
     288        ELSE                   ;   zthick(:,:) = 0._wp       ;   zhtc3(:,:) = 0._wp 
     289        ENDIF 
     290        ! integration down to ilevel 
     291        DO jk = 1, ilevel 
     292           zthick(:,:) = zthick(:,:) + e3t_n(:,:,jk) 
     293           zhtc3  (:,:) = zhtc3  (:,:) + e3t_n(:,:,jk) * tsn(:,:,jk,jp_tem) * tmask(:,:,jk) 
     294        END DO 
     295        ! deepest layer 
     296        zthick(:,:) = 300. - zthick(:,:)   !   remaining thickness to reach 300m 
     297        DO jj = 1, jpj 
     298           DO ji = 1, jpi 
     299              zhtc3(ji,jj) = zhtc3(ji,jj) + tsn(ji,jj,ilevel+1,jp_tem)                  & 
     300                   &                      * MIN( e3t_n(ji,jj,ilevel+1), zthick(ji,jj) ) * tmask(ji,jj,ilevel+1) 
     301           END DO 
     302        END DO 
     303        ! from temperature to heat contain 
     304        zcoef = rau0 * rcp 
     305        zhtc3(:,:) = zcoef * zhtc3(:,:) 
     306        CALL iom_put( "hc300", zhtc3*tmask(:,:,1) )      ! first 300m heat content 
     307     ENDIF 
     308     ! 
     309   END SUBROUTINE dia_hth 
     310 
     311 
     312   SUBROUTINE dia_hth_init 
     313      !!--------------------------------------------------------------------------- 
     314      !!                  ***  ROUTINE dia_hth_init  *** 
     315      !! 
     316      !! ** Purpose: Initialization for ML and thermocline depths 
     317      !! 
     318      !! ** Action : If any upper ocean diagnostic required by xml file, set in dia_hth 
     319      !!--------------------------------------------------------------------------- 
     320       ! 
     321      IF(lwp) THEN 
     322         WRITE(numout,*) 
     323         WRITE(numout,*) 'dia_hth_init : heat and salt budgets diagnostics' 
     324         WRITE(numout,*) '~~~~~~~~~~~~ ' 
     325      ENDIF 
     326      ll_diahth = iom_use("mlddzt") .OR. iom_use("mldr0_3") .OR. iom_use("mldr0_1") .OR.  & 
     327           & iom_use("mld_dt02") .OR. iom_use("topthdep") .OR. iom_use("mldr10_3") .OR.  & 
     328           & iom_use("pycndep") .OR. iom_use("tinv") .OR. iom_use("depti").OR. & 
     329           & iom_use("20d") .OR. iom_use("28d") .OR. iom_use("hc300") 
     330      IF(lwp) THEN 
     331         WRITE(numout,*) '   output upper ocean diagnostics (T) or not (F)       ll_diahth = ', ll_diahth 
     332      ENDIF 
    53333      ! 
    54       ALLOCATE( hth(jpi,jpj), hd20(jpi,jpj), hd28(jpi,jpj), htc3(jpi,jpj), STAT=dia_hth_alloc ) 
    55       ! 
    56       CALL mpp_sum ( 'diahth', dia_hth_alloc ) 
    57       IF(dia_hth_alloc /= 0)   CALL ctl_stop( 'STOP', 'dia_hth_alloc: failed to allocate arrays.' ) 
    58       ! 
    59    END FUNCTION dia_hth_alloc 
    60  
    61  
    62    SUBROUTINE dia_hth( kt ) 
    63       !!--------------------------------------------------------------------- 
    64       !!                  ***  ROUTINE dia_hth  *** 
    65       !! 
    66       !! ** Purpose : Computes 
    67       !!      the mixing layer depth (turbocline): avt = 5.e-4 
    68       !!      the depth of strongest vertical temperature gradient 
    69       !!      the mixed layer depth with density     criteria: rho = rho(10m or surf) + 0.03(or 0.01) 
    70       !!      the mixed layer depth with temperature criteria: abs( tn - tn(10m) ) = 0.2        
    71       !!      the top of the thermochine: tn = tn(10m) - ztem2  
    72       !!      the pycnocline depth with density criteria equivalent to a temperature variation  
    73       !!                rho = rho10m + (dr/dT)(T,S,10m)*(-0.2 degC)  
    74       !!      the barrier layer thickness 
    75       !!      the maximal verical inversion of temperature and its depth max( 0, max of tn - tn(10m) ) 
    76       !!      the depth of the 20 degree isotherm (linear interpolation) 
    77       !!      the depth of the 28 degree isotherm (linear interpolation) 
    78       !!      the heat content of first 300 m 
    79       !! 
    80       !! ** Method :  
    81       !!------------------------------------------------------------------- 
    82       INTEGER, INTENT( in ) ::   kt      ! ocean time-step index 
    83       !! 
    84       INTEGER                          ::   ji, jj, jk            ! dummy loop arguments 
    85       INTEGER                          ::   iid, ilevel           ! temporary integers 
    86       INTEGER, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ik20, ik28  ! levels 
    87       REAL(wp)                         ::   zavt5 = 5.e-4_wp      ! Kz criterion for the turbocline depth 
    88       REAL(wp)                         ::   zrho3 = 0.03_wp       ! density     criterion for mixed layer depth 
    89       REAL(wp)                         ::   zrho1 = 0.01_wp       ! density     criterion for mixed layer depth 
    90       REAL(wp)                         ::   ztem2 = 0.2_wp        ! temperature criterion for mixed layer depth 
    91       REAL(wp)                         ::   zthick_0, zcoef       ! temporary scalars 
    92       REAL(wp)                         ::   zztmp, zzdep          ! temporary scalars inside do loop 
    93       REAL(wp)                         ::   zu, zv, zw, zut, zvt  ! temporary workspace 
    94       REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   zabs2      ! MLD: abs( tn - tn(10m) ) = ztem2  
    95       REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ztm2       ! Top of thermocline: tn = tn(10m) - ztem2      
    96       REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   zrho10_3   ! MLD: rho = rho10m + zrho3       
    97       REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   zpycn      ! pycnocline: rho = rho10m + (dr/dT)(T,S,10m)*(-0.2 degC) 
    98       REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ztinv      ! max of temperature inversion 
    99       REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   zdepinv    ! depth of temperature inversion 
    100       REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   zrho0_3    ! MLD rho = rho(surf) = 0.03 
    101       REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   zrho0_1    ! MLD rho = rho(surf) = 0.01 
    102       REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   zmaxdzT    ! max of dT/dz 
    103       REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   zthick     ! vertical integration thickness  
    104       REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   zdelr      ! delta rho equivalent to deltaT = 0.2 
    105       !!---------------------------------------------------------------------- 
    106       IF( ln_timing )   CALL timing_start('dia_hth') 
    107  
    108       IF( kt == nit000 ) THEN 
    109          !                                      ! allocate dia_hth array 
    110          IF( dia_hth_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'dia_hth : unable to allocate standard arrays' ) 
    111  
    112          IF(.NOT. ALLOCATED(ik20) ) THEN 
    113             ALLOCATE(ik20(jpi,jpj), ik28(jpi,jpj), & 
    114                &      zabs2(jpi,jpj),   & 
    115                &      ztm2(jpi,jpj),    & 
    116                &      zrho10_3(jpi,jpj),& 
    117                &      zpycn(jpi,jpj),   & 
    118                &      ztinv(jpi,jpj),   & 
    119                &      zdepinv(jpi,jpj), & 
    120                &      zrho0_3(jpi,jpj), & 
    121                &      zrho0_1(jpi,jpj), & 
    122                &      zmaxdzT(jpi,jpj), & 
    123                &      zthick(jpi,jpj),  & 
    124                &      zdelr(jpi,jpj), STAT=ji) 
    125             CALL mpp_sum('diahth', ji) 
    126             IF( ji /= 0 )   CALL ctl_stop( 'STOP', 'dia_hth : unable to allocate standard ocean arrays' ) 
    127          END IF 
    128  
    129          IF(lwp) WRITE(numout,*) 
    130          IF(lwp) WRITE(numout,*) 'dia_hth : diagnostics of the thermocline depth' 
    131          IF(lwp) WRITE(numout,*) '~~~~~~~ ' 
    132          IF(lwp) WRITE(numout,*) 
    133       ENDIF 
    134  
    135       ! initialization 
    136       ztinv  (:,:) = 0._wp   
    137       zdepinv(:,:) = 0._wp   
    138       zmaxdzT(:,:) = 0._wp   
    139       DO jj = 1, jpj 
    140          DO ji = 1, jpi 
    141             zztmp = gdepw_n(ji,jj,mbkt(ji,jj)+1)  
    142             hth     (ji,jj) = zztmp 
    143             zabs2   (ji,jj) = zztmp 
    144             ztm2    (ji,jj) = zztmp 
    145             zrho10_3(ji,jj) = zztmp 
    146             zpycn   (ji,jj) = zztmp 
    147         END DO 
    148       END DO 
    149       IF( nla10 > 1 ) THEN  
    150          DO jj = 1, jpj 
    151             DO ji = 1, jpi 
    152                zztmp = gdepw_n(ji,jj,mbkt(ji,jj)+1)  
    153                zrho0_3(ji,jj) = zztmp 
    154                zrho0_1(ji,jj) = zztmp 
    155             END DO 
    156          END DO 
    157       ENDIF 
    158        
    159       ! Preliminary computation 
    160       ! computation of zdelr = (dr/dT)(T,S,10m)*(-0.2 degC) 
    161       DO jj = 1, jpj 
    162          DO ji = 1, jpi 
    163             IF( tmask(ji,jj,nla10) == 1. ) THEN 
    164                zu  =  1779.50 + 11.250 * tsn(ji,jj,nla10,jp_tem) - 3.80   * tsn(ji,jj,nla10,jp_sal)                             & 
    165                   &                                              - 0.0745 * tsn(ji,jj,nla10,jp_tem) * tsn(ji,jj,nla10,jp_tem)   & 
    166                   &                                              - 0.0100 * tsn(ji,jj,nla10,jp_tem) * tsn(ji,jj,nla10,jp_sal) 
    167                zv  =  5891.00 + 38.000 * tsn(ji,jj,nla10,jp_tem) + 3.00   * tsn(ji,jj,nla10,jp_sal)                             & 
    168                   &                                              - 0.3750 * tsn(ji,jj,nla10,jp_tem) * tsn(ji,jj,nla10,jp_tem) 
    169                zut =    11.25 -  0.149 * tsn(ji,jj,nla10,jp_tem) - 0.01   * tsn(ji,jj,nla10,jp_sal) 
    170                zvt =    38.00 -  0.750 * tsn(ji,jj,nla10,jp_tem) 
    171                zw  = (zu + 0.698*zv) * (zu + 0.698*zv) 
    172                zdelr(ji,jj) = ztem2 * (1000.*(zut*zv - zvt*zu)/zw) 
    173             ELSE 
    174                zdelr(ji,jj) = 0._wp 
    175             ENDIF 
    176          END DO 
    177       END DO 
    178  
    179       ! ------------------------------------------------------------- ! 
    180       ! thermocline depth: strongest vertical gradient of temperature ! 
    181       ! turbocline depth (mixing layer depth): avt = zavt5            ! 
    182       ! MLD: rho = rho(1) + zrho3                                     ! 
    183       ! MLD: rho = rho(1) + zrho1                                     ! 
    184       ! ------------------------------------------------------------- ! 
    185       DO jk = jpkm1, 2, -1   ! loop from bottom to 2 
    186          DO jj = 1, jpj 
    187             DO ji = 1, jpi 
    188                ! 
    189                zzdep = gdepw_n(ji,jj,jk) 
    190                zztmp = ( tsn(ji,jj,jk-1,jp_tem) - tsn(ji,jj,jk,jp_tem) ) / zzdep * tmask(ji,jj,jk)   ! vertical gradient of temperature (dT/dz) 
    191                zzdep = zzdep * tmask(ji,jj,1) 
    192  
    193                IF( zztmp > zmaxdzT(ji,jj) ) THEN                         
    194                   zmaxdzT(ji,jj) = zztmp   ;   hth    (ji,jj) = zzdep                ! max and depth of dT/dz 
    195                ENDIF 
    196                 
    197                IF( nla10 > 1 ) THEN  
    198                   zztmp = rhop(ji,jj,jk) - rhop(ji,jj,1)                             ! delta rho(1) 
    199                   IF( zztmp > zrho3 )          zrho0_3(ji,jj) = zzdep                ! > 0.03 
    200                   IF( zztmp > zrho1 )          zrho0_1(ji,jj) = zzdep                ! > 0.01 
    201                ENDIF 
    202  
    203             END DO 
    204          END DO 
    205       END DO 
    206        
    207       CALL iom_put( "mlddzt", hth )            ! depth of the thermocline 
    208       IF( nla10 > 1 ) THEN  
    209          CALL iom_put( "mldr0_3", zrho0_3 )   ! MLD delta rho(surf) = 0.03 
    210          CALL iom_put( "mldr0_1", zrho0_1 )   ! MLD delta rho(surf) = 0.01 
    211       ENDIF 
    212  
    213       ! ------------------------------------------------------------- ! 
    214       ! MLD: abs( tn - tn(10m) ) = ztem2                              ! 
    215       ! Top of thermocline: tn = tn(10m) - ztem2                      ! 
    216       ! MLD: rho = rho10m + zrho3                                     ! 
    217       ! pycnocline: rho = rho10m + (dr/dT)(T,S,10m)*(-0.2 degC)       ! 
    218       ! temperature inversion: max( 0, max of tn - tn(10m) )          ! 
    219       ! depth of temperature inversion                                ! 
    220       ! ------------------------------------------------------------- ! 
    221       DO jk = jpkm1, nlb10, -1   ! loop from bottom to nlb10 
    222          DO jj = 1, jpj 
    223             DO ji = 1, jpi 
    224                ! 
    225                zzdep = gdepw_n(ji,jj,jk) * tmask(ji,jj,1) 
    226                ! 
    227                zztmp = tsn(ji,jj,nla10,jp_tem) - tsn(ji,jj,jk,jp_tem)  ! - delta T(10m) 
    228                IF( ABS(zztmp) > ztem2 )      zabs2   (ji,jj) = zzdep   ! abs > 0.2 
    229                IF(     zztmp  > ztem2 )      ztm2    (ji,jj) = zzdep   ! > 0.2 
    230                zztmp = -zztmp                                          ! delta T(10m) 
    231                IF( zztmp >  ztinv(ji,jj) ) THEN                        ! temperature inversion 
    232                   ztinv(ji,jj) = zztmp   ;   zdepinv (ji,jj) = zzdep   ! max value and depth 
    233                ENDIF 
    234  
    235                zztmp = rhop(ji,jj,jk) - rhop(ji,jj,nla10)              ! delta rho(10m) 
    236                IF( zztmp > zrho3        )    zrho10_3(ji,jj) = zzdep   ! > 0.03 
    237                IF( zztmp > zdelr(ji,jj) )    zpycn   (ji,jj) = zzdep   ! > equi. delta T(10m) - 0.2 
    238                ! 
    239             END DO 
    240          END DO 
    241       END DO 
    242  
    243       CALL iom_put( "mld_dt02", zabs2        )   ! MLD abs(delta t) - 0.2 
    244       CALL iom_put( "topthdep", ztm2         )   ! T(10) - 0.2 
    245       CALL iom_put( "mldr10_3", zrho10_3     )   ! MLD delta rho(10m) = 0.03 
    246       CALL iom_put( "pycndep" , zpycn        )   ! MLD delta rho equi. delta T(10m) = 0.2 
    247       CALL iom_put( "tinv"    , ztinv        )   ! max. temp. inv. (t10 ref)  
    248       CALL iom_put( "depti"   , zdepinv      )   ! depth of max. temp. inv. (t10 ref)  
    249  
    250  
    251       ! ----------------------------------- ! 
    252       ! search deepest level above 20C/28C  ! 
    253       ! ----------------------------------- ! 
    254       ik20(:,:) = 1 
    255       ik28(:,:) = 1 
    256       DO jk = 1, jpkm1   ! beware temperature is not always decreasing with depth => loop from top to bottom 
    257          DO jj = 1, jpj 
    258             DO ji = 1, jpi 
    259                zztmp = tsn(ji,jj,jk,jp_tem) 
    260                IF( zztmp >= 20. )   ik20(ji,jj) = jk 
    261                IF( zztmp >= 28. )   ik28(ji,jj) = jk 
    262             END DO 
    263          END DO 
    264       END DO 
    265  
    266       ! --------------------------- ! 
    267       !  Depth of 20C/28C isotherm  ! 
    268       ! --------------------------- ! 
    269       DO jj = 1, jpj 
    270          DO ji = 1, jpi 
    271             ! 
    272             zzdep = gdepw_n(ji,jj,mbkt(ji,jj)+1)       ! depth of the oean bottom 
    273             ! 
    274             iid = ik20(ji,jj) 
    275             IF( iid /= 1 ) THEN  
    276                zztmp =      gdept_n(ji,jj,iid  )   &                     ! linear interpolation 
    277                   &  + (    gdept_n(ji,jj,iid+1) - gdept_n(ji,jj,iid)                       )   & 
    278                   &  * ( 20.*tmask(ji,jj,iid+1) - tsn(ji,jj,iid,jp_tem)                       )   & 
    279                   &  / ( tsn(ji,jj,iid+1,jp_tem) - tsn(ji,jj,iid,jp_tem) + (1.-tmask(ji,jj,1)) ) 
    280                hd20(ji,jj) = MIN( zztmp , zzdep) * tmask(ji,jj,1)       ! bound by the ocean depth 
    281             ELSE  
    282                hd20(ji,jj) = 0._wp 
    283             ENDIF 
    284             ! 
    285             iid = ik28(ji,jj) 
    286             IF( iid /= 1 ) THEN  
    287                zztmp =      gdept_n(ji,jj,iid  )   &                     ! linear interpolation 
    288                   &  + (    gdept_n(ji,jj,iid+1) - gdept_n(ji,jj,iid)                       )   & 
    289                   &  * ( 28.*tmask(ji,jj,iid+1) -    tsn(ji,jj,iid,jp_tem)                       )   & 
    290                   &  / (  tsn(ji,jj,iid+1,jp_tem) -    tsn(ji,jj,iid,jp_tem) + (1.-tmask(ji,jj,1)) ) 
    291                hd28(ji,jj) = MIN( zztmp , zzdep ) * tmask(ji,jj,1)      ! bound by the ocean depth 
    292             ELSE  
    293                hd28(ji,jj) = 0._wp 
    294             ENDIF 
    295  
    296          END DO 
    297       END DO 
    298       CALL iom_put( "20d", hd20 )   ! depth of the 20 isotherm 
    299       CALL iom_put( "28d", hd28 )   ! depth of the 28 isotherm 
    300  
    301       ! ----------------------------- ! 
    302       !  Heat content of first 300 m  ! 
    303       ! ----------------------------- ! 
    304  
    305       ! find ilevel with (ilevel+1) the deepest W-level above 300m (we assume we can use e3t_1d to do this search...) 
    306       ilevel   = 0 
    307       zthick_0 = 0._wp 
    308       DO jk = 1, jpkm1                       
    309          zthick_0 = zthick_0 + e3t_1d(jk) 
    310          IF( zthick_0 < 300. )   ilevel = jk 
    311       END DO 
    312       ! surface boundary condition 
    313       IF( ln_linssh ) THEN   ;   zthick(:,:) = sshn(:,:)   ;   htc3(:,:) = tsn(:,:,1,jp_tem) * sshn(:,:) * tmask(:,:,1)   
    314       ELSE                   ;   zthick(:,:) = 0._wp       ;   htc3(:,:) = 0._wp                                    
    315       ENDIF 
    316       ! integration down to ilevel 
    317       DO jk = 1, ilevel 
    318          zthick(:,:) = zthick(:,:) + e3t_n(:,:,jk) 
    319          htc3  (:,:) = htc3  (:,:) + e3t_n(:,:,jk) * tsn(:,:,jk,jp_tem) * tmask(:,:,jk) 
    320       END DO 
    321       ! deepest layer 
    322       zthick(:,:) = 300. - zthick(:,:)   !   remaining thickness to reach 300m 
    323       DO jj = 1, jpj 
    324          DO ji = 1, jpi 
    325             htc3(ji,jj) = htc3(ji,jj) + tsn(ji,jj,ilevel+1,jp_tem)                  & 
    326                &                      * MIN( e3t_n(ji,jj,ilevel+1), zthick(ji,jj) ) * tmask(ji,jj,ilevel+1) 
    327          END DO 
    328       END DO 
    329       ! from temperature to heat contain 
    330       zcoef = rau0 * rcp 
    331       htc3(:,:) = zcoef * htc3(:,:) 
    332       CALL iom_put( "hc300", htc3 )      ! first 300m heat content 
    333       ! 
    334       IF( ln_timing )   CALL timing_stop('dia_hth') 
    335       ! 
    336    END SUBROUTINE dia_hth 
    337  
    338 #else 
    339    !!---------------------------------------------------------------------- 
    340    !!   Default option :                                       Empty module 
    341    !!---------------------------------------------------------------------- 
    342    LOGICAL , PUBLIC, PARAMETER ::   lk_diahth = .FALSE.  !: thermocline-20d depths flag 
    343 CONTAINS 
    344    SUBROUTINE dia_hth( kt )         ! Empty routine 
    345       IMPLICIT NONE 
    346       INTEGER, INTENT( in ) :: kt 
    347       WRITE(*,*) 'dia_hth: You should not have seen this print! error?', kt 
    348    END SUBROUTINE dia_hth 
    349 #endif 
    350  
    351    !!====================================================================== 
     334   END SUBROUTINE dia_hth_init 
    352335END MODULE diahth 
  • NEMO/branches/UKMO/NEMO_4.0_FKOSM/src/OCE/step.F90

    r10888 r12516  
    101101      IF( kstp == nit000 ) THEN                       ! initialize IOM context (must be done after nemo_init for AGRIF+XIOS+OASIS) 
    102102                             CALL iom_init(      cxios_context          )  ! for model grid (including passible AGRIF zoom) 
    103          IF( ln_crs      )   CALL iom_init( TRIM(cxios_context)//"_crs" )  ! for coarse grid 
     103                             IF( ln_crs      )   CALL iom_init( TRIM(cxios_context)//"_crs" )  ! for coarse grid 
     104                             CALL dia_hth_init    ! extra ML depth diagnostics, thermocline depths 
    104105      ENDIF 
    105106      IF( kstp /= nit000 )   CALL day( kstp )         ! Calendar (day was already called at nit000 in day_init) 
     
    214215      IF( lk_floats  )   CALL flo_stp ( kstp )        ! drifting Floats 
    215216      IF( ln_diacfl  )   CALL dia_cfl ( kstp )        ! Courant number diagnostics 
    216       IF( lk_diahth  )   CALL dia_hth ( kstp )        ! Thermocline depth (20 degres isotherm depth) 
     217      IF( ll_diahth  )   CALL dia_hth ( kstp )        ! Thermocline depth (20 degres isotherm depth) 
    217218      IF( lk_diadct  )   CALL dia_dct ( kstp )        ! Transports 
    218219                         CALL dia_ar5 ( kstp )        ! ar5 diag 
Note: See TracChangeset for help on using the changeset viewer.