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 13710 for NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/OCE/DIA – NEMO

Ignore:
Timestamp:
2020-11-02T10:56:42+01:00 (4 years ago)
Author:
emanuelaclementi
Message:

branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves: merge with trunk@13708, see #2155 and #2339

Location:
NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves
Files:
9 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves

    • 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_r12970_AGRIF_CMEMS      ext/AGRIF 
        66^/vendors/FCM@HEAD            ext/FCM 
        77^/vendors/IOIPSL@HEAD         ext/IOIPSL 
        88 
        99# SETTE 
        10 ^/utils/CI/sette@HEAD         sette 
         10^/utils/CI/sette@13559        sette 
  • NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/OCE/DIA/diaar5.F90

    r12630 r13710  
    3939   !! * Substitutions 
    4040#  include "do_loop_substitute.h90" 
     41#  include "domzgr_substitute.h90" 
    4142   !!---------------------------------------------------------------------- 
    4243   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
     
    7677      ! 
    7778      REAL(wp), ALLOCATABLE, DIMENSION(:,:)     :: zarea_ssh , zbotpres       ! 2D workspace  
    78       REAL(wp), ALLOCATABLE, DIMENSION(:,:)     :: zpe, z2d                   ! 2D workspace  
    79       REAL(wp), ALLOCATABLE, DIMENSION(:,:,:)   :: zrhd , zrhop, ztpot   ! 3D workspace 
     79      REAL(wp), ALLOCATABLE, DIMENSION(:,:)     :: z2d, zpe                   ! 2D workspace  
     80      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:)   :: z3d, zrhd, ztpot, zgdept   ! 3D workspace (zgdept: needed to use the substitute) 
    8081      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) :: ztsn                       ! 4D workspace 
    8182 
     
    8788      IF( l_ar5 ) THEN  
    8889         ALLOCATE( zarea_ssh(jpi,jpj), zbotpres(jpi,jpj), z2d(jpi,jpj) ) 
    89          ALLOCATE( zrhd(jpi,jpj,jpk) , zrhop(jpi,jpj,jpk) ) 
     90         ALLOCATE( zrhd(jpi,jpj,jpk) ) 
    9091         ALLOCATE( ztsn(jpi,jpj,jpk,jpts) ) 
    9192         zarea_ssh(:,:) = e1e2t(:,:) * ssh(:,:,Kmm) 
     
    101102            zrhd(:,:,jk) = e1e2t(:,:) * e3t(:,:,jk,Kmm) * tmask(:,:,jk) 
    102103         END DO 
     104         DO jk = 1, jpk 
     105            z3d(:,:,jk) =  rho0 * e3t(:,:,jk,Kmm) * tmask(:,:,jk) 
     106         END DO  
    103107         CALL iom_put( 'volcello'  , zrhd(:,:,:)  )  ! WARNING not consistent with CMIP DR where volcello is at ca. 2000 
    104          CALL iom_put( 'masscello' , rho0 * e3t(:,:,:,Kmm) * tmask(:,:,:) )  ! ocean mass 
     108         CALL iom_put( 'masscello' , z3d (:,:,:) )   ! ocean mass 
    105109      ENDIF  
    106110      ! 
    107111      IF( iom_use( 'e3tb' ) )  THEN    ! bottom layer thickness 
    108          DO_2D_11_11 
     112         DO_2D( 1, 1, 1, 1 ) 
    109113            ikb = mbkt(ji,jj) 
    110114            z2d(ji,jj) = e3t(ji,jj,ikb,Kmm) 
     
    128132         ztsn(:,:,:,jp_tem) = ts(:,:,:,jp_tem,Kmm)                    ! thermosteric ssh 
    129133         ztsn(:,:,:,jp_sal) = sn0(:,:,:) 
    130          CALL eos( ztsn, zrhd, gdept(:,:,:,Kmm) )                       ! now in situ density using initial salinity 
     134         ALLOCATE( zgdept(jpi,jpj,jpk) ) 
     135         DO jk = 1, jpk 
     136            zgdept(:,:,jk) = gdept(:,:,jk,Kmm) 
     137         END DO 
     138         CALL eos( ztsn, zrhd, zgdept)                       ! now in situ density using initial salinity 
    131139         ! 
    132140         zbotpres(:,:) = 0._wp                        ! no atmospheric surface pressure, levitating sea-ice 
     
    136144         IF( ln_linssh ) THEN 
    137145            IF( ln_isfcav ) THEN 
    138                DO ji = 1, jpi 
    139                   DO jj = 1, jpj 
    140                      iks = mikt(ji,jj) 
    141                      zbotpres(ji,jj) = zbotpres(ji,jj) + ssh(ji,jj,Kmm) * zrhd(ji,jj,iks) + riceload(ji,jj) 
    142                   END DO 
    143                END DO 
     146               DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 
     147                  iks = mikt(ji,jj) 
     148                  zbotpres(ji,jj) = zbotpres(ji,jj) + ssh(ji,jj,Kmm) * zrhd(ji,jj,iks) + riceload(ji,jj) 
     149               END_2D 
    144150            ELSE 
    145151               zbotpres(:,:) = zbotpres(:,:) + ssh(:,:,Kmm) * zrhd(:,:,1) 
     
    155161       
    156162         !                                         ! steric sea surface height 
    157          CALL eos( ts(:,:,:,:,Kmm), zrhd, zrhop, gdept(:,:,:,Kmm) )                 ! now in situ and potential density 
    158          zrhop(:,:,jpk) = 0._wp 
    159          CALL iom_put( 'rhop', zrhop ) 
    160          ! 
    161163         zbotpres(:,:) = 0._wp                        ! no atmospheric surface pressure, levitating sea-ice 
    162164         DO jk = 1, jpkm1 
    163             zbotpres(:,:) = zbotpres(:,:) + e3t(:,:,jk,Kmm) * zrhd(:,:,jk) 
     165            zbotpres(:,:) = zbotpres(:,:) + e3t(:,:,jk,Kmm) * rhd(:,:,jk) 
    164166         END DO 
    165167         IF( ln_linssh ) THEN 
     
    168170                  DO jj = 1,jpj 
    169171                     iks = mikt(ji,jj) 
    170                      zbotpres(ji,jj) = zbotpres(ji,jj) + ssh(ji,jj,Kmm) * zrhd(ji,jj,iks) + riceload(ji,jj) 
     172                     zbotpres(ji,jj) = zbotpres(ji,jj) + ssh(ji,jj,Kmm) * rhd(ji,jj,iks) + riceload(ji,jj) 
    171173                  END DO 
    172174               END DO 
    173175            ELSE 
    174                zbotpres(:,:) = zbotpres(:,:) + ssh(:,:,Kmm) * zrhd(:,:,1) 
     176               zbotpres(:,:) = zbotpres(:,:) + ssh(:,:,Kmm) * rhd(:,:,1) 
    175177            END IF 
    176178         END IF 
     
    184186         CALL iom_put( 'botpres', zbotpres ) 
    185187         ! 
     188         DEALLOCATE( zgdept ) 
     189         ! 
    186190      ENDIF 
    187191 
     
    189193          !                                         ! Mean density anomalie, temperature and salinity 
    190194          ztsn(:,:,:,:) = 0._wp                    ! ztsn(:,:,1,jp_tem/sal) is used here as 2D Workspace for temperature & salinity 
    191           DO_3D_11_11( 1, jpkm1 ) 
     195          DO_3D( 1, 1, 1, 1, 1, jpkm1 ) 
    192196             zztmp = e1e2t(ji,jj) * e3t(ji,jj,jk,Kmm) 
    193197             ztsn(ji,jj,1,jp_tem) = ztsn(ji,jj,1,jp_tem) + zztmp * ts(ji,jj,jk,jp_tem,Kmm) 
     
    249253             IF( iom_use( 'tosmint_pot') ) THEN 
    250254               z2d(:,:) = 0._wp 
    251                DO_3D_11_11( 1, jpkm1 ) 
     255               DO_3D( 1, 1, 1, 1, 1, jpkm1 ) 
    252256                  z2d(ji,jj) = z2d(ji,jj) + rho0 * e3t(ji,jj,jk,Kmm) *  ztpot(ji,jj,jk) 
    253257               END_3D 
     
    270274         zpe(:,:) = 0._wp 
    271275         IF( ln_zdfddm ) THEN 
    272             DO_3D_11_11( 2, jpk ) 
     276            DO_3D( 1, 1, 1, 1, 2, jpk ) 
    273277               IF( rn2(ji,jj,jk) > 0._wp ) THEN 
    274278                  zrw = ( gdept(ji,jj,jk,Kmm) - gdepw(ji,jj,jk,Kmm) ) / e3w(ji,jj,jk,Kmm) 
     
    283287            END_3D 
    284288          ELSE 
    285             DO_3D_11_11( 1, jpk ) 
     289            DO_3D( 1, 1, 1, 1, 1, jpk ) 
    286290               zpe(ji,jj) = zpe(ji,jj) + avt(ji,jj,jk) * MIN(0._wp,rn2(ji,jj,jk)) * rho0 * e3w(ji,jj,jk,Kmm) 
    287291            END_3D 
     
    293297      IF( l_ar5 ) THEN 
    294298        DEALLOCATE( zarea_ssh , zbotpres, z2d ) 
    295         DEALLOCATE( zrhd      , zrhop    ) 
    296299        DEALLOCATE( ztsn                 ) 
    297300      ENDIF 
     
    319322     
    320323      z2d(:,:) = puflx(:,:,1)  
    321       DO_3D_00_00( 1, jpkm1 ) 
     324      DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 
    322325         z2d(ji,jj) = z2d(ji,jj) + puflx(ji,jj,jk)  
    323326      END_3D 
    324        CALL lbc_lnk( 'diaar5', z2d, 'U', -1. ) 
     327       CALL lbc_lnk( 'diaar5', z2d, 'U', -1.0_wp ) 
    325328       IF( cptr == 'adv' ) THEN 
    326329          IF( ktra == jp_tem ) CALL iom_put( 'uadv_heattr' , rho0_rcp * z2d )  ! advective heat transport in i-direction 
     
    333336       ! 
    334337       z2d(:,:) = pvflx(:,:,1)  
    335        DO_3D_00_00( 1, jpkm1 ) 
     338       DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 
    336339          z2d(ji,jj) = z2d(ji,jj) + pvflx(ji,jj,jk)  
    337340       END_3D 
    338        CALL lbc_lnk( 'diaar5', z2d, 'V', -1. ) 
     341       CALL lbc_lnk( 'diaar5', z2d, 'V', -1.0_wp ) 
    339342       IF( cptr == 'adv' ) THEN 
    340343          IF( ktra == jp_tem ) CALL iom_put( 'vadv_heattr' , rho0_rcp * z2d )  ! advective heat transport in j-direction 
     
    367370      IF(   iom_use( 'voltot'  ) .OR. iom_use( 'sshtot'    )  .OR. iom_use( 'sshdyn' )  .OR.  &  
    368371         &  iom_use( 'masstot' ) .OR. iom_use( 'temptot'   )  .OR. iom_use( 'saltot' ) .OR.  &     
    369          &  iom_use( 'botpres' ) .OR. iom_use( 'sshthster' )  .OR. iom_use( 'sshsteric' )  ) L_ar5 = .TRUE. 
     372         &  iom_use( 'botpres' ) .OR. iom_use( 'sshthster' )  .OR. iom_use( 'sshsteric' ) .OR. & 
     373         &  iom_use( 'rhop' )  ) L_ar5 = .TRUE. 
    370374   
    371375      IF( l_ar5 ) THEN 
     
    379383         zvol0 (:,:) = 0._wp 
    380384         thick0(:,:) = 0._wp 
    381          DO_3D_11_11( 1, jpkm1 ) 
     385         DO_3D( 1, 1, 1, 1, 1, jpkm1 )   ! interpolation of salinity at the last ocean level (i.e. the partial step) 
    382386            idep = tmask(ji,jj,jk) * e3t_0(ji,jj,jk) 
    383387            zvol0 (ji,jj) = zvol0 (ji,jj) +  idep * e1e2t(ji,jj) 
     
    390394            ALLOCATE( zsaldta(jpi,jpj,jpk,jpts) ) 
    391395            CALL iom_open ( 'sali_ref_clim_monthly', inum ) 
    392             CALL iom_get  ( inum, jpdom_data, 'vosaline' , zsaldta(:,:,:,1), 1  ) 
    393             CALL iom_get  ( inum, jpdom_data, 'vosaline' , zsaldta(:,:,:,2), 12 ) 
     396            CALL iom_get  ( inum, jpdom_global, 'vosaline' , zsaldta(:,:,:,1), 1  ) 
     397            CALL iom_get  ( inum, jpdom_global, 'vosaline' , zsaldta(:,:,:,2), 12 ) 
    394398            CALL iom_close( inum ) 
    395399 
     
    397401            sn0(:,:,:) = sn0(:,:,:) * tmask(:,:,:) 
    398402            IF( ln_zps ) THEN               ! z-coord. partial steps 
    399                DO_2D_11_11 
     403               DO_2D( 1, 1, 1, 1 )          ! interpolation of salinity at the last ocean level (i.e. the partial step) 
    400404                  ik = mbkt(ji,jj) 
    401405                  IF( ik > 1 ) THEN 
  • NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/OCE/DIA/diacfl.F90

    r12489 r13710  
    3434   !! * Substitutions 
    3535#  include "do_loop_substitute.h90" 
     36#  include "domzgr_substitute.h90" 
    3637   !!---------------------------------------------------------------------- 
    3738   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
     
    5556      INTEGER , DIMENSION(3)           ::   iloc_u , iloc_v , iloc_w , iloc  ! workspace 
    5657      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zCu_cfl, zCv_cfl, zCw_cfl        ! workspace 
     58      LOGICAL , DIMENSION(jpi,jpj,jpk) ::   llmsk 
    5759      !!---------------------------------------------------------------------- 
    5860      ! 
    5961      IF( ln_timing )   CALL timing_start('dia_cfl') 
    6062      ! 
    61       DO_3D_11_11( 1, jpk ) 
     63      llmsk(   1:Nis1,:,:) = .FALSE.                                              ! exclude halos from the checked region 
     64      llmsk(Nie1: jpi,:,:) = .FALSE. 
     65      llmsk(:,   1:Njs1,:) = .FALSE. 
     66      llmsk(:,Nje1: jpj,:) = .FALSE. 
     67      ! 
     68      DO_3D( 0, 0, 0, 0, 1, jpk )      ! calculate Courant numbers 
    6269         zCu_cfl(ji,jj,jk) = ABS( uu(ji,jj,jk,Kmm) ) * rDt / e1u  (ji,jj)      ! for i-direction 
    6370         zCv_cfl(ji,jj,jk) = ABS( vv(ji,jj,jk,Kmm) ) * rDt / e2v  (ji,jj)      ! for j-direction 
    64          zCw_cfl(ji,jj,jk) = ABS( ww(ji,jj,jk) ) * rDt / e3w(ji,jj,jk,Kmm)   ! for k-direction 
     71         zCw_cfl(ji,jj,jk) = ABS( ww(ji,jj,jk) ) * rDt / e3w(ji,jj,jk,Kmm)     ! for k-direction 
    6572      END_3D 
    6673      ! 
    6774      ! write outputs 
    68       IF( iom_use('cfl_cu') )   CALL iom_put( 'cfl_cu', MAXVAL( zCu_cfl, dim=3 ) ) 
    69       IF( iom_use('cfl_cv') )   CALL iom_put( 'cfl_cv', MAXVAL( zCv_cfl, dim=3 ) ) 
    70       IF( iom_use('cfl_cw') )   CALL iom_put( 'cfl_cw', MAXVAL( zCw_cfl, dim=3 ) ) 
     75      IF( iom_use('cfl_cu') ) THEN 
     76         llmsk(Nis0:Nie0,Njs0:Nje0,:) = umask(Nis0:Nie0,Njs0:Nje0,:) == 1._wp        ! define only the inner domain 
     77         CALL iom_put( 'cfl_cu', MAXVAL( zCu_cfl, mask = llmsk, dim=3 ) ) 
     78      ENDIF 
     79      IF( iom_use('cfl_cv') ) THEN 
     80         llmsk(Nis0:Nie0,Njs0:Nje0,:) = vmask(Nis0:Nie0,Njs0:Nje0,:) == 1._wp        ! define only the inner domain 
     81         CALL iom_put( 'cfl_cv', MAXVAL( zCv_cfl, mask = llmsk, dim=3 ) ) 
     82      ENDIF 
     83      IF( iom_use('cfl_cw') ) THEN 
     84         llmsk(Nis0:Nie0,Njs0:Nje0,:) = wmask(Nis0:Nie0,Njs0:Nje0,:) == 1._wp        ! define only the inner domain 
     85         CALL iom_put( 'cfl_cw', MAXVAL( zCw_cfl, mask = llmsk, dim=3 ) ) 
     86      ENDIF 
    7187 
    7288      !                    ! calculate maximum values and locations 
    73       IF( lk_mpp ) THEN 
    74          CALL mpp_maxloc( 'diacfl', zCu_cfl, umask, zCu_max, iloc_u ) 
    75          CALL mpp_maxloc( 'diacfl', zCv_cfl, vmask, zCv_max, iloc_v ) 
    76          CALL mpp_maxloc( 'diacfl', zCw_cfl, wmask, zCw_max, iloc_w ) 
    77       ELSE 
    78          iloc = MAXLOC( ABS( zcu_cfl(:,:,:) ) ) 
    79          iloc_u(1) = iloc(1) + nimpp - 1 
    80          iloc_u(2) = iloc(2) + njmpp - 1 
    81          iloc_u(3) = iloc(3) 
    82          zCu_max = zCu_cfl(iloc(1),iloc(2),iloc(3)) 
    83          ! 
    84          iloc = MAXLOC( ABS( zcv_cfl(:,:,:) ) ) 
    85          iloc_v(1) = iloc(1) + nimpp - 1 
    86          iloc_v(2) = iloc(2) + njmpp - 1 
    87          iloc_v(3) = iloc(3) 
    88          zCv_max = zCv_cfl(iloc(1),iloc(2),iloc(3)) 
    89          ! 
    90          iloc = MAXLOC( ABS( zcw_cfl(:,:,:) ) ) 
    91          iloc_w(1) = iloc(1) + nimpp - 1 
    92          iloc_w(2) = iloc(2) + njmpp - 1 
    93          iloc_w(3) = iloc(3) 
    94          zCw_max = zCw_cfl(iloc(1),iloc(2),iloc(3)) 
    95       ENDIF 
     89      llmsk(Nis0:Nie0,Njs0:Nje0,:) = umask(Nis0:Nie0,Njs0:Nje0,:) == 1._wp        ! define only the inner domain 
     90      CALL mpp_maxloc( 'diacfl', zCu_cfl, llmsk, zCu_max, iloc_u ) 
     91      llmsk(Nis0:Nie0,Njs0:Nje0,:) = vmask(Nis0:Nie0,Njs0:Nje0,:) == 1._wp        ! define only the inner domain 
     92      CALL mpp_maxloc( 'diacfl', zCv_cfl, llmsk, zCv_max, iloc_v ) 
     93      llmsk(Nis0:Nie0,Njs0:Nje0,:) = wmask(Nis0:Nie0,Njs0:Nje0,:) == 1._wp        ! define only the inner domain 
     94      CALL mpp_maxloc( 'diacfl', zCw_cfl, llmsk, zCw_max, iloc_w ) 
    9695      ! 
    97       !                    ! write out to file 
    98       IF( lwp ) THEN 
     96      IF( lwp ) THEN       ! write out to file 
    9997         WRITE(numcfl,FMT='(2x,i6,3x,a6,4x,f7.4,1x,i4,1x,i4,1x,i4)') kt, 'Max Cu', zCu_max, iloc_u(1), iloc_u(2), iloc_u(3) 
    10098         WRITE(numcfl,FMT='(11x,     a6,4x,f7.4,1x,i4,1x,i4,1x,i4)')     'Max Cv', zCv_max, iloc_v(1), iloc_v(2), iloc_v(3) 
  • NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/OCE/DIA/diadct.F90

    r12489 r13710  
    1111   !!            3.4  ! 09/2011 (C Bricaud) 
    1212   !!---------------------------------------------------------------------- 
    13    !! does not work with agrif 
    1413#if ! defined key_agrif 
     14   !!                        ==>>  CAUTION: does not work with agrif 
    1515   !!---------------------------------------------------------------------- 
    1616   !!   dia_dct      :  Compute the transport through a sec. 
     
    6666   TYPE SECTION 
    6767      CHARACTER(len=60)                            :: name              ! name of the sec 
    68       LOGICAL                                      :: llstrpond         ! true if you want the computation of salt and 
    69                                                                        ! heat transports 
     68      LOGICAL                                      :: llstrpond         ! true if you want the computation of salt and heat transports 
    7069      LOGICAL                                      :: ll_ice_section    ! ice surface and ice volume computation 
    7170      LOGICAL                                      :: ll_date_line      ! = T if the section crosses the date-line 
     
    7473      INTEGER, DIMENSION(nb_point_max)             :: direction         ! vector direction of the point in the section 
    7574      CHARACTER(len=40),DIMENSION(nb_class_max)    :: classname         ! characteristics of the class 
    76       REAL(wp), DIMENSION(nb_class_max)            :: zsigi           ,&! in-situ   density classes    (99 if you don't want) 
    77                                                       zsigp           ,&! potential density classes    (99 if you don't want) 
    78                                                       zsal            ,&! salinity classes   (99 if you don't want) 
    79                                                       ztem            ,&! temperature classes(99 if you don't want) 
    80                                                       zlay              ! level classes      (99 if you don't want) 
     75      REAL(wp), DIMENSION(nb_class_max)            :: zsigi             ! in-situ   density classes    (99 if you don't want) 
     76      REAL(wp), DIMENSION(nb_class_max)            :: zsigp             ! potential density classes    (99 if you don't want) 
     77      REAL(wp), DIMENSION(nb_class_max)            :: zsal              ! salinity classes   (99 if you don't want) 
     78      REAL(wp), DIMENSION(nb_class_max)            :: ztem              ! temperature classes(99 if you don't want) 
     79      REAL(wp), DIMENSION(nb_class_max)            :: zlay              ! level classes      (99 if you don't want) 
    8180      REAL(wp), DIMENSION(nb_type_class,nb_class_max)  :: transport     ! transport output 
    8281      REAL(wp)                                         :: slopeSection  ! slope of the section 
     
    9089   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:)   ::  transports_2d   
    9190 
     91 
     92   !! * Substitutions 
     93#  include "domzgr_substitute.h90" 
    9294   !!---------------------------------------------------------------------- 
    9395   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
     
    9597   !! Software governed by the CeCILL license (see ./LICENSE) 
    9698   !!---------------------------------------------------------------------- 
     99 
    97100CONTAINS 
    98101  
     
    409412              ijloc=ijglo-njmpp+1   !  " 
    410413 
    411               !verify if the point is on the local domain:(1,nlei)*(1,nlej) 
    412               IF( iiloc >= 1 .AND. iiloc <= nlei .AND. & 
    413                   ijloc >= 1 .AND. ijloc <= nlej       )THEN 
     414              !verify if the point is on the local domain:(1,Nie0)*(1,Nje0) 
     415              IF( iiloc >= 1 .AND. iiloc <= Nie0 .AND. & 
     416                  ijloc >= 1 .AND. ijloc <= Nje0       )THEN 
    414417                 iptloc = iptloc + 1                                                 ! count local points 
    415418                 secs(jsec)%listPoint(iptloc) = POINT_SECTION(mi0(iiglo),mj0(ijglo)) ! store local coordinates 
     
    516519  
    517520     !which coordinate shall we verify ? 
    518      IF      ( cdind=='I' )THEN   ; itest=nlei ; iind=1 
    519      ELSE IF ( cdind=='J' )THEN   ; itest=nlej ; iind=2 
     521     IF      ( cdind=='I' )THEN   ; itest=Nie0 ; iind=1 
     522     ELSE IF ( cdind=='J' )THEN   ; itest=Nje0 ; iind=2 
    520523     ELSE    ; CALL ctl_stop("removepoints :Wrong value for cdind")  
    521524     ENDIF 
     
    11191122  !!    |               |                  |       interpolation between ptab(I,J,K) and ptab(I,J,K+1) 
    11201123  !!    |               |                  |       zbis =  
    1121   !!    |               |                  |      [ e3w(I+1,J,K)*ptab(I,J,K) + ( e3w(I,J,K) - e3w(I+1,J,K) ) * ptab(I,J,K-1) ] 
    1122   !!    |               |                  |      /[ e3w(I+1,J,K) + e3w(I,J,K) - e3w(I+1,J,K) ]  
     1124  !!    |               |                  |      [ e3w_n(I+1,J,K,NOW)*ptab(I,J,K) + ( e3w_n(I,J,K,NOW) - e3w_n(I+1,J,K,NOW) ) * ptab(I,J,K-1) ] 
     1125  !!    |               |                  |     /[ e3w_n(I+1,J,K,NOW)             +   e3w_n(I,J,K,NOW) - e3w_n(I+1,J,K,NOW) ]  
    11231126  !!    |               |                  |  
    11241127  !!    |               |                  |    2. Horizontal interpolation: compute value at U/V point 
     
    12131216 
    12141217     ze3t  = e3t(ii2,ij2,kk,Kmm) - e3t(ii1,ij1,kk,Kmm)  
    1215      zwgt1 = ( e3w(ii2,ij2,kk,Kmm) - e3w(ii1,ij1,kk,Kmm) ) / e3w(ii2,ij2,kk,Kmm) 
    1216      zwgt2 = ( e3w(ii1,ij1,kk,Kmm) - e3w(ii2,ij2,kk,Kmm) ) / e3w(ii1,ij1,kk,Kmm) 
     1218     zwgt1 = ( e3w(ii2,ij2,kk,Kmm) - e3w(ii1,ij1,kk,Kmm) )   & 
     1219        &    / e3w(ii2,ij2,kk,Kmm) 
     1220     zwgt2 = ( e3w(ii1,ij1,kk,Kmm) - e3w(ii2,ij2,kk,Kmm) )   & 
     1221        &    / e3w(ii1,ij1,kk,Kmm) 
    12171222 
    12181223     IF(kk .NE. 1)THEN 
  • NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/OCE/DIA/diahsb.F90

    r12489 r13710  
    5050   REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   tmask_ini 
    5151 
     52   !! * Substitutions 
     53#  include "domzgr_substitute.h90" 
    5254   !!---------------------------------------------------------------------- 
    5355   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
     
    156158      ! 
    157159      DO jk = 1, jpkm1           ! volume variation (calculated with scale factors) 
    158          zwrk(:,:,jk) = surf(:,:)*e3t(:,:,jk,Kmm)*tmask(:,:,jk) - surf_ini(:,:)*e3t_ini(:,:,jk)*tmask_ini(:,:,jk) 
     160         zwrk(:,:,jk) =   surf    (:,:) * e3t    (:,:,jk,Kmm)*tmask    (:,:,jk)   & 
     161            &           - surf_ini(:,:) * e3t_ini(:,:,jk    )*tmask_ini(:,:,jk) 
    159162      END DO 
    160163      zdiff_v2 = glob_sum_full( 'diahsb', zwrk(:,:,:) )     ! glob_sum_full needed as tmask and tmask_ini could be different 
    161164      DO jk = 1, jpkm1           ! heat content variation 
    162          zwrk(:,:,jk) = ( surf(:,:)*e3t(:,:,jk,Kmm)*ts(:,:,jk,jp_tem,Kmm) - surf_ini(:,:)*hc_loc_ini(:,:,jk) ) 
     165         zwrk(:,:,jk) = ( surf    (:,:) * e3t(:,:,jk,Kmm)*ts(:,:,jk,jp_tem,Kmm)   & 
     166            &           - surf_ini(:,:) *         hc_loc_ini(:,:,jk) ) 
    163167      END DO 
    164168      zdiff_hc = glob_sum_full( 'diahsb', zwrk(:,:,:) ) 
    165169      DO jk = 1, jpkm1           ! salt content variation 
    166          zwrk(:,:,jk) = ( surf(:,:)*e3t(:,:,jk,Kmm)*ts(:,:,jk,jp_sal,Kmm) - surf_ini(:,:)*sc_loc_ini(:,:,jk) ) 
     170         zwrk(:,:,jk) = ( surf    (:,:) * e3t(:,:,jk,Kmm)*ts(:,:,jk,jp_sal,Kmm)   & 
     171            &           - surf_ini(:,:) *         sc_loc_ini(:,:,jk) ) 
    167172      END DO 
    168173      zdiff_sc = glob_sum_full( 'diahsb', zwrk(:,:,:) ) 
     
    269274               CALL iom_get( numror, 'frc_wn_s', frc_wn_s, ldxios = lrxios ) 
    270275            ENDIF 
    271             CALL iom_get( numror, jpdom_autoglo, 'surf_ini'  , surf_ini  , ldxios = lrxios ) ! ice sheet coupling 
    272             CALL iom_get( numror, jpdom_autoglo, 'ssh_ini'   , ssh_ini   , ldxios = lrxios ) 
    273             CALL iom_get( numror, jpdom_autoglo, 'e3t_ini'   , e3t_ini   , ldxios = lrxios ) 
    274             CALL iom_get( numror, jpdom_autoglo, 'tmask_ini' , tmask_ini , ldxios = lrxios ) 
    275             CALL iom_get( numror, jpdom_autoglo, 'hc_loc_ini', hc_loc_ini, ldxios = lrxios ) 
    276             CALL iom_get( numror, jpdom_autoglo, 'sc_loc_ini', sc_loc_ini, ldxios = lrxios ) 
     276            CALL iom_get( numror, jpdom_auto, 'surf_ini'  , surf_ini  , ldxios = lrxios ) ! ice sheet coupling 
     277            CALL iom_get( numror, jpdom_auto, 'ssh_ini'   , ssh_ini   , ldxios = lrxios ) 
     278            CALL iom_get( numror, jpdom_auto, 'e3t_ini'   , e3t_ini   , ldxios = lrxios ) 
     279            CALL iom_get( numror, jpdom_auto, 'tmask_ini' , tmask_ini , ldxios = lrxios ) 
     280            CALL iom_get( numror, jpdom_auto, 'hc_loc_ini', hc_loc_ini, ldxios = lrxios ) 
     281            CALL iom_get( numror, jpdom_auto, 'sc_loc_ini', sc_loc_ini, ldxios = lrxios ) 
    277282            IF( ln_linssh ) THEN 
    278                CALL iom_get( numror, jpdom_autoglo, 'ssh_hc_loc_ini', ssh_hc_loc_ini, ldxios = lrxios ) 
    279                CALL iom_get( numror, jpdom_autoglo, 'ssh_sc_loc_ini', ssh_sc_loc_ini, ldxios = lrxios ) 
     283               CALL iom_get( numror, jpdom_auto, 'ssh_hc_loc_ini', ssh_hc_loc_ini, ldxios = lrxios ) 
     284               CALL iom_get( numror, jpdom_auto, 'ssh_sc_loc_ini', ssh_sc_loc_ini, ldxios = lrxios ) 
    280285            ENDIF 
    281286         ELSE 
  • NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/OCE/DIA/diahth.F90

    r12489 r13710  
    4242   !! * Substitutions 
    4343#  include "do_loop_substitute.h90" 
     44#  include "domzgr_substitute.h90" 
    4445   !!---------------------------------------------------------------------- 
    4546   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
     
    129130            zdepinv(:,:) = 0._wp   
    130131            zmaxdzT(:,:) = 0._wp   
    131             DO_2D_11_11 
     132            DO_2D( 1, 1, 1, 1 ) 
    132133               zztmp = gdepw(ji,jj,mbkt(ji,jj)+1,Kmm)  
    133134               hth     (ji,jj) = zztmp 
     
    138139            END_2D 
    139140            IF( nla10 > 1 ) THEN  
    140                DO_2D_11_11 
     141               DO_2D( 1, 1, 1, 1 ) 
    141142                  zztmp = gdepw(ji,jj,mbkt(ji,jj)+1,Kmm)  
    142143                  zrho0_3(ji,jj) = zztmp 
     
    147148            ! Preliminary computation 
    148149            ! computation of zdelr = (dr/dT)(T,S,10m)*(-0.2 degC) 
    149             DO_2D_11_11 
     150            DO_2D( 1, 1, 1, 1 ) 
    150151               IF( tmask(ji,jj,nla10) == 1. ) THEN 
    151152                  zu  =  1779.50 + 11.250 * ts(ji,jj,nla10,jp_tem,Kmm) - 3.80   * ts(ji,jj,nla10,jp_sal,Kmm)  & 
     
    169170            ! MLD: rho = rho(1) + zrho1                                     ! 
    170171            ! ------------------------------------------------------------- ! 
    171             DO_3DS_11_11( jpkm1, 2, -1 ) 
     172            DO_3DS( 1, 1, 1, 1, jpkm1, 2, -1 )   ! loop from bottom to 2 
    172173               ! 
    173174               zzdep = gdepw(ji,jj,jk,Kmm) 
     
    206207            ! depth of temperature inversion                                ! 
    207208            ! ------------------------------------------------------------- ! 
    208             DO_3DS_11_11( jpkm1, nlb10, -1 ) 
     209            DO_3DS( 1, 1, 1, 1, jpkm1, nlb10, -1 )   ! loop from bottom to nlb10 
    209210               ! 
    210211               zzdep = gdepw(ji,jj,jk,Kmm) * tmask(ji,jj,1) 
     
    304305      ! --------------------------------------- ! 
    305306      iktem(:,:) = 1 
    306       DO_3D_11_11( 1, jpkm1 ) 
     307      DO_3D( 1, 1, 1, 1, 1, jpkm1 )   ! beware temperature is not always decreasing with depth => loop from top to bottom 
    307308         zztmp = ts(ji,jj,jk,jp_tem,Kmm) 
    308309         IF( zztmp >= ptem )   iktem(ji,jj) = jk 
     
    312313      !  Depth of ptem isotherm         ! 
    313314      ! ------------------------------- ! 
    314       DO_2D_11_11 
     315      DO_2D( 1, 1, 1, 1 ) 
    315316         ! 
    316317         zzdep = gdepw(ji,jj,mbkt(ji,jj)+1,Kmm)       ! depth of the ocean bottom 
     
    350351      ! 
    351352      ilevel(:,:) = 1 
    352       DO_3D_11_11( 2, jpkm1 ) 
     353      DO_3D( 1, 1, 1, 1, 2, jpkm1 ) 
    353354         IF( ( gdept(ji,jj,jk,Kmm) < pdep ) .AND. ( tmask(ji,jj,jk) == 1 ) ) THEN 
    354355             ilevel(ji,jj) = jk 
     
    358359      END_3D 
    359360      ! 
    360       DO_2D_11_11 
     361      DO_2D( 1, 1, 1, 1 ) 
    361362         ik = ilevel(ji,jj) 
    362363         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         phtc(ji,jj)   = phtc(ji,jj)    & 
     365            &           + pt (ji,jj,ik+1) * MIN( e3t(ji,jj,ik+1,Kmm), zthick(ji,jj) ) & 
    364366                                                       * tmask(ji,jj,ik+1) 
    365367      END_2D 
  • NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/OCE/DIA/diamlr.F90

    r12642 r13710  
    44   !! Management of the IOM context for multiple-linear-regression analysis 
    55   !!====================================================================== 
    6    !! History :       !  2019  (S. Mueller) 
     6   !! History :  4.0  !  2019  (S. Mueller)   Original code 
    77   !!---------------------------------------------------------------------- 
    88 
    99   USE par_oce        , ONLY :   wp, jpi, jpj 
    1010   USE phycst         , ONLY :   rpi 
     11   USE dom_oce        , ONLY :   adatrj 
     12   USE tide_mod 
     13   ! 
    1114   USE in_out_manager , ONLY :   lwp, numout, ln_timing 
    1215   USE iom            , ONLY :   iom_put, iom_use, iom_update_file_name 
    13    USE dom_oce        , ONLY :   adatrj 
    1416   USE timing         , ONLY :   timing_start, timing_stop 
    1517#if defined key_iomput 
    1618   USE xios 
    1719#endif 
    18    USE tide_mod 
    1920 
    2021   IMPLICIT NONE 
    2122   PRIVATE 
    2223 
    23    LOGICAL, PUBLIC ::   lk_diamlr = .FALSE. 
     24   LOGICAL, PUBLIC ::   lk_diamlr = .FALSE.   !:         ===>>>   NOT a DOCTOR norm name :  use l_diamlr 
     25   !                                                              lk_  is used only for logical controlled by a CPP key 
    2426 
    2527   PUBLIC ::   dia_mlr_init, dia_mlr_iom_init, dia_mlr 
     
    4244      !! 
    4345      !!---------------------------------------------------------------------- 
    44  
     46      ! 
    4547      lk_diamlr = .TRUE. 
    46  
     48      ! 
    4749      IF(lwp) THEN 
    4850         WRITE(numout, *) 
     
    5052         WRITE(numout, *) '~~~~~~~~~~~~   multiple-linear-regression analysis' 
    5153      END IF 
    52  
     54      ! 
    5355   END SUBROUTINE dia_mlr_init 
     56 
    5457 
    5558   SUBROUTINE dia_mlr_iom_init 
     
    397400   END SUBROUTINE dia_mlr_iom_init 
    398401 
     402 
    399403   SUBROUTINE dia_mlr 
    400404      !!---------------------------------------------------------------------- 
     
    404408      !! 
    405409      !!---------------------------------------------------------------------- 
    406  
    407410      REAL(wp), DIMENSION(jpi,jpj) ::   zadatrj2d 
     411      !!---------------------------------------------------------------------- 
    408412 
    409413      IF( ln_timing )   CALL timing_start('dia_mlr') 
     
    412416      ! (value of adatrj converted to time in units of seconds) 
    413417      ! 
    414       ! A 2-dimensional field of constant value is sent, and subsequently used 
    415       ! directly or transformed to a scalar or a constant 3-dimensional field as 
    416       ! required. 
     418      ! A 2-dimensional field of constant value is sent, and subsequently used directly  
     419      ! or transformed to a scalar or a constant 3-dimensional field as required. 
    417420      zadatrj2d(:,:) = adatrj*86400.0_wp 
    418421      IF ( iom_use('diamlr_time') ) CALL iom_put('diamlr_time', zadatrj2d) 
    419        
     422      ! 
    420423      IF( ln_timing )   CALL timing_stop('dia_mlr') 
    421  
     424      ! 
    422425   END SUBROUTINE dia_mlr 
    423426 
     427   !!====================================================================== 
    424428END MODULE diamlr 
  • NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/OCE/DIA/diaptr.F90

    r12489 r13710  
    3636   END INTERFACE 
    3737 
    38    PUBLIC   ptr_sj         ! call by tra_ldf & tra_adv routines 
    39    PUBLIC   ptr_sjk        !  
    40    PUBLIC   dia_ptr_init   ! call in memogcm 
    4138   PUBLIC   dia_ptr        ! call in step module 
    4239   PUBLIC   dia_ptr_hst    ! called from tra_ldf/tra_adv routines 
    4340 
    44    !                                  !!** namelist  namptr  ** 
    4541   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   hstr_adv, hstr_ldf, hstr_eiv   !: Heat/Salt TRansports(adv, diff, Bolus.) 
    4642   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   hstr_ove, hstr_btr, hstr_vtr   !: heat Salt TRansports(overturn, baro, merional) 
    4743 
    48    LOGICAL , PUBLIC ::   l_diaptr        !: tracers  trend flag (set from namelist in trdini) 
    49    INTEGER, PARAMETER, PUBLIC ::   nptr = 5  ! (glo, atl, pac, ind, ipc) 
     44   LOGICAL, PUBLIC ::   l_diaptr       !: tracers  trend flag 
    5045 
    5146   REAL(wp) ::   rc_sv    = 1.e-6_wp   ! conversion from m3/s to Sverdrup 
     
    5954   REAL(wp), TARGET, ALLOCATABLE, SAVE, DIMENSION(:,:) :: p_fval2d 
    6055 
    61    LOGICAL ::   ll_init = .TRUE.        !: tracers  trend flag (set from namelist in trdini) 
     56   LOGICAL ::   ll_init = .TRUE.        !: tracers  trend flag 
     57    
    6258   !! * Substitutions 
    6359#  include "do_loop_substitute.h90" 
     60#  include "domzgr_substitute.h90" 
    6461   !!---------------------------------------------------------------------- 
    6562   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
     
    8683      ! 
    8784      !overturning calculation 
    88       REAL(wp), DIMENSION(jpj,jpk,nptr) :: sjk, r1_sjk, v_msf  ! i-mean i-k-surface and its inverse 
    89       REAL(wp), DIMENSION(jpj,jpk,nptr) :: zt_jk, zs_jk ! i-mean T and S, j-Stream-Function 
    90  
    91       REAL(wp), DIMENSION(jpi,jpj,jpk,nptr)  :: z4d1, z4d2 
    92       REAL(wp), DIMENSION(jpi,jpj,nptr)      :: z3dtr ! i-mean T and S, j-Stream-Function 
     85      REAL(wp), DIMENSION(:,:,:  ), ALLOCATABLE ::  sjk, r1_sjk, v_msf  ! i-mean i-k-surface and its inverse 
     86      REAL(wp), DIMENSION(:,:,:  ), ALLOCATABLE ::   zt_jk, zs_jk        ! i-mean T and S, j-Stream-Function 
     87 
     88      REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE ::  z4d1, z4d2 
     89      REAL(wp), DIMENSION(:,:,:  ), ALLOCATABLE ::   z3dtr 
    9390      !!---------------------------------------------------------------------- 
    9491      ! 
    9592      IF( ln_timing )   CALL timing_start('dia_ptr') 
    9693 
    97       IF( kt == nit000 .AND. ll_init )   CALL dia_ptr_init 
    98       ! 
    99       IF( .NOT. l_diaptr )   RETURN 
    100  
     94      IF( kt == nit000 .AND. ll_init )   CALL dia_ptr_init   ! -> will define l_diaptr and nbasin 
     95      ! 
     96      IF( .NOT. l_diaptr ) THEN 
     97         IF( ln_timing ) CALL timing_stop('dia_ptr') 
     98         RETURN 
     99      ENDIF 
     100      ! 
     101      ALLOCATE( z3dtr(jpi,jpj,nbasin) ) 
     102      ! 
    101103      IF( PRESENT( pvtr ) ) THEN 
    102104         IF( iom_use( 'zomsf' ) ) THEN    ! effective MSF 
    103             DO jn = 1, nptr                                    ! by sub-basins 
     105            ALLOCATE( z4d1(jpi,jpj,jpk,nbasin) ) 
     106            DO jn = 1, nbasin                                    ! by sub-basins 
    104107               z4d1(1,:,:,jn) =  ptr_sjk( pvtr(:,:,:), btmsk34(:,:,jn) )  ! zonal cumulative effective transport excluding closed seas 
    105108               DO jk = jpkm1, 1, -1  
     
    111114            END DO 
    112115            CALL iom_put( 'zomsf', z4d1 * rc_sv ) 
     116            DEALLOCATE( z4d1 ) 
    113117         ENDIF 
    114118         IF(  iom_use( 'sopstove' ) .OR. iom_use( 'sophtove' ) .OR.   & 
     
    117121            zmask(:,:,:) = 0._wp 
    118122            zts(:,:,:,:) = 0._wp 
    119             DO_3D_10_11( 1, jpkm1 ) 
     123            DO_3D( 1, 0, 1, 1, 1, jpkm1 ) 
    120124               zvfc = e1v(ji,jj) * e3v(ji,jj,jk,Kmm) 
    121125               zmask(ji,jj,jk)      = vmask(ji,jj,jk)      * zvfc 
     
    125129         ENDIF 
    126130         IF( iom_use( 'sopstove' ) .OR. iom_use( 'sophtove' ) ) THEN 
    127             DO jn = 1, nptr 
     131            DO jn = 1, nbasin 
     132               ALLOCATE( sjk(jpj,jpk,nbasin), r1_sjk(jpj,jpk,nbasin), v_msf(jpj,jpk,nbasin),   & 
     133                  &                          zt_jk(jpj,jpk,nbasin), zs_jk(jpj,jpk,nbasin) ) 
    128134               sjk(:,:,jn) = ptr_sjk( zmask(:,:,:), btmsk(:,:,jn) ) 
    129135               r1_sjk(:,:,jn) = 0._wp 
     
    135141               hstr_ove(:,jp_tem,jn) = SUM( v_msf(:,:,jn)*zt_jk(:,:,jn), 2 ) 
    136142               hstr_ove(:,jp_sal,jn) = SUM( v_msf(:,:,jn)*zs_jk(:,:,jn), 2 ) 
     143               DEALLOCATE( sjk, r1_sjk, v_msf, zt_jk, zs_jk ) 
    137144               ! 
    138145            ENDDO 
    139             DO jn = 1, nptr 
     146            DO jn = 1, nbasin 
    140147               z3dtr(1,:,jn) = hstr_ove(:,jp_tem,jn) * rc_pwatt  !  (conversion in PW) 
    141148               DO ji = 1, jpi 
     
    144151            ENDDO 
    145152            CALL iom_put( 'sophtove', z3dtr ) 
    146             DO jn = 1, nptr 
     153            DO jn = 1, nbasin 
    147154               z3dtr(1,:,jn) = hstr_ove(:,jp_sal,jn) * rc_ggram !  (conversion in Gg) 
    148155               DO ji = 1, jpi 
     
    155162         IF( iom_use( 'sopstbtr' ) .OR. iom_use( 'sophtbtr' ) ) THEN 
    156163            ! Calculate barotropic heat and salt transport here  
    157             DO jn = 1, nptr 
     164            DO jn = 1, nbasin 
     165               ALLOCATE( sjk(jpj,1,nbasin), r1_sjk(jpj,1,nbasin) ) 
    158166               sjk(:,1,jn) = ptr_sj( zmask(:,:,:), btmsk(:,:,jn) ) 
    159167               r1_sjk(:,1,jn) = 0._wp 
     
    165173               hstr_btr(:,jp_tem,jn) = zvsum(:) * ztsum(:) * r1_sjk(:,1,jn) 
    166174               hstr_btr(:,jp_sal,jn) = zvsum(:) * zssum(:) * r1_sjk(:,1,jn) 
     175               DEALLOCATE( sjk, r1_sjk ) 
    167176               ! 
    168177            ENDDO 
    169             DO jn = 1, nptr 
     178            DO jn = 1, nbasin 
    170179               z3dtr(1,:,jn) = hstr_btr(:,jp_tem,jn) * rc_pwatt  !  (conversion in PW) 
    171180               DO ji = 1, jpi 
     
    174183            ENDDO 
    175184            CALL iom_put( 'sophtbtr', z3dtr ) 
    176             DO jn = 1, nptr 
     185            DO jn = 1, nbasin 
    177186               z3dtr(1,:,jn) = hstr_btr(:,jp_sal,jn) * rc_ggram !  (conversion in Gg) 
    178187               DO ji = 1, jpi 
     
    188197         zts(:,:,:,:) = 0._wp 
    189198         IF( iom_use( 'zotem' ) .OR. iom_use( 'zosal' ) .OR. iom_use( 'zosrf' )  ) THEN    ! i-mean i-k-surface  
    190             DO_3D_11_11( 1, jpkm1 ) 
     199            ALLOCATE( z4d1(jpi,jpj,jpk,nbasin), z4d2(jpi,jpj,jpk,nbasin) ) 
     200            DO_3D( 1, 1, 1, 1, 1, jpkm1 ) 
    191201               zsfc = e1t(ji,jj) * e3t(ji,jj,jk,Kmm) 
    192202               zmask(ji,jj,jk)      = tmask(ji,jj,jk)      * zsfc 
     
    195205            END_3D 
    196206            ! 
    197             DO jn = 1, nptr 
     207            DO jn = 1, nbasin 
    198208               zmask(1,:,:) = ptr_sjk( zmask(:,:,:), btmsk(:,:,jn) ) 
     209               DO ji = 1, jpi 
     210                  zmask(ji,:,:) = zmask(1,:,:) 
     211               ENDDO 
    199212               z4d1(:,:,:,jn) = zmask(:,:,:) 
    200213            ENDDO 
    201214            CALL iom_put( 'zosrf', z4d1 ) 
    202215            ! 
    203             DO jn = 1, nptr 
     216            DO jn = 1, nbasin 
    204217               z4d2(1,:,:,jn) = ptr_sjk( zts(:,:,:,jp_tem), btmsk(:,:,jn) ) & 
    205218                  &            / MAX( z4d1(1,:,:,jn), 10.e-15 ) 
     
    210223            CALL iom_put( 'zotem', z4d2 ) 
    211224            ! 
    212             DO jn = 1, nptr 
     225            DO jn = 1, nbasin 
    213226               z4d2(1,:,:,jn) = ptr_sjk( zts(:,:,:,jp_sal), btmsk(:,:,jn) ) & 
    214227                  &            / MAX( z4d1(1,:,:,jn), 10.e-15 ) 
     
    218231            ENDDO 
    219232            CALL iom_put( 'zosal', z4d2 ) 
     233            DEALLOCATE( z4d1, z4d2 ) 
    220234            ! 
    221235         ENDIF 
     
    224238         IF( iom_use( 'sophtadv' ) .OR. iom_use( 'sopstadv' ) ) THEN   
    225239            !  
    226             DO jn = 1, nptr 
     240            DO jn = 1, nbasin 
    227241               z3dtr(1,:,jn) = hstr_adv(:,jp_tem,jn) * rc_pwatt  !  (conversion in PW) 
    228242               DO ji = 1, jpi 
     
    231245            ENDDO 
    232246            CALL iom_put( 'sophtadv', z3dtr ) 
    233             DO jn = 1, nptr 
     247            DO jn = 1, nbasin 
    234248               z3dtr(1,:,jn) = hstr_adv(:,jp_sal,jn) * rc_ggram !  (conversion in Gg) 
    235249               DO ji = 1, jpi 
     
    242256         IF( iom_use( 'sophtldf' ) .OR. iom_use( 'sopstldf' ) ) THEN   
    243257            !  
    244             DO jn = 1, nptr 
     258            DO jn = 1, nbasin 
    245259               z3dtr(1,:,jn) = hstr_ldf(:,jp_tem,jn) * rc_pwatt  !  (conversion in PW) 
    246260               DO ji = 1, jpi 
     
    249263            ENDDO 
    250264            CALL iom_put( 'sophtldf', z3dtr ) 
    251             DO jn = 1, nptr 
     265            DO jn = 1, nbasin 
    252266               z3dtr(1,:,jn) = hstr_ldf(:,jp_sal,jn) * rc_ggram !  (conversion in Gg) 
    253267               DO ji = 1, jpi 
     
    260274         IF( iom_use( 'sophteiv' ) .OR. iom_use( 'sopsteiv' ) ) THEN   
    261275            !  
    262             DO jn = 1, nptr 
     276            DO jn = 1, nbasin 
    263277               z3dtr(1,:,jn) = hstr_eiv(:,jp_tem,jn) * rc_pwatt  !  (conversion in PW) 
    264278               DO ji = 1, jpi 
     
    267281            ENDDO 
    268282            CALL iom_put( 'sophteiv', z3dtr ) 
    269             DO jn = 1, nptr 
     283            DO jn = 1, nbasin 
    270284               z3dtr(1,:,jn) = hstr_eiv(:,jp_sal,jn) * rc_ggram !  (conversion in Gg) 
    271285               DO ji = 1, jpi 
     
    278292         IF( iom_use( 'sopstvtr' ) .OR. iom_use( 'sophtvtr' ) ) THEN 
    279293            zts(:,:,:,:) = 0._wp 
    280             DO_3D_10_11( 1, jpkm1 ) 
     294            DO_3D( 1, 0, 1, 1, 1, jpkm1 ) 
    281295               zvfc = e1v(ji,jj) * e3v(ji,jj,jk,Kmm) 
    282296               zts(ji,jj,jk,jp_tem) = (ts(ji,jj,jk,jp_tem,Kmm)+ts(ji,jj+1,jk,jp_tem,Kmm)) * 0.5 * zvfc  !Tracers averaged onto V grid 
     
    285299             CALL dia_ptr_hst( jp_tem, 'vtr', zts(:,:,:,jp_tem) ) 
    286300             CALL dia_ptr_hst( jp_sal, 'vtr', zts(:,:,:,jp_sal) ) 
    287              DO jn = 1, nptr 
     301             DO jn = 1, nbasin 
    288302                z3dtr(1,:,jn) = hstr_vtr(:,jp_tem,jn) * rc_pwatt  !  (conversion in PW) 
    289303                DO ji = 1, jpi 
     
    292306             ENDDO 
    293307             CALL iom_put( 'sophtvtr', z3dtr ) 
    294              DO jn = 1, nptr 
     308             DO jn = 1, nbasin 
    295309               z3dtr(1,:,jn) = hstr_vtr(:,jp_sal,jn) * rc_ggram !  (conversion in Gg) 
    296310               DO ji = 1, jpi 
     
    309323      ENDIF 
    310324      ! 
     325      DEALLOCATE( z3dtr ) 
     326      ! 
    311327      IF( ln_timing )   CALL timing_stop('dia_ptr') 
    312328      ! 
     
    318334      !!                  ***  ROUTINE dia_ptr_init  *** 
    319335      !!                    
    320       !! ** Purpose :   Initialization, namelist read 
     336      !! ** Purpose :   Initialization 
    321337      !!---------------------------------------------------------------------- 
    322338      INTEGER ::  inum, jn           ! local integers 
     
    324340      REAL(wp), DIMENSION(jpi,jpj) :: zmsk 
    325341      !!---------------------------------------------------------------------- 
    326  
    327       l_diaptr = .FALSE. 
    328       IF(   iom_use( 'zomsf'    ) .OR. iom_use( 'zotem'    ) .OR. iom_use( 'zosal'    ) .OR.  & 
    329          &  iom_use( 'zosrf'    ) .OR. iom_use( 'sopstove' ) .OR. iom_use( 'sophtove' ) .OR.  & 
    330          &  iom_use( 'sopstbtr' ) .OR. iom_use( 'sophtbtr' ) .OR. iom_use( 'sophtadv' ) .OR.  & 
    331          &  iom_use( 'sopstadv' ) .OR. iom_use( 'sophtldf' ) .OR. iom_use( 'sopstldf' ) .OR.  &  
    332          &  iom_use( 'sophteiv' ) .OR. iom_use( 'sopsteiv' ) .OR. iom_use( 'sopstvtr' ) .OR.  & 
    333          &  iom_use( 'sophtvtr' ) .OR. iom_use( 'uocetr_vsum_cumul' ) )  l_diaptr  = .TRUE. 
    334  
     342       
     343      ! l_diaptr is defined with iom_use 
     344      !   --> dia_ptr_init must be done after the call to iom_init 
     345      !   --> cannot be .TRUE. without cpp key: key_iom -->  nbasin define by iom_init is initialized 
     346      l_diaptr = iom_use( 'zomsf'    ) .OR. iom_use( 'zotem'    ) .OR. iom_use( 'zosal'    ) .OR.  & 
     347         &       iom_use( 'zosrf'    ) .OR. iom_use( 'sopstove' ) .OR. iom_use( 'sophtove' ) .OR.  & 
     348         &       iom_use( 'sopstbtr' ) .OR. iom_use( 'sophtbtr' ) .OR. iom_use( 'sophtadv' ) .OR.  & 
     349         &       iom_use( 'sopstadv' ) .OR. iom_use( 'sophtldf' ) .OR. iom_use( 'sopstldf' ) .OR.  &  
     350         &       iom_use( 'sophteiv' ) .OR. iom_use( 'sopsteiv' ) .OR. iom_use( 'sopstvtr' ) .OR.  & 
     351         &       iom_use( 'sophtvtr' ) .OR. iom_use( 'uocetr_vsum_cumul' )  
    335352  
    336353      IF(lwp) THEN                     ! Control print 
     
    338355         WRITE(numout,*) 'dia_ptr_init : poleward transport and msf initialization' 
    339356         WRITE(numout,*) '~~~~~~~~~~~~' 
    340          WRITE(numout,*) '   Namelist namptr : set ptr parameters' 
    341357         WRITE(numout,*) '      Poleward heat & salt transport (T) or not (F)      l_diaptr  = ', l_diaptr 
    342358      ENDIF 
     
    345361         ! 
    346362         IF( dia_ptr_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'dia_ptr_init : unable to allocate arrays' ) 
    347  
     363         ! 
    348364         rc_pwatt = rc_pwatt * rho0_rcp          ! conversion from K.s-1 to PetaWatt 
    349365         rc_ggram = rc_ggram * rho0              ! conversion from m3/s to Gg/s 
     
    352368 
    353369         btmsk(:,:,1) = tmask_i(:,:)                  
    354          CALL iom_open( 'subbasins', inum,  ldstop = .FALSE.  ) 
    355          CALL iom_get( inum, jpdom_data, 'atlmsk', btmsk(:,:,2) )   ! Atlantic basin 
    356          CALL iom_get( inum, jpdom_data, 'pacmsk', btmsk(:,:,3) )   ! Pacific  basin 
    357          CALL iom_get( inum, jpdom_data, 'indmsk', btmsk(:,:,4) )   ! Indian   basin 
    358          CALL iom_close( inum ) 
    359          btmsk(:,:,5) = MAX ( btmsk(:,:,3), btmsk(:,:,4) )          ! Indo-Pacific basin 
    360          DO jn = 2, nptr 
    361             btmsk(:,:,jn) = btmsk(:,:,jn) * tmask_i(:,:)               ! interior domain only 
     370         IF( nbasin == 5 ) THEN   ! nbasin has been initialized in iom_init to define the axis "basin" 
     371            CALL iom_open( 'subbasins', inum ) 
     372            CALL iom_get( inum, jpdom_global, 'atlmsk', btmsk(:,:,2) )   ! Atlantic basin 
     373            CALL iom_get( inum, jpdom_global, 'pacmsk', btmsk(:,:,3) )   ! Pacific  basin 
     374            CALL iom_get( inum, jpdom_global, 'indmsk', btmsk(:,:,4) )   ! Indian   basin 
     375            CALL iom_close( inum ) 
     376            btmsk(:,:,5) = MAX ( btmsk(:,:,3), btmsk(:,:,4) )            ! Indo-Pacific basin 
     377         ENDIF 
     378         DO jn = 2, nbasin 
     379            btmsk(:,:,jn) = btmsk(:,:,jn) * tmask_i(:,:)                 ! interior domain only 
    362380         END DO 
    363381         ! JD : modification so that overturning streamfunction is available in Atlantic at 34S to compare with observations 
     
    368386         END WHERE 
    369387         btmsk34(:,:,1) = btmsk(:,:,1)                  
    370          DO jn = 2, nptr 
    371             btmsk34(:,:,jn) = btmsk(:,:,jn) * zmsk(:,:)               ! interior domain only 
     388         DO jn = 2, nbasin 
     389            btmsk34(:,:,jn) = btmsk(:,:,jn) * zmsk(:,:)                  ! interior domain only 
    372390         ENDDO 
    373391 
     
    403421      IF( cptr == 'adv' ) THEN 
    404422         IF( ktra == jp_tem )  THEN 
    405              DO jn = 1, nptr 
     423             DO jn = 1, nbasin 
    406424                hstr_adv(:,jp_tem,jn) = ptr_sj( pvflx(:,:,:), btmsk(:,:,jn) ) 
    407425             ENDDO 
    408426         ENDIF 
    409427         IF( ktra == jp_sal )  THEN 
    410              DO jn = 1, nptr 
     428             DO jn = 1, nbasin 
    411429                hstr_adv(:,jp_sal,jn) = ptr_sj( pvflx(:,:,:), btmsk(:,:,jn) ) 
    412430             ENDDO 
     
    416434      IF( cptr == 'ldf' ) THEN 
    417435         IF( ktra == jp_tem )  THEN 
    418              DO jn = 1, nptr 
     436             DO jn = 1, nbasin 
    419437                hstr_ldf(:,jp_tem,jn) = ptr_sj( pvflx(:,:,:), btmsk(:,:,jn) ) 
    420438             ENDDO 
    421439         ENDIF 
    422440         IF( ktra == jp_sal )  THEN 
    423              DO jn = 1, nptr 
     441             DO jn = 1, nbasin 
    424442                hstr_ldf(:,jp_sal,jn) = ptr_sj( pvflx(:,:,:), btmsk(:,:,jn) ) 
    425443             ENDDO 
     
    429447      IF( cptr == 'eiv' ) THEN 
    430448         IF( ktra == jp_tem )  THEN 
    431              DO jn = 1, nptr 
     449             DO jn = 1, nbasin 
    432450                hstr_eiv(:,jp_tem,jn) = ptr_sj( pvflx(:,:,:), btmsk(:,:,jn) ) 
    433451             ENDDO 
    434452         ENDIF 
    435453         IF( ktra == jp_sal )  THEN 
    436              DO jn = 1, nptr 
     454             DO jn = 1, nbasin 
    437455                hstr_eiv(:,jp_sal,jn) = ptr_sj( pvflx(:,:,:), btmsk(:,:,jn) ) 
    438456             ENDDO 
     
    442460      IF( cptr == 'vtr' ) THEN 
    443461         IF( ktra == jp_tem )  THEN 
    444              DO jn = 1, nptr 
     462             DO jn = 1, nbasin 
    445463                hstr_vtr(:,jp_tem,jn) = ptr_sj( pvflx(:,:,:), btmsk(:,:,jn) ) 
    446464             ENDDO 
    447465         ENDIF 
    448466         IF( ktra == jp_sal )  THEN 
    449              DO jn = 1, nptr 
     467             DO jn = 1, nbasin 
    450468                hstr_vtr(:,jp_sal,jn) = ptr_sj( pvflx(:,:,:), btmsk(:,:,jn) ) 
    451469             ENDDO 
     
    465483      ierr(:) = 0 
    466484      ! 
     485      ! nbasin has been initialized in iom_init to define the axis "basin" 
     486      ! 
    467487      IF( .NOT. ALLOCATED( btmsk ) ) THEN 
    468          ALLOCATE( btmsk(jpi,jpj,nptr)    , btmsk34(jpi,jpj,nptr),   & 
    469             &      hstr_adv(jpj,jpts,nptr), hstr_eiv(jpj,jpts,nptr), & 
    470             &      hstr_ove(jpj,jpts,nptr), hstr_btr(jpj,jpts,nptr), & 
    471             &      hstr_ldf(jpj,jpts,nptr), hstr_vtr(jpj,jpts,nptr), STAT=ierr(1)  ) 
     488         ALLOCATE( btmsk(jpi,jpj,nbasin)    , btmsk34(jpi,jpj,nbasin),   & 
     489            &      hstr_adv(jpj,jpts,nbasin), hstr_eiv(jpj,jpts,nbasin), & 
     490            &      hstr_ove(jpj,jpts,nbasin), hstr_btr(jpj,jpts,nbasin), & 
     491            &      hstr_ldf(jpj,jpts,nbasin), hstr_vtr(jpj,jpts,nbasin), STAT=ierr(1)  ) 
    472492            ! 
    473493         ALLOCATE( p_fval1d(jpj), p_fval2d(jpj,jpk), Stat=ierr(2)) 
     
    503523      ijpj = jpj 
    504524      p_fval(:) = 0._wp 
    505       DO_3D_00_00( 1, jpkm1 ) 
     525      DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 
    506526         p_fval(jj) = p_fval(jj) + pvflx(ji,jj,jk) * pmsk(ji,jj) * tmask_i(ji,jj) 
    507527      END_3D 
     
    536556      ijpj = jpj 
    537557      p_fval(:) = 0._wp 
    538       DO_2D_00_00 
     558      DO_2D( 0, 0, 0, 0 ) 
    539559         p_fval(jj) = p_fval(jj) + pvflx(ji,jj) * pmsk(ji,jj) * tmask_i(ji,jj) 
    540560      END_2D 
     
    565585      p_fval(:,:) = 0._wp 
    566586      DO jc = 1, jpnj ! looping over all processors in j axis 
    567          DO_2D_00_00 
     587         DO_2D( 0, 0, 0, 0 ) 
    568588            p_fval(ji,jj) = p_fval(ji,jj-1) + pva(ji,jj) * tmask_i(ji,jj) 
    569589         END_2D 
    570          CALL lbc_lnk( 'diaptr', p_fval, 'U', -1. ) 
     590         CALL lbc_lnk( 'diaptr', p_fval, 'U', -1.0_wp ) 
    571591      END DO 
    572592      !  
     
    604624      p_fval(:,:) = 0._wp 
    605625      ! 
    606       DO_3D_00_00( 1, jpkm1 ) 
     626      DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 
    607627         p_fval(jj,jk) = p_fval(jj,jk) + pta(ji,jj,jk) * pmsk(ji,jj) * tmask_i(ji,jj) 
    608628      END_3D 
  • NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/OCE/DIA/diawri.F90

    r12649 r13710  
    8585   !! * Substitutions 
    8686#  include "do_loop_substitute.h90" 
     87#  include "domzgr_substitute.h90" 
    8788   !!---------------------------------------------------------------------- 
    8889   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
     
    117118      INTEGER ::   ji, jj, jk       ! dummy loop indices 
    118119      INTEGER ::   ikbot            ! local integer 
     120      REAL(wp)::   ze3 
    119121      REAL(wp)::   zztmp , zztmpx   ! local scalar 
    120122      REAL(wp)::   zztmp2, zztmpy   !   -      - 
     
    136138      CALL iom_put("e3v_0", e3v_0(:,:,:) ) 
    137139      ! 
    138       CALL iom_put( "e3t" , e3t(:,:,:,Kmm) ) 
    139       CALL iom_put( "e3u" , e3u(:,:,:,Kmm) ) 
    140       CALL iom_put( "e3v" , e3v(:,:,:,Kmm) ) 
    141       CALL iom_put( "e3w" , e3w(:,:,:,Kmm) ) 
    142       IF( iom_use("e3tdef") )   & 
    143          CALL iom_put( "e3tdef"  , ( ( e3t(:,:,:,Kmm) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100 * tmask(:,:,:) ) ** 2 ) 
    144  
    145       IF( ll_wd ) THEN 
    146          CALL iom_put( "ssh" , (ssh(:,:,Kmm)+ssh_ref)*tmask(:,:,1) )   ! sea surface height (brought back to the reference used for wetting and drying) 
     140      IF ( iom_use("e3t") .OR. iom_use("e3tdef") ) THEN  ! time-varying e3t 
     141         DO jk = 1, jpk 
     142            z3d(:,:,jk) =  e3t(:,:,jk,Kmm) 
     143         END DO 
     144         CALL iom_put( "e3t"     ,     z3d(:,:,:) ) 
     145         CALL iom_put( "e3tdef"  , ( ( z3d(:,:,:) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100 * tmask(:,:,:) ) ** 2 )  
     146      ENDIF  
     147      IF ( iom_use("e3u") ) THEN                         ! time-varying e3u 
     148         DO jk = 1, jpk 
     149            z3d(:,:,jk) =  e3u(:,:,jk,Kmm) 
     150         END DO  
     151         CALL iom_put( "e3u" , z3d(:,:,:) ) 
     152      ENDIF 
     153      IF ( iom_use("e3v") ) THEN                         ! time-varying e3v 
     154         DO jk = 1, jpk 
     155            z3d(:,:,jk) =  e3v(:,:,jk,Kmm) 
     156         END DO  
     157         CALL iom_put( "e3v" , z3d(:,:,:) ) 
     158      ENDIF 
     159      IF ( iom_use("e3w") ) THEN                         ! time-varying e3w 
     160         DO jk = 1, jpk 
     161            z3d(:,:,jk) =  e3w(:,:,jk,Kmm) 
     162         END DO  
     163         CALL iom_put( "e3w" , z3d(:,:,:) ) 
     164      ENDIF 
     165 
     166      IF( ll_wd ) THEN                                   ! sea surface height (brought back to the reference used for wetting and drying) 
     167         CALL iom_put( "ssh" , (ssh(:,:,Kmm)+ssh_ref)*tmask(:,:,1) ) 
    147168      ELSE 
    148169         CALL iom_put( "ssh" , ssh(:,:,Kmm) )              ! sea surface height 
     
    155176      CALL iom_put(  "sst", ts(:,:,1,jp_tem,Kmm) )    ! surface temperature 
    156177      IF ( iom_use("sbt") ) THEN 
    157          DO_2D_11_11 
     178         DO_2D( 0, 0, 0, 0 ) 
    158179            ikbot = mbkt(ji,jj) 
    159180            z2d(ji,jj) = ts(ji,jj,ikbot,jp_tem,Kmm) 
     
    165186      CALL iom_put(  "sss", ts(:,:,1,jp_sal,Kmm) )    ! surface salinity 
    166187      IF ( iom_use("sbs") ) THEN 
    167          DO_2D_11_11 
     188         DO_2D( 0, 0, 0, 0 ) 
    168189            ikbot = mbkt(ji,jj) 
    169190            z2d(ji,jj) = ts(ji,jj,ikbot,jp_sal,Kmm) 
     
    172193      ENDIF 
    173194 
     195#if ! defined key_qco 
     196      CALL iom_put( "rhop", rhop(:,:,:) )          ! 3D potential density (sigma0) 
     197#endif 
     198 
    174199      IF ( iom_use("taubot") ) THEN                ! bottom stress 
    175200         zztmp = rho0 * 0.25 
    176201         z2d(:,:) = 0._wp 
    177          DO_2D_00_00 
     202         DO_2D( 0, 0, 0, 0 ) 
    178203            zztmp2 = (  ( rCdU_bot(ji+1,jj)+rCdU_bot(ji  ,jj) ) * uu(ji  ,jj,mbku(ji  ,jj),Kmm)  )**2   & 
    179204               &   + (  ( rCdU_bot(ji  ,jj)+rCdU_bot(ji-1,jj) ) * uu(ji-1,jj,mbku(ji-1,jj),Kmm)  )**2   & 
     
    183208            ! 
    184209         END_2D 
    185          CALL lbc_lnk( 'diawri', z2d, 'T', 1. ) 
    186210         CALL iom_put( "taubot", z2d )            
    187211      ENDIF 
     
    190214      CALL iom_put(  "ssu", uu(:,:,1,Kmm) )            ! surface i-current 
    191215      IF ( iom_use("sbu") ) THEN 
    192          DO_2D_11_11 
     216         DO_2D( 0, 0, 0, 0 ) 
    193217            ikbot = mbku(ji,jj) 
    194218            z2d(ji,jj) = uu(ji,jj,ikbot,Kmm) 
     
    200224      CALL iom_put(  "ssv", vv(:,:,1,Kmm) )            ! surface j-current 
    201225      IF ( iom_use("sbv") ) THEN 
    202          DO_2D_11_11 
     226         DO_2D( 0, 0, 0, 0 ) 
    203227            ikbot = mbkv(ji,jj) 
    204228            z2d(ji,jj) = vv(ji,jj,ikbot,Kmm) 
     
    208232 
    209233      IF( ln_zad_Aimp ) ww = ww + wi               ! Recombine explicit and implicit parts of vertical velocity for diagnostic output 
    210       ! 
    211234      CALL iom_put( "woce", ww )                   ! vertical velocity 
     235 
    212236      IF( iom_use('w_masstr') .OR. iom_use('w_masstr2') ) THEN   ! vertical mass transport & its square value 
    213237         ! Caution: in the VVL case, it only correponds to the baroclinic mass transport. 
     
    229253      IF( iom_use('logavs') )   CALL iom_put( "logavs", LOG( MAX( 1.e-20_wp, avs(:,:,:) ) ) ) 
    230254 
     255      IF ( iom_use("socegrad") .OR. iom_use("socegrad2") ) THEN 
     256         z3d(:,:,jpk) = 0. 
     257         DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 
     258            zztmp  = ts(ji,jj,jk,jp_sal,Kmm) 
     259            zztmpx = (ts(ji+1,jj,jk,jp_sal,Kmm) - zztmp) * r1_e1u(ji,jj) + (zztmp - ts(ji-1,jj  ,jk,jp_sal,Kmm)) * r1_e1u(ji-1,jj) 
     260            zztmpy = (ts(ji,jj+1,jk,jp_sal,Kmm) - zztmp) * r1_e2v(ji,jj) + (zztmp - ts(ji  ,jj-1,jk,jp_sal,Kmm)) * r1_e2v(ji,jj-1) 
     261            z3d(ji,jj,jk) = 0.25 * ( zztmpx * zztmpx + zztmpy * zztmpy )   & 
     262               &                 * umask(ji,jj,jk) * umask(ji-1,jj,jk) * vmask(ji,jj,jk) * umask(ji,jj-1,jk) 
     263         END_3D 
     264         CALL iom_put( "socegrad2",  z3d )          ! square of module of sal gradient 
     265         DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 
     266            z3d(ji,jj,jk) = SQRT( z3d(ji,jj,jk) ) 
     267         END_3D 
     268         CALL iom_put( "socegrad" ,  z3d )          ! module of sal gradient 
     269      ENDIF 
     270          
    231271      IF ( iom_use("sstgrad") .OR. iom_use("sstgrad2") ) THEN 
    232          DO_2D_00_00 
     272         DO_2D( 0, 0, 0, 0 )                                 ! sst gradient 
    233273            zztmp  = ts(ji,jj,1,jp_tem,Kmm) 
    234274            zztmpx = ( ts(ji+1,jj,1,jp_tem,Kmm) - zztmp ) * r1_e1u(ji,jj) + ( zztmp - ts(ji-1,jj  ,1,jp_tem,Kmm) ) * r1_e1u(ji-1,jj) 
     
    237277               &              * umask(ji,jj,1) * umask(ji-1,jj,1) * vmask(ji,jj,1) * umask(ji,jj-1,1) 
    238278         END_2D 
    239          CALL lbc_lnk( 'diawri', z2d, 'T', 1. ) 
    240279         CALL iom_put( "sstgrad2",  z2d )          ! square of module of sst gradient 
    241          z2d(:,:) = SQRT( z2d(:,:) ) 
     280         DO_2D( 0, 0, 0, 0 ) 
     281            z2d(ji,jj) = SQRT( z2d(ji,jj) ) 
     282         END_2D 
    242283         CALL iom_put( "sstgrad" ,  z2d )          ! module of sst gradient 
    243284      ENDIF 
     
    246287      IF( iom_use("heatc") ) THEN 
    247288         z2d(:,:)  = 0._wp  
    248          DO_3D_11_11( 1, jpkm1 ) 
     289         DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 
    249290            z2d(ji,jj) = z2d(ji,jj) + e3t(ji,jj,jk,Kmm) * ts(ji,jj,jk,jp_tem,Kmm) * tmask(ji,jj,jk) 
    250291         END_3D 
     
    254295      IF( iom_use("saltc") ) THEN 
    255296         z2d(:,:)  = 0._wp  
    256          DO_3D_11_11( 1, jpkm1 ) 
     297         DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 
    257298            z2d(ji,jj) = z2d(ji,jj) + e3t(ji,jj,jk,Kmm) * ts(ji,jj,jk,jp_sal,Kmm) * tmask(ji,jj,jk) 
    258299         END_3D 
     
    260301      ENDIF 
    261302      ! 
    262       IF ( iom_use("eken") ) THEN 
     303      IF( iom_use("salt2c") ) THEN 
     304         z2d(:,:)  = 0._wp  
     305         DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 
     306            z2d(ji,jj) = z2d(ji,jj) + e3t(ji,jj,jk,Kmm) * ts(ji,jj,jk,jp_sal,Kmm) * ts(ji,jj,jk,jp_sal,Kmm) * tmask(ji,jj,jk) 
     307         END_3D 
     308         CALL iom_put( "salt2c", rho0 * z2d )          ! vertically integrated salt content (PSU*kg/m2) 
     309      ENDIF 
     310      ! 
     311      IF ( iom_use("ke") .OR. iom_use("ke_int") ) THEN 
    263312         z3d(:,:,jpk) = 0._wp  
    264          DO_3D_00_00( 1, jpkm1 ) 
    265             zztmp  = 0.25_wp * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm) 
    266             z3d(ji,jj,jk) = zztmp * (  uu(ji-1,jj,jk,Kmm)**2 * e2u(ji-1,jj) * e3u(ji-1,jj,jk,Kmm)   & 
    267                &                     + uu(ji  ,jj,jk,Kmm)**2 * e2u(ji  ,jj) * e3u(ji  ,jj,jk,Kmm)   & 
    268                &                     + vv(ji,jj-1,jk,Kmm)**2 * e1v(ji,jj-1) * e3v(ji,jj-1,jk,Kmm)   & 
    269                &                     + vv(ji,jj  ,jk,Kmm)**2 * e1v(ji,jj  ) * e3v(ji,jj  ,jk,Kmm)   ) 
    270          END_3D 
    271          CALL lbc_lnk( 'diawri', z3d, 'T', 1. ) 
    272          CALL iom_put( "eken", z3d )                 ! kinetic energy 
     313         DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 
     314            zztmpx = 0.5 * ( uu(ji-1,jj  ,jk,Kmm) + uu(ji,jj,jk,Kmm) ) 
     315            zztmpy = 0.5 * ( vv(ji  ,jj-1,jk,Kmm) + vv(ji,jj,jk,Kmm) ) 
     316            z3d(ji,jj,jk) = 0.5 * ( zztmpx*zztmpx + zztmpy*zztmpy ) 
     317         END_3D 
     318         CALL iom_put( "ke", z3d )                 ! kinetic energy 
     319 
     320         z2d(:,:)  = 0._wp  
     321         DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 
     322            z2d(ji,jj) = z2d(ji,jj) + e3t(ji,jj,jk,Kmm) * z3d(ji,jj,jk) * e1e2t(ji,jj) * tmask(ji,jj,jk) 
     323         END_3D 
     324         CALL iom_put( "ke_int", z2d )   ! vertically integrated kinetic energy 
    273325      ENDIF 
    274326      ! 
    275327      CALL iom_put( "hdiv", hdiv )                  ! Horizontal divergence 
     328 
     329      IF ( iom_use("relvor") .OR. iom_use("absvor") .OR. iom_use("potvor") ) THEN 
     330          
     331         z3d(:,:,jpk) = 0._wp  
     332         DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 
     333            z3d(ji,jj,jk) = (   e2v(ji+1,jj  ) * vv(ji+1,jj  ,jk,Kmm) - e2v(ji,jj) * vv(ji,jj,jk,Kmm)    & 
     334               &              - e1u(ji  ,jj+1) * uu(ji  ,jj+1,jk,Kmm) + e1u(ji,jj) * uu(ji,jj,jk,Kmm)  ) * r1_e1e2f(ji,jj) 
     335         END_3D 
     336         CALL iom_put( "relvor", z3d )                  ! relative vorticity 
     337 
     338         DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 
     339            z3d(ji,jj,jk) = ff_f(ji,jj) + z3d(ji,jj,jk)  
     340         END_3D 
     341         CALL iom_put( "absvor", z3d )                  ! absolute vorticity 
     342 
     343         DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 
     344            ze3  = (  e3t(ji,jj+1,jk,Kmm)*tmask(ji,jj+1,jk) + e3t(ji+1,jj+1,jk,Kmm)*tmask(ji+1,jj+1,jk)   & 
     345               &    + e3t(ji,jj  ,jk,Kmm)*tmask(ji,jj  ,jk) + e3t(ji+1,jj  ,jk,Kmm)*tmask(ji+1,jj  ,jk)  ) 
     346            IF( ze3 /= 0._wp ) THEN   ;   ze3 = 4._wp / ze3 
     347            ELSE                      ;   ze3 = 0._wp 
     348            ENDIF 
     349            z3d(ji,jj,jk) = ze3 * z3d(ji,jj,jk)  
     350         END_3D 
     351         CALL iom_put( "potvor", z3d )                  ! potential vorticity 
     352 
     353      ENDIF 
    276354      ! 
    277355      IF( iom_use("u_masstr") .OR. iom_use("u_masstr_vint") .OR. iom_use("u_heattr") .OR. iom_use("u_salttr") ) THEN 
     
    288366      IF( iom_use("u_heattr") ) THEN 
    289367         z2d(:,:) = 0._wp  
    290          DO_3D_00_00( 1, jpkm1 ) 
     368         DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 
    291369            z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( ts(ji,jj,jk,jp_tem,Kmm) + ts(ji+1,jj,jk,jp_tem,Kmm) ) 
    292370         END_3D 
    293          CALL lbc_lnk( 'diawri', z2d, 'U', -1. ) 
    294371         CALL iom_put( "u_heattr", 0.5*rcp * z2d )    ! heat transport in i-direction 
    295372      ENDIF 
     
    297374      IF( iom_use("u_salttr") ) THEN 
    298375         z2d(:,:) = 0.e0  
    299          DO_3D_00_00( 1, jpkm1 ) 
     376         DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 
    300377            z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( ts(ji,jj,jk,jp_sal,Kmm) + ts(ji+1,jj,jk,jp_sal,Kmm) ) 
    301378         END_3D 
    302          CALL lbc_lnk( 'diawri', z2d, 'U', -1. ) 
    303379         CALL iom_put( "u_salttr", 0.5 * z2d )        ! heat transport in i-direction 
    304380      ENDIF 
     
    315391      IF( iom_use("v_heattr") ) THEN 
    316392         z2d(:,:) = 0.e0  
    317          DO_3D_00_00( 1, jpkm1 ) 
     393         DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 
    318394            z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( ts(ji,jj,jk,jp_tem,Kmm) + ts(ji,jj+1,jk,jp_tem,Kmm) ) 
    319395         END_3D 
    320          CALL lbc_lnk( 'diawri', z2d, 'V', -1. ) 
    321396         CALL iom_put( "v_heattr", 0.5*rcp * z2d )    !  heat transport in j-direction 
    322397      ENDIF 
     
    324399      IF( iom_use("v_salttr") ) THEN 
    325400         z2d(:,:) = 0._wp  
    326          DO_3D_00_00( 1, jpkm1 ) 
     401         DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 
    327402            z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( ts(ji,jj,jk,jp_sal,Kmm) + ts(ji,jj+1,jk,jp_sal,Kmm) ) 
    328403         END_3D 
    329          CALL lbc_lnk( 'diawri', z2d, 'V', -1. ) 
    330404         CALL iom_put( "v_salttr", 0.5 * z2d )        !  heat transport in j-direction 
    331405      ENDIF 
     
    333407      IF( iom_use("tosmint") ) THEN 
    334408         z2d(:,:) = 0._wp 
    335          DO_3D_00_00( 1, jpkm1 ) 
     409         DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 
    336410            z2d(ji,jj) = z2d(ji,jj) + e3t(ji,jj,jk,Kmm) *  ts(ji,jj,jk,jp_tem,Kmm) 
    337411         END_3D 
    338          CALL lbc_lnk( 'diawri', z2d, 'T', -1. ) 
    339412         CALL iom_put( "tosmint", rho0 * z2d )        ! Vertical integral of temperature 
    340413      ENDIF 
    341414      IF( iom_use("somint") ) THEN 
    342415         z2d(:,:)=0._wp 
    343          DO_3D_00_00( 1, jpkm1 ) 
     416         DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 
    344417            z2d(ji,jj) = z2d(ji,jj) + e3t(ji,jj,jk,Kmm) * ts(ji,jj,jk,jp_sal,Kmm) 
    345418         END_3D 
    346          CALL lbc_lnk( 'diawri', z2d, 'T', -1. ) 
    347419         CALL iom_put( "somint", rho0 * z2d )         ! Vertical integral of salinity 
    348420      ENDIF 
     
    415487      ! 
    416488      REAL(wp), DIMENSION(jpi,jpj)   :: zw2d       ! 2D workspace 
    417       REAL(wp), DIMENSION(jpi,jpj,jpk) :: zw3d       ! 3D workspace 
     489      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zw3d, ze3t, zgdept       ! 3D workspace 
    418490      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: zw3d_abl   ! ABL 3D workspace 
    419491      !!---------------------------------------------------------------------- 
     
    447519 
    448520      ! Define indices of the horizontal output zoom and vertical limit storage 
    449       iimi = 1      ;      iima = jpi 
    450       ijmi = 1      ;      ijma = jpj 
     521      iimi = Nis0   ;   iima = Nie0 
     522      ijmi = Njs0   ;   ijma = Nje0 
    451523      ipk = jpk 
    452524      IF(ln_abl) ipka = jpkam1 
     
    455527      it = kt 
    456528      itmod = kt - nit000 + 1 
     529 
     530      ! store e3t for subsitute 
     531      DO jk = 1, jpk 
     532         ze3t  (:,:,jk) =  e3t  (:,:,jk,Kmm) 
     533         zgdept(:,:,jk) =  gdept(:,:,jk,Kmm) 
     534      END DO 
    457535 
    458536 
     
    569647         DEALLOCATE(zw3d_abl) 
    570648         ENDIF 
     649         ! 
    571650 
    572651         ! Declare all the output fields as NETCDF variables 
     
    578657            &          jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout ) 
    579658         IF(  .NOT.ln_linssh  ) THEN 
    580             CALL histdef( nid_T, "vovvle3t", "Level thickness"                    , "m"      ,&  ! e3t(:,:,:,Kmm) 
     659            CALL histdef( nid_T, "vovvle3t", "Level thickness"                    , "m"      ,&  ! e3t n 
    581660            &             jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout ) 
    582             CALL histdef( nid_T, "vovvldep", "T point depth"                      , "m"      ,&  ! e3t(:,:,:,Kmm) 
     661            CALL histdef( nid_T, "vovvldep", "T point depth"                      , "m"      ,&  ! e3t n 
    583662            &             jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout ) 
    584             CALL histdef( nid_T, "vovvldef", "Squared level deformation"          , "%^2"    ,&  ! e3t(:,:,:,Kmm) 
     663            CALL histdef( nid_T, "vovvldef", "Squared level deformation"          , "%^2"    ,&  ! e3t n 
    585664            &             jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout ) 
    586665         ENDIF 
     
    766845 
    767846      IF( .NOT.ln_linssh ) THEN 
    768          CALL histwrite( nid_T, "votemper", it, ts(:,:,:,jp_tem,Kmm) * e3t(:,:,:,Kmm) , ndim_T , ndex_T  )   ! heat content 
    769          CALL histwrite( nid_T, "vosaline", it, ts(:,:,:,jp_sal,Kmm) * e3t(:,:,:,Kmm) , ndim_T , ndex_T  )   ! salt content 
    770          CALL histwrite( nid_T, "sosstsst", it, ts(:,:,1,jp_tem,Kmm) * e3t(:,:,1,Kmm) , ndim_hT, ndex_hT )   ! sea surface heat content 
    771          CALL histwrite( nid_T, "sosaline", it, ts(:,:,1,jp_sal,Kmm) * e3t(:,:,1,Kmm) , ndim_hT, ndex_hT )   ! sea surface salinity content 
     847         CALL histwrite( nid_T, "votemper", it, ts(:,:,:,jp_tem,Kmm) * ze3t(:,:,:) , ndim_T , ndex_T  )   ! heat content 
     848         CALL histwrite( nid_T, "vosaline", it, ts(:,:,:,jp_sal,Kmm) * ze3t(:,:,:) , ndim_T , ndex_T  )   ! salt content 
     849         CALL histwrite( nid_T, "sosstsst", it, ts(:,:,1,jp_tem,Kmm) * ze3t(:,:,1) , ndim_hT, ndex_hT )   ! sea surface heat content 
     850         CALL histwrite( nid_T, "sosaline", it, ts(:,:,1,jp_sal,Kmm) * ze3t(:,:,1) , ndim_hT, ndex_hT )   ! sea surface salinity content 
    772851      ELSE 
    773852         CALL histwrite( nid_T, "votemper", it, ts(:,:,:,jp_tem,Kmm) , ndim_T , ndex_T  )   ! temperature 
     
    777856      ENDIF 
    778857      IF( .NOT.ln_linssh ) THEN 
    779          zw3d(:,:,:) = ( ( e3t(:,:,:,Kmm) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100 * tmask(:,:,:) ) ** 2 
    780          CALL histwrite( nid_T, "vovvle3t", it, e3t (:,:,:,Kmm) , ndim_T , ndex_T  )   ! level thickness 
    781          CALL histwrite( nid_T, "vovvldep", it, gdept(:,:,:,Kmm) , ndim_T , ndex_T  )   ! t-point depth 
     858         zw3d(:,:,:) = ( ( ze3t(:,:,:) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100 * tmask(:,:,:) ) ** 2 
     859         CALL histwrite( nid_T, "vovvle3t", it, ze3t (:,:,:)    , ndim_T , ndex_T  )   ! level thickness 
     860         CALL histwrite( nid_T, "vovvldep", it, zgdept , ndim_T , ndex_T  )   ! t-point depth  
    782861         CALL histwrite( nid_T, "vovvldef", it, zw3d             , ndim_T , ndex_T  )   ! level thickness deformation 
    783862      ENDIF 
     
    918997      !! 
    919998      INTEGER :: inum, jk 
     999      REAL(wp), DIMENSION(jpi,jpj,jpk) :: ze3t, zgdept      ! 3D workspace !!st patch to use substitution 
    9201000      !!---------------------------------------------------------------------- 
    9211001      !  
    922       IF(lwp) WRITE(numout,*) 
    923       IF(lwp) WRITE(numout,*) 'dia_wri_state : single instantaneous ocean state' 
    924       IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~   and forcing fields file created ' 
    925       IF(lwp) WRITE(numout,*) '                and named :', cdfile_name, '...nc' 
     1002      IF(lwp) THEN 
     1003         WRITE(numout,*) 
     1004         WRITE(numout,*) 'dia_wri_state : single instantaneous ocean state' 
     1005         WRITE(numout,*) '~~~~~~~~~~~~~   and forcing fields file created ' 
     1006         WRITE(numout,*) '                and named :', cdfile_name, '...nc' 
     1007      ENDIF  
     1008      ! 
     1009      DO jk = 1, jpk 
     1010         ze3t(:,:,jk) =  e3t(:,:,jk,Kmm) 
     1011         zgdept(:,:,jk) =  gdept(:,:,jk,Kmm) 
     1012      END DO 
    9261013      ! 
    9271014      CALL iom_open( TRIM(cdfile_name), inum, ldwrt = .TRUE. ) 
     
    9381025      ENDIF 
    9391026      CALL iom_rstput( 0, 0, inum, 'risfdep', risfdep            )    ! now k-velocity 
    940       CALL iom_rstput( 0, 0, inum, 'ht'     , ht                 )    ! now water column height 
     1027      CALL iom_rstput( 0, 0, inum, 'ht'     , ht(:,:)            )    ! now water column height 
    9411028      ! 
    9421029      IF ( ln_isf ) THEN 
     
    9751062      CALL iom_rstput( 0, 0, inum, 'sometauy', vtau              )    ! j-wind stress 
    9761063      IF(  .NOT.ln_linssh  ) THEN              
    977          CALL iom_rstput( 0, 0, inum, 'vovvldep', gdept(:,:,:,Kmm)        )    !  T-cell depth  
    978          CALL iom_rstput( 0, 0, inum, 'vovvle3t', e3t(:,:,:,Kmm)          )    !  T-cell thickness   
     1064         CALL iom_rstput( 0, 0, inum, 'vovvldep', zgdept        )    !  T-cell depth  
     1065         CALL iom_rstput( 0, 0, inum, 'vovvle3t', ze3t          )    !  T-cell thickness   
    9791066      END IF 
    9801067      IF( ln_wave .AND. ln_sdw ) THEN 
     
    9981085         CALL iom_close( inum ) 
    9991086      ENDIF 
     1087      ! 
    10001088#endif 
    1001  
    10021089   END SUBROUTINE dia_wri_state 
    10031090 
Note: See TracChangeset for help on using the changeset viewer.