Changeset 11100


Ignore:
Timestamp:
2019-06-11T16:07:19+02:00 (17 months ago)
Author:
agn
Message:

Properly undid premature check into release-4.0

File:
1 edited

Legend:

Unmodified
Added
Removed
  • NEMO/releases/release-4.0/src/OCE/DIA/diahth.F90

    r11098 r11100  
    3030 
    3131   PUBLIC   dia_hth       ! routine called by step.F90 
     32   PUBLIC   dia_hth_alloc ! routine called by nemogcm.F90 
    3233 
    3334   LOGICAL , PUBLIC, PARAMETER ::   lk_diahth = .TRUE.    !: thermocline-20d depths flag 
    3435    
     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] 
     41 
    3542   !!---------------------------------------------------------------------- 
    3643   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
     
    4047CONTAINS 
    4148 
     49   FUNCTION dia_hth_alloc() 
     50      !!--------------------------------------------------------------------- 
     51      INTEGER :: dia_hth_alloc 
     52      !!--------------------------------------------------------------------- 
     53      ! 
     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 
    4261 
    4362   SUBROUTINE dia_hth( kt ) 
    44      !!--------------------------------------------------------------------- 
    45      !!                  ***  ROUTINE dia_hth  *** 
    46      !! 
    47      !! ** Purpose : Computes 
    48      !!      the mixing layer depth (turbocline): avt = 5.e-4 
    49      !!      the depth of strongest vertical temperature gradient 
    50      !!      the mixed layer depth with density     criteria: rho = rho(10m or surf) + 0.03(or 0.01) 
    51      !!      the mixed layer depth with temperature criteria: abs( tn - tn(10m) ) = 0.2        
    52      !!      the top of the thermochine: tn = tn(10m) - ztem2  
    53      !!      the pycnocline depth with density criteria equivalent to a temperature variation  
    54      !!                rho = rho10m + (dr/dT)(T,S,10m)*(-0.2 degC)  
    55      !!      the barrier layer thickness 
    56      !!      the maximal verical inversion of temperature and its depth max( 0, max of tn - tn(10m) ) 
    57      !!      the depth of the 20 degree isotherm (linear interpolation) 
    58      !!      the depth of the 28 degree isotherm (linear interpolation) 
    59      !!      the heat content of first 300 m 
    60      !! 
    61      !! ** Method :  
    62      !!------------------------------------------------------------------- 
    63      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index 
    64      !! 
    65      INTEGER                          ::   ji, jj, jk            ! dummy loop arguments 
    66      INTEGER                          ::   iid, ilevel           ! temporary integers 
    67      INTEGER, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ik20, ik28  ! levels 
    68      REAL(wp)                         ::   zavt5 = 5.e-4_wp      ! Kz criterion for the turbocline depth 
    69      REAL(wp)                         ::   zrho3 = 0.03_wp       ! density     criterion for mixed layer depth 
    70      REAL(wp)                         ::   zrho1 = 0.01_wp       ! density     criterion for mixed layer depth 
    71      REAL(wp)                         ::   ztem2 = 0.2_wp        ! temperature criterion for mixed layer depth 
    72      REAL(wp)                         ::   zthick_0, zcoef       ! temporary scalars 
    73      REAL(wp)                         ::   zztmp, zzdep          ! temporary scalars inside do loop 
    74      REAL(wp)                         ::   zu, zv, zw, zut, zvt  ! temporary workspace 
    75      REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   zabs2      ! MLD: abs( tn - tn(10m) ) = ztem2  
    76      REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ztm2       ! Top of thermocline: tn = tn(10m) - ztem2      
    77      REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   zrho10_3   ! MLD: rho = rho10m + zrho3       
    78      REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   zpycn      ! pycnocline: rho = rho10m + (dr/dT)(T,S,10m)*(-0.2 degC) 
    79      REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ztinv      ! max of temperature inversion 
    80      REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   zdepinv    ! depth of temperature inversion 
    81      REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   zrho0_3    ! MLD rho = rho(surf) = 0.03 
    82      REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   zrho0_1    ! MLD rho = rho(surf) = 0.01 
    83      REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   zmaxdzT    ! max of dT/dz 
    84      REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   zthick     ! vertical integration thickness  
    85      REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   zdelr      ! delta rho equivalent to deltaT = 0.2 
    86    ! note: following variables should move to local variables once iom_put is always used  
    87      REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   zhth    !: depth of the max vertical temperature gradient [m] 
    88      REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   zhd20   !: depth of 20 C isotherm                         [m] 
    89      REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   zhd28   !: depth of 28 C isotherm                         [m] 
    90      REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   zhtc3   !: heat content of first 300 m                    [W] 
    91  
    92      !!---------------------------------------------------------------------- 
    93      IF( ln_timing )   CALL timing_start('dia_hth') 
    94  
    95      IF( kt == nit000 ) THEN 
    96         !                                      ! allocate dia_hth array 
    97         IF(.NOT. ALLOCATED(ik20) ) THEN 
    98            ALLOCATE(ik20(jpi,jpj), ik28(jpi,jpj), & 
    99                 &      zhth(jpi,jpj),   & 
    100                 &      zhd20(jpi,jpj),   & 
    101                 &      zhd28(jpi,jpj),   & 
    102                 &      zhtc3(jpi,jpj),   & 
    103                 &      zabs2(jpi,jpj),   & 
    104                 &      ztm2(jpi,jpj),    & 
    105                 &      zrho10_3(jpi,jpj),& 
    106                 &      zpycn(jpi,jpj),   & 
    107                 &      ztinv(jpi,jpj),   & 
    108                 &      zdepinv(jpi,jpj), & 
    109                 &      zrho0_3(jpi,jpj), & 
    110                 &      zrho0_1(jpi,jpj), & 
    111                 &      zmaxdzT(jpi,jpj), & 
    112                 &      zthick(jpi,jpj),  & 
    113                 &      zdelr(jpi,jpj), STAT=ji) 
    114            IF( lk_mpp  )   CALL mpp_sum(ji) 
    115            IF( ji /= 0 )   CALL ctl_stop( 'STOP', 'dia_hth : unable to allocate standard ocean arrays' ) 
    116         END IF 
    117  
    118         IF(lwp) WRITE(numout,*) 
    119         IF(lwp) WRITE(numout,*) 'dia_hth : diagnostics of the thermocline depth' 
    120         IF(lwp) WRITE(numout,*) '~~~~~~~ ' 
    121         IF(lwp) WRITE(numout,*) 
    122      ENDIF 
    123  
    124      IF(iom_use("mlddzt").OR.iom_use("mldr0_3").OR.iom_use("mldr0_1").OR.iom_use("mld_dt02") & 
    125           & .OR.iom_use("topthdep").OR.iom_use("mldr10_3").OR.iom_use("pycndep").OR.iom_use("tinv").OR.iom_use("depti")) THEN 
    126         ! initialization 
    127         ztinv  (:,:) = 0._wp   
    128         zdepinv(:,:) = 0._wp   
    129         zmaxdzT(:,:) = 0._wp   
    130         DO jj = 1, jpj 
    131            DO ji = 1, jpi 
    132               zztmp = gdepw_n(ji,jj,mbkt(ji,jj)+1)  
    133               zhth     (ji,jj) = zztmp 
    134               zabs2   (ji,jj) = zztmp 
    135               ztm2    (ji,jj) = zztmp 
    136               zrho10_3(ji,jj) = zztmp 
    137               zpycn   (ji,jj) = zztmp 
    138            END DO 
     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 
    139147        END DO 
    140         IF( nla10 > 1 ) THEN  
    141            DO jj = 1, jpj 
    142               DO ji = 1, jpi 
    143                  zztmp = gdepw_n(ji,jj,mbkt(ji,jj)+1)  
    144                  zrho0_3(ji,jj) = zztmp 
    145                  zrho0_1(ji,jj) = zztmp 
    146               END DO 
    147            END DO 
    148         ENDIF 
    149      ENDIF 
    150  
    151      IF (iom_use("mlddzt").OR.iom_use("mldr0_3").OR.iom_use("mldr0_1")) THEN 
    152         ! ------------------------------------------------------------- ! 
    153         ! thermocline depth: strongest vertical gradient of temperature ! 
    154         ! turbocline depth (mixing layer depth): avt = zavt5            ! 
    155         ! MLD: rho = rho(1) + zrho3                                     ! 
    156         ! MLD: rho = rho(1) + zrho1                                     ! 
    157         ! ------------------------------------------------------------- ! 
    158         DO jk = jpkm1, 2, -1   ! loop from bottom to 2 
    159            DO jj = 1, jpj 
    160               DO ji = 1, jpi 
    161                  ! 
    162                  zzdep = gdepw_n(ji,jj,jk) 
    163                  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) 
    164                  zzdep = zzdep * tmask(ji,jj,1) 
    165  
    166                  IF( zztmp > zmaxdzT(ji,jj) ) THEN                         
    167                     zmaxdzT(ji,jj) = zztmp   ;   zhth    (ji,jj) = zzdep                ! max and depth of dT/dz 
    168                  ENDIF 
    169  
    170                  IF( nla10 > 1 ) THEN  
    171                     zztmp = rhop(ji,jj,jk) - rhop(ji,jj,1)                             ! delta rho(1) 
    172                     IF( zztmp > zrho3 )          zrho0_3(ji,jj) = zzdep                ! > 0.03 
    173                     IF( zztmp > zrho1 )          zrho0_1(ji,jj) = zzdep                ! > 0.01 
    174                  ENDIF 
    175  
    176               END DO 
    177            END DO 
    178         END DO 
    179  
    180         IF (iom_use("mlddzt")) CALL iom_put( "mlddzt", zhth*tmask(:,:,1) )            ! depth of the thermocline 
    181         IF( nla10 > 1 ) THEN  
    182            IF (iom_use("mldr0_3")) CALL iom_put( "mldr0_3", zrho0_3*tmask(:,:,1) )   ! MLD delta rho(surf) = 0.03 
    183            IF (iom_use("mldr0_1")) CALL iom_put( "mldr0_1", zrho0_1*tmask(:,:,1) )   ! MLD delta rho(surf) = 0.01 
    184         ENDIF 
    185      ENDIF 
    186  
    187      IF (iom_use("mld_dt02").OR.iom_use("topthdep").OR.iom_use("mldr10_3").OR. & 
    188           & iom_use("pycndep").OR.iom_use("tinv").OR.iom_use("depti")) THEN 
    189         IF (iom_use("pycndep")) THEN 
    190            ! Preliminary computation 
    191            ! computation of zdelr = (dr/dT)(T,S,10m)*(-0.2 degC) 
    192            DO jj = 1, jpj 
    193               DO ji = 1, jpi 
    194                  IF( tmask(ji,jj,nla10) == 1. ) THEN 
    195                     zu  =  1779.50 + 11.250 * tsn(ji,jj,nla10,jp_tem) - 3.80   * tsn(ji,jj,nla10,jp_sal)                             & 
    196                          &                                              - 0.0745 * tsn(ji,jj,nla10,jp_tem) * tsn(ji,jj,nla10,jp_tem)   & 
    197                          &                                              - 0.0100 * tsn(ji,jj,nla10,jp_tem) * tsn(ji,jj,nla10,jp_sal) 
    198                     zv  =  5891.00 + 38.000 * tsn(ji,jj,nla10,jp_tem) + 3.00   * tsn(ji,jj,nla10,jp_sal)                             & 
    199                          &                                              - 0.3750 * tsn(ji,jj,nla10,jp_tem) * tsn(ji,jj,nla10,jp_tem) 
    200                     zut =    11.25 -  0.149 * tsn(ji,jj,nla10,jp_tem) - 0.01   * tsn(ji,jj,nla10,jp_sal) 
    201                     zvt =    38.00 -  0.750 * tsn(ji,jj,nla10,jp_tem) 
    202                     zw  = (zu + 0.698*zv) * (zu + 0.698*zv) 
    203                     zdelr(ji,jj) = ztem2 * (1000.*(zut*zv - zvt*zu)/zw) 
    204                  ELSE 
    205                     zdelr(ji,jj) = 0._wp 
    206                  ENDIF 
    207               END DO 
    208            END DO 
    209         ENDIF 
    210  
    211         ! ------------------------------------------------------------- ! 
    212         ! MLD: abs( tn - tn(10m) ) = ztem2                              ! 
    213         ! Top of thermocline: tn = tn(10m) - ztem2                      ! 
    214         ! MLD: rho = rho10m + zrho3                                     ! 
    215         ! pycnocline: rho = rho10m + (dr/dT)(T,S,10m)*(-0.2 degC)       ! 
    216         ! temperature inversion: max( 0, max of tn - tn(10m) )          ! 
    217         ! depth of temperature inversion                                ! 
    218         ! ------------------------------------------------------------- ! 
    219         DO jk = jpkm1, nlb10, -1   ! loop from bottom to nlb10 
    220            DO jj = 1, jpj 
    221               DO ji = 1, jpi 
    222                  ! 
    223                  zzdep = gdepw_n(ji,jj,jk) * tmask(ji,jj,1) 
    224                  ! 
    225                  zztmp = tsn(ji,jj,nla10,jp_tem) - tsn(ji,jj,jk,jp_tem)  ! - delta T(10m) 
    226                  IF( ABS(zztmp) > ztem2 )      zabs2   (ji,jj) = zzdep   ! abs > 0.2 
    227                  IF(     zztmp  > ztem2 )      ztm2    (ji,jj) = zzdep   ! > 0.2 
    228                  zztmp = -zztmp                                          ! delta T(10m) 
    229                  IF( zztmp >  ztinv(ji,jj) ) THEN                        ! temperature inversion 
    230                     ztinv(ji,jj) = zztmp   ;   zdepinv (ji,jj) = zzdep   ! max value and depth 
    231                  ENDIF 
    232  
    233                  zztmp = rhop(ji,jj,jk) - rhop(ji,jj,nla10)              ! delta rho(10m) 
    234                  IF( zztmp > zrho3        )    zrho10_3(ji,jj) = zzdep   ! > 0.03 
    235                  IF( zztmp > zdelr(ji,jj) )    zpycn   (ji,jj) = zzdep   ! > equi. delta T(10m) - 0.2 
    236                  ! 
    237               END DO 
    238            END DO 
    239         END DO 
    240  
    241         IF (iom_use("mld_dt02")) CALL iom_put( "mld_dt02", zabs2*tmask(:,:,1)        )   ! MLD abs(delta t) - 0.2 
    242         IF (iom_use("topthdep")) CALL iom_put( "topthdep", ztm2*tmask(:,:,1)         )   ! T(10) - 0.2 
    243         IF (iom_use("mldr10_3")) CALL iom_put( "mldr10_3", zrho10_3*tmask(:,:,1)     )   ! MLD delta rho(10m) = 0.03 
    244         IF (iom_use("pycndep")) CALL iom_put( "pycndep" , zpycn*tmask(:,:,1)        )   ! MLD delta rho equi. delta T(10m) = 0.2 
    245         IF (iom_use("tinv")) CALL iom_put( "tinv"    , ztinv*tmask(:,:,1)        )   ! max. temp. inv. (t10 ref)  
    246         IF (iom_use("depti")) CALL iom_put( "depti"   , zdepinv*tmask(:,:,1)      )   ! depth of max. temp. inv. (t10 ref)  
    247      ENDIF 
    248  
    249      IF(iom_use("20d").OR.iom_use("28d")) THEN 
    250         ! ----------------------------------- ! 
    251         ! search deepest level above 20C/28C  ! 
    252         ! ----------------------------------- ! 
    253         ik20(:,:) = 1 
    254         ik28(:,:) = 1 
    255         DO jk = 1, jpkm1   ! beware temperature is not always decreasing with depth => loop from top to bottom 
    256            DO jj = 1, jpj 
    257               DO ji = 1, jpi 
    258                  zztmp = tsn(ji,jj,jk,jp_tem) 
    259                  IF( zztmp >= 20. )   ik20(ji,jj) = jk 
    260                  IF( zztmp >= 28. )   ik28(ji,jj) = jk 
    261               END DO 
    262            END DO 
    263         END DO 
    264  
    265         ! --------------------------- ! 
    266         !  Depth of 20C/28C isotherm  ! 
    267         ! --------------------------- ! 
    268         DO jj = 1, jpj 
    269            DO ji = 1, jpi 
    270               ! 
    271               zzdep = gdepw_n(ji,jj,mbkt(ji,jj)+1)       ! depth of the oean bottom 
    272               ! 
    273               iid = ik20(ji,jj) 
    274               IF( iid /= 1 ) THEN  
    275                  zztmp =      gdept_n(ji,jj,iid  )   &                     ! linear interpolation 
    276                       &  + (    gdept_n(ji,jj,iid+1) - gdept_n(ji,jj,iid)                       )   & 
    277                       &  * ( 20.*tmask(ji,jj,iid+1) - tsn(ji,jj,iid,jp_tem)                       )   & 
    278                       &  / ( tsn(ji,jj,iid+1,jp_tem) - tsn(ji,jj,iid,jp_tem) + (1.-tmask(ji,jj,1)) ) 
    279                  zhd20(ji,jj) = MIN( zztmp , zzdep) * tmask(ji,jj,1)       ! bound by the ocean depth 
    280               ELSE  
    281                  zhd20(ji,jj) = 0._wp 
    282               ENDIF 
    283               ! 
    284               iid = ik28(ji,jj) 
    285               IF( iid /= 1 ) THEN  
    286                  zztmp =      gdept_n(ji,jj,iid  )   &                     ! linear interpolation 
    287                       &  + (    gdept_n(ji,jj,iid+1) - gdept_n(ji,jj,iid)                       )   & 
    288                       &  * ( 28.*tmask(ji,jj,iid+1) -    tsn(ji,jj,iid,jp_tem)                       )   & 
    289                       &  / (  tsn(ji,jj,iid+1,jp_tem) -    tsn(ji,jj,iid,jp_tem) + (1.-tmask(ji,jj,1)) ) 
    290                  zhd28(ji,jj) = MIN( zztmp , zzdep ) * tmask(ji,jj,1)      ! bound by the ocean depth 
    291               ELSE  
    292                  zhd28(ji,jj) = 0._wp 
    293               ENDIF 
    294  
    295            END DO 
    296         END DO 
    297         CALL iom_put( "20d", zhd20 )   ! depth of the 20 isotherm 
    298         CALL iom_put( "28d", zhd28 )   ! depth of the 28 isotherm 
    299      ENDIF 
    300  
    301      ! ----------------------------- ! 
    302      !  Heat content of first 300 m  ! 
    303      ! ----------------------------- ! 
    304      IF (iom_use("hc300")) THEN 
    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(:,:)   ;   zhtc3(:,:) = tsn(:,:,1,jp_tem) * sshn(:,:) * tmask(:,:,1)   
    314         ELSE                   ;   zthick(:,:) = 0._wp       ;   zhtc3(:,:) = 0._wp                                    
    315         ENDIF 
    316         ! integration down to ilevel 
    317         DO jk = 1, ilevel 
    318            zthick(:,:) = zthick(:,:) + e3t_n(:,:,jk) 
    319            zhtc3  (:,:) = zhtc3  (:,:) + 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               zhtc3(ji,jj) = zhtc3(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         zhtc3(:,:) = zcoef * zhtc3(:,:) 
    332         CALL iom_put( "hc300", zhtc3*tmask(:,:,1) )      ! first 300m heat content 
    333      ENDIF 
    334      ! 
     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      ! 
    335334      IF( ln_timing )   CALL timing_stop('dia_hth') 
    336      ! 
     335      ! 
    337336   END SUBROUTINE dia_hth 
    338337 
Note: See TracChangeset for help on using the changeset viewer.