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 12377 for NEMO/trunk/src/OCE/DIA/diahth.F90 – NEMO

Ignore:
Timestamp:
2020-02-12T15:39:06+01:00 (4 years ago)
Author:
acc
Message:

The big one. Merging all 2019 developments from the option 1 branch back onto the trunk.

This changeset reproduces 2019/dev_r11943_MERGE_2019 on the trunk using a 2-URL merge
onto a working copy of the trunk. I.e.:

svn merge --ignore-ancestry \

svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/trunk \
svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/branches/2019/dev_r11943_MERGE_2019 ./

The --ignore-ancestry flag avoids problems that may otherwise arise from the fact that
the merge history been trunk and branch may have been applied in a different order but
care has been taken before this step to ensure that all applicable fixes and updates
are present in the merge branch.

The trunk state just before this step has been branched to releases/release-4.0-HEAD
and that branch has been immediately tagged as releases/release-4.0.2. Any fixes
or additions in response to tickets on 4.0, 4.0.1 or 4.0.2 should be done on
releases/release-4.0-HEAD. From now on future 'point' releases (e.g. 4.0.2) will
remain unchanged with periodic releases as needs demand. Note release-4.0-HEAD is a
transitional naming convention. Future full releases, say 4.2, will have a release-4.2
branch which fulfills this role and the first point release (e.g. 4.2.0) will be made
immediately following the release branch creation.

2020 developments can be started from any trunk revision later than this one.

Location:
NEMO/trunk
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • NEMO/trunk

    • Property svn:externals
      •  

        old new  
        33^/utils/build/mk@HEAD         mk 
        44^/utils/tools@HEAD            tools 
        5 ^/vendors/AGRIF/dev@HEAD      ext/AGRIF 
         5^/vendors/AGRIF/dev_r11615_ENHANCE-04_namelists_as_internalfiles_agrif@HEAD      ext/AGRIF 
        66^/vendors/FCM@HEAD            ext/FCM 
        77^/vendors/IOIPSL@HEAD         ext/IOIPSL 
  • NEMO/trunk/src/OCE/DIA/diahth.F90

    r12276 r12377  
    4040 
    4141 
     42   !! * Substitutions 
     43#  include "do_loop_substitute.h90" 
    4244   !!---------------------------------------------------------------------- 
    4345   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
     
    6163 
    6264 
    63    SUBROUTINE dia_hth( kt ) 
     65   SUBROUTINE dia_hth( kt, Kmm ) 
    6466      !!--------------------------------------------------------------------- 
    6567      !!                  ***  ROUTINE dia_hth  *** 
     
    8284      !!------------------------------------------------------------------- 
    8385      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index 
     86      INTEGER, INTENT( in ) ::   Kmm     ! ocean time level index 
    8487      !! 
    8588      INTEGER                      ::   ji, jj, jk            ! dummy loop arguments 
     
    126129            zdepinv(:,:) = 0._wp   
    127130            zmaxdzT(:,:) = 0._wp   
    128             DO jj = 1, jpj 
    129                DO ji = 1, jpi 
    130                   zztmp = gdepw_n(ji,jj,mbkt(ji,jj)+1)  
    131                   hth     (ji,jj) = zztmp 
    132                   zabs2   (ji,jj) = zztmp 
    133                   ztm2    (ji,jj) = zztmp 
    134                   zrho10_3(ji,jj) = zztmp 
    135                   zpycn   (ji,jj) = zztmp 
    136                  END DO 
    137             END DO 
     131            DO_2D_11_11 
     132               zztmp = gdepw(ji,jj,mbkt(ji,jj)+1,Kmm)  
     133               hth     (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_2D 
    138139            IF( nla10 > 1 ) THEN  
    139                DO jj = 1, jpj 
    140                   DO ji = 1, jpi 
    141                      zztmp = gdepw_n(ji,jj,mbkt(ji,jj)+1)  
    142                      zrho0_3(ji,jj) = zztmp 
    143                      zrho0_1(ji,jj) = zztmp 
    144                   END DO 
    145                END DO 
     140               DO_2D_11_11 
     141                  zztmp = gdepw(ji,jj,mbkt(ji,jj)+1,Kmm)  
     142                  zrho0_3(ji,jj) = zztmp 
     143                  zrho0_1(ji,jj) = zztmp 
     144               END_2D 
    146145            ENDIF 
    147146       
    148147            ! Preliminary computation 
    149148            ! computation of zdelr = (dr/dT)(T,S,10m)*(-0.2 degC) 
    150             DO jj = 1, jpj 
    151                DO ji = 1, jpi 
    152                   IF( tmask(ji,jj,nla10) == 1. ) THEN 
    153                      zu  =  1779.50 + 11.250 * tsn(ji,jj,nla10,jp_tem) - 3.80   * tsn(ji,jj,nla10,jp_sal)  & 
    154                         &           - 0.0745 * tsn(ji,jj,nla10,jp_tem) * tsn(ji,jj,nla10,jp_tem)   & 
    155                         &           - 0.0100 * tsn(ji,jj,nla10,jp_tem) * tsn(ji,jj,nla10,jp_sal) 
    156                      zv  =  5891.00 + 38.000 * tsn(ji,jj,nla10,jp_tem) + 3.00   * tsn(ji,jj,nla10,jp_sal)  & 
    157                         &           - 0.3750 * tsn(ji,jj,nla10,jp_tem) * tsn(ji,jj,nla10,jp_tem) 
    158                      zut =    11.25 -  0.149 * tsn(ji,jj,nla10,jp_tem) - 0.01   * tsn(ji,jj,nla10,jp_sal) 
    159                      zvt =    38.00 -  0.750 * tsn(ji,jj,nla10,jp_tem) 
    160                      zw  = (zu + 0.698*zv) * (zu + 0.698*zv) 
    161                      zdelr(ji,jj) = ztem2 * (1000.*(zut*zv - zvt*zu)/zw) 
    162                   ELSE 
    163                      zdelr(ji,jj) = 0._wp 
    164                   ENDIF 
    165                END DO 
    166             END DO 
     149            DO_2D_11_11 
     150               IF( tmask(ji,jj,nla10) == 1. ) THEN 
     151                  zu  =  1779.50 + 11.250 * ts(ji,jj,nla10,jp_tem,Kmm) - 3.80   * ts(ji,jj,nla10,jp_sal,Kmm)  & 
     152                     &           - 0.0745 * ts(ji,jj,nla10,jp_tem,Kmm) * ts(ji,jj,nla10,jp_tem,Kmm)   & 
     153                     &           - 0.0100 * ts(ji,jj,nla10,jp_tem,Kmm) * ts(ji,jj,nla10,jp_sal,Kmm) 
     154                  zv  =  5891.00 + 38.000 * ts(ji,jj,nla10,jp_tem,Kmm) + 3.00   * ts(ji,jj,nla10,jp_sal,Kmm)  & 
     155                     &           - 0.3750 * ts(ji,jj,nla10,jp_tem,Kmm) * ts(ji,jj,nla10,jp_tem,Kmm) 
     156                  zut =    11.25 -  0.149 * ts(ji,jj,nla10,jp_tem,Kmm) - 0.01   * ts(ji,jj,nla10,jp_sal,Kmm) 
     157                  zvt =    38.00 -  0.750 * ts(ji,jj,nla10,jp_tem,Kmm) 
     158                  zw  = (zu + 0.698*zv) * (zu + 0.698*zv) 
     159                  zdelr(ji,jj) = ztem2 * (1000.*(zut*zv - zvt*zu)/zw) 
     160               ELSE 
     161                  zdelr(ji,jj) = 0._wp 
     162               ENDIF 
     163            END_2D 
    167164 
    168165            ! ------------------------------------------------------------- ! 
     
    172169            ! MLD: rho = rho(1) + zrho1                                     ! 
    173170            ! ------------------------------------------------------------- ! 
    174             DO jk = jpkm1, 2, -1   ! loop from bottom to 2 
    175                DO jj = 1, jpj 
    176                   DO ji = 1, jpi 
    177                      ! 
    178                      zzdep = gdepw_n(ji,jj,jk) 
    179                      zztmp = ( tsn(ji,jj,jk-1,jp_tem) - tsn(ji,jj,jk,jp_tem) ) & 
    180                             & / zzdep * tmask(ji,jj,jk)   ! vertical gradient of temperature (dT/dz) 
    181                      zzdep = zzdep * tmask(ji,jj,1) 
    182  
    183                      IF( zztmp > zmaxdzT(ji,jj) ) THEN                         
    184                          zmaxdzT(ji,jj) = zztmp    
    185                          hth    (ji,jj) = zzdep                ! max and depth of dT/dz 
    186                      ENDIF 
    187                 
    188                      IF( nla10 > 1 ) THEN  
    189                         zztmp = rhop(ji,jj,jk) - rhop(ji,jj,1)                       ! delta rho(1) 
    190                         IF( zztmp > zrho3 )   zrho0_3(ji,jj) = zzdep                ! > 0.03 
    191                         IF( zztmp > zrho1 )   zrho0_1(ji,jj) = zzdep                ! > 0.01 
    192                      ENDIF 
    193                   END DO 
    194                END DO 
    195             END DO 
     171            DO_3DS_11_11( jpkm1, 2, -1 ) 
     172               ! 
     173               zzdep = gdepw(ji,jj,jk,Kmm) 
     174               zztmp = ( ts(ji,jj,jk-1,jp_tem,Kmm) - ts(ji,jj,jk,jp_tem,Kmm) ) & 
     175                      & / zzdep * tmask(ji,jj,jk)   ! vertical gradient of temperature (dT/dz) 
     176               zzdep = zzdep * tmask(ji,jj,1) 
     177 
     178               IF( zztmp > zmaxdzT(ji,jj) ) THEN                         
     179                   zmaxdzT(ji,jj) = zztmp    
     180                   hth    (ji,jj) = zzdep                ! max and depth of dT/dz 
     181               ENDIF 
     182          
     183               IF( nla10 > 1 ) THEN  
     184                  zztmp = rhop(ji,jj,jk) - rhop(ji,jj,1)                       ! delta rho(1) 
     185                  IF( zztmp > zrho3 )   zrho0_3(ji,jj) = zzdep                ! > 0.03 
     186                  IF( zztmp > zrho1 )   zrho0_1(ji,jj) = zzdep                ! > 0.01 
     187               ENDIF 
     188            END_3D 
    196189          
    197190            CALL iom_put( 'mlddzt', hth )            ! depth of the thermocline 
     
    213206            ! depth of temperature inversion                                ! 
    214207            ! ------------------------------------------------------------- ! 
    215             DO jk = jpkm1, nlb10, -1   ! loop from bottom to nlb10 
    216                DO jj = 1, jpj 
    217                   DO ji = 1, jpi 
    218                      ! 
    219                      zzdep = gdepw_n(ji,jj,jk) * tmask(ji,jj,1) 
    220                      ! 
    221                      zztmp = tsn(ji,jj,nla10,jp_tem) - tsn(ji,jj,jk,jp_tem)  ! - delta T(10m) 
    222                      IF( ABS(zztmp) > ztem2 )      zabs2   (ji,jj) = zzdep   ! abs > 0.2 
    223                      IF(     zztmp  > ztem2 )      ztm2    (ji,jj) = zzdep   ! > 0.2 
    224                      zztmp = -zztmp                                          ! delta T(10m) 
    225                      IF( zztmp >  ztinv(ji,jj) ) THEN                        ! temperature inversion 
    226                         ztinv(ji,jj) = zztmp    
    227                         zdepinv (ji,jj) = zzdep   ! max value and depth 
    228                      ENDIF 
    229  
    230                      zztmp = rhop(ji,jj,jk) - rhop(ji,jj,nla10)              ! delta rho(10m) 
    231                      IF( zztmp > zrho3        )    zrho10_3(ji,jj) = zzdep   ! > 0.03 
    232                      IF( zztmp > zdelr(ji,jj) )    zpycn   (ji,jj) = zzdep   ! > equi. delta T(10m) - 0.2 
    233                      ! 
    234                   END DO 
    235                END DO 
    236             END DO 
     208            DO_3DS_11_11( jpkm1, nlb10, -1 ) 
     209               ! 
     210               zzdep = gdepw(ji,jj,jk,Kmm) * tmask(ji,jj,1) 
     211               ! 
     212               zztmp = ts(ji,jj,nla10,jp_tem,Kmm) - ts(ji,jj,jk,jp_tem,Kmm)  ! - delta T(10m) 
     213               IF( ABS(zztmp) > ztem2 )      zabs2   (ji,jj) = zzdep   ! abs > 0.2 
     214               IF(     zztmp  > ztem2 )      ztm2    (ji,jj) = zzdep   ! > 0.2 
     215               zztmp = -zztmp                                          ! delta T(10m) 
     216               IF( zztmp >  ztinv(ji,jj) ) THEN                        ! temperature inversion 
     217                  ztinv(ji,jj) = zztmp    
     218                  zdepinv (ji,jj) = zzdep   ! max value and depth 
     219               ENDIF 
     220 
     221               zztmp = rhop(ji,jj,jk) - rhop(ji,jj,nla10)              ! delta rho(10m) 
     222               IF( zztmp > zrho3        )    zrho10_3(ji,jj) = zzdep   ! > 0.03 
     223               IF( zztmp > zdelr(ji,jj) )    zpycn   (ji,jj) = zzdep   ! > equi. delta T(10m) - 0.2 
     224               ! 
     225            END_3D 
    237226 
    238227            CALL iom_put( 'mld_dt02', zabs2    )   ! MLD abs(delta t) - 0.2 
     
    250239         IF( iom_use ('20d') ) THEN  ! depth of the 20 isotherm 
    251240            ztem2 = 20. 
    252             CALL dia_hth_dep( ztem2, hd20 )   
     241            CALL dia_hth_dep( Kmm, ztem2, hd20 )   
    253242            CALL iom_put( '20d', hd20 )     
    254243         ENDIF 
     
    256245         IF( iom_use ('26d') ) THEN  ! depth of the 26 isotherm 
    257246            ztem2 = 26. 
    258             CALL dia_hth_dep( ztem2, hd26 )   
     247            CALL dia_hth_dep( Kmm, ztem2, hd26 )   
    259248            CALL iom_put( '26d', hd26 )     
    260249         ENDIF 
     
    262251         IF( iom_use ('28d') ) THEN  ! depth of the 28 isotherm 
    263252            ztem2 = 28. 
    264             CALL dia_hth_dep( ztem2, hd28 )   
     253            CALL dia_hth_dep( Kmm, ztem2, hd28 )   
    265254            CALL iom_put( '28d', hd28 )     
    266255         ENDIF 
     
    271260         IF( iom_use ('hc300') ) THEN   
    272261            zzdep = 300. 
    273             CALL  dia_hth_htc( zzdep, tsn(:,:,:,jp_tem), htc3 ) 
     262            CALL  dia_hth_htc( Kmm, zzdep, ts(:,:,:,jp_tem,Kmm), htc3 ) 
    274263            CALL iom_put( 'hc300', rau0_rcp * htc3 )  ! vertically integrated heat content (J/m2) 
    275264         ENDIF 
     
    280269         IF( iom_use ('hc700') ) THEN   
    281270            zzdep = 700. 
    282             CALL  dia_hth_htc( zzdep, tsn(:,:,:,jp_tem), htc7 ) 
     271            CALL  dia_hth_htc( Kmm, zzdep, ts(:,:,:,jp_tem,Kmm), htc7 ) 
    283272            CALL iom_put( 'hc700', rau0_rcp * htc7 )  ! vertically integrated heat content (J/m2) 
    284273   
     
    290279         IF( iom_use ('hc2000') ) THEN   
    291280            zzdep = 2000. 
    292             CALL  dia_hth_htc( zzdep, tsn(:,:,:,jp_tem), htc20 ) 
     281            CALL  dia_hth_htc( Kmm, zzdep, ts(:,:,:,jp_tem,Kmm), htc20 ) 
    293282            CALL iom_put( 'hc2000', rau0_rcp * htc20 )  ! vertically integrated heat content (J/m2)   
    294283         ENDIF 
     
    301290   END SUBROUTINE dia_hth 
    302291 
    303    SUBROUTINE dia_hth_dep( ptem, pdept ) 
    304       ! 
     292   SUBROUTINE dia_hth_dep( Kmm, ptem, pdept ) 
     293      ! 
     294      INTEGER , INTENT(in) :: Kmm      ! ocean time level index 
    305295      REAL(wp), INTENT(in) :: ptem 
    306296      REAL(wp), DIMENSION(jpi,jpj), INTENT(out) :: pdept      
     
    314304      ! --------------------------------------- ! 
    315305      iktem(:,:) = 1 
    316       DO jk = 1, jpkm1   ! beware temperature is not always decreasing with depth => loop from top to bottom 
    317          DO jj = 1, jpj 
    318             DO ji = 1, jpi 
    319                zztmp = tsn(ji,jj,jk,jp_tem) 
    320                IF( zztmp >= ptem )   iktem(ji,jj) = jk 
    321             END DO 
    322          END DO 
    323       END DO 
     306      DO_3D_11_11( 1, jpkm1 ) 
     307         zztmp = ts(ji,jj,jk,jp_tem,Kmm) 
     308         IF( zztmp >= ptem )   iktem(ji,jj) = jk 
     309      END_3D 
    324310 
    325311      ! ------------------------------- ! 
    326312      !  Depth of ptem isotherm         ! 
    327313      ! ------------------------------- ! 
    328       DO jj = 1, jpj 
    329          DO ji = 1, jpi 
    330             ! 
    331             zzdep = gdepw_n(ji,jj,mbkt(ji,jj)+1)       ! depth of the ocean bottom 
    332             ! 
    333             iid = iktem(ji,jj) 
    334             IF( iid /= 1 ) THEN  
    335                 zztmp =     gdept_n(ji,jj,iid  )   &                     ! linear interpolation 
    336                   &  + (    gdept_n(ji,jj,iid+1) - gdept_n(ji,jj,iid)                       )   & 
    337                   &  * ( 20.*tmask(ji,jj,iid+1) - tsn(ji,jj,iid,jp_tem)                       )   & 
    338                   &  / ( tsn(ji,jj,iid+1,jp_tem) - tsn(ji,jj,iid,jp_tem) + (1.-tmask(ji,jj,1)) ) 
    339                pdept(ji,jj) = MIN( zztmp , zzdep) * tmask(ji,jj,1)       ! bound by the ocean depth 
    340             ELSE  
    341                pdept(ji,jj) = 0._wp 
    342             ENDIF 
    343          END DO 
    344       END DO 
     314      DO_2D_11_11 
     315         ! 
     316         zzdep = gdepw(ji,jj,mbkt(ji,jj)+1,Kmm)       ! depth of the ocean bottom 
     317         ! 
     318         iid = iktem(ji,jj) 
     319         IF( iid /= 1 ) THEN  
     320             zztmp =     gdept(ji,jj,iid  ,Kmm)   &                     ! linear interpolation 
     321               &  + (    gdept(ji,jj,iid+1,Kmm) - gdept(ji,jj,iid,Kmm)                       )   & 
     322               &  * ( 20.*tmask(ji,jj,iid+1) - ts(ji,jj,iid,jp_tem,Kmm)                       )   & 
     323               &  / ( ts(ji,jj,iid+1,jp_tem,Kmm) - ts(ji,jj,iid,jp_tem,Kmm) + (1.-tmask(ji,jj,1)) ) 
     324            pdept(ji,jj) = MIN( zztmp , zzdep) * tmask(ji,jj,1)       ! bound by the ocean depth 
     325         ELSE  
     326            pdept(ji,jj) = 0._wp 
     327         ENDIF 
     328      END_2D 
    345329      ! 
    346330   END SUBROUTINE dia_hth_dep 
    347331 
    348332 
    349    SUBROUTINE dia_hth_htc( pdep, ptn, phtc ) 
    350       ! 
     333   SUBROUTINE dia_hth_htc( Kmm, pdep, pt, phtc ) 
     334      ! 
     335      INTEGER , INTENT(in) :: Kmm      ! ocean time level index 
    351336      REAL(wp), INTENT(in) :: pdep     ! depth over the heat content 
    352       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in) :: ptn    
     337      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in) :: pt    
    353338      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: phtc   
    354339      ! 
     
    361346       
    362347      IF( .NOT. ln_linssh ) THEN   ;   zthick(:,:) = 0._wp       ;   phtc(:,:) = 0._wp                                    
    363       ELSE                         ;   zthick(:,:) = sshn(:,:)   ;   phtc(:,:) = ptn(:,:,1) * sshn(:,:) * tmask(:,:,1)    
     348      ELSE                         ;   zthick(:,:) = ssh(:,:,Kmm)   ;   phtc(:,:) = pt(:,:,1) * ssh(:,:,Kmm) * tmask(:,:,1)    
    364349      ENDIF 
    365350      ! 
    366351      ilevel(:,:) = 1 
    367       DO jk = 2, jpkm1 
    368          DO jj = 1, jpj 
    369             DO ji = 1, jpi 
    370                IF( ( gdept_n(ji,jj,jk) < pdep ) .AND. ( tmask(ji,jj,jk) == 1 ) ) THEN 
    371                    ilevel(ji,jj) = jk 
    372                    zthick(ji,jj) = zthick(ji,jj) + e3t_n(ji,jj,jk) 
    373                    phtc  (ji,jj) = phtc  (ji,jj) + e3t_n(ji,jj,jk) * ptn(ji,jj,jk) 
    374                ENDIF 
    375             ENDDO 
    376          ENDDO 
    377       ENDDO 
    378       ! 
    379       DO jj = 1, jpj 
    380          DO ji = 1, jpi 
    381             ik = ilevel(ji,jj) 
    382             zthick(ji,jj) = pdep - zthick(ji,jj)   !   remaining thickness to reach depht pdep 
    383             phtc(ji,jj)   = phtc(ji,jj) + ptn(ji,jj,ik+1) * MIN( e3t_n(ji,jj,ik+1), zthick(ji,jj) ) & 
    384                                                           * tmask(ji,jj,ik+1) 
    385          END DO 
    386       ENDDO 
     352      DO_3D_11_11( 2, jpkm1 ) 
     353         IF( ( gdept(ji,jj,jk,Kmm) < pdep ) .AND. ( tmask(ji,jj,jk) == 1 ) ) THEN 
     354             ilevel(ji,jj) = jk 
     355             zthick(ji,jj) = zthick(ji,jj) + e3t(ji,jj,jk,Kmm) 
     356             phtc  (ji,jj) = phtc  (ji,jj) + e3t(ji,jj,jk,Kmm) * pt(ji,jj,jk) 
     357         ENDIF 
     358      END_3D 
     359      ! 
     360      DO_2D_11_11 
     361         ik = ilevel(ji,jj) 
     362         zthick(ji,jj) = pdep - zthick(ji,jj)   !   remaining thickness to reach depht pdep 
     363         phtc(ji,jj)   = phtc(ji,jj) + pt(ji,jj,ik+1) * MIN( e3t(ji,jj,ik+1,Kmm), zthick(ji,jj) ) & 
     364                                                       * tmask(ji,jj,ik+1) 
     365      END_2D 
    387366      ! 
    388367      ! 
Note: See TracChangeset for help on using the changeset viewer.