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

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

Ignore:
Timestamp:
2009-08-04T16:56:59+02:00 (15 years ago)
Author:
smasson
Message:

update diahth and zdfmxl, see ticket:468

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/NEMO/OPA_SRC/DIA/diahth.F90

    r1551 r1577  
    99   !!                 !  1999-07  (E. Guilyardi)  hd28 + heat content  
    1010   !!            8.5  !  2002-06  (G. Madec)  F90: Free form and module 
    11    !!   NEMO     3.2  !  2009-07  (S. Masson) hc300 bugfix + cleaning 
     11   !!   NEMO     3.2  !  2009-07  (S. Masson) hc300 bugfix + cleaning + add new diag 
    1212   !!---------------------------------------------------------------------- 
    1313 
     
    1616   !!   'key_diahth' :                              thermocline depth diag. 
    1717   !!---------------------------------------------------------------------- 
    18    !!   dia_hth      : Compute diagnostics associated with the thermocline 
     18   !!   dia_hth      : Compute varius diagnostics associated with the mixed layer 
    1919   !!---------------------------------------------------------------------- 
    2020   !! * Modules used 
    2121   USE oce             ! ocean dynamics and tracers 
    2222   USE dom_oce         ! ocean space and time domain 
     23   USE zdf_oce         ! ocean vertical physics 
    2324   USE phycst          ! physical constants 
    2425   USE in_out_manager  ! I/O manager 
     
    3233 
    3334   !! * Shared module variables 
    34    LOGICAL , PUBLIC, PARAMETER ::   lk_diahth = .TRUE.   !: thermocline-20d depths flag 
    35    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   &  !: 
    36       hth  ,      &  !: depth of the max vertical temperature gradient (m) 
    37       hd20 ,      &  !: depth of 20 C isotherm (m) 
    38       hd28 ,      &  !: depth of 28 C isotherm (m) 
    39       htc3           !: heat content of first 300 m 
     35   LOGICAL , PUBLIC, PARAMETER          ::   lk_diahth = .TRUE.   !: thermocline-20d depths flag 
     36   ! note: following variables should move to local variables once iom_put is always used  
     37   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   hmld                 !: mixing layer depth (turbocline)                [m] 
     38   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   hth                  !: depth of the max vertical temperature gradient [m] 
     39   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   hd20                 !: depth of 20 C isotherm                         [m] 
     40   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   hd28                 !: depth of 28 C isotherm                         [m] 
     41   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   htc3                 !: heat content of first 300 m                    [W] 
    4042 
    4143   !! * Substitutions 
     
    5355      !!                  ***  ROUTINE dia_hth  *** 
    5456      !! 
    55       !! ** Purpose : 
    56       !!      Computes the depth of strongest vertical temperature gradient 
    57       !!      Computes the depth of the 20 degree isotherm 
    58       !!      Computes the depth of the 28 degree isotherm 
    59       !!      Computes the heat content of first 300 m 
     57      !! ** Purpose : Computes 
     58      !!      the mixing layer depth (turbocline): avt = 5.e-4 
     59      !!      the depth of strongest vertical temperature gradient 
     60      !!      the mixed layer depth with density     criteria: rho = rho(10m or surf) + 0.03(or 0.01) 
     61      !!      the mixed layer depth with temperature criteria: abs( tn - tn(10m) ) = 0.2        
     62      !!      the top of the thermochine: tn = tn(10m) - ztem2  
     63      !!      the pycnocline depth with density criteria equivalent to a temperature variation  
     64      !!                rho = rho10m + (dr/dT)(T,S,10m)*(-0.2 degC)  
     65      !!      the barrier layer thickness 
     66      !!      the maximal verical inversion of temperature and its depth max( 0, max of tn - tn(10m) ) 
     67      !!      the depth of the 20 degree isotherm (linear interpolation) 
     68      !!      the depth of the 28 degree isotherm (linear interpolation) 
     69      !!      the heat content of first 300 m 
    6070      !! 
    6171      !! ** Method :  
     
    6676      INTEGER                          ::   ji, jj, jk            ! dummy loop arguments 
    6777      INTEGER                          ::   iid, iif, ilevel      ! temporary integers 
    68       INTEGER, DIMENSION(jpi,jpj)      ::   ikc                   ! levels 
    69       REAL(wp)                         ::   zd, zthick_0, zcoef   ! temporary scalars 
    70       REAL(wp), DIMENSION(jpi,jpj)     ::   zthick 
    71       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zdzt 
     78      INTEGER, DIMENSION(jpi,jpj)      ::   ik20, ik28            ! levels 
     79      REAL(wp)                         ::   zavt5 = 5.e-4_wp      ! Kz criterion for the turbocline depth 
     80      REAL(wp)                         ::   zrho3 = 0.03_wp       ! density     criterion for mixed layer depth 
     81      REAL(wp)                         ::   zrho1 = 0.01_wp       ! density     criterion for mixed layer depth 
     82      REAL(wp)                         ::   ztem2 = 0.2_wp        ! temperature criterion for mixed layer depth 
     83      REAL(wp)                         ::   zthick_0, zcoef       ! temporary scalars 
     84      REAL(wp)                         ::   zztmp, zzdep          ! temporary scalars inside do loop 
     85      REAL(wp)                         ::   zu, zv, zw, zut, zvt  ! temporary workspace 
     86      REAL(wp), DIMENSION(jpi,jpj)     ::   zabs2                 ! MLD: abs( tn - tn(10m) ) = ztem2  
     87      REAL(wp), DIMENSION(jpi,jpj)     ::   ztm2                  ! Top of thermocline: tn = tn(10m) - ztem2      
     88      REAL(wp), DIMENSION(jpi,jpj)     ::   zrho10_3              ! MLD: rho = rho10m + zrho3       
     89      REAL(wp), DIMENSION(jpi,jpj)     ::   zpycn                 ! pycnocline: rho = rho10m + (dr/dT)(T,S,10m)*(-0.2 degC) 
     90      REAL(wp), DIMENSION(jpi,jpj)     ::   ztinv                 ! max of temperature inversion 
     91      REAL(wp), DIMENSION(jpi,jpj)     ::   zdepinv               ! depth of temperature inversion 
     92      REAL(wp), DIMENSION(jpi,jpj)     ::   zrho0_3               ! MLD rho = rho(surf) = 0.03 
     93      REAL(wp), DIMENSION(jpi,jpj)     ::   zrho0_1               ! MLD rho = rho(surf) = 0.01 
     94      REAL(wp), DIMENSION(jpi,jpj)     ::   zmaxdzT               ! max of dT/dz 
     95      REAL(wp), DIMENSION(jpi,jpj)     ::   zthick                ! vertical integration thickness  
     96      REAL(wp), DIMENSION(jpi,jpj)     ::   zdelr                 ! delta rho equivalent to deltaT = 0.2 
    7297      !!---------------------------------------------------------------------- 
    7398 
     
    79104      ENDIF 
    80105 
    81       ! -------------------------- ! 
    82       !  Depth of the thermocline  ! 
    83       ! -------------------------- ! 
    84       ! The depth of the thermocline is defined as the depth of the  
    85       ! strongest vertical temperature gradient 
    86       zdzt(:,:,1) = 0.e0 
    87       DO jk = 2, jpk                      ! vertical gradient of temperature 
    88          zdzt(:,:,jk) = ( tn(:,:,jk-1) - tn(:,:,jk) ) / fse3w(:,:,jk) * tmask(:,:,jk) 
    89       END DO 
     106      ! initialization 
     107      ztinv  (:,:) = 0.e0_wp   
     108      zdepinv(:,:) = 0.e0_wp   
     109      zmaxdzT(:,:) = 0.e0_wp   
    90110      DO jj = 1, jpj 
    91111         DO ji = 1, jpi 
    92             ilevel = MAXLOC( zdzt(ji,jj,:), dim= 1 )      ! level of maximum vertical temperature gradient 
    93             hth(ji,jj) = fsdepw(ji,jj,ilevel)             ! depth of the thermocline 
    94          END DO          
    95       END DO 
    96  
    97       CALL iom_put( "thermod", hth )      ! depth of the thermocline 
    98  
    99       ! ----------------------- ! 
    100       !  Depth of 20C isotherm  ! 
    101       ! ----------------------- ! 
     112            zztmp = bathy(ji,jj) 
     113            hth     (ji,jj) = zztmp 
     114            hmld    (ji,jj) = zztmp 
     115            zabs2   (ji,jj) = zztmp 
     116            ztm2    (ji,jj) = zztmp 
     117            zrho10_3(ji,jj) = zztmp 
     118            zpycn   (ji,jj) = zztmp 
     119        END DO 
     120      END DO 
     121      IF( nla10 > 1 ) THEN  
     122         DO jj = 1, jpj 
     123            DO ji = 1, jpi 
     124               zztmp = bathy(ji,jj) 
     125               zrho0_3(ji,jj) = zztmp 
     126               zrho0_1(ji,jj) = zztmp 
     127            END DO 
     128         END DO 
     129      ENDIF 
    102130       
    103       ! search last level above 20C (beware temperature is not always decreasing with depth) 
    104       ikc(:,:) = MAXLOC( fsdept(:,:,1:jpkm1), mask = tn(:,:,1:jpkm1) >= 20., dim = 3 )   !   ikc = 0 if mask always equal to false 
    105       ! Depth of 20C isotherm, linear interpolation 
     131      ! Preliminary computation 
     132      ! computation of zdelr = (dr/dT)(T,S,10m)*(-0.2 degC) 
     133      DO jj=1, jpj 
     134         DO ji=1, jpi 
     135            IF( tmask(ji,jj,nla10) == 1. ) THEN 
     136               zu  =  1779.50 + 11.250*tn(ji,jj,nla10) - 3.80*sn(ji,jj,nla10) - 0.0745*tn(ji,jj,nla10)*tn(ji,jj,nla10)   & 
     137                  &                                                           - 0.0100*tn(ji,jj,nla10)*sn(ji,jj,nla10) 
     138               zv  =  5891.00 + 38.000*tn(ji,jj,nla10) + 3.00*sn(ji,jj,nla10) - 0.3750*tn(ji,jj,nla10)*tn(ji,jj,nla10) 
     139               zut =    11.25 -  0.149*tn(ji,jj,nla10) - 0.01*sn(ji,jj,nla10) 
     140               zvt =    38.00 -  0.750*tn(ji,jj,nla10) 
     141               zw  = (zu + 0.698*zv) * (zu + 0.698*zv) 
     142               zdelr(ji,jj) = ztem2 * (1000.*(zut*zv - zvt*zu)/zw) 
     143            ELSE 
     144               zdelr(ji,jj) = 0.e0 
     145            ENDIF 
     146         END DO 
     147      END DO 
     148 
     149      ! ------------------------------------------------------------- ! 
     150      ! thermocline depth: strongest vertical gradient of temperature ! 
     151      ! turbocline depth (mixing layer depth): avt = zavt5            ! 
     152      ! MLD: rho = rho(1) + zrho3                                     ! 
     153      ! MLD: rho = rho(1) + zrho1                                     ! 
     154      ! ------------------------------------------------------------- ! 
     155      DO jk = jpkm1, 2, -1   ! loop from bottom to 2 
     156         DO jj = 1, jpj 
     157            DO ji = 1, jpi 
     158 
     159               zzdep = fsdepw(ji,jj,jk) 
     160               zztmp = ( tn(ji,jj,jk-1) - tn(ji,jj,jk) ) / zzdep * tmask(ji,jj,jk)   ! vertical gradient of temperature (dT/dz) 
     161               zzdep = zzdep * tmask(ji,jj,1) 
     162 
     163               IF( zztmp > zmaxdzT(ji,jj) ) THEN                         
     164                  zmaxdzT(ji,jj) = zztmp   ;   hth    (ji,jj) = zzdep                ! max and depth of dT/dz 
     165               ENDIF 
     166                
     167               IF( avt (ji,jj,jk) < zavt5 )    hmld   (ji,jj) = zzdep                ! avt < zavt5 
     168 
     169               IF( nla10 > 1 ) THEN  
     170                  zztmp = rhop(ji,jj,jk) - rhop(ji,jj,1)                             ! delta rho(1) 
     171                  IF( zztmp > zrho3 )          zrho0_3(ji,jj) = zzdep                ! > 0.03 
     172                  IF( zztmp > zrho1 )          zrho0_1(ji,jj) = zzdep                ! > 0.01 
     173               ENDIF 
     174 
     175            END DO 
     176         END DO 
     177      END DO 
     178       
     179      CALL iom_put( "mlddzt", hth )            ! depth of the thermocline 
     180      CALL iom_put( "mldkz5", hmld )           ! turbocline depth 
     181      IF( nla10 > 1 ) THEN  
     182         CALL iom_put( "mldr0_3", zrho0_3 )   ! MLD delta rho(surf) = 0.03 
     183         CALL iom_put( "mldr0_1", zrho0_1 )   ! MLD delta rho(surf) = 0.01 
     184      ENDIF 
     185 
     186      ! ------------------------------------------------------------- ! 
     187      ! MLD: abs( tn - tn(10m) ) = ztem2                              ! 
     188      ! Top of thermocline: tn = tn(10m) - ztem2                      ! 
     189      ! MLD: rho = rho10m + zrho3                                     ! 
     190      ! pycnocline: rho = rho10m + (dr/dT)(T,S,10m)*(-0.2 degC)       ! 
     191      ! temperature inversion: max( 0, max of tn - tn(10m) )          ! 
     192      ! depth of temperature inversion                                ! 
     193      ! ------------------------------------------------------------- ! 
     194      DO jk = jpkm1, nlb10, -1   ! loop from bottom to nlb10 
     195         DO jj = 1, jpj 
     196            DO ji = 1, jpi 
     197 
     198               zzdep = fsdepw(ji,jj,jk) * tmask(ji,jj,1) 
     199 
     200               zztmp = tn(ji,jj,nla10) - tn(ji,jj,jk)                  ! - delta T(10m) 
     201               IF( ABS(zztmp) > ztem2 )      zabs2   (ji,jj) = zzdep   ! abs > 0.2 
     202               IF(     zztmp  > ztem2 )      ztm2    (ji,jj) = zzdep   ! > 0.2 
     203               zztmp = -zztmp                                          ! delta T(10m) 
     204               IF( zztmp >  ztinv(ji,jj) ) THEN                        ! temperature inversion 
     205                  ztinv(ji,jj) = zztmp   ;   zdepinv (ji,jj) = zzdep   ! max value and depth 
     206               ENDIF 
     207 
     208               zztmp = rhop(ji,jj,jk) - rhop(ji,jj,nla10)              ! delta rho(10m) 
     209               IF( zztmp > zrho3        )    zrho10_3(ji,jj) = zzdep   ! > 0.03 
     210               IF( zztmp > zdelr(ji,jj) )    zpycn   (ji,jj) = zzdep   ! > equi. delta T(10m) - 0.2 
     211 
     212            END DO 
     213         END DO 
     214      END DO 
     215 
     216      CALL iom_put( "mld|dt|" , zabs2        )   ! MLD abs(delta t) - 0.2 
     217      CALL iom_put( "topthdep", ztm2         )   ! T(10) - 0.2 
     218      CALL iom_put( "mldr10_3", zrho10_3     )   ! MLD delta rho(10m) = 0.03 
     219      CALL iom_put( "pycndep" , zpycn        )   ! MLD delta rho equi. delta T(10m) = 0.2 
     220      CALL iom_put( "BLT"     , ztm2 - zpycn )   ! Barrier Layer Thickness 
     221      CALL iom_put( "tinv"    , ztinv        )   ! max. temp. inv. (t10 ref)  
     222      CALL iom_put( "depti"   , zdepinv      )   ! depth of max. temp. inv. (t10 ref)  
     223 
     224 
     225      ! ----------------------------------- ! 
     226      ! search deepest level above 20C/28C  ! 
     227      ! ----------------------------------- ! 
     228      ik20(:,:) = 1 
     229      ik28(:,:) = 1 
     230      DO jk = 1, jpkm1   ! beware temperature is not always decreasing with depth => loop from top to bottom 
     231         DO jj = 1, jpj 
     232            DO ji = 1, jpi 
     233               zztmp = tn(ji,jj,jk) 
     234               IF( zztmp >= 20. )   ik20(ji,jj) = jk 
     235               IF( zztmp >= 28. )   ik28(ji,jj) = jk 
     236            END DO 
     237         END DO 
     238      END DO 
     239 
     240      ! --------------------------- ! 
     241      !  Depth of 20C/28C isotherm  ! 
     242      ! --------------------------- ! 
    106243      DO jj = 1, jpj 
    107244         DO ji = 1, jpi 
    108             iid = MAX( 1, ikc(ji,jj) )    ! make sure that iid /= 0 (these points are masked later) 
    109             zd = fsdept(ji,jj,iid) + (    fsdept(ji,jj,iid+1) - fsdept(ji,jj,iid)                           )   & 
    110                &                   * ( 20.*tmask(ji,jj,iid+1) -     tn(ji,jj,iid)                           )   & 
    111                &                   / (        tn(ji,jj,iid+1) -     tn(ji,jj,iid) + ( 1. - tmask(ji,jj,1) ) ) 
    112             ! bound by the ocean depth, minimum value, first T-point depth 
     245 
    113246            iif = mbathy(ji,jj) 
    114             hd20(ji,jj) = MIN( zd*tmask(ji,jj,1), fsdepw(ji,jj,iif) ) 
    115          END DO 
    116       END DO 
    117       WHERE(ikc == 0 )   hd20 = 0.e0      ! set to 0 points where tn is always < 20 (not very good for temporal mean...) 
    118       CALL iom_put( "20d", hd20 )         ! depth of the 20 isotherm 
    119  
    120       ! ----------------------- ! 
    121       !  Depth of 28C isotherm  !  
    122       ! ----------------------- ! 
    123        
    124       ! search last level above 28C (beware temperature is not always decreasing with depth) 
    125       ikc(:,:) = MAXLOC( fsdept(:,:,1:jpkm1), mask = tn(:,:,1:jpkm1) >= 28., dim = 3 )   !   ikc = 0 if mask always equal to false 
    126       ! Depth of 28C isotherm, linear interpolation 
    127       DO jj = 1, jpj 
    128          DO ji = 1, jpi 
    129             iid = MAX( 1, ikc(ji,jj) )    ! make sure that iid /= 0 (these points are masked later) 
    130             zd = fsdept(ji,jj,iid) + (    fsdept(ji,jj,iid+1) - fsdept(ji,jj,iid)                           )   & 
    131                &                   * ( 28.*tmask(ji,jj,iid+1) -     tn(ji,jj,iid)                           )   & 
    132                &                   / (        tn(ji,jj,iid+1) -     tn(ji,jj,iid) + ( 1. - tmask(ji,jj,1) ) ) 
    133             ! bound by the ocean depth, minimum value, first T-point depth 
    134             iif = mbathy(ji,jj) 
    135             hd28(ji,jj) = MIN( zd*tmask(ji,jj,1), fsdepw(ji,jj,iif) ) 
    136          END DO 
    137       END DO 
    138       WHERE(ikc == 0 )   hd28 = 0.e0      ! set to 0 points where tn is always < 28 (not very good for temporal mean...) 
    139       CALL iom_put( "28d", hd28 )         ! depth of the 28 isotherm 
     247            zzdep = fsdepw(ji,jj,iif) 
     248 
     249            iid = ik20(ji,jj) 
     250            IF( iid /= 1 ) THEN  
     251               ! linear interpolation 
     252               zztmp =      fsdept(ji,jj,iid  )   & 
     253                  &  + (    fsdept(ji,jj,iid+1) - fsdept(ji,jj,iid)                       )   & 
     254                  &  * ( 20.*tmask(ji,jj,iid+1) -     tn(ji,jj,iid)                       )   & 
     255                  &  / (        tn(ji,jj,iid+1) -     tn(ji,jj,iid) + (1.-tmask(ji,jj,1)) ) 
     256               ! bound by the ocean depth, minimum value, first T-point depth 
     257               hd20(ji,jj) = MIN( zztmp*tmask(ji,jj,1), zzdep) 
     258            ELSE  
     259               hd20(ji,jj)=0. 
     260            ENDIF 
     261 
     262            iid = ik28(ji,jj) 
     263            IF( iid /= 1 ) THEN  
     264               ! linear interpolation 
     265               zztmp =      fsdept(ji,jj,iid  )   & 
     266                  &  + (    fsdept(ji,jj,iid+1) - fsdept(ji,jj,iid)                       )   & 
     267                  &  * ( 28.*tmask(ji,jj,iid+1) -     tn(ji,jj,iid)                       )   & 
     268                  &  / (        tn(ji,jj,iid+1) -     tn(ji,jj,iid) + (1.-tmask(ji,jj,1)) ) 
     269               ! bound by the ocean depth, minimum value, first T-point depth 
     270               hd28(ji,jj) = MIN( zztmp*tmask(ji,jj,1), zzdep ) 
     271            ELSE  
     272               hd28(ji,jj) = 0. 
     273            ENDIF 
     274 
     275         END DO 
     276      END DO 
     277      CALL iom_put( "20d", hd20 )   ! depth of the 20 isotherm 
     278      CALL iom_put( "28d", hd28 )   ! depth of the 28 isotherm 
    140279 
    141280      ! ----------------------------- ! 
     
    145284      ! find ilevel with (ilevel+1) the deepest W-level above 300m (we assume we can use e3t_0 to do this search...) 
    146285      ilevel = 0 
    147       zthick_0 = 0.e0 
     286      zthick_0 = 0.e0_wp 
    148287      DO jk = 1, jpkm1                       
    149288         zthick_0 = zthick_0 + e3t_0(jk) 
     
    151290      END DO 
    152291      ! surface boundary condition 
    153       IF( lk_vvl ) THEN   ;   zthick(:,:) = 0.e0        ;   htc3(:,:) = 0.e0                                      
     292      IF( lk_vvl ) THEN   ;   zthick(:,:) = 0.e0_wp     ;   htc3(:,:) = 0.e0_wp                                    
    154293      ELSE                ;   zthick(:,:) = sshn(:,:)   ;   htc3(:,:) = tn(:,:,jk) * sshn(:,:) * tmask(:,:,jk)    
    155294      ENDIF 
     
    169308      zcoef = rau0 * rcp 
    170309      htc3(:,:) = zcoef * htc3(:,:) 
    171       CALL iom_put( "hc300", htc3 )      ! first 300m heaat content 
     310      CALL iom_put( "hc300", htc3 )      ! first 300m heat content 
    172311 
    173312 
Note: See TracChangeset for help on using the changeset viewer.