Ignore:
Timestamp:
2019-06-11T14:58:36+02:00 (18 months ago)
Author:
agn
Message:

new diahth to work without key_diahth: changes in diahth.F90

File:
1 edited

Legend:

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

    r11090 r11091  
    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 
    1212   !!---------------------------------------------------------------------- 
    13 #if   defined key_diahth 
    1413   !!---------------------------------------------------------------------- 
    1514   !!   'key_diahth' :                              thermocline depth diag. 
     
    2423   USE lib_mpp         ! MPP library 
    2524   USE iom             ! I/O library 
    26    USE timing          ! preformance summary 
    2725 
    2826   IMPLICIT NONE 
     
    3028 
    3129   PUBLIC   dia_hth       ! routine called by step.F90 
    32  
    33    LOGICAL , PUBLIC, PARAMETER ::   lk_diahth = .TRUE.    !: thermocline-20d depths flag 
    34     
     30   PUBLIC   dia_hth_init  ! routine called by nemogcm.F90 
     31 
     32   LOGICAL, PUBLIC ::   ln_diahth   !: Compute further diagnostics of ML and thermocline depth 
     33 
    3534   !!---------------------------------------------------------------------- 
    3635   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
     
    4948     !!      the depth of strongest vertical temperature gradient 
    5049     !!      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)  
     50     !!      the mixed layer depth with temperature criteria: abs( tn - tn(10m) ) = 0.2 
     51     !!      the top of the thermochine: tn = tn(10m) - ztem2 
     52     !!      the pycnocline depth with density criteria equivalent to a temperature variation 
     53     !!                rho = rho10m + (dr/dT)(T,S,10m)*(-0.2 degC) 
    5554     !!      the barrier layer thickness 
    5655     !!      the maximal verical inversion of temperature and its depth max( 0, max of tn - tn(10m) ) 
     
    5958     !!      the heat content of first 300 m 
    6059     !! 
    61      !! ** Method :  
     60     !! ** Method : 
    6261     !!------------------------------------------------------------------- 
    6362     INTEGER, INTENT( in ) ::   kt      ! ocean time-step index 
     
    6564     INTEGER                          ::   ji, jj, jk            ! dummy loop arguments 
    6665     INTEGER                          ::   iid, ilevel           ! temporary integers 
    67      INTEGER, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ik20, ik28  ! levels 
     66     INTEGER, DIMENSION(jpi,jpj) ::   ik20, ik28  ! levels 
    6867     REAL(wp)                         ::   zavt5 = 5.e-4_wp      ! Kz criterion for the turbocline depth 
    6968     REAL(wp)                         ::   zrho3 = 0.03_wp       ! density     criterion for mixed layer depth 
     
    7372     REAL(wp)                         ::   zztmp, zzdep          ! temporary scalars inside do loop 
    7473     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 
    139         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 
     74     REAL(wp), DIMENSION(jpi,jpj) ::   zabs2      ! MLD: abs( tn - tn(10m) ) = ztem2 
     75     REAL(wp), DIMENSION(jpi,jpj) ::   ztm2       ! Top of thermocline: tn = tn(10m) - ztem2 
     76     REAL(wp), DIMENSION(jpi,jpj) ::   zrho10_3   ! MLD: rho = rho10m + zrho3 
     77     REAL(wp), DIMENSION(jpi,jpj) ::   zpycn      ! pycnocline: rho = rho10m + (dr/dT)(T,S,10m)*(-0.2 degC) 
     78     REAL(wp), DIMENSION(jpi,jpj) ::   ztinv      ! max of temperature inversion 
     79     REAL(wp), DIMENSION(jpi,jpj) ::   zdepinv    ! depth of temperature inversion 
     80     REAL(wp), DIMENSION(jpi,jpj) ::   zrho0_3    ! MLD rho = rho(surf) = 0.03 
     81     REAL(wp), DIMENSION(jpi,jpj) ::   zrho0_1    ! MLD rho = rho(surf) = 0.01 
     82     REAL(wp), DIMENSION(jpi,jpj) ::   zmaxdzT    ! max of dT/dz 
     83     REAL(wp), DIMENSION(jpi,jpj) ::   zthick     ! vertical integration thickness 
     84     REAL(wp), DIMENSION(jpi,jpj) ::   zdelr      ! delta rho equivalent to deltaT = 0.2 
     85   ! note: following variables should move to local variables once iom_put is always used 
     86     REAL(wp), DIMENSION(jpi,jpj) ::   zhth    !: depth of the max vertical temperature gradient [m] 
     87     REAL(wp), DIMENSION(jpi,jpj) ::   zhd20   !: depth of 20 C isotherm                         [m] 
     88     REAL(wp), DIMENSION(jpi,jpj) ::   zhd28   !: depth of 28 C isotherm                         [m] 
     89     REAL(wp), DIMENSION(jpi,jpj) ::   zhtc3   !: heat content of first 300 m                    [W] 
    15090 
    15191     IF (iom_use("mlddzt").OR.iom_use("mldr0_3").OR.iom_use("mldr0_1")) THEN 
     
    15696        ! MLD: rho = rho(1) + zrho1                                     ! 
    15797        ! ------------------------------------------------------------- ! 
     98        zmaxdzT(:,:) = 0._wp 
     99        IF( nla10 > 1 ) THEN 
     100           DO jj = 1, jpj 
     101              DO ji = 1, jpi 
     102                 zztmp = gdepw_n(ji,jj,mbkt(ji,jj)+1) 
     103                 zrho0_3(ji,jj) = zztmp 
     104                 zrho0_1(ji,jj) = zztmp 
     105                 zhth(ji,jj) = zztmp 
     106              END DO 
     107           END DO 
     108        ELSE IF (iom_use("mlddzt")) THEN 
     109           DO jj = 1, jpj 
     110              DO ji = 1, jpi 
     111                 zztmp = gdepw_n(ji,jj,mbkt(ji,jj)+1) 
     112                 zhth(ji,jj) = zztmp 
     113              END DO 
     114           END DO 
     115        ELSE 
     116           zhth(:,:) = 0._wp 
     117        ENDIF 
     118 
    158119        DO jk = jpkm1, 2, -1   ! loop from bottom to 2 
    159120           DO jj = 1, jpj 
     
    164125                 zzdep = zzdep * tmask(ji,jj,1) 
    165126 
    166                  IF( zztmp > zmaxdzT(ji,jj) ) THEN                         
     127                 IF( zztmp > zmaxdzT(ji,jj) ) THEN 
    167128                    zmaxdzT(ji,jj) = zztmp   ;   zhth    (ji,jj) = zzdep                ! max and depth of dT/dz 
    168129                 ENDIF 
    169130 
    170                  IF( nla10 > 1 ) THEN  
     131                 IF( nla10 > 1 ) THEN 
    171132                    zztmp = rhop(ji,jj,jk) - rhop(ji,jj,1)                             ! delta rho(1) 
    172133                    IF( zztmp > zrho3 )          zrho0_3(ji,jj) = zzdep                ! > 0.03 
     
    179140 
    180141        IF (iom_use("mlddzt")) CALL iom_put( "mlddzt", zhth*tmask(:,:,1) )            ! depth of the thermocline 
    181         IF( nla10 > 1 ) THEN  
     142        IF( nla10 > 1 ) THEN 
    182143           IF (iom_use("mldr0_3")) CALL iom_put( "mldr0_3", zrho0_3*tmask(:,:,1) )   ! MLD delta rho(surf) = 0.03 
    183144           IF (iom_use("mldr0_1")) CALL iom_put( "mldr0_1", zrho0_1*tmask(:,:,1) )   ! MLD delta rho(surf) = 0.01 
     
    187148     IF (iom_use("mld_dt02").OR.iom_use("topthdep").OR.iom_use("mldr10_3").OR. & 
    188149          & iom_use("pycndep").OR.iom_use("tinv").OR.iom_use("depti")) THEN 
     150        DO jj = 1, jpj 
     151           DO ji = 1, jpi 
     152              zztmp = gdepw_n(ji,jj,mbkt(ji,jj)+1) 
     153              zabs2   (ji,jj) = zztmp 
     154              ztm2    (ji,jj) = zztmp 
     155              zrho10_3(ji,jj) = zztmp 
     156              zpycn   (ji,jj) = zztmp 
     157           END DO 
     158        END DO 
     159        ztinv  (:,:) = 0._wp 
     160        zdepinv(:,:) = 0._wp 
     161 
    189162        IF (iom_use("pycndep")) THEN 
    190163           ! Preliminary computation 
     
    207180              END DO 
    208181           END DO 
     182        ELSE 
     183           zdelr(:,:) = 0._wp 
    209184        ENDIF 
    210185 
     
    243218        IF (iom_use("mldr10_3")) CALL iom_put( "mldr10_3", zrho10_3*tmask(:,:,1)     )   ! MLD delta rho(10m) = 0.03 
    244219        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)  
     220        IF (iom_use("tinv")) CALL iom_put( "tinv"    , ztinv*tmask(:,:,1)        )   ! max. temp. inv. (t10 ref) 
     221        IF (iom_use("depti")) CALL iom_put( "depti"   , zdepinv*tmask(:,:,1)      )   ! depth of max. temp. inv. (t10 ref) 
    247222     ENDIF 
    248223 
     
    272247              ! 
    273248              iid = ik20(ji,jj) 
    274               IF( iid /= 1 ) THEN  
     249              IF( iid /= 1 ) THEN 
    275250                 zztmp =      gdept_n(ji,jj,iid  )   &                     ! linear interpolation 
    276251                      &  + (    gdept_n(ji,jj,iid+1) - gdept_n(ji,jj,iid)                       )   & 
     
    278253                      &  / ( tsn(ji,jj,iid+1,jp_tem) - tsn(ji,jj,iid,jp_tem) + (1.-tmask(ji,jj,1)) ) 
    279254                 zhd20(ji,jj) = MIN( zztmp , zzdep) * tmask(ji,jj,1)       ! bound by the ocean depth 
    280               ELSE  
     255              ELSE 
    281256                 zhd20(ji,jj) = 0._wp 
    282257              ENDIF 
    283258              ! 
    284259              iid = ik28(ji,jj) 
    285               IF( iid /= 1 ) THEN  
     260              IF( iid /= 1 ) THEN 
    286261                 zztmp =      gdept_n(ji,jj,iid  )   &                     ! linear interpolation 
    287262                      &  + (    gdept_n(ji,jj,iid+1) - gdept_n(ji,jj,iid)                       )   & 
     
    289264                      &  / (  tsn(ji,jj,iid+1,jp_tem) -    tsn(ji,jj,iid,jp_tem) + (1.-tmask(ji,jj,1)) ) 
    290265                 zhd28(ji,jj) = MIN( zztmp , zzdep ) * tmask(ji,jj,1)      ! bound by the ocean depth 
    291               ELSE  
     266              ELSE 
    292267                 zhd28(ji,jj) = 0._wp 
    293268              ENDIF 
     
    306281        ilevel   = 0 
    307282        zthick_0 = 0._wp 
    308         DO jk = 1, jpkm1                       
     283        DO jk = 1, jpkm1 
    309284           zthick_0 = zthick_0 + e3t_1d(jk) 
    310285           IF( zthick_0 < 300. )   ilevel = jk 
    311286        END DO 
    312287        ! 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                                    
     288        IF( ln_linssh ) THEN   ;   zthick(:,:) = sshn(:,:)   ;   zhtc3(:,:) = tsn(:,:,1,jp_tem) * sshn(:,:) * tmask(:,:,1) 
     289        ELSE                   ;   zthick(:,:) = 0._wp       ;   zhtc3(:,:) = 0._wp 
    315290        ENDIF 
    316291        ! integration down to ilevel 
     
    333308     ENDIF 
    334309     ! 
    335       IF( ln_timing )   CALL timing_stop('dia_hth') 
    336      ! 
    337310   END SUBROUTINE dia_hth 
    338311 
    339 #else 
    340    !!---------------------------------------------------------------------- 
    341    !!   Default option :                                       Empty module 
    342    !!---------------------------------------------------------------------- 
    343    LOGICAL , PUBLIC, PARAMETER ::   lk_diahth = .FALSE.  !: thermocline-20d depths flag 
    344 CONTAINS 
    345    SUBROUTINE dia_hth( kt )         ! Empty routine 
    346       IMPLICIT NONE 
    347       INTEGER, INTENT( in ) :: kt 
    348       WRITE(*,*) 'dia_hth: You should not have seen this print! error?', kt 
    349    END SUBROUTINE dia_hth 
    350 #endif 
    351  
    352    !!====================================================================== 
     312 
     313   SUBROUTINE dia_hth_init 
     314      !!--------------------------------------------------------------------------- 
     315      !!                  ***  ROUTINE dia_hth_init  *** 
     316      !! 
     317      !! ** Purpose: Initialization for ML and thermocline depths 
     318      !! 
     319      !! ** Action : Read namelist flag to permit or not extra Ml depths and thermocline depths 
     320      !!--------------------------------------------------------------------------- 
     321      INTEGER ::   ierror, ios   ! local integer 
     322      !! 
     323      NAMELIST/namhth/ ln_diahth 
     324      !!--------------------------------------------------------------------------- 
     325       ! 
     326      IF(lwp) THEN 
     327         WRITE(numout,*) 
     328         WRITE(numout,*) 'dia_hth_init : heat and salt budgets diagnostics' 
     329         WRITE(numout,*) '~~~~~~~~~~~~ ' 
     330      ENDIF 
     331      REWIND( numnam_ref )              ! Namelist namhth in reference namelist 
     332      READ  ( numnam_ref, namhth, IOSTAT = ios, ERR = 901) 
     333901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namhth in reference namelist', lwp ) 
     334      REWIND( numnam_cfg )              ! Namelist namhth in configuration namelist 
     335      READ  ( numnam_cfg, namhth, IOSTAT = ios, ERR = 902 ) 
     336902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namhth in configuration namelist', lwp ) 
     337      IF(lwm) WRITE( numond, namhth ) 
     338 
     339      IF(lwp) THEN 
     340         WRITE(numout,*) '   Namelist  namhth :' 
     341         WRITE(numout,*) '      check the heat and salt budgets (T) or not (F)       ln_diahth = ', ln_diahth 
     342      ENDIF 
     343      ! 
     344   END SUBROUTINE dia_hth_init 
    353345END MODULE diahth 
Note: See TracChangeset for help on using the changeset viewer.