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 7367 – NEMO

Changeset 7367


Ignore:
Timestamp:
2016-11-29T11:52:31+01:00 (8 years ago)
Author:
deazer
Message:

Contains merged code for CO5 reference.

Location:
branches/UKMO/CO5_package_branch/NEMOGCM/NEMO/OPA_SRC
Files:
14 added
94 edited

Legend:

Unmodified
Added
Removed
  • branches/UKMO/CO5_package_branch/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90

    r7363 r7367  
    1010   !!   NEMO     3.3  ! 2010-05  (D. Lea)  Update to work with NEMO v3.2 
    1111   !!             -   ! 2010-05  (D. Lea)  add calc_month_len routine based on day_init  
     12   !!            3.4  ! 2012-10  (A. Weaver and K. Mogensen) Fix for direct initialization 
    1213   !!---------------------------------------------------------------------- 
    1314 
     
    2021   !!   dyn_asm_inc  : Apply the dynamic (u and v) increments 
    2122   !!   ssh_asm_inc  : Apply the SSH increment 
     23   !!   seaice_asm_inc  : Apply the seaice increment 
    2224   !!---------------------------------------------------------------------- 
    2325   USE wrk_nemo         ! Memory Allocation 
     
    2527   USE dom_oce          ! Ocean space and time domain 
    2628   USE oce              ! Dynamics and active tracers defined in memory 
    27    USE divcur           ! Horizontal divergence and relative vorticity 
    2829   USE ldfdyn_oce       ! ocean dynamics: lateral physics 
    2930   USE eosbn2           ! Equation of state - in situ and potential density 
     
    3334   USE c1d              ! 1D initialization 
    3435   USE in_out_manager   ! I/O manager 
    35    USE lib_mpp           ! MPP library 
     36   USE lib_mpp          ! MPP library 
     37#if defined key_lim3 || defined key_lim2 || defined key_cice 
     38#if defined key_lim3  
     39   USE ice_3, ONLY : &               ! LIM Ice model variables  () 
     40      & frld, pfrld, hicif, hsnif, phicif 
     41   USE sbc_oce, ONLY : & 
     42      & fr_i                         ! ice fraction 
     43#endif 
     44#if defined key_lim2 
     45   USE ice_2, ONLY : &               ! LIM Ice model variables 
     46      & frld, pfrld, hicif, hsnif, phicif 
     47   USE sbc_oce, ONLY : & 
     48      & fr_i                         ! ice fraction 
     49#endif 
     50#if defined key_cice 
     51   USE sbc_oce, ONLY : & 
     52      & fr_i                          ! ice fraction 
     53   USE sbc_ice, ONLY : &             ! CICE Ice model variables 
     54      & naicet, ndaice_da, nfresh_da, nfsalt_da, nTf  
     55   USE ice_constants, only: Lfresh, rhoi,rhos        ! for updating ice and snow enthalphy 
     56!   USE ice_therm_itd, only: hfrazilmin               ! thickness at new ice points 
     57   USE ice_domain_size, only: ncat,ntilyr,ntslyr 
     58#endif 
     59   USE phycst, ONLY : &              ! Physical Ice variables 
     60      & soce, sice, rhoic, rhosn, rday 
     61#endif 
     62   USE sbc_oce          ! Surface boundary condition variables. 
     63    
     64   USE eosbn2, only: tfreez  
     65   
     66   USE zdfmxl, ONLY :  &   
     67   &  hmld_tref,       &    
     68#if defined key_karaml 
     69   &  hmld_kara,       & 
     70#endif    
     71   &  hmld,            &  
     72   &  hmlp 
     73  
     74#if defined key_bdy  
     75   USE bdy_oce, ONLY: bdytmask   
     76#endif   
     77   USE histcom 
    3678 
    3779   IMPLICIT NONE 
     
    4385   PUBLIC   dyn_asm_inc    !: Apply the dynamic (u and v) increments 
    4486   PUBLIC   ssh_asm_inc    !: Apply the SSH increment 
     87   PUBLIC   seaice_asm_inc !: Apply the seaice increment 
    4588 
    4689#if defined key_asminc 
     
    5699   LOGICAL, PUBLIC :: ln_dyninc = .FALSE. !: No dynamics (u and v) assimilation increments 
    57100   LOGICAL, PUBLIC :: ln_sshinc = .FALSE. !: No sea surface height assimilation increment 
     101   LOGICAL, PUBLIC :: ln_seaiceinc = .FALSE. !: No sea ice concentration increment 
    58102   LOGICAL, PUBLIC :: ln_salfix = .FALSE. !: Apply minimum salinity check 
     103   LOGICAL, PUBLIC :: ln_temnofreeze = .FALSE. !: Don't allow the temperature to drop below freezing 
    59104   INTEGER, PUBLIC :: nn_divdmp = 0       !: Apply divergence damping filter nn_divdmp times 
    60105 
     
    78123 
    79124   REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   ssh_bkg, ssh_bkginc   ! Background sea surface height and its increment 
    80  
     125   REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   seaice_bkginc         ! Increment to the background sea ice conc 
     126    
     127  INTEGER :: mld_choice = 4     !: choice of mld criteria to use                   
     128                             !: 1) turbocline depth   
     129                             !: 2) surface to 0.001 kg/m^3 change   
     130                             !: 3) Kara MLD   
     131                             !: 4) Temperature criteria.                                 
     132                                  
    81133   !! * Substitutions 
    82134#  include "domzgr_substitute.h90" 
     
    122174      REAL(wp) :: zdate_bkg    ! Date in background state file for DI 
    123175      REAL(wp) :: zdate_inc    ! Time axis in increments file 
    124  
     176       
     177      REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: &  
     178          &       t_bkginc_2d  !file for reading in 2D   
     179                               !temperature increments  
     180      REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: &  
     181          &       z_mld     !Mixed layer depth  
     182           
    125183      REAL(wp), POINTER, DIMENSION(:,:) :: hdiv 
     184             
    126185      !! 
    127186      NAMELIST/nam_asminc/ ln_bkgwri, ln_trjwri,                           & 
     
    130189         &                 nitbkg, nitdin, nitiaustr, nitiaufin, niaufn,   & 
    131190         &                 nittrjfrq, ln_salfix, salfixmin,                & 
    132          &                 nn_divdmp 
     191         &                 nn_divdmp, mld_choice 
    133192      !!---------------------------------------------------------------------- 
    134193 
     
    143202      ln_dyninc = .FALSE. 
    144203      ln_sshinc = .FALSE. 
     204      ln_seaiceinc = .FALSE. 
    145205      ln_asmdin = .FALSE. 
    146206      ln_asmiau = .TRUE. 
    147207      ln_salfix = .FALSE. 
     208      ln_temnofreeze = .FALSE. 
    148209      salfixmin = -9999 
    149210      nitbkg    = 0 
     
    156217      REWIND ( numnam ) 
    157218      READ   ( numnam, nam_asminc ) 
    158  
     219      
     220      ! Set the data time for diagnostics to the end of the IAU period   
     221      ! and multiply by the timestep to get seconds from start of run  
     222      data_time = rdt * nitiaufin  
     223        
     224      IF( ln_sco .AND. (ln_sshinc .OR. ln_seaiceinc .OR. ln_asmdin &  
     225         &              .OR. ln_dyninc ) )THEN  
     226         CALL ctl_warn("Only SST assimilation currently supported in "//&  
     227         &              "s-coordinates")  
     228         ln_sshinc = .FALSE.  
     229         ln_seaiceinc = .FALSE.  
     230         ln_asmdin = .FALSE.  
     231         ln_dyninc = .FALSE.  
     232      ENDIF  
     233       
    159234      ! Control print 
    160235      IF(lwp) THEN 
     
    169244         WRITE(numout,*) '      Logical switch for applying SSH increments               ln_sshinc = ', ln_sshinc 
    170245         WRITE(numout,*) '      Logical switch for Direct Initialization (DI)            ln_asmdin = ', ln_asmdin 
     246         WRITE(numout,*) '      Logical switch for applying sea ice increments        ln_seaiceinc = ', ln_seaiceinc 
    171247         WRITE(numout,*) '      Logical switch for Incremental Analysis Updating (IAU)   ln_asmiau = ', ln_asmiau 
    172248         WRITE(numout,*) '      Timestep of background in [0,nitend-nit000-1]            nitbkg    = ', nitbkg 
     
    235311 
    236312      IF (      ( ( .NOT. ln_asmdin ).AND.( .NOT. ln_asmiau ) ) & 
    237            .AND.( ( ln_trainc ).OR.( ln_dyninc ).OR.( ln_sshinc ) ) ) & 
    238          & CALL ctl_stop( ' One or more of ln_trainc, ln_dyninc and ln_sshinc is set to .true.', & 
     313           .AND.( ( ln_trainc ).OR.( ln_dyninc ).OR.( ln_sshinc ) .OR. ( ln_seaiceinc) )) & 
     314         & CALL ctl_stop( ' One or more of ln_trainc, ln_dyninc, ln_sshinc and ln_seaiceinc is set to .true.', & 
    239315         &                ' but ln_asmdin and ln_asmiau are both set to .false. :', & 
    240316         &                ' Inconsistent options') 
     
    248324         &                ' Type IAU weighting function is invalid') 
    249325 
    250       IF ( ( .NOT. ln_trainc ).AND.( .NOT. ln_dyninc ).AND.( .NOT. ln_sshinc ) & 
     326      IF ( ( .NOT. ln_trainc ).AND.( .NOT. ln_dyninc ).AND.( .NOT. ln_sshinc ).AND.( .NOT. ln_seaiceinc ) & 
    251327         &                     )  & 
    252          & CALL ctl_warn( ' ln_trainc, ln_dyninc and ln_sshinc are set to .false. :', & 
     328         & CALL ctl_warn( ' ln_trainc, ln_dyninc, ln_sshinc and ln_seaiceinc are set to .false. :', & 
    253329         &                ' The assimilation increments are not applied') 
    254330 
     
    353429      ALLOCATE( v_bkginc(jpi,jpj,jpk) ) 
    354430      ALLOCATE( ssh_bkginc(jpi,jpj)   ) 
     431      ALLOCATE( seaice_bkginc(jpi,jpj)) 
    355432#if defined key_asminc 
    356433      ALLOCATE( ssh_iau(jpi,jpj)      ) 
     
    361438      v_bkginc(:,:,:) = 0.0 
    362439      ssh_bkginc(:,:) = 0.0 
     440      seaice_bkginc(:,:) = 0.0 
    363441#if defined key_asminc 
    364442      ssh_iau(:,:)    = 0.0 
    365443#endif 
    366       IF ( ( ln_trainc ).OR.( ln_dyninc ).OR.( ln_sshinc ) ) THEN 
     444      IF ( ( ln_trainc ).OR.( ln_dyninc ).OR.( ln_sshinc ).OR.( ln_seaiceinc ) ) THEN 
    367445 
    368446         !-------------------------------------------------------------------- 
     
    397475 
    398476         IF ( ln_trainc ) THEN    
    399             CALL iom_get( inum, jpdom_autoglo, 'bckint', t_bkginc, 1 ) 
    400             CALL iom_get( inum, jpdom_autoglo, 'bckins', s_bkginc, 1 ) 
    401             ! Apply the masks 
    402             t_bkginc(:,:,:) = t_bkginc(:,:,:) * tmask(:,:,:) 
    403             s_bkginc(:,:,:) = s_bkginc(:,:,:) * tmask(:,:,:) 
    404             ! Set missing increments to 0.0 rather than 1e+20 
    405             ! to allow for differences in masks 
    406             WHERE( ABS( t_bkginc(:,:,:) ) > 1.0e+10 ) t_bkginc(:,:,:) = 0.0 
    407             WHERE( ABS( s_bkginc(:,:,:) ) > 1.0e+10 ) s_bkginc(:,:,:) = 0.0 
     477             
     478            IF (ln_sco) THEN  
     479                 
     480               ALLOCATE(z_mld(jpi,jpj))  
     481               SELECT CASE(mld_choice)  
     482               CASE(1)  
     483                  z_mld = hmld  
     484               CASE(2)  
     485                  z_mld = hmlp  
     486               CASE(3)  
     487#if defined key_karaml 
     488                  z_mld = hmld_kara 
     489#endif 
     490                  CALL ctl_stop("Kara mixed layer not defined in current version of NEMO")  ! JW: Safty feature, should be removed 
     491                                                                                            ! once the kara mixed layer is availible  
     492               CASE(4)  
     493                  z_mld = hmld_tref  
     494               END SELECT        
     495                       
     496               ALLOCATE( t_bkginc_2d(jpi,jpj) )  
     497               CALL iom_get( inum, jpdom_autoglo, 'bckinsurft', t_bkginc_2d, 1)  
     498#if defined key_bdy                 
     499               DO jk = 1,jpkm1  
     500                  WHERE( z_mld(:,:) > fsdepw(:,:,jk) )  
     501                     t_bkginc(:,:,jk) = t_bkginc_2d(:,:) * bdytmask(:,:)  
     502                  ELSEWHERE  
     503                     t_bkginc(:,:,jk) = 0.  
     504                  ENDWHERE  
     505               ENDDO  
     506#else  
     507               t_bkginc(:,:,:) = 0.  
     508#endif                 
     509               s_bkginc(:,:,:) = 0.  
     510                 
     511               !DEALLOCATE temporary arrays  
     512               DEALLOCATE(z_mld, t_bkginc_2d)  
     513             
     514            ELSE  
     515                
     516               CALL iom_get( inum, jpdom_autoglo, 'bckint', t_bkginc, 1 ) 
     517               CALL iom_get( inum, jpdom_autoglo, 'bckins', s_bkginc, 1 ) 
     518               ! Apply the masks 
     519               t_bkginc(:,:,:) = t_bkginc(:,:,:) * tmask(:,:,:) 
     520               s_bkginc(:,:,:) = s_bkginc(:,:,:) * tmask(:,:,:) 
     521               ! Set missing increments to 0.0 rather than 1e+20 
     522               ! to allow for differences in masks 
     523               WHERE( ABS( t_bkginc(:,:,:) ) > 1.0e+10 ) t_bkginc(:,:,:) = 0.0 
     524               WHERE( ABS( s_bkginc(:,:,:) ) > 1.0e+10 ) s_bkginc(:,:,:) = 0.0 
     525          
     526            ENDIF 
     527          
    408528         ENDIF 
    409529 
     
    429549         ENDIF 
    430550 
     551         IF ( ln_seaiceinc ) THEN 
     552            CALL iom_get( inum, jpdom_autoglo, 'bckinseaice', seaice_bkginc, 1 ) 
     553            ! Apply the masks 
     554            seaice_bkginc(:,:) = seaice_bkginc(:,:) * tmask(:,:,1) 
     555            ! Set missing increments to 0.0 rather than 1e+20 
     556            ! to allow for differences in masks 
     557            WHERE( ABS( seaice_bkginc(:,:) ) > 1.0e+10 ) seaice_bkginc(:,:) = 0.0 
     558         ENDIF 
     559 
    431560         CALL iom_close( inum ) 
    432561  
     
    437566      !----------------------------------------------------------------------- 
    438567 
    439  
    440568      IF ( ln_dyninc .AND. nn_divdmp > 0 ) THEN 
    441569 
    442       CALL wrk_alloc(jpi,jpj,hdiv)  
    443  
    444        DO  jt = 1, nn_divdmp 
    445  
    446            DO jk = 1, jpkm1 
    447  
    448                   hdiv(:,:) = 0._wp 
    449  
    450             DO jj = 2, jpjm1 
    451                DO ji = fs_2, fs_jpim1   ! vector opt. 
    452                   hdiv(ji,jj) =   & 
    453                      (  e2u(ji  ,jj)*fse3u(ji  ,jj,jk) * u_bkginc(ji  ,jj,jk)       & 
    454                       - e2u(ji-1,jj)*fse3u(ji-1,jj,jk) * u_bkginc(ji-1,jj,jk)       & 
    455                       + e1v(ji,jj  )*fse3v(ji,jj  ,jk) * v_bkginc(ji,jj  ,jk)       & 
    456                       - e1v(ji,jj-1)*fse3v(ji,jj-1,jk) * v_bkginc(ji,jj-1,jk)  )    & 
    457                       / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 
     570         CALL wrk_alloc(jpi,jpj,hdiv)  
     571 
     572         DO  jt = 1, nn_divdmp 
     573 
     574            DO jk = 1, jpkm1 
     575 
     576               hdiv(:,:) = 0._wp 
     577 
     578               DO jj = 2, jpjm1 
     579                  DO ji = fs_2, fs_jpim1   ! vector opt. 
     580                     hdiv(ji,jj) =   & 
     581                        (  e2u(ji  ,jj  ) * fse3u(ji  ,jj  ,jk) * u_bkginc(ji  ,jj  ,jk)     & 
     582                         - e2u(ji-1,jj  ) * fse3u(ji-1,jj  ,jk) * u_bkginc(ji-1,jj  ,jk)     & 
     583                         + e1v(ji  ,jj  ) * fse3v(ji  ,jj  ,jk) * v_bkginc(ji  ,jj  ,jk)     & 
     584                         - e1v(ji  ,jj-1) * fse3v(ji  ,jj-1,jk) * v_bkginc(ji  ,jj-1,jk)  )  & 
     585                         / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 
     586                  END DO 
    458587               END DO 
     588 
     589               CALL lbc_lnk( hdiv, 'T', 1. )   ! lateral boundary cond. (no sign change) 
     590 
     591               DO jj = 2, jpjm1 
     592                  DO ji = fs_2, fs_jpim1   ! vector opt. 
     593                     u_bkginc(ji,jj,jk) = u_bkginc(ji,jj,jk) + 0.2_wp * ( e1t(ji+1,jj)*e2t(ji+1,jj) * hdiv(ji+1,jj)   & 
     594                                                                        - e1t(ji  ,jj)*e2t(ji  ,jj) * hdiv(ji  ,jj) ) & 
     595                                                                      / e1u(ji,jj) * umask(ji,jj,jk)  
     596                     v_bkginc(ji,jj,jk) = v_bkginc(ji,jj,jk) + 0.2_wp * ( e1t(ji,jj+1)*e2t(ji,jj+1) * hdiv(ji,jj+1)   & 
     597                                                                        - e1t(ji,jj  )*e2t(ji,jj  ) * hdiv(ji,jj  ) ) & 
     598                                                                      / e2v(ji,jj) * vmask(ji,jj,jk)  
     599                  END DO 
     600               END DO 
     601 
    459602            END DO 
    460603 
    461             CALL lbc_lnk( hdiv, 'T', 1. )   ! lateral boundary cond. (no sign change) 
    462  
    463             DO jj = 2, jpjm1 
    464                DO ji = fs_2, fs_jpim1   ! vector opt. 
    465                   u_bkginc(ji,jj,jk) = u_bkginc(ji,jj,jk) + 0.2 * ( e1t(ji+1,jj)*e2t(ji+1,jj) * hdiv(ji+1,jj)   & 
    466                                                                   - e1t(ji  ,jj)*e2t(ji  ,jj) * hdiv(ji  ,jj) ) & 
    467                                                                 / e1u(ji,jj) * umask(ji,jj,jk)  
    468                   v_bkginc(ji,jj,jk) = v_bkginc(ji,jj,jk) + 0.2 * ( e1t(ji,jj+1)*e2t(ji,jj+1) * hdiv(ji,jj+1)   & 
    469                                                                   - e1t(ji,jj  )*e2t(ji,jj  ) * hdiv(ji,jj  ) ) & 
    470                                                                 / e2v(ji,jj) * vmask(ji,jj,jk)  
    471                END DO 
    472             END DO 
    473  
    474            END DO 
    475  
    476        END DO 
    477  
    478        CALL wrk_dealloc(jpi,jpj,hdiv)  
     604         END DO 
     605 
     606         CALL wrk_dealloc(jpi,jpj,hdiv)  
    479607 
    480608      ENDIF 
     
    506634         CALL iom_open( c_asmdin, inum ) 
    507635 
    508          CALL iom_get( inum, 'zdate', zdate_bkg )  
     636         CALL iom_get( inum, 'rdastp', zdate_bkg )  
    509637         
    510638         IF(lwp) THEN 
     
    662790      INTEGER :: it 
    663791      REAL(wp) :: zincwgt  ! IAU weight for current time step 
    664       !!---------------------------------------------------------------------- 
     792      REAL (wp), DIMENSION(jpi,jpj,jpk) :: fzptnz ! 3d freezing point values 
     793      !!---------------------------------------------------------------------- 
     794 
     795      ! freezing point calculation taken from oc_fz_pt (but calculated for all depths)  
     796      ! used to prevent the applied increments taking the temperature below the local freezing point  
     797 
     798      ! Note:  For NEMO/CICE this will be identical to nTf (for the surface), but defined at the now point.  
     799  
     800      DO jk=1, jpkm1  
     801         fzptnz (:,:,jk) = tfreez(tsn(:,:,jk,jp_sal )) - 7.53e-4_wp * fsdepw(:,:,jk)  
     802      ENDDO  
    665803 
    666804      IF ( ln_asmiau ) THEN 
     
    684822            ! Update the tracer tendencies 
    685823            DO jk = 1, jpkm1 
    686                tsa(:,:,jk,jp_tem) = tsa(:,:,jk,jp_tem) + t_bkginc(:,:,jk) * zincwgt   
    687                tsa(:,:,jk,jp_sal) = tsa(:,:,jk,jp_sal) + s_bkginc(:,:,jk) * zincwgt 
     824               IF (ln_temnofreeze) THEN 
     825                  ! Do not apply negative increments if the temperature will fall below freezing 
     826                  WHERE(t_bkginc(:,:,jk) > 0.0_wp .OR. & 
     827                     &   tsn(:,:,jk,jp_tem) + tsa(:,:,jk,jp_tem) + t_bkginc(:,:,jk) * wgtiau(it) > fzptnz(:,:,jk) )  
     828                     tsa(:,:,jk,jp_tem) = tsa(:,:,jk,jp_tem) + t_bkginc(:,:,jk) * zincwgt   
     829                  END WHERE 
     830               ELSE 
     831                  tsa(:,:,jk,jp_tem) = tsa(:,:,jk,jp_tem) + t_bkginc(:,:,jk) * zincwgt   
     832               ENDIF 
     833               IF (ln_salfix) THEN 
     834                  ! Do not apply negative increments if the salinity will fall below a specified 
     835                  ! minimum value salfixmin 
     836                  WHERE(s_bkginc(:,:,jk) > 0.0_wp .OR. & 
     837                     &   tsn(:,:,jk,jp_sal) + tsa(:,:,jk,jp_sal) + s_bkginc(:,:,jk) * wgtiau(it) > salfixmin )  
     838                     tsa(:,:,jk,jp_sal) = tsa(:,:,jk,jp_sal) + s_bkginc(:,:,jk) * zincwgt 
     839                  END WHERE 
     840               ELSE 
     841                  tsa(:,:,jk,jp_sal) = tsa(:,:,jk,jp_sal) + s_bkginc(:,:,jk) * zincwgt 
     842               ENDIF 
    688843            END DO 
    689  
    690             ! Salinity fix 
    691             IF (ln_salfix) THEN 
    692                DO jk = 1, jpkm1 
    693                   DO jj = 1, jpj 
    694                      DO ji= 1, jpi 
    695                         tsa(ji,jj,jk,jp_sal) = MAX( tsa(ji,jj,jk,jp_sal), salfixmin ) 
    696                      END DO 
    697                   END DO 
    698                END DO 
    699             ENDIF 
    700844 
    701845         ENDIF 
     
    718862 
    719863            ! Initialize the now fields with the background + increment 
    720             tsn(:,:,:,jp_tem) = t_bkg(:,:,:) + t_bkginc(:,:,:)    
    721             tsn(:,:,:,jp_sal) = s_bkg(:,:,:) + s_bkginc(:,:,:)    
    722  
    723             ! Optional salinity fix 
     864            IF (ln_temnofreeze) THEN 
     865               ! Do not apply negative increments if the temperature will fall below freezing 
     866               WHERE(t_bkginc(:,:,:) > 0.0_wp .OR. & 
     867                  &   tsn(:,:,:,jp_tem) + t_bkginc(:,:,:) > fzptnz(:,:,:) )  
     868                  tsn(:,:,:,jp_tem) = t_bkg(:,:,:) + t_bkginc(:,:,:)    
     869               END WHERE 
     870            ELSE 
     871               tsn(:,:,:,jp_tem) = t_bkg(:,:,:) + t_bkginc(:,:,:)    
     872            ENDIF 
    724873            IF (ln_salfix) THEN 
    725                DO jk = 1, jpkm1 
    726                   DO jj = 1, jpj 
    727                      DO ji= 1, jpi 
    728                         tsn(ji,jj,jk,jp_sal) = MAX( tsn(ji,jj,jk,jp_sal), salfixmin ) 
    729                      END DO 
    730                   END DO 
    731                END DO 
     874               ! Do not apply negative increments if the salinity will fall below a specified 
     875               ! minimum value salfixmin 
     876               WHERE(s_bkginc(:,:,:) > 0.0_wp .OR. & 
     877                  &   tsn(:,:,:,jp_sal) + s_bkginc(:,:,:) > salfixmin )  
     878                  tsn(:,:,:,jp_sal) = s_bkg(:,:,:) + s_bkginc(:,:,:)    
     879               END WHERE 
     880            ELSE 
     881               tsn(:,:,:,jp_sal) = s_bkg(:,:,:) + s_bkginc(:,:,:)    
    732882            ENDIF 
    733883 
    734             tsb(:,:,:,:) = tsn(:,:,:,:)                        ! Update before fields 
     884            tsb(:,:,:,:) = tsn(:,:,:,:)               ! Update before fields 
    735885 
    736886            CALL eos( tsb, rhd, rhop )                ! Before potential and in situ densities 
    737887          
    738888            IF( ln_zps .AND. .NOT. lk_c1d ) & 
    739                &  CALL zps_hde( nit000, jpts, tsb,   &  ! Partial steps: before horizontal derivative 
    740                &                gtsu, gtsv, rhd,        &  ! of T, S, rd at the bottom ocean level 
     889               &  CALL zps_hde( nit000, jpts, tsb, &  ! Partial steps: before horizontal derivative 
     890               &                gtsu, gtsv, rhd,   &  ! of T, S, rd at the bottom ocean level 
    741891               &                gru , grv ) 
     892 
     893#if defined key_zdfkpp 
     894            CALL eos( tsn, rhd )                      ! Compute rhd 
     895#endif 
    742896 
    743897            DEALLOCATE( t_bkginc ) 
     
    748902         !   
    749903      ENDIF 
     904      ! Perhaps the following call should be in step 
     905      IF   ( ln_seaiceinc  )   CALL seaice_asm_inc ( kt )   ! apply sea ice concentration increment 
    750906      ! 
    751907   END SUBROUTINE tra_asm_inc 
     
    817973            vb(:,:,:) = vn(:,:,:) 
    818974  
    819             CALL div_cur( kt )            ! Compute divergence and curl for now fields 
    820  
    821             rotb (:,:,:) = rotn (:,:,:)   ! Update before fields 
    822             hdivb(:,:,:) = hdivn(:,:,:) 
    823  
    824975            DEALLOCATE( u_bkg    ) 
    825976            DEALLOCATE( v_bkg    ) 
     
    846997      ! 
    847998      INTEGER :: it 
     999      INTEGER :: jk 
    8481000      REAL(wp) :: zincwgt  ! IAU weight for current time step 
    8491001      !!---------------------------------------------------------------------- 
     
    8911043            sshn(:,:) = ssh_bkg(:,:) + ssh_bkginc(:,:)   
    8921044 
    893             sshb(:,:) = sshn(:,:)         ! Update before fields 
     1045            ! Update before fields 
     1046            sshb(:,:) = sshn(:,:) 
    8941047 
    8951048            DEALLOCATE( ssh_bkg    ) 
     
    9021055   END SUBROUTINE ssh_asm_inc 
    9031056 
     1057   SUBROUTINE seaice_asm_inc( kt, kindic ) 
     1058      !!---------------------------------------------------------------------- 
     1059      !!                    ***  ROUTINE seaice_asm_inc  *** 
     1060      !!           
     1061      !! ** Purpose : Apply the sea ice assimilation increment. 
     1062      !! 
     1063      !! ** Method  : Direct initialization or Incremental Analysis Updating. 
     1064      !! 
     1065      !! ** Action  :  
     1066      !! 
     1067      !! History : 
     1068      !!        !  07-2011  (D. Lea)  Initial version based on ssh_asm_inc 
     1069      !!---------------------------------------------------------------------- 
     1070 
     1071      IMPLICIT NONE 
     1072 
     1073      !! * Arguments 
     1074      INTEGER, INTENT(IN) :: kt   ! Current time step 
     1075      INTEGER, OPTIONAL, INTENT(IN) :: kindic ! flag for disabling the deallocation 
     1076 
     1077      !! * Local declarations 
     1078      INTEGER :: it 
     1079      REAL(wp) :: zincwgt  ! IAU weight for current time step 
     1080 
     1081#if defined key_lim3 || defined key_lim2 
     1082      REAL(wp), DIMENSION(jpi,jpj) :: zofrld, zohicif, zseaicendg, zhicifinc  ! LIM 
     1083      REAL(wp) :: zhicifmin=0.5_wp      ! ice minimum depth in metres 
     1084 
     1085#endif 
     1086 
     1087 
     1088      IF ( ln_asmiau ) THEN 
     1089 
     1090         !-------------------------------------------------------------------- 
     1091         ! Incremental Analysis Updating 
     1092         !-------------------------------------------------------------------- 
     1093 
     1094         IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN 
     1095 
     1096            it = kt - nit000 + 1 
     1097            zincwgt = wgtiau(it)      ! IAU weight for the current time step  
     1098            ! note this is not a tendency so should not be divided by rdt (as with the tracer and other increments) 
     1099 
     1100            IF(lwp) THEN 
     1101               WRITE(numout,*)  
     1102               WRITE(numout,*) 'seaice_asm_inc : sea ice conc IAU at time step = ', & 
     1103                  &  kt,' with IAU weight = ', wgtiau(it) 
     1104               WRITE(numout,*) '~~~~~~~~~~~~' 
     1105            ENDIF 
     1106 
     1107#if defined key_lim3 || defined key_lim2 
     1108 
     1109            zofrld(:,:)=frld(:,:) 
     1110            zohicif(:,:)=hicif(:,:) 
     1111 
     1112            frld = MIN( MAX( frld(:,:) - seaice_bkginc(:,:) * zincwgt, 0.0_wp), 1.0_wp) 
     1113            pfrld = MIN( MAX( pfrld(:,:) - seaice_bkginc(:,:) * zincwgt, 0.0_wp), 1.0_wp) 
     1114            fr_i(:,:) = 1.0_wp - frld(:,:)        ! adjust ice fraction 
     1115 
     1116            zseaicendg(:,:)=zofrld(:,:) - frld(:,:)         ! find out actual sea ice nudge applied 
     1117 
     1118            ! Nudge sea ice depth to bring it up to a required minimum depth 
     1119 
     1120            WHERE( zseaicendg(:,:) > 0.0_wp .AND. hicif(:,:) < zhicifmin )  
     1121               zhicifinc(:,:) = (zhicifmin - hicif(:,:)) * zincwgt     
     1122            ELSEWHERE 
     1123               zhicifinc(:,:) = 0.0_wp 
     1124            END WHERE 
     1125 
     1126! nudge ice depth 
     1127            hicif(:,:)=hicif(:,:) + zhicifinc(:,:) 
     1128            phicif(:,:)=phicif(:,:) + zhicifinc(:,:)        
     1129 
     1130! seaice salinity balancing (to add) 
     1131 
     1132#endif 
     1133 
     1134#if defined key_cice 
     1135 
     1136! Pass ice increment tendency into CICE 
     1137            ndaice_da(:,:) = seaice_bkginc(:,:) * zincwgt / rdt 
     1138 
     1139#endif 
     1140 
     1141            IF ( kt == nitiaufin_r ) THEN 
     1142               DEALLOCATE( seaice_bkginc ) 
     1143            ENDIF 
     1144 
     1145         ELSE 
     1146 
     1147#if defined key_cice 
     1148 
     1149! Zero ice increment tendency into CICE 
     1150            ndaice_da(:,:) = 0.0_wp 
     1151 
     1152#endif 
     1153 
     1154         ENDIF 
     1155 
     1156      ELSEIF ( ln_asmdin ) THEN 
     1157 
     1158         !-------------------------------------------------------------------- 
     1159         ! Direct Initialization 
     1160         !-------------------------------------------------------------------- 
     1161 
     1162         IF ( kt == nitdin_r ) THEN 
     1163 
     1164            neuler = 0                    ! Force Euler forward step 
     1165 
     1166#if defined key_lim3 || defined key_lim2 
     1167 
     1168            zofrld(:,:)=frld(:,:) 
     1169            zohicif(:,:)=hicif(:,:) 
     1170  
     1171            ! Initialize the now fields the background + increment 
     1172 
     1173            frld(:,:) = MIN( MAX( frld(:,:) - seaice_bkginc(:,:), 0.0_wp), 1.0_wp) 
     1174            pfrld(:,:) = frld(:,:)  
     1175            fr_i(:,:) = 1.0_wp - frld(:,:)        ! adjust ice fraction 
     1176 
     1177            zseaicendg(:,:)=zofrld(:,:) - frld(:,:)         ! find out actual sea ice nudge applied 
     1178 
     1179            ! Nudge sea ice depth to bring it up to a required minimum depth 
     1180 
     1181            WHERE( zseaicendg(:,:) > 0.0_wp .AND. hicif(:,:) < zhicifmin )  
     1182               zhicifinc(:,:) = (zhicifmin - hicif(:,:)) * zincwgt     
     1183            ELSEWHERE 
     1184               zhicifinc(:,:) = 0.0_wp 
     1185            END WHERE 
     1186 
     1187! nudge ice depth 
     1188            hicif(:,:)=hicif(:,:) + zhicifinc(:,:) 
     1189            phicif(:,:)=phicif(:,:)        
     1190 
     1191! seaice salinity balancing (to add) 
     1192   
     1193#endif 
     1194  
     1195#if defined key_cice 
     1196 
     1197! Pass ice increment tendency into CICE - is this correct? 
     1198           ndaice_da(:,:) = seaice_bkginc(:,:) / rdt 
     1199 
     1200#endif 
     1201           IF ( .NOT. PRESENT(kindic) ) THEN 
     1202              DEALLOCATE( seaice_bkginc ) 
     1203           END IF 
     1204 
     1205         ELSE 
     1206 
     1207#if defined key_cice 
     1208 
     1209! Zero ice increment tendency into CICE  
     1210            ndaice_da(:,:) = 0.0_wp 
     1211 
     1212#endif 
     1213          
     1214         ENDIF 
     1215 
     1216!#if defined key_lim3 || defined key_lim2 || defined key_cice 
     1217! 
     1218!            IF (ln_seaicebal ) THEN        
     1219!             !! balancing salinity increments 
     1220!             !! simple case from limflx.F90 (doesn't include a mass flux) 
     1221!             !! assumption is that as ice concentration is reduced or increased 
     1222!             !! the snow and ice depths remain constant 
     1223!             !! note that snow is being created where ice concentration is being increased 
     1224!             !! - could be more sophisticated and 
     1225!             !! not do this (but would need to alter h_snow) 
     1226! 
     1227!             usave(:,:,:)=sb(:,:,:)   ! use array as a temporary store 
     1228! 
     1229!             DO jj = 1, jpj 
     1230!               DO ji = 1, jpi  
     1231!           ! calculate change in ice and snow mass per unit area 
     1232!           ! positive values imply adding salt to the ocean (results from ice formation) 
     1233!           ! fwf : ice formation and melting 
     1234! 
     1235!                 zfons = ( -nfresh_da(ji,jj)*soce + nfsalt_da(ji,jj) )*rdt 
     1236! 
     1237!           ! change salinity down to mixed layer depth 
     1238!                 mld=hmld_kara(ji,jj) 
     1239! 
     1240!           ! prevent small mld 
     1241!           ! less than 10m can cause salinity instability  
     1242!                 IF (mld < 10) mld=10 
     1243! 
     1244!           ! set to bottom of a level  
     1245!                 DO jk = jpk-1, 2, -1 
     1246!                   IF ((mld > gdepw(ji,jj,jk)) .and. (mld < gdepw(ji,jj,jk+1))) THEN  
     1247!                     mld=gdepw(ji,jj,jk+1) 
     1248!                     jkmax=jk 
     1249!                   ENDIF 
     1250!                 ENDDO 
     1251! 
     1252!            ! avoid applying salinity balancing in shallow water or on land 
     1253!            !  
     1254! 
     1255!            ! dsal_ocn (psu kg m^-2) / (kg m^-3 * m) 
     1256! 
     1257!                 dsal_ocn=0.0_wp 
     1258!                 sal_thresh=5.0_wp        ! minimum salinity threshold for salinity balancing 
     1259! 
     1260!                 if (tmask(ji,jj,1) > 0 .AND. tmask(ji,jj,jkmax) > 0 ) & 
     1261!                              dsal_ocn = zfons / (rhop(ji,jj,1) * mld) 
     1262! 
     1263!           ! put increments in for levels in the mixed layer 
     1264!           ! but prevent salinity below a threshold value  
     1265! 
     1266!                   DO jk = 1, jkmax               
     1267! 
     1268!                     IF (dsal_ocn > 0.0_wp .or. sb(ji,jj,jk)+dsal_ocn > sal_thresh) THEN  
     1269!                           sb(ji,jj,jk) = sb(ji,jj,jk) + dsal_ocn 
     1270!                           sn(ji,jj,jk) = sn(ji,jj,jk) + dsal_ocn 
     1271!                     ENDIF 
     1272! 
     1273!                   ENDDO 
     1274! 
     1275!      !            !  salt exchanges at the ice/ocean interface 
     1276!      !            zpmess         = zfons / rdt_ice    ! rdt_ice is ice timestep 
     1277!      ! 
     1278!      !! Adjust fsalt. A +ve fsalt means adding salt to ocean 
     1279!      !!           fsalt(ji,jj) =  fsalt(ji,jj) + zpmess     ! adjust fsalt   
     1280!      !!                
     1281!      !!           emps(ji,jj) = emps(ji,jj) + zpmess        ! or adjust emps (see icestp1d)  
     1282!      !!                                                     ! E-P (kg m-2 s-2) 
     1283!      !            emp(ji,jj) = emp(ji,jj) + zpmess          ! E-P (kg m-2 s-2) 
     1284!               ENDDO !ji 
     1285!             ENDDO !jj! 
     1286! 
     1287!            ENDIF !ln_seaicebal 
     1288! 
     1289!#endif 
     1290 
     1291 
     1292      ENDIF 
     1293 
     1294   END SUBROUTINE seaice_asm_inc 
    9041295   !!====================================================================== 
    9051296END MODULE asminc 
  • branches/UKMO/CO5_package_branch/NEMOGCM/NEMO/OPA_SRC/ASM/asmtrj.F90

    r7363 r7367  
    1414   !!                 ! 2009-06 (F. Vigilant) asm_trj_wri: special case when kt=nit000-1 
    1515   !!                 ! 2009-07 (F. Vigilant) asm_trj_wri: add computation of eiv at restart 
     16    
     17   !!                 ! 2012-11 (A. Weaver) Save avt_bkg for mixing layer computation, remove en_bkg 
    1618   !!---------------------------------------------------------------------- 
    1719 
     
    3537   USE zdfmxl             ! Mixed layer depth 
    3638   USE dom_oce, ONLY :   ndastp 
    37    USE sol_oce, ONLY :   gcx   ! Solver variables defined in memory 
    3839   USE in_out_manager     ! I/O manager 
    3940   USE iom                ! I/O module 
     
    4344   USE ldfeiv             ! eddy induced velocity coef.      (ldf_eiv routine) 
    4445#endif 
    45  
     46#if defined key_lim2 
     47   USE ice_2 
     48#endif 
     49#if defined key_lim3 
     50   USE ice 
     51#endif 
    4652   IMPLICIT NONE 
    4753   PRIVATE 
     
    110116            CALL iom_rstput( kt, nitbkg_r, inum, 'tn'     , tsn(:,:,:,jp_tem) ) 
    111117            CALL iom_rstput( kt, nitbkg_r, inum, 'sn'     , tsn(:,:,:,jp_sal) ) 
     118            CALL iom_rstput( kt, nitbkg_r, inum, 'avt'    , avt               )             
    112119            CALL iom_rstput( kt, nitbkg_r, inum, 'sshn'   , sshn              ) 
    113 #if defined key_zdftke 
    114             CALL iom_rstput( kt, nitbkg_r, inum, 'en'     , en                ) 
    115 #endif 
    116             CALL iom_rstput( kt, nitbkg_r, inum, 'gcx'    , gcx               ) 
    117120            ! 
    118121            CALL iom_close( inum ) 
     
    148151            CALL iom_rstput( kt, nitdin_r, inum, 'tn'     , tsn(:,:,:,jp_tem) ) 
    149152            CALL iom_rstput( kt, nitdin_r, inum, 'sn'     , tsn(:,:,:,jp_sal) ) 
     153            CALL iom_rstput( kt, nitbkg_r, inum, 'avt'    , avt               )             
    150154            CALL iom_rstput( kt, nitdin_r, inum, 'sshn'   , sshn              ) 
     155#if defined key_lim2 || defined key_lim3 
     156            IF(( nn_ice == 2 ) .OR. ( nn_ice == 3 )) THEN 
     157               CALL iom_rstput( kt, nitdin_r, inum, 'iceconc', 1.0 - frld(:,:)   ) 
     158            ENDIF 
     159#endif 
    151160            ! 
    152161            CALL iom_close( inum ) 
     
    222231         CALL iom_rstput( it, it, inum, 'avt'   , avt    ) 
    223232#if defined key_ldfslp 
    224          CALL iom_rstput( it, it, inum, 'uslp'  , uslp   ) 
    225          CALL iom_rstput( it, it, inum, 'vslp'  , vslp   ) 
    226          CALL iom_rstput( it, it, inum, 'wslpi' , wslpi  ) 
    227          CALL iom_rstput( it, it, inum, 'wslpj' , wslpj  ) 
     233         CALL iom_rstput( it, it, inum, 'uslp_hor'  , uslp_hor   ) 
     234         CALL iom_rstput( it, it, inum, 'vslp_hor'  , vslp_hor   ) 
     235         CALL iom_rstput( it, it, inum, 'wslpi_hor' , wslpi_hor  ) 
     236         CALL iom_rstput( it, it, inum, 'wslpj_hor' , wslpj_hor  ) 
     237         CALL iom_rstput( it, it, inum, 'uslp_iso'  , uslp_iso   ) 
     238         CALL iom_rstput( it, it, inum, 'vslp_iso'  , vslp_iso   ) 
     239         CALL iom_rstput( it, it, inum, 'wslpi_iso' , wslpi_iso  ) 
     240         CALL iom_rstput( it, it, inum, 'wslpj_iso' , wslpj_iso  ) 
    228241#endif 
    229242#if defined key_zdfddm 
  • branches/UKMO/CO5_package_branch/NEMOGCM/NEMO/OPA_SRC/BDY/bdy_oce.F90

    r7363 r7367  
    5757   LOGICAL                    ::   ln_mask_file             !: =T read bdymask from file 
    5858   LOGICAL                    ::   ln_vol                   !: =T volume correction              
     59   LOGICAL, DIMENSION(jp_bdy) ::   ln_sponge                !: =T use sponge layer  
    5960   ! 
    6061   INTEGER                    ::   nb_bdy                   !: number of open boundary sets 
     
    6263   INTEGER                    ::   nn_volctl                !: = 0 the total volume will have the variability of the surface Flux E-P  
    6364   !                                                        !  = 1 the volume will be constant during all the integration. 
     65   REAL(wp)                   ::   rn_sponge                !: multiplier of diffusion for sponge layer 
     66 
    6467   INTEGER, DIMENSION(jp_bdy) ::   nn_dyn2d                 ! Choice of boundary condition for barotropic variables (U,V,SSH) 
    6568   INTEGER, DIMENSION(jp_bdy) ::   nn_dyn2d_dta           !: = 0 use the initial state as bdy dta ;  
     
    8386   !! Global variables 
    8487   !!---------------------------------------------------------------------- 
    85    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   bdytmask   !: Mask defining computational domain at T-points 
    86    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   bdyumask   !: Mask defining computational domain at U-points 
    87    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   bdyvmask   !: Mask defining computational domain at V-points 
     88   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   bdytmask      !: Mask defining computational domain at T-points 
     89   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   bdyumask      !: Mask defining computational domain at U-points 
     90   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   bdyvmask      !: Mask defining computational domain at V-points 
     91   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   sponge_factor !: Multiplier for diffusion for sponge layer 
    8892 
    8993   REAL(wp)                                    ::   bdysurftot !: Lateral surface of unstructured open boundary 
     
    120124      ! 
    121125      ALLOCATE( bdytmask(jpi,jpj) , bdyumask(jpi,jpj), bdyvmask(jpi,jpj),                    &   
    122          &      STAT=bdy_oce_alloc ) 
     126         &      sponge_factor(jpi,jpj), STAT=bdy_oce_alloc ) 
    123127         ! 
    124128      IF( lk_mpp             )   CALL mpp_sum ( bdy_oce_alloc ) 
  • branches/UKMO/CO5_package_branch/NEMOGCM/NEMO/OPA_SRC/BDY/bdydta.F90

    r7363 r7367  
    3232   USE ice_2 
    3333#endif 
     34   USE sbcapr 
    3435 
    3536   IMPLICIT NONE 
     
    238239               ENDIF 
    239240            ENDIF 
    240             jstart = jend+1 
     241            jstart = jstart + nb_bdy_fld(ib_bdy) 
    241242 
    242243            ! If full velocities in boundary data then split into barotropic and baroclinic data 
     
    281282         END IF ! nn_dta(ib_bdy) = 1 
    282283      END DO  ! ib_bdy 
     284 
     285      IF ( ln_apr_obc ) THEN 
     286         DO ib_bdy = 1, nb_bdy 
     287            IF (nn_tra(ib_bdy).NE.4)THEN 
     288               igrd = 1                      ! meridional velocity 
     289               DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd) 
     290                  ii   = idx_bdy(ib_bdy)%nbi(ib,igrd) 
     291                  ij   = idx_bdy(ib_bdy)%nbj(ib,igrd) 
     292                  dta_bdy(ib_bdy)%ssh(ib) = dta_bdy(ib_bdy)%ssh(ib) + ssh_ib(ii,ij) 
     293               ENDDO 
     294            ENDIF 
     295         ENDDO 
     296      ENDIF 
    283297 
    284298      IF( nn_timing == 1 ) CALL timing_stop('bdy_dta') 
     
    317331      TYPE(FLD_N) ::   bn_frld, bn_hicif, bn_hsnif      ! 
    318332#endif 
    319       NAMELIST/nambdy_dta/ cn_dir, bn_tem, bn_sal, bn_u3d, bn_v3d, bn_ssh, bn_u2d, bn_v2d  
    320 #if defined key_lim2 
    321       NAMELIST/nambdy_dta/ bn_frld, bn_hicif, bn_hsnif 
    322 #endif 
    323       NAMELIST/nambdy_dta/ ln_full_vel 
     333      NAMELIST/nambdy_dta_1/ cn_dir, bn_tem, bn_sal, bn_u3d, bn_v3d, bn_ssh, bn_u2d, bn_v2d  
     334      NAMELIST/nambdy_dta_2/ cn_dir, bn_tem, bn_sal, bn_u3d, bn_v3d, bn_ssh, bn_u2d, bn_v2d  
     335#if defined key_lim2 
     336      NAMELIST/nambdy_dta_1/ bn_frld, bn_hicif, bn_hsnif 
     337      NAMELIST/nambdy_dta_2/ bn_frld, bn_hicif, bn_hsnif 
     338#endif 
     339      NAMELIST/nambdy_dta_1/ ln_full_vel 
     340      NAMELIST/nambdy_dta_2/ ln_full_vel 
    324341      !!--------------------------------------------------------------------------- 
    325342 
     
    403420 
    404421            ! Important NOT to rewind here. 
    405             READ( numnam, nambdy_dta ) 
     422            if ( ib_bdy == 1 )  READ( numnam, nambdy_dta_1 ) 
     423            if ( ib_bdy == 2 )  READ( numnam, nambdy_dta_2 ) 
    406424 
    407425            cn_dir_array(ib_bdy) = cn_dir 
     
    554572            ! Recalculate field counts 
    555573            !------------------------- 
    556             nb_bdy_fld_sum = 0 
    557574            IF( ib_bdy .eq. 1 ) THEN  
     575               nb_bdy_fld_sum = 0 
    558576               nb_bdy_fld(ib_bdy) = jfld 
    559577               nb_bdy_fld_sum     = jfld               
  • branches/UKMO/CO5_package_branch/NEMOGCM/NEMO/OPA_SRC/BDY/bdyice_lim2.F90

    r7363 r7367  
    5353            CYCLE 
    5454         CASE(jp_frs) 
    55             CALL bdy_ice_frs( idx_bdy(ib_bdy), dta_idx(ib_bdy) ) 
     55            CALL bdy_ice_frs( idx_bdy(ib_bdy), dta_bdy(ib_bdy) ) 
    5656         CASE DEFAULT 
    5757            CALL ctl_stop( 'bdy_ice_lim_2 : unrecognised option for open boundaries for ice fields' ) 
  • branches/UKMO/CO5_package_branch/NEMOGCM/NEMO/OPA_SRC/BDY/bdyini.F90

    r7363 r7367  
    8383         &             nn_ice_lim2, nn_ice_lim2_dta,                       & 
    8484#endif 
    85          &             ln_vol, nn_volctl, nn_rimwidth 
     85         &             ln_vol, nn_volctl, ln_sponge, rn_sponge, nn_rimwidth 
    8686      !! 
    8787      NAMELIST/nambdy_index/ nbdysege, jpieob, jpjedt, jpjeft,             & 
     
    127127      ln_vol            = .false. 
    128128      nn_volctl         = -1  ! uninitialised flag 
     129      ln_sponge(:)      = .false.  
     130      rn_sponge         = 0.0  
    129131      nn_rimwidth(:)    = -1  ! uninitialised flag 
    130132 
     
    224226        IF(lwp) WRITE(numout,*) 
    225227 
     228        IF( ln_sponge(ib_bdy) ) THEN                     ! check sponge layer choice  
     229          IF(lwp) WRITE(numout,*) 'Sponge layer applied at open boundaries'  
     230          IF(lwp) WRITE(numout,*) 'Multiplier for diffusion in sponge layer : ', rn_sponge  
     231          IF(lwp) WRITE(numout,*)  
     232        ELSE  
     233          IF(lwp) WRITE(numout,*) 'No volume correction applied at open boundaries'  
     234          IF(lwp) WRITE(numout,*)  
     235        ENDIF  
     236  
    226237      ENDDO 
    227238 
    228      IF( ln_vol ) THEN                     ! check volume conservation (nn_volctl value) 
    229        IF(lwp) WRITE(numout,*) 'Volume correction applied at open boundaries' 
    230        IF(lwp) WRITE(numout,*) 
    231        SELECT CASE ( nn_volctl ) 
    232          CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '      The total volume will be constant' 
    233          CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '      The total volume will vary according to the surface E-P flux' 
    234          CASE DEFAULT   ;   CALL ctl_stop( 'nn_volctl must be 0 or 1' ) 
    235        END SELECT 
    236        IF(lwp) WRITE(numout,*) 
    237      ELSE 
    238        IF(lwp) WRITE(numout,*) 'No volume correction applied at open boundaries' 
    239        IF(lwp) WRITE(numout,*) 
    240      ENDIF 
     239      IF( ln_vol ) THEN                     ! check volume conservation (nn_volctl value) 
     240        IF(lwp) WRITE(numout,*) 'Volume correction applied at open boundaries' 
     241        IF(lwp) WRITE(numout,*) 
     242        SELECT CASE ( nn_volctl ) 
     243          CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '      The total volume will be constant' 
     244          CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '      The total volume will vary according to the surface E-P flux' 
     245          CASE DEFAULT   ;   CALL ctl_stop( 'nn_volctl must be 0 or 1' ) 
     246        END SELECT 
     247        IF(lwp) WRITE(numout,*) 
     248      ELSE 
     249        IF(lwp) WRITE(numout,*) 'No volume correction applied at open boundaries' 
     250        IF(lwp) WRITE(numout,*) 
     251      ENDIF 
     252 
     253      sponge_factor(:,:) = 1.0  
    241254 
    242255      ! ------------------------------------------------- 
     
    247260      ! --------------------------------------------- 
    248261      REWIND( numnam )                     
     262      jpbdta = 1 
    249263      DO ib_bdy = 1, nb_bdy 
    250264 
    251          jpbdta = 1 
    252265         IF( .NOT. ln_coords_file(ib_bdy) ) THEN ! Work out size of global arrays from namelist parameters 
    253266  
     
    317330 
    318331            nblendta(:,ib_bdy) = nblendta(:,ib_bdy) * nn_rimwidth(ib_bdy) 
    319             jpbdta = MAXVAL(nblendta(:,ib_bdy))                
     332            jpbdta = MAX( jpbdta, MAXVAL(nblendta(:,ib_bdy)) ) 
    320333 
    321334 
     
    324337 
    325338            CALL iom_open( cn_coords_file(ib_bdy), inum ) 
    326             jpbdta = 1 
     339 
    327340            DO igrd = 1, jpbgrd 
    328341               id_dummy = iom_varid( inum, 'nbi'//cgrid(igrd), kdimsz=kdimsz )   
     
    330343               jpbdta = MAX(jpbdta, kdimsz(1)) 
    331344            ENDDO 
     345            CALL iom_close( inum ) 
    332346 
    333347         ENDIF  
     
    507521         ELSE            ! Read global index arrays from boundary coordinates file. 
    508522 
     523            CALL iom_open( cn_coords_file(ib_bdy), inum ) 
    509524            DO igrd = 1, jpbgrd 
    510525               CALL iom_get( inum, jpdom_unknown, 'nbi'//cgrid(igrd), dta_global(1:nblendta(igrd,ib_bdy),:,1) ) 
     
    616631         END DO  
    617632 
    618       ENDDO 
     633         ! Compute multiplier for diffusion for sponge layer  
     634         ! -------------------------------------------------  
     635         IF( ln_sponge(ib_bdy) ) THEN  
     636            igrd=1 
     637            DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd)  
     638               nbi => idx_bdy(ib_bdy)%nbi(ib,igrd)  
     639               nbj => idx_bdy(ib_bdy)%nbj(ib,igrd)  
     640               nbr => idx_bdy(ib_bdy)%nbr(ib,igrd)  
     641               sponge_factor(nbi,nbj) = 1.0 + (rn_sponge-1.0) * ( 1.- TANH( FLOAT( nbr - 1 ) *0.5 ) )  
     642            END DO  
     643         ENDIF  
     644 
     645      ENDDO  ! ib_bdy 
    619646 
    620647      ! ------------------------------------------------------ 
     
    773800            DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd) 
    774801               nbi => idx_bdy(ib_bdy)%nbi(ib,igrd) 
    775                nbj => idx_bdy(ib_bdy)%nbi(ib,igrd) 
     802               nbj => idx_bdy(ib_bdy)%nbj(ib,igrd) 
    776803               flagu => idx_bdy(ib_bdy)%flagu(ib) 
    777804               bdysurftot = bdysurftot + hu     (nbi  , nbj)                           & 
     
    786813            DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd) 
    787814               nbi => idx_bdy(ib_bdy)%nbi(ib,igrd) 
    788                nbj => idx_bdy(ib_bdy)%nbi(ib,igrd) 
     815               nbj => idx_bdy(ib_bdy)%nbj(ib,igrd) 
    789816               flagv => idx_bdy(ib_bdy)%flagv(ib) 
    790817               bdysurftot = bdysurftot + hv     (nbi, nbj  )                           & 
  • branches/UKMO/CO5_package_branch/NEMOGCM/NEMO/OPA_SRC/DIA/diadct.F90

    r3294 r7367  
    3838  USE dianam          ! build name of file 
    3939  USE lib_mpp         ! distributed memory computing library 
    40 #if defined key_lim2 || defined key_lim3 
    41   USE ice 
     40#if defined key_lim2 
     41  USE ice_2 
     42#endif 
     43#if defined key_lim3 
     44  USE ice_3 
    4245#endif 
    4346  USE domvvl 
     
    4952 
    5053  !! * Routine accessibility 
    51   PUBLIC   dia_dct     ! routine called by step.F90 
    52   PUBLIC   dia_dct_init! routine called by opa.F90 
     54  PUBLIC   dia_dct      ! routine called by step.F90 
     55  PUBLIC   dia_dct_init ! routine called by opa.F90 
     56  PUBLIC   diadct_alloc ! routine called by nemo_init in nemogcm.F90 
    5357  PRIVATE  readsec 
    5458  PRIVATE  removepoints 
    5559  PRIVATE  transport 
    5660  PRIVATE  dia_dct_wri 
     61  PRIVATE  dia_dct_wri_NOOS 
    5762 
    5863#include "domzgr_substitute.h90" 
     
    6570  INTEGER :: nn_dctwri   = 1     ! Frequency of output 
    6671  INTEGER :: nn_secdebug = 0     ! Number of the section to debug 
     72  INTEGER :: nn_dct_h    = 1     ! Frequency of computation for NOOS hourly files 
     73  INTEGER :: nn_dctwri_h = 1     ! Frequency of output for NOOS hourly files 
    6774    
    68   INTEGER, PARAMETER :: nb_class_max  = 10 
    69   INTEGER, PARAMETER :: nb_sec_max    = 150 
    70   INTEGER, PARAMETER :: nb_point_max  = 2000 
    71   INTEGER, PARAMETER :: nb_type_class = 14 
     75  INTEGER, PARAMETER :: nb_class_max  = 11   ! maximum number of classes, i.e. depth levels or density classes 
     76  INTEGER, PARAMETER :: nb_sec_max    = 30   ! maximum number of sections 
     77  INTEGER, PARAMETER :: nb_point_max  = 375  ! maximum number of points in a single section 
     78  INTEGER, PARAMETER :: nb_type       = 14   ! types of calculations, i.e. pos transport, neg transport, heat transport, salt transport 
     79  INTEGER, PARAMETER :: nb_3d_vars    = 5 
     80  INTEGER, PARAMETER :: nb_2d_vars    = 2 
    7281  INTEGER            :: nb_sec  
    7382 
     
    8291  TYPE SECTION 
    8392     CHARACTER(len=60)                            :: name              ! name of the sec 
    84      LOGICAL                                      :: llstrpond         ! true if you want the computation of salt and 
    85                                                                        ! heat transports 
     93     LOGICAL                                      :: llstrpond         ! true if you want the computation of salinity and heat transports 
    8694     LOGICAL                                      :: ll_ice_section    ! ice surface and ice volume computation 
    8795     LOGICAL                                      :: ll_date_line      ! = T if the section crosses the date-line 
     
    95103                                                     ztem            ,&! temperature classes(99 if you don't want) 
    96104                                                     zlay              ! level classes      (99 if you don't want) 
    97      REAL(wp), DIMENSION(nb_type_class,nb_class_max)  :: transport     ! transport output 
     105     REAL(wp), DIMENSION(nb_type,nb_class_max)        :: transport     ! transport output 
     106     REAL(wp), DIMENSION(nb_type,nb_class_max)        :: transport_h   ! transport output 
    98107     REAL(wp)                                         :: slopeSection  ! slope of the section 
    99108     INTEGER                                          :: nb_point      ! number of points in the section 
     
    102111 
    103112  TYPE(SECTION),DIMENSION(nb_sec_max) :: secs ! Array of sections 
    104   
     113 
     114  REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) ::  transports_3d 
     115  REAL(wp), ALLOCATABLE, DIMENSION(:,:,:)   ::  transports_2d 
     116  REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) ::  transports_3d_h 
     117  REAL(wp), ALLOCATABLE, DIMENSION(:,:,:)   ::  transports_2d_h 
     118  REAL(wp), ALLOCATABLE, DIMENSION(:,:,:)   ::  z_hr_output 
    105119  
    106120CONTAINS 
     121 
     122  INTEGER FUNCTION diadct_alloc() 
     123     !!---------------------------------------------------------------------- 
     124     !!                   ***  FUNCTION diadct_alloc  *** 
     125     !!---------------------------------------------------------------------- 
     126     INTEGER :: ierr(2) 
     127     !!---------------------------------------------------------------------- 
     128     ! 
     129     ALLOCATE(transports_3d(nb_3d_vars,nb_sec_max,nb_point_max,jpk),   STAT=ierr(1) ) 
     130     ALLOCATE(transports_2d(nb_2d_vars,nb_sec_max,nb_point_max)    ,   STAT=ierr(2) ) 
     131     ALLOCATE(transports_3d_h(nb_3d_vars,nb_sec_max,nb_point_max,jpk), STAT=ierr(3) ) 
     132     ALLOCATE(transports_2d_h(nb_2d_vars,nb_sec_max,nb_point_max)    , STAT=ierr(4) ) 
     133     ALLOCATE(z_hr_output(nb_sec_max,24,nb_class_max)                , STAT=ierr(5) ) 
     134        ! 
     135     diadct_alloc = MAXVAL( ierr ) 
     136     IF( diadct_alloc /= 0 )   CALL ctl_warn('diadct_alloc: failed to allocate arrays') 
     137     ! 
     138  END FUNCTION diadct_alloc 
     139 
    107140 
    108141  SUBROUTINE dia_dct_init 
     
    110143     !!               ***  ROUTINE diadct  ***   
    111144     !! 
    112      !!  ** Purpose: Read the namelist parametres 
     145     !!  ** Purpose: Read the namelist parameters 
    113146     !!              Open output files 
    114147     !! 
     
    121154     REWIND( numnam ) 
    122155     READ  ( numnam, namdct ) 
     156 
     157     IF( ln_NOOS ) THEN 
     158        nn_dct=3600./rdt         ! hard coded for NOOS transects, to give 25 hour means  
     159        nn_dctwri=86400./rdt 
     160        nn_dct_h=1       ! hard coded for NOOS transects, to give hourly data 
     161        nn_dctwri_h=3600./rdt 
     162     ENDIF 
    123163 
    124164     IF( lwp ) THEN 
     
    126166        WRITE(numout,*) "diadct_init: compute transports through sections " 
    127167        WRITE(numout,*) "~~~~~~~~~~~~~~~~~~~~~" 
    128         WRITE(numout,*) "       Frequency of computation: nn_dct    = ",nn_dct 
    129         WRITE(numout,*) "       Frequency of write:       nn_dctwri = ",nn_dctwri 
     168        IF( ln_NOOS ) THEN 
     169           WRITE(numout,*) "       Frequency of computation hard coded to be every hour: nn_dct    = ",nn_dct 
     170           WRITE(numout,*) "       Frequency of write hard coded to average 25 instantaneous hour values: nn_dctwri = ",nn_dctwri 
     171           WRITE(numout,*) "       Frequency of hourly computation hard coded to be every timestep: nn_dct_h  = ",nn_dct_h 
     172           WRITE(numout,*) "       Frequency of hourly write hard coded to every hour: nn_dctwri_h = ",nn_dctwri_h 
     173        ELSE 
     174           WRITE(numout,*) "       Frequency of computation: nn_dct    = ",nn_dct 
     175           WRITE(numout,*) "       Frequency of write:       nn_dctwri = ",nn_dctwri 
     176        ENDIF 
    130177 
    131178        IF      ( nn_secdebug .GE. 1 .AND. nn_secdebug .LE. nb_sec_max )THEN 
     
    146193     !open output file 
    147194     IF( lwp ) THEN 
    148         CALL ctl_opn( numdct_vol,  'volume_transport', 'NEW', 'FORMATTED', 'SEQUENTIAL', -1, numout,  .FALSE. ) 
    149         CALL ctl_opn( numdct_heat, 'heat_transport'  , 'NEW', 'FORMATTED', 'SEQUENTIAL', -1, numout,  .FALSE. ) 
    150         CALL ctl_opn( numdct_salt, 'salt_transport'  , 'NEW', 'FORMATTED', 'SEQUENTIAL', -1, numout,  .FALSE. ) 
     195        IF( ln_NOOS ) THEN 
     196           CALL ctl_opn( numdct_NOOS  ,'NOOS_transport'  , 'NEW', 'FORMATTED', 'SEQUENTIAL', -1, numout,  .FALSE. ) 
     197           CALL ctl_opn( numdct_NOOS_h,'NOOS_transport_h', 'NEW', 'FORMATTED', 'SEQUENTIAL', -1, numout,  .FALSE. ) 
     198        ELSE 
     199           CALL ctl_opn( numdct_vol , 'volume_transport', 'NEW', 'FORMATTED', 'SEQUENTIAL', -1, numout,  .FALSE. ) 
     200           CALL ctl_opn( numdct_temp, 'heat_transport'  , 'NEW', 'FORMATTED', 'SEQUENTIAL', -1, numout,  .FALSE. ) 
     201           CALL ctl_opn( numdct_sal , 'salt_transport'  , 'NEW', 'FORMATTED', 'SEQUENTIAL', -1, numout,  .FALSE. ) 
     202        ENDIF 
    151203     ENDIF 
     204 
     205     ! Initialise arrays to zero 
     206     transports_3d(:,:,:,:)  =0._wp 
     207     transports_2d(:,:,:)    =0._wp 
     208     transports_3d_h(:,:,:,:)=0._wp 
     209     transports_2d_h(:,:,:)  =0._wp 
     210     z_hr_output(:,:,:)      =0._wp 
    152211 
    153212     IF( nn_timing == 1 )   CALL timing_stop('dia_dct_init') 
     
    160219     !!               ***  ROUTINE diadct  ***   
    161220     !! 
    162      !!  ** Purpose: Compute sections tranport and write it in numdct file 
     221     !!  Purpose :: Compute section transports and write it in numdct files 
     222     !!   
     223     !!  Method  :: All arrays initialised to zero in dct_init 
     224     !!             Each nn_dct time step call subroutine 'transports' for 
     225     !!               each section to sum the transports. 
     226     !!             Each nn_dctwri time step: 
     227     !!               Divide the arrays by the number of summations to gain 
     228     !!               an average value 
     229     !!               Call dia_dct_sum to sum relevant grid boxes to obtain 
     230     !!               totals for each class (density, depth, temp or sal) 
     231     !!               Call dia_dct_wri to write the transports into file 
     232     !!               Reinitialise all relevant arrays to zero 
    163233     !!--------------------------------------------------------------------- 
    164234     !! * Arguments 
     
    167237     !! * Local variables 
    168238     INTEGER             :: jsec,            &! loop on sections 
    169                             iost,            &! error for opening fileout 
    170                             itotal            ! nb_sec_max*nb_type_class*nb_class_max 
     239                            itotal            ! nb_sec_max*nb_type*nb_class_max 
    171240     LOGICAL             :: lldebug =.FALSE.  ! debug a section   
    172      CHARACTER(len=160)  :: clfileout         ! fileout name 
    173241 
    174242      
    175      INTEGER , DIMENSION(1)             :: ish   ! tmp array for mpp_sum 
    176      INTEGER , DIMENSION(3)             :: ish2  !   " 
    177      REAL(wp), POINTER, DIMENSION(:)    :: zwork !   "  
    178      REAL(wp), POINTER, DIMENSION(:,:,:):: zsum  !   " 
     243     INTEGER , DIMENSION(1)             :: ish      ! tmp array for mpp_sum 
     244     INTEGER , DIMENSION(3)             :: ish2     !   " 
     245     REAL(wp), POINTER, DIMENSION(:)    :: zwork    !   "  
     246     REAL(wp), POINTER, DIMENSION(:,:,:):: zsum     !   " 
    179247 
    180248     !!---------------------------------------------------------------------     
     
    182250 
    183251     IF( lk_mpp )THEN 
    184         itotal = nb_sec_max*nb_type_class*nb_class_max 
    185         CALL wrk_alloc( itotal                                , zwork )  
    186         CALL wrk_alloc( nb_sec_max,nb_type_class,nb_class_max , zsum  ) 
     252        itotal = nb_sec_max*nb_type*nb_class_max 
     253        CALL wrk_alloc( itotal                          , zwork )  
     254        CALL wrk_alloc( nb_sec_max,nb_type,nb_class_max , zsum  ) 
    187255     ENDIF     
    188256  
     257     zwork(:) = 0.0 
     258     zsum(:,:,:) = 0.0 
     259 
    189260     IF( lwp .AND. kt==nit000+nn_dct-1 ) THEN 
    190261         WRITE(numout,*) " " 
     
    194265     ENDIF 
    195266 
    196   
    197      ! Compute transport and write only at nn_dctwri 
    198      IF( MOD(kt,nn_dct)==0 ) THEN  
     267     IF ( MOD(kt,nn_dct)==0 .or. &               ! compute transport every nn_dct time steps 
     268         (ln_NOOS .and. kt==nn_it000 ) )  THEN   ! also include first time step when calculating NOOS 25 hour averages 
    199269 
    200270        DO jsec=1,nb_sec 
    201271 
    202            !debug this section computing ? 
    203272           lldebug=.FALSE. 
    204273           IF( (jsec==nn_secdebug .OR. nn_secdebug==-1) .AND.  kt==nit000+nn_dct-1 .AND. lwp ) lldebug=.TRUE.  
    205274 
    206275           !Compute transport through section   
    207            CALL transport(secs(jsec),lldebug)  
     276           CALL transport(secs(jsec),lldebug,jsec)  
    208277 
    209278        ENDDO 
     
    211280        IF( MOD(kt,nn_dctwri)==0 )THEN 
    212281 
    213            IF( lwp .AND. kt==nit000+nn_dctwri-1 )WRITE(numout,*)"      diadct: write at kt = ",kt          
     282           IF( lwp .AND. kt==nit000+nn_dctwri-1 )WRITE(numout,*)"      diadct: average and write at kt = ",kt          
    214283   
     284           !! divide arrays by nn_dctwri/nn_dct to obtain average 
     285           transports_3d(:,:,:,:)=transports_3d(:,:,:,:)/(nn_dctwri/nn_dct) 
     286           transports_2d(:,:,:)  =transports_2d(:,:,:)  /(nn_dctwri/nn_dct) 
     287 
     288           ! Sum over each class 
     289           DO jsec=1,nb_sec 
     290              CALL dia_dct_sum(secs(jsec),jsec) 
     291           ENDDO 
     292  
    215293           !Sum on all procs  
    216294           IF( lk_mpp )THEN 
    217               ish(1)  =  nb_sec_max*nb_type_class*nb_class_max  
    218               ish2    = (/nb_sec_max,nb_type_class,nb_class_max/) 
     295              zsum(:,:,:)=0.0_wp 
     296              ish(1)  =  nb_sec_max*nb_type*nb_class_max  
     297              ish2    = (/nb_sec_max,nb_type,nb_class_max/) 
    219298              DO jsec=1,nb_sec ; zsum(jsec,:,:) = secs(jsec)%transport(:,:) ; ENDDO 
    220299              zwork(:)= RESHAPE(zsum(:,:,:), ish ) 
     
    227306           DO jsec=1,nb_sec 
    228307 
    229               IF( lwp )CALL dia_dct_wri(kt,jsec,secs(jsec)) 
     308              IF( lwp .and. .not. ln_NOOS )CALL dia_dct_wri(kt,jsec,secs(jsec)) 
     309              IF( lwp .and.       ln_NOOS )CALL dia_dct_wri_NOOS(kt,jsec,secs(jsec))   ! use NOOS specific formatting 
    230310             
    231311              !nullify transports values after writing 
     312              transports_3d(:,jsec,:,:)=0.0 
     313              transports_2d(:,jsec,:)=0.0 
    232314              secs(jsec)%transport(:,:)=0.   
     315              IF ( ln_NOOS ) CALL transport(secs(jsec),lldebug,jsec)  ! reinitialise for next 25 hour instantaneous average (overlapping values) 
    233316 
    234317           ENDDO 
     
    238321     ENDIF 
    239322 
     323     IF ( MOD(kt,nn_dct_h)==0 ) THEN            ! compute transport every nn_dct_h time steps 
     324 
     325        DO jsec=1,nb_sec 
     326 
     327           lldebug=.FALSE. 
     328           IF( (jsec==nn_secdebug .OR. nn_secdebug==-1) .AND.  kt==nit000+nn_dct_h-1 .AND. lwp ) lldebug=.TRUE.  
     329 
     330           !Compute transport through section   
     331           CALL transport_h(secs(jsec),lldebug,jsec)  
     332 
     333        ENDDO 
     334              
     335        IF( MOD(kt,nn_dctwri_h)==0 )THEN 
     336 
     337           IF( lwp .AND. kt==nit000+nn_dctwri_h-1 )WRITE(numout,*)"      diadct: average and write hourly files at kt = ",kt          
     338   
     339           !! divide arrays by nn_dctwri/nn_dct to obtain average 
     340           transports_3d_h(:,:,:,:)=transports_3d_h(:,:,:,:)/(nn_dctwri_h/nn_dct_h) 
     341           transports_2d_h(:,:,:)  =transports_2d_h(:,:,:)  /(nn_dctwri_h/nn_dct_h) 
     342 
     343           ! Sum over each class 
     344           DO jsec=1,nb_sec 
     345              CALL dia_dct_sum_h(secs(jsec),jsec) 
     346           ENDDO 
     347  
     348           !Sum on all procs  
     349          IF( lk_mpp )THEN 
     350              ish(1)  =  nb_sec_max*nb_type*nb_class_max  
     351              ish2    = (/nb_sec_max,nb_type,nb_class_max/) 
     352              DO jsec=1,nb_sec ; zsum(jsec,:,:) = secs(jsec)%transport_h(:,:) ; ENDDO 
     353              zwork(:)= RESHAPE(zsum(:,:,:), ish ) 
     354              CALL mpp_sum(zwork, ish(1)) 
     355              zsum(:,:,:)= RESHAPE(zwork,ish2) 
     356              DO jsec=1,nb_sec ; secs(jsec)%transport_h(:,:) = zsum(jsec,:,:) ; ENDDO 
     357           ENDIF 
     358 
     359           !Write the transport 
     360           DO jsec=1,nb_sec 
     361 
     362              IF( lwp .and.       ln_NOOS )CALL dia_dct_wri_NOOS_h(kt/nn_dctwri_h,jsec,secs(jsec))   ! use NOOS specific formatting 
     363             
     364              !nullify transports values after writing 
     365              transports_3d_h(:,jsec,:,:)=0.0 
     366              transports_2d_h(:,jsec,:)=0.0 
     367              secs(jsec)%transport_h(:,:)=0.   
     368              IF ( ln_NOOS ) CALL transport_h(secs(jsec),lldebug,jsec)  ! reinitialise for next 25 hour instantaneous average (overlapping values) 
     369 
     370           ENDDO 
     371 
     372        ENDIF  
     373 
     374     ENDIF     
     375 
    240376     IF( lk_mpp )THEN 
    241         itotal = nb_sec_max*nb_type_class*nb_class_max 
    242         CALL wrk_dealloc( itotal                                , zwork )  
    243         CALL wrk_dealloc( nb_sec_max,nb_type_class,nb_class_max , zsum  ) 
     377        itotal = nb_sec_max*nb_type*nb_class_max 
     378        CALL wrk_dealloc( itotal                          , zwork )  
     379        CALL wrk_dealloc( nb_sec_max,nb_type,nb_class_max , zsum  ) 
    244380     ENDIF     
    245381 
     
    299435        secs(jsec)%zlay         = 99._wp          
    300436        secs(jsec)%transport    =  0._wp   ; secs(jsec)%nb_point       = 0 
     437        secs(jsec)%transport_h  =  0._wp   ; secs(jsec)%nb_point       = 0 
    301438 
    302439        !read section's number / name / computing choices / classes / slopeSection / points number 
     
    331468 
    332469            WRITE(numout,*)       "   Section name :                       ",TRIM(secs(jsec)%name) 
    333             WRITE(numout,*)       "      Compute heat and salt transport ? ",secs(jsec)%llstrpond 
     470            WRITE(numout,*)       "      Compute temperature and salinity transports ? ",secs(jsec)%llstrpond 
    334471            WRITE(numout,*)       "      Compute ice transport ?           ",secs(jsec)%ll_ice_section 
    335472            WRITE(numout,*)       "      Section crosses date-line ?       ",secs(jsec)%ll_date_line 
     
    362499              WRITE(numout,*)"      List of points in global domain:" 
    363500              DO jpt=1,iptglo 
    364                  WRITE(numout,*)'        # I J ',jpt,coordtemp(jpt) 
     501                 WRITE(numout,*)'        # I J ',jpt,coordtemp(jpt),directemp(jpt) 
    365502              ENDDO                   
    366503           ENDIF 
     
    403540 
    404541              IF(jsec==nn_secdebug .AND. secs(jsec)%nb_point .NE. 0)THEN 
    405               WRITE(narea+200,*)'avant secs(jsec)%nb_point iptloc ',secs(jsec)%nb_point,iptloc 
    406542              DO jpt = 1,iptloc 
    407543                 iiglo = secs(jsec)%listPoint(jpt)%I + jpizoom - 1 + nimpp - 1 
    408544                 ijglo = secs(jsec)%listPoint(jpt)%J + jpjzoom - 1 + njmpp - 1 
    409                  WRITE(narea+200,*)'avant # I J : ',iiglo,ijglo 
    410545              ENDDO 
    411546              ENDIF 
     
    421556           ENDIF 
    422557           IF(jsec==nn_secdebug .AND. secs(jsec)%nb_point .NE. 0)THEN 
    423               WRITE(narea+200,*)'apres secs(jsec)%nb_point iptloc ',secs(jsec)%nb_point,iptloc 
    424558              DO jpt = 1,secs(jsec)%nb_point 
    425559                 iiglo = secs(jsec)%listPoint(jpt)%I + jpizoom - 1 + nimpp - 1 
    426560                 ijglo = secs(jsec)%listPoint(jpt)%J + jpjzoom - 1 + njmpp - 1 
    427                  WRITE(narea+200,*)'apres # I J : ',iiglo,ijglo 
    428561              ENDDO 
    429562           ENDIF 
     
    534667     CALL wrk_dealloc(    nb_point_max, idirec ) 
    535668     CALL wrk_dealloc( 2, nb_point_max, icoord ) 
     669 
    536670  END SUBROUTINE removepoints 
    537671 
    538   SUBROUTINE transport(sec,ld_debug) 
     672  SUBROUTINE transport(sec,ld_debug,jsec) 
    539673     !!------------------------------------------------------------------------------------------- 
    540674     !!                     ***  ROUTINE transport  *** 
    541675     !! 
    542      !!  ** Purpose : Compute the transport through a section 
    543      !! 
    544      !!  ** Method  :Transport through a given section is equal to the sum of transports 
    545      !!              computed on each proc. 
    546      !!              On each proc,transport is equal to the sum of transport computed through 
    547      !!               segments linking each point of sec%listPoint  with the next one.    
    548      !! 
    549      !!              !BE carefull :           
    550      !!              one section is a sum of segments 
    551      !!              one segment is defined by 2 consectuve points in sec%listPoint 
    552      !!              all points of sec%listPoint are positioned on the F-point of the cell.  
     676     !!  Purpose ::  Compute the transport for each point in a section 
     677     !! 
     678     !!  Method  ::  Loop over each segment, and each vertical level and add the transport 
     679     !!              Be aware :           
     680     !!              One section is a sum of segments 
     681     !!              One segment is defined by 2 consecutive points in sec%listPoint 
     682     !!              All points of sec%listPoint are positioned on the F-point of the cell 
    553683     !!  
    554      !!              There are several loops:                  
    555      !!              loop on the density/temperature/salinity/level classes 
     684     !!              There are two loops:                  
    556685     !!              loop on the segment between 2 nodes 
    557686     !!              loop on the level jk 
    558      !!              test on the density/temperature/salinity/level 
    559      !! 
    560      !! ** Output: sec%transport: volume/mass/ice/heat/salt transport in the 2 directions 
    561      !! 
     687     !! 
     688     !! ** Output: Arrays containing the volume,density,salinity,temperature etc 
     689     !!            transports for each point in a section, summed over each nn_dct. 
    562690     !! 
    563691     !!------------------------------------------------------------------------------------------- 
     
    565693     TYPE(SECTION),INTENT(INOUT) :: sec 
    566694     LOGICAL      ,INTENT(IN)    :: ld_debug 
     695     INTEGER      ,INTENT(IN)    :: jsec        ! numeric identifier of section 
    567696     
    568697     !! * Local variables 
    569      INTEGER             :: jk,jseg,jclass,   &!loop on level/segment/classes  
    570                             isgnu  , isgnv     ! 
    571      INTEGER :: ii, ij ! local integer 
    572      REAL(wp):: zumid        , zvmid        ,&!U/V velocity on a cell segment 
    573                 zumid_ice    , zvmid_ice    ,&!U/V ice velocity 
    574                 zTnorm                      ,&!transport of velocity through one cell's sides 
    575                 ztransp1     , ztransp2     ,&!total        transport in directions 1 and 2 
    576                 ztemp1       , ztemp2       ,&!temperature  transport     " 
    577                 zrhoi1       , zrhoi2       ,&!mass         transport     " 
    578                 zrhop1       , zrhop2       ,&!mass         transport     " 
    579                 zsal1        , zsal2        ,&!salinity     transport     " 
    580                 zice_vol_pos , zice_vol_neg ,&!volume  ice  transport     " 
    581                 zice_surf_pos, zice_surf_neg  !surface ice  transport     " 
     698     INTEGER             :: jk,jseg,jclass,    &!loop on level/segment/classes  
     699                            isgnu  , isgnv      ! 
     700     REAL(wp):: zumid        , zvmid        ,  &!U/V velocity on a cell segment 
     701                zumid_ice    , zvmid_ice    ,  &!U/V ice velocity 
     702                zTnorm                          !transport of velocity through one cell's sides 
    582703     REAL(wp):: ztn, zsn, zrhoi, zrhop, zsshn, zfsdep ! temperature/salinity/ssh/potential density /depth at u/v point 
    583704 
    584705     TYPE(POINT_SECTION) :: k 
    585      REAL(wp), POINTER, DIMENSION(:,:):: zsum ! 2D work array 
    586706     !!-------------------------------------------------------- 
    587      CALL wrk_alloc( nb_type_class , nb_class_max , zsum   ) 
    588707 
    589708     IF( ld_debug )WRITE(numout,*)'      Compute transport' 
    590  
    591      !----------------! 
    592      ! INITIALIZATION ! 
    593      !----------------! 
    594      zsum    = 0._wp 
    595      zice_surf_neg = 0._wp ; zice_surf_pos = 0._wp 
    596      zice_vol_pos  = 0._wp ; zice_vol_neg  = 0._wp 
    597709 
    598710     !---------------------------! 
     
    626738        ELSE                                ; isgnv =  1 
    627739        ENDIF 
    628  
    629         IF( ld_debug )write(numout,*)"isgnu isgnv ",isgnu,isgnv 
     740        IF( sec%slopeSection .GE. 9999. )     isgnv =  1 
     741 
     742        IF( ld_debug )write(numout,*)"sec%slopeSection isgnu isgnv ",sec%slopeSection,isgnu,isgnv 
    630743 
    631744        !--------------------------------------! 
     
    670783           END SELECT 
    671784 
    672            !------------------------------- 
    673            !  LOOP ON THE DENSITY CLASSES | 
    674            !------------------------------- 
    675            !The computation is made for each density class 
    676            DO jclass=1,MAX(1,sec%nb_class-1) 
    677  
    678               ztransp1=0._wp ; zrhoi1=0._wp ; zrhop1=0._wp ; ztemp1=0._wp ;zsal1=0._wp 
    679               ztransp2=0._wp ; zrhoi2=0._wp ; zrhop2=0._wp ; ztemp2=0._wp ;zsal2=0._wp 
    680      
    681               !---------------------------| 
    682               !     LOOP ON THE LEVEL     | 
    683               !---------------------------| 
    684               !Sum of the transport on the vertical  
    685               DO jk=1,jpk 
    686                      
    687  
    688                  ! compute temparature, salinity, insitu & potential density, ssh and depth at U/V point 
    689                  SELECT CASE( sec%direction(jseg) ) 
    690                  CASE(0,1) 
    691                     ztn   = interp(k%I,k%J,jk,'V',tsn(:,:,:,jp_tem) ) 
    692                     zsn   = interp(k%I,k%J,jk,'V',tsn(:,:,:,jp_sal) ) 
    693                     zrhop = interp(k%I,k%J,jk,'V',rhop) 
    694                     zrhoi = interp(k%I,k%J,jk,'V',rhd*rau0+rau0) 
    695                     zsshn =  0.5*( sshn(k%I,k%J)    + sshn(k%I,k%J+1)    ) * vmask(k%I,k%J,1) 
    696                  CASE(2,3) 
    697                     ztn   = interp(k%I,k%J,jk,'U',tsn(:,:,:,jp_tem) ) 
    698                     zsn   = interp(k%I,k%J,jk,'U',tsn(:,:,:,jp_sal) ) 
    699                     zrhop = interp(k%I,k%J,jk,'U',rhop) 
    700                     zrhoi = interp(k%I,k%J,jk,'U',rhd*rau0+rau0) 
    701                     zsshn =  0.5*( sshn(k%I,k%J)    + sshn(k%I+1,k%J)    ) * umask(k%I,k%J,1)  
    702                  END SELECT 
    703  
    704                  zfsdep= gdept(k%I,k%J,jk) 
    705   
    706                  !----------------------------------------------! 
    707                  !TEST ON THE DENSITY/SALINITY/TEMPERATURE/LEVEL!  
    708                  !----------------------------------------------! 
    709   
    710                  IF ( (    ((( zrhop .GE. (sec%zsigp(jclass)+1000.  )) .AND.    & 
    711                            (   zrhop .LE. (sec%zsigp(jclass+1)+1000. ))) .OR.    & 
    712                            ( sec%zsigp(jclass) .EQ. 99.)) .AND.                 & 
    713                            ((( zrhoi .GE. (sec%zsigi(jclass) + 1000.  )) .AND.    & 
    714                            (   zrhoi .LE. (sec%zsigi(jclass+1)+1000. ))) .OR.    & 
    715                            ( sec%zsigi(jclass) .EQ. 99.)) .AND.                 & 
    716                            ((( zsn .GT. sec%zsal(jclass)) .AND.                & 
    717                            (   zsn .LE. sec%zsal(jclass+1))) .OR.              & 
    718                            ( sec%zsal(jclass) .EQ. 99.)) .AND.                 & 
    719                            ((( ztn .GE. sec%ztem(jclass)) .AND.                & 
    720                            (   ztn .LE. sec%ztem(jclass+1))) .OR.              & 
    721                            ( sec%ztem(jclass) .EQ.99.)) .AND.                  & 
    722                            ((( zfsdep .GE. sec%zlay(jclass)) .AND.            & 
    723                            (   zfsdep .LE. sec%zlay(jclass+1))) .OR.          & 
    724                            ( sec%zlay(jclass) .EQ. 99. ))))   THEN 
    725  
    726  
    727                     !compute velocity with the correct direction 
    728                     SELECT CASE( sec%direction(jseg) ) 
    729                     CASE(0,1)   
    730                        zumid=0. 
    731                        zvmid=isgnv*vn(k%I,k%J,jk)*vmask(k%I,k%J,jk) 
    732                     CASE(2,3) 
    733                        zumid=isgnu*un(k%I,k%J,jk)*umask(k%I,k%J,jk) 
    734                        zvmid=0. 
    735                     END SELECT 
    736  
    737                     !velocity* cell's length * cell's thickness 
    738                     zTnorm=zumid*e2u(k%I,k%J)*  fse3u(k%I,k%J,jk)+     & 
    739                            zvmid*e1v(k%I,k%J)*  fse3v(k%I,k%J,jk) 
     785           !---------------------------| 
     786           !     LOOP ON THE LEVEL     | 
     787           !---------------------------| 
     788           !Sum of the transport on the vertical  
     789           DO jk=1,mbathy(k%I,k%J) 
     790 
     791              ! compute temparature, salinity, insitu & potential density, ssh and depth at U/V point 
     792              SELECT CASE( sec%direction(jseg) ) 
     793              CASE(0,1) 
     794                 ztn   = interp(k%I,k%J,jk,'V',tsn(:,:,:,jp_tem) ) 
     795                 zsn   = interp(k%I,k%J,jk,'V',tsn(:,:,:,jp_sal) ) 
     796                 zrhop = interp(k%I,k%J,jk,'V',rhop) 
     797                 zrhoi = interp(k%I,k%J,jk,'V',rhd*rau0+rau0) 
     798                 zsshn =  0.5*( sshn(k%I,k%J)    + sshn(k%I,k%J+1)    ) * vmask(k%I,k%J,1) 
     799              CASE(2,3) 
     800                 ztn   = interp(k%I,k%J,jk,'U',tsn(:,:,:,jp_tem) ) 
     801                 zsn   = interp(k%I,k%J,jk,'U',tsn(:,:,:,jp_sal) ) 
     802                 zrhop = interp(k%I,k%J,jk,'U',rhop) 
     803                 zrhoi = interp(k%I,k%J,jk,'U',rhd*rau0+rau0) 
     804                 zsshn =  0.5*( sshn(k%I,k%J)    + sshn(k%I+1,k%J)    ) * umask(k%I,k%J,1)  
     805              END SELECT 
     806 
     807              zfsdep= gdept(k%I,k%J,jk) 
     808  
     809              !compute velocity with the correct direction 
     810              SELECT CASE( sec%direction(jseg) ) 
     811              CASE(0,1)   
     812                 zumid=0. 
     813                 zvmid=isgnv*vn(k%I,k%J,jk)*vmask(k%I,k%J,jk) 
     814              CASE(2,3) 
     815                 zumid=isgnu*un(k%I,k%J,jk)*umask(k%I,k%J,jk) 
     816                 zvmid=0. 
     817              END SELECT 
     818 
     819              !zTnorm=transport through one cell; 
     820              !velocity* cell's length * cell's thickness 
     821              zTnorm=zumid*e2u(k%I,k%J)*  fse3u(k%I,k%J,jk)+     & 
     822                     zvmid*e1v(k%I,k%J)*  fse3v(k%I,k%J,jk) 
    740823 
    741824#if ! defined key_vvl 
    742                     !add transport due to free surface 
    743                     IF( jk==1 )THEN 
    744                        zTnorm = zTnorm + zumid* e2u(k%I,k%J) * zsshn * umask(k%I,k%J,jk) + & 
    745                                          zvmid* e1v(k%I,k%J) * zsshn * vmask(k%I,k%J,jk) 
    746                     ENDIF 
     825              !add transport due to free surface 
     826              IF( jk==1 )THEN 
     827                 zTnorm = zTnorm + zumid* e2u(k%I,k%J) * zsshn * umask(k%I,k%J,jk) + & 
     828                                   zvmid* e1v(k%I,k%J) * zsshn * vmask(k%I,k%J,jk) 
     829              ENDIF 
    747830#endif 
    748                     !COMPUTE TRANSPORT  
    749                     !zTnorm=transport through one cell for one class 
    750                     !ztransp1 or ztransp2=transport through one cell i 
    751                     !                     for one class for one direction 
    752                     IF( zTnorm .GE. 0 )THEN 
    753  
    754                        ztransp1=zTnorm+ztransp1 
    755   
    756                        IF ( sec%llstrpond ) THEN 
    757                           ztemp1 = ztemp1  + zTnorm * ztn  
    758                           zsal1  = zsal1   + zTnorm * zsn 
    759                           zrhoi1 = zrhoi1  + zTnorm * zrhoi 
    760                           zrhop1 = zrhop1  + zTnorm * zrhop 
    761                        ENDIF 
    762  
    763                     ELSE 
    764  
    765                        ztransp2=(zTnorm)+ztransp2 
    766  
    767                        IF ( sec%llstrpond ) THEN 
    768                           ztemp2 = ztemp2  + zTnorm * ztn  
    769                           zsal2  = zsal2   + zTnorm * zsn 
    770                           zrhoi2 = zrhoi2  + zTnorm * zrhoi 
    771                           zrhop2 = zrhop2  + zTnorm * zrhop 
    772                        ENDIF 
    773                     ENDIF 
    774   
    775              
    776                  ENDIF ! end of density test 
    777               ENDDO!end of loop on the level 
    778  
    779               !ZSUM=TRANSPORT FOR EACH CLASSES FOR THE  DIRECTIONS 
    780               !--------------------------------------------------- 
    781               zsum(1,jclass)     = zsum(1,jclass)+ztransp1 
    782               zsum(2,jclass)     = zsum(2,jclass)+ztransp2 
    783               IF( sec%llstrpond )THEN 
    784                  zsum(3 ,jclass) = zsum( 3,jclass)+zrhoi1 
    785                  zsum(4 ,jclass) = zsum( 4,jclass)+zrhoi2 
    786                  zsum(5 ,jclass) = zsum( 5,jclass)+zrhop1 
    787                  zsum(6 ,jclass) = zsum( 6,jclass)+zrhop2 
    788                  zsum(7 ,jclass) = zsum( 7,jclass)+ztemp1 
    789                  zsum(8 ,jclass) = zsum( 8,jclass)+ztemp2 
    790                  zsum(9 ,jclass) = zsum( 9,jclass)+zsal1 
    791                  zsum(10,jclass) = zsum(10,jclass)+zsal2 
     831              !COMPUTE TRANSPORT  
     832 
     833              transports_3d(1,jsec,jseg,jk) = transports_3d(1,jsec,jseg,jk) + zTnorm 
     834  
     835              IF ( sec%llstrpond ) THEN 
     836                 transports_3d(2,jsec,jseg,jk) = transports_3d(2,jsec,jseg,jk)  + zTnorm * zrhoi 
     837                 transports_3d(3,jsec,jseg,jk) = transports_3d(3,jsec,jseg,jk)  + zTnorm * zrhop 
     838                 transports_3d(4,jsec,jseg,jk) = transports_3d(4,jsec,jseg,jk)  + zTnorm * 4.e+3_wp * (ztn+273.15) * 1026._wp 
     839                 transports_3d(5,jsec,jseg,jk) = transports_3d(5,jsec,jseg,jk)  + zTnorm * 0.001 * zsn * 1026._wp 
    792840              ENDIF 
    793     
    794            ENDDO !end of loop on the density classes 
     841 
     842           ENDDO !end of loop on the level 
    795843 
    796844#if defined key_lim2 || defined key_lim3 
     
    816864              zTnorm=zumid_ice*e2u(k%I,k%J)+zvmid_ice*e1v(k%I,k%J) 
    817865    
    818               IF( zTnorm .GE. 0)THEN 
    819                  zice_vol_pos = (zTnorm)*   & 
    820                                       (1.0 - frld(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J))  & 
    821                                      *(hsnif(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J) +  & 
    822                                        hicif(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J)) & 
    823                                       +zice_vol_pos 
    824                  zice_surf_pos = (zTnorm)*   & 
    825                                        (1.0 -  frld(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J))  & 
    826                                       +zice_surf_pos 
    827               ELSE 
    828                  zice_vol_neg=(zTnorm)*   & 
     866              transports_2d(1,jsec,jseg) = transports_2d(1,jsec,jseg) + (zTnorm)*   & 
    829867                                   (1.0 - frld(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J))  & 
    830868                                  *(hsnif(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J) +  & 
    831869                                    hicif(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J)) & 
    832                                   +zice_vol_neg 
    833                  zice_surf_neg=(zTnorm)*   & 
    834                                     (1.0 - frld(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J))  & 
    835                                      +zice_surf_neg 
    836               ENDIF 
    837     
    838               zsum(11,1) = zsum(11,1)+zice_vol_pos 
    839               zsum(12,1) = zsum(12,1)+zice_vol_neg 
    840               zsum(13,1) = zsum(13,1)+zice_surf_pos 
    841               zsum(14,1) = zsum(14,1)+zice_surf_neg 
     870                                   +zice_vol_pos 
     871              transports_2d(2,jsec,jseg) = transports_2d(2,jsec,jseg) + (zTnorm)*   & 
     872                                    (1.0 -  frld(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J))  & 
     873                                   +zice_surf_pos 
    842874    
    843875           ENDIF !end of ice case 
     
    846878        ENDDO !end of loop on the segment 
    847879 
    848  
    849      ELSE  !if sec%nb_point =0 
    850         zsum(1:2,:)=0. 
    851         IF (sec%llstrpond) zsum(3:10,:)=0. 
    852         zsum( 11:14,:)=0. 
    853880     ENDIF   !end of sec%nb_point =0 case 
    854  
    855      !-------------------------------| 
    856      !FINISH COMPUTING TRANSPORTS    | 
    857      !-------------------------------| 
    858      DO jclass=1,MAX(1,sec%nb_class-1) 
    859         sec%transport(1,jclass)=sec%transport(1,jclass)+zsum(1,jclass)*1.E-6 
    860         sec%transport(2,jclass)=sec%transport(2,jclass)+zsum(2,jclass)*1.E-6 
    861         IF( sec%llstrpond ) THEN 
    862            IF( zsum(1,jclass) .NE. 0._wp ) THEN 
    863               sec%transport( 3,jclass) = sec%transport( 3,jclass) + zsum( 3,jclass)/zsum(1,jclass) 
    864               sec%transport( 5,jclass) = sec%transport( 5,jclass) + zsum( 5,jclass)/zsum(1,jclass) 
    865               sec%transport( 7,jclass) = sec%transport( 7,jclass) + zsum( 7,jclass) 
    866               sec%transport( 9,jclass) = sec%transport( 9,jclass) + zsum( 9,jclass) 
    867            ENDIF 
    868            IF( zsum(2,jclass) .NE. 0._wp )THEN 
    869               sec%transport( 4,jclass) = sec%transport( 4,jclass) + zsum( 4,jclass)/zsum(2,jclass) 
    870               sec%transport( 6,jclass) = sec%transport( 6,jclass) + zsum( 6,jclass)/zsum(2,jclass) 
    871               sec%transport( 8,jclass) = sec%transport( 8,jclass) + zsum( 8,jclass) 
    872               sec%transport(10,jclass) = sec%transport(10,jclass) + zsum(10,jclass) 
    873            ENDIF 
    874         ELSE 
    875            sec%transport( 3,jclass) = 0._wp 
    876            sec%transport( 4,jclass) = 0._wp 
    877            sec%transport( 5,jclass) = 0._wp 
    878            sec%transport( 6,jclass) = 0._wp 
    879            sec%transport( 7,jclass) = 0._wp 
    880            sec%transport( 8,jclass) = 0._wp 
    881            sec%transport(10,jclass) = 0._wp 
    882         ENDIF 
    883      ENDDO    
    884  
    885      IF( sec%ll_ice_section ) THEN 
    886         sec%transport( 9,1)=sec%transport( 9,1)+zsum( 9,1)*1.E-6 
    887         sec%transport(10,1)=sec%transport(10,1)+zsum(10,1)*1.E-6 
    888         sec%transport(11,1)=sec%transport(11,1)+zsum(11,1)*1.E-6 
    889         sec%transport(12,1)=sec%transport(12,1)+zsum(12,1)*1.E-6 
    890      ENDIF 
    891  
    892      CALL wrk_dealloc( nb_type_class , nb_class_max , zsum   ) 
    893881     ! 
    894882  END SUBROUTINE transport 
    895    
     883 
     884  SUBROUTINE transport_h(sec,ld_debug,jsec) 
     885     !!------------------------------------------------------------------------------------------- 
     886     !!                     ***  ROUTINE hourly_transport  *** 
     887     !! 
     888     !!              Exactly as routine transport but for data summed at 
     889     !!              each time step and averaged each hour 
     890     !! 
     891     !!  Purpose ::  Compute the transport for each point in a section 
     892     !! 
     893     !!  Method  ::  Loop over each segment, and each vertical level and add the transport 
     894     !!              Be aware :           
     895     !!              One section is a sum of segments 
     896     !!              One segment is defined by 2 consecutive points in sec%listPoint 
     897     !!              All points of sec%listPoint are positioned on the F-point of the cell 
     898     !!  
     899     !!              There are two loops:                  
     900     !!              loop on the segment between 2 nodes 
     901     !!              loop on the level jk 
     902     !! 
     903     !! ** Output: Arrays containing the volume,density,salinity,temperature etc 
     904     !!            transports for each point in a section, summed over each nn_dct. 
     905     !! 
     906     !!------------------------------------------------------------------------------------------- 
     907     !! * Arguments 
     908     TYPE(SECTION),INTENT(INOUT) :: sec 
     909     LOGICAL      ,INTENT(IN)    :: ld_debug 
     910     INTEGER      ,INTENT(IN)    :: jsec        ! numeric identifier of section 
     911     
     912     !! * Local variables 
     913     INTEGER             :: jk,jseg,jclass,    &!loop on level/segment/classes  
     914                            isgnu  , isgnv      ! 
     915     REAL(wp):: zumid        , zvmid        ,  &!U/V velocity on a cell segment 
     916                zumid_ice    , zvmid_ice    ,  &!U/V ice velocity 
     917                zTnorm                          !transport of velocity through one cell's sides 
     918     REAL(wp):: ztn, zsn, zrhoi, zrhop, zsshn, zfsdep ! temperature/salinity/ssh/potential density /depth at u/v point 
     919 
     920     TYPE(POINT_SECTION) :: k 
     921     !!-------------------------------------------------------- 
     922 
     923     IF( ld_debug )WRITE(numout,*)'      Compute transport' 
     924 
     925     !---------------------------! 
     926     !  COMPUTE TRANSPORT        ! 
     927     !---------------------------! 
     928     IF(sec%nb_point .NE. 0)THEN    
     929 
     930        !---------------------------------------------------------------------------------------------------- 
     931        !Compute sign for velocities: 
     932        ! 
     933        !convention: 
     934        !   non horizontal section: direction + is toward left hand of section 
     935        !       horizontal section: direction + is toward north of section 
     936        ! 
     937        ! 
     938        !       slopeSection < 0     slopeSection > 0       slopeSection=inf            slopeSection=0 
     939        !       ----------------      -----------------     ---------------             -------------- 
     940        ! 
     941        !   isgnv=1         direction +       
     942        !  ______         _____             ______                                                    
     943        !        |           //|            |                  |                         direction +    
     944        !        | isgnu=1  // |            |isgnu=1           |isgnu=1                     /|\ 
     945        !        |_______  //         ______|    \\            | ---\                        | 
     946        !               |             | isgnv=-1  \\ |         | ---/ direction +       ____________ 
     947        !               |             |          __\\|         |                     
     948        !               |             |     direction +        |                      isgnv=1                                  
     949        !                                                       
     950        !---------------------------------------------------------------------------------------------------- 
     951        isgnu = 1 
     952        IF( sec%slopeSection .GT. 0 ) THEN  ; isgnv = -1  
     953        ELSE                                ; isgnv =  1 
     954        ENDIF 
     955 
     956        IF( ld_debug )write(numout,*)"isgnu isgnv ",isgnu,isgnv 
     957 
     958        !--------------------------------------! 
     959        ! LOOP ON THE SEGMENT BETWEEN 2 NODES  ! 
     960        !--------------------------------------! 
     961        DO jseg=1,MAX(sec%nb_point-1,0) 
     962               
     963           !------------------------------------------------------------------------------------------- 
     964           ! Select the appropriate coordinate for computing the velocity of the segment 
     965           ! 
     966           !                      CASE(0)                                    Case (2) 
     967           !                      -------                                    -------- 
     968           !  listPoint(jseg)                 listPoint(jseg+1)       listPoint(jseg)  F(i,j)       
     969           !      F(i,j)----------V(i+1,j)-------F(i+1,j)                               | 
     970           !                                                                            | 
     971           !                                                                            | 
     972           !                                                                            | 
     973           !                      Case (3)                                            U(i,j) 
     974           !                      --------                                              | 
     975           !                                                                            | 
     976           !  listPoint(jseg+1) F(i,j+1)                                                | 
     977           !                        |                                                   | 
     978           !                        |                                                   | 
     979           !                        |                                 listPoint(jseg+1) F(i,j-1) 
     980           !                        |                                             
     981           !                        |                                             
     982           !                     U(i,j+1)                                             
     983           !                        |                                       Case(1)      
     984           !                        |                                       ------       
     985           !                        |                                             
     986           !                        |                 listPoint(jseg+1)             listPoint(jseg)                            
     987           !                        |                 F(i-1,j)-----------V(i,j) -------f(jseg)                            
     988           ! listPoint(jseg)     F(i,j) 
     989           !  
     990           !------------------------------------------------------------------------------------------- 
     991 
     992           SELECT CASE( sec%direction(jseg) ) 
     993           CASE(0)  ;   k = sec%listPoint(jseg) 
     994           CASE(1)  ;   k = POINT_SECTION(sec%listPoint(jseg)%I+1,sec%listPoint(jseg)%J) 
     995           CASE(2)  ;   k = sec%listPoint(jseg) 
     996           CASE(3)  ;   k = POINT_SECTION(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J+1) 
     997           END SELECT 
     998 
     999           !---------------------------| 
     1000           !     LOOP ON THE LEVEL     | 
     1001           !---------------------------| 
     1002           !Sum of the transport on the vertical  
     1003           DO jk=1,mbathy(k%I,k%J) 
     1004 
     1005              ! compute temparature, salinity, insitu & potential density, ssh and depth at U/V point 
     1006              SELECT CASE( sec%direction(jseg) ) 
     1007              CASE(0,1) 
     1008                 ztn   = interp(k%I,k%J,jk,'V',tsn(:,:,:,jp_tem) ) 
     1009                 zsn   = interp(k%I,k%J,jk,'V',tsn(:,:,:,jp_sal) ) 
     1010                 zrhop = interp(k%I,k%J,jk,'V',rhop) 
     1011                 zrhoi = interp(k%I,k%J,jk,'V',rhd*rau0+rau0) 
     1012                 zsshn =  0.5*( sshn(k%I,k%J)    + sshn(k%I,k%J+1)    ) * vmask(k%I,k%J,1) 
     1013              CASE(2,3) 
     1014                 ztn   = interp(k%I,k%J,jk,'U',tsn(:,:,:,jp_tem) ) 
     1015                 zsn   = interp(k%I,k%J,jk,'U',tsn(:,:,:,jp_sal) ) 
     1016                 zrhop = interp(k%I,k%J,jk,'U',rhop) 
     1017                 zrhoi = interp(k%I,k%J,jk,'U',rhd*rau0+rau0) 
     1018                 zsshn =  0.5*( sshn(k%I,k%J)    + sshn(k%I+1,k%J)    ) * umask(k%I,k%J,1)  
     1019              END SELECT 
     1020 
     1021              zfsdep= gdept(k%I,k%J,jk) 
     1022  
     1023              !compute velocity with the correct direction 
     1024              SELECT CASE( sec%direction(jseg) ) 
     1025              CASE(0,1)   
     1026                 zumid=0. 
     1027                 zvmid=isgnv*vn(k%I,k%J,jk)*vmask(k%I,k%J,jk) 
     1028              CASE(2,3) 
     1029                 zumid=isgnu*un(k%I,k%J,jk)*umask(k%I,k%J,jk) 
     1030                 zvmid=0. 
     1031              END SELECT 
     1032 
     1033              !zTnorm=transport through one cell; 
     1034              !velocity* cell's length * cell's thickness 
     1035              zTnorm=zumid*e2u(k%I,k%J)*  fse3u(k%I,k%J,jk)+     & 
     1036                     zvmid*e1v(k%I,k%J)*  fse3v(k%I,k%J,jk) 
     1037 
     1038#if ! defined key_vvl 
     1039              !add transport due to free surface 
     1040              IF( jk==1 )THEN 
     1041                 zTnorm = zTnorm + zumid* e2u(k%I,k%J) * zsshn * umask(k%I,k%J,jk) + & 
     1042                                   zvmid* e1v(k%I,k%J) * zsshn * vmask(k%I,k%J,jk) 
     1043              ENDIF 
     1044#endif 
     1045              !COMPUTE TRANSPORT  
     1046 
     1047              transports_3d_h(1,jsec,jseg,jk) = transports_3d_h(1,jsec,jseg,jk) + zTnorm 
     1048  
     1049              IF ( sec%llstrpond ) THEN 
     1050                 transports_3d_h(2,jsec,jseg,jk) = transports_3d_h(2,jsec,jseg,jk)  + zTnorm * zrhoi 
     1051                 transports_3d_h(3,jsec,jseg,jk) = transports_3d_h(3,jsec,jseg,jk)  + zTnorm * zrhop 
     1052                 transports_3d_h(4,jsec,jseg,jk) = transports_3d_h(4,jsec,jseg,jk)  + zTnorm * 4.e+3_wp * (ztn+273.15) * 1026._wp 
     1053                 transports_3d_h(5,jsec,jseg,jk) = transports_3d_h(5,jsec,jseg,jk)  + zTnorm * 0.001 * zsn * 1026._wp 
     1054              ENDIF 
     1055 
     1056           ENDDO !end of loop on the level 
     1057 
     1058#if defined key_lim2 || defined key_lim3 
     1059 
     1060           !ICE CASE     
     1061           !------------ 
     1062           IF( sec%ll_ice_section )THEN 
     1063              SELECT CASE (sec%direction(jseg)) 
     1064              CASE(0) 
     1065                 zumid_ice = 0 
     1066                 zvmid_ice =  isgnv*0.5*(v_ice(k%I,k%J+1)+v_ice(k%I+1,k%J+1)) 
     1067              CASE(1) 
     1068                 zumid_ice = 0 
     1069                 zvmid_ice =  isgnv*0.5*(v_ice(k%I,k%J+1)+v_ice(k%I+1,k%J+1)) 
     1070              CASE(2) 
     1071                 zvmid_ice = 0 
     1072                 zumid_ice =  isgnu*0.5*(u_ice(k%I+1,k%J)+u_ice(k%I+1,k%J+1)) 
     1073              CASE(3) 
     1074                 zvmid_ice = 0 
     1075                 zumid_ice =  isgnu*0.5*(u_ice(k%I+1,k%J)+u_ice(k%I+1,k%J+1)) 
     1076              END SELECT 
     1077    
     1078              zTnorm=zumid_ice*e2u(k%I,k%J)+zvmid_ice*e1v(k%I,k%J) 
     1079    
     1080              transports_2d_h(1,jsec,jseg) = transports_2d_h(1,jsec,jseg) + (zTnorm)*   & 
     1081                                   (1.0 - frld(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J))  & 
     1082                                  *(hsnif(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J) +  & 
     1083                                    hicif(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J)) & 
     1084                                   +zice_vol_pos 
     1085              transports_2d_h(2,jsec,jseg) = transports_2d_h(2,jsec,jseg) + (zTnorm)*   & 
     1086                                    (1.0 -  frld(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J))  & 
     1087                                   +zice_surf_pos 
     1088    
     1089           ENDIF !end of ice case 
     1090#endif 
     1091  
     1092        ENDDO !end of loop on the segment 
     1093 
     1094     ENDIF   !end of sec%nb_point =0 case 
     1095     ! 
     1096  END SUBROUTINE transport_h 
     1097  
     1098  SUBROUTINE dia_dct_sum(sec,jsec) 
     1099     !!------------------------------------------------------------- 
     1100     !! Purpose: Average the transport over nn_dctwri time steps  
     1101     !! and sum over the density/salinity/temperature/depth classes 
     1102     !! 
     1103     !! Method:  
     1104     !!           Sum over relevant grid cells to obtain values 
     1105     !!           for each 
     1106     !!              There are several loops:                  
     1107     !!              loop on the segment between 2 nodes 
     1108     !!              loop on the level jk 
     1109     !!              loop on the density/temperature/salinity/level classes 
     1110     !!              test on the density/temperature/salinity/level 
     1111     !! 
     1112     !!  ** Method  :Transport through a given section is equal to the sum of transports 
     1113     !!              computed on each proc. 
     1114     !!              On each proc,transport is equal to the sum of transport computed through 
     1115     !!               segments linking each point of sec%listPoint  with the next one.    
     1116     !! 
     1117     !!------------------------------------------------------------- 
     1118     !! * arguments 
     1119     TYPE(SECTION),INTENT(INOUT) :: sec 
     1120     INTEGER      ,INTENT(IN)    :: jsec        ! numeric identifier of section 
     1121 
     1122     TYPE(POINT_SECTION) :: k 
     1123     INTEGER  :: jk,jseg,jclass                        !loop on level/segment/classes  
     1124     REAL(wp) :: ztn, zsn, zrhoi, zrhop, zsshn, zfsdep ! temperature/salinity/ssh/potential density /depth at u/v point 
     1125     !!------------------------------------------------------------- 
     1126 
     1127     !! Sum the relevant segments to obtain values for each class 
     1128     IF(sec%nb_point .NE. 0)THEN    
     1129 
     1130        !--------------------------------------! 
     1131        ! LOOP ON THE SEGMENT BETWEEN 2 NODES  ! 
     1132        !--------------------------------------! 
     1133        DO jseg=1,MAX(sec%nb_point-1,0) 
     1134            
     1135           !------------------------------------------------------------------------------------------- 
     1136           ! Select the appropriate coordinate for computing the velocity of the segment 
     1137           ! 
     1138           !                      CASE(0)                                    Case (2) 
     1139           !                      -------                                    -------- 
     1140           !  listPoint(jseg)                 listPoint(jseg+1)       listPoint(jseg)  F(i,j)       
     1141           !      F(i,j)----------V(i+1,j)-------F(i+1,j)                               | 
     1142           !                                                                            | 
     1143           !                                                                            | 
     1144           !                                                                            | 
     1145           !                      Case (3)                                            U(i,j) 
     1146           !                      --------                                              | 
     1147           !                                                                            | 
     1148           !  listPoint(jseg+1) F(i,j+1)                                                | 
     1149           !                        |                                                   | 
     1150           !                        |                                                   | 
     1151           !                        |                                 listPoint(jseg+1) F(i,j-1) 
     1152           !                        |                                             
     1153           !                        |                                             
     1154           !                     U(i,j+1)                                             
     1155           !                        |                                       Case(1)      
     1156           !                        |                                       ------       
     1157           !                        |                                             
     1158           !                        |                 listPoint(jseg+1)             listPoint(jseg)                            
     1159           !                        |                 F(i-1,j)-----------V(i,j) -------f(jseg)                            
     1160           ! listPoint(jseg)     F(i,j) 
     1161           !  
     1162           !------------------------------------------------------------------------------------------- 
     1163 
     1164           SELECT CASE( sec%direction(jseg) ) 
     1165           CASE(0)  ;   k = sec%listPoint(jseg) 
     1166           CASE(1)  ;   k = POINT_SECTION(sec%listPoint(jseg)%I+1,sec%listPoint(jseg)%J) 
     1167           CASE(2)  ;   k = sec%listPoint(jseg) 
     1168           CASE(3)  ;   k = POINT_SECTION(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J+1) 
     1169           END SELECT 
     1170 
     1171           !---------------------------| 
     1172           !     LOOP ON THE LEVEL     | 
     1173           !---------------------------| 
     1174           !Sum of the transport on the vertical  
     1175           DO jk=1,mbathy(k%I,k%J) 
     1176 
     1177              ! compute temperature, salinity, insitu & potential density, ssh and depth at U/V point 
     1178              SELECT CASE( sec%direction(jseg) ) 
     1179              CASE(0,1) 
     1180                 ztn   = interp(k%I,k%J,jk,'V',tsn(:,:,:,jp_tem) ) 
     1181                 zsn   = interp(k%I,k%J,jk,'V',tsn(:,:,:,jp_sal) ) 
     1182                 zrhop = interp(k%I,k%J,jk,'V',rhop) 
     1183                 zrhoi = interp(k%I,k%J,jk,'V',rhd*rau0+rau0) 
     1184                 zsshn =  0.5*( sshn(k%I,k%J)    + sshn(k%I,k%J+1)    ) * vmask(k%I,k%J,1) 
     1185              CASE(2,3) 
     1186                 ztn   = interp(k%I,k%J,jk,'U',tsn(:,:,:,jp_tem) ) 
     1187                 zsn   = interp(k%I,k%J,jk,'U',tsn(:,:,:,jp_sal) ) 
     1188                 zrhop = interp(k%I,k%J,jk,'U',rhop) 
     1189                 zrhoi = interp(k%I,k%J,jk,'U',rhd*rau0+rau0) 
     1190                 zsshn =  0.5*( sshn(k%I,k%J)    + sshn(k%I+1,k%J)    ) * umask(k%I,k%J,1)  
     1191              END SELECT 
     1192 
     1193              zfsdep= gdept(k%I,k%J,jk) 
     1194  
     1195              !------------------------------- 
     1196              !  LOOP ON THE DENSITY CLASSES | 
     1197              !------------------------------- 
     1198              !The computation is made for each density/heat/salt/... class 
     1199              DO jclass=1,MAX(1,sec%nb_class-1) 
     1200 
     1201                 !----------------------------------------------! 
     1202                 !TEST ON THE DENSITY/SALINITY/TEMPERATURE/LEVEL!  
     1203                 !----------------------------------------------! 
     1204  
     1205                 IF ( (                                                    & 
     1206                    ((( zrhop .GE. (sec%zsigp(jclass)+1000.  )) .AND.      & 
     1207                    (   zrhop .LE. (sec%zsigp(jclass+1)+1000. ))) .OR.     & 
     1208                    ( sec%zsigp(jclass) .EQ. 99.)) .AND.                   & 
     1209 
     1210                    ((( zrhoi .GE. (sec%zsigi(jclass) + 1000.  )) .AND.    & 
     1211                    (   zrhoi .LE. (sec%zsigi(jclass+1)+1000. ))) .OR.     & 
     1212                    ( sec%zsigi(jclass) .EQ. 99.)) .AND.                   & 
     1213 
     1214                    ((( zsn .GT. sec%zsal(jclass)) .AND.                   & 
     1215                    (   zsn .LE. sec%zsal(jclass+1))) .OR.                 & 
     1216                    ( sec%zsal(jclass) .EQ. 99.)) .AND.                    & 
     1217 
     1218                    ((( ztn .GE. sec%ztem(jclass)) .AND.                   & 
     1219                    (   ztn .LE. sec%ztem(jclass+1))) .OR.                 & 
     1220                    ( sec%ztem(jclass) .EQ.99.)) .AND.                     & 
     1221 
     1222                    ((( zfsdep .GE. sec%zlay(jclass)) .AND.                & 
     1223                    (   zfsdep .LE. sec%zlay(jclass+1))) .OR.              & 
     1224                    ( sec%zlay(jclass) .EQ. 99. ))                         & 
     1225                                                                   ))   THEN 
     1226 
     1227                    !SUM THE TRANSPORTS FOR EACH CLASSES FOR THE POSITIVE AND NEGATIVE DIRECTIONS 
     1228                    !---------------------------------------------------------------------------- 
     1229                    IF (transports_3d(1,jsec,jseg,jk) .GE. 0.0) THEN  
     1230                       sec%transport(1,jclass) = sec%transport(1,jclass)+transports_3d(1,jsec,jseg,jk) 
     1231                    ELSE 
     1232                       sec%transport(2,jclass) = sec%transport(2,jclass)+transports_3d(1,jsec,jseg,jk) 
     1233                    ENDIF 
     1234                    IF( sec%llstrpond )THEN 
     1235 
     1236                       IF( transports_3d(1,jsec,jseg,jk) .NE. 0._wp ) THEN 
     1237 
     1238                          IF (transports_3d(2,jsec,jseg,jk)/transports_3d(1,jsec,jseg,jk) .GE. 0.0 ) THEN 
     1239                             sec%transport(3,jclass) = sec%transport(3,jclass)+transports_3d(2,jsec,jseg,jk)/transports_3d(1,jsec,jseg,jk) 
     1240                          ELSE 
     1241                             sec%transport(4,jclass) = sec%transport(4,jclass)+transports_3d(2,jsec,jseg,jk)/transports_3d(1,jsec,jseg,jk) 
     1242                          ENDIF 
     1243 
     1244                          IF ( transports_3d(3,jsec,jseg,jk)/transports_3d(1,jsec,jseg,jk) .GE. 0.0 ) THEN 
     1245                             sec%transport(5,jclass) = sec%transport(5,jclass)+transports_3d(3,jsec,jseg,jk)/transports_3d(1,jsec,jseg,jk) 
     1246                          ELSE 
     1247                             sec%transport(6,jclass) = sec%transport(6,jclass)+transports_3d(3,jsec,jseg,jk)/transports_3d(1,jsec,jseg,jk) 
     1248                          ENDIF 
     1249 
     1250                       ENDIF 
     1251 
     1252                       IF ( transports_3d(4,jsec,jseg,jk) .GE. 0.0 ) THEN 
     1253                          sec%transport(7,jclass) = sec%transport(7,jclass)+transports_3d(4,jsec,jseg,jk) 
     1254                       ELSE 
     1255                          sec%transport(8,jclass) = sec%transport(8,jclass)+transports_3d(4,jsec,jseg,jk) 
     1256                       ENDIF 
     1257 
     1258                       IF ( transports_3d(5,jsec,jseg,jk) .GE. 0.0 ) THEN 
     1259                          sec%transport( 9,jclass) = sec%transport( 9,jclass)+transports_3d(5,jsec,jseg,jk) 
     1260                       ELSE 
     1261                          sec%transport(10,jclass) = sec%transport(10,jclass)+transports_3d(5,jsec,jseg,jk) 
     1262                       ENDIF 
     1263 
     1264                    ELSE 
     1265                       sec%transport( 3,jclass) = 0._wp 
     1266                       sec%transport( 4,jclass) = 0._wp 
     1267                       sec%transport( 5,jclass) = 0._wp 
     1268                       sec%transport( 6,jclass) = 0._wp 
     1269                       sec%transport( 7,jclass) = 0._wp 
     1270                       sec%transport( 8,jclass) = 0._wp 
     1271                       sec%transport( 9,jclass) = 0._wp 
     1272                       sec%transport(10,jclass) = 0._wp 
     1273                    ENDIF 
     1274 
     1275                 ENDIF ! end of test if point is in class 
     1276    
     1277              ENDDO ! end of loop on the classes 
     1278 
     1279           ENDDO ! loop over jk 
     1280 
     1281#if defined key_lim2 || defined key_lim3 
     1282 
     1283           !ICE CASE     
     1284           IF( sec%ll_ice_section )THEN 
     1285 
     1286              IF ( transports_2d(1,jsec,jseg) .GE. 0.0 ) THEN 
     1287                 sec%transport(11,1) = sec%transport(11,1)+transports_2d(1,jsec,jseg) 
     1288              ELSE 
     1289                 sec%transport(12,1) = sec%transport(12,1)+transports_2d(1,jsec,jseg) 
     1290              ENDIF 
     1291 
     1292              IF ( transports_2d(3,jsec,jseg) .GE. 0.0 ) THEN 
     1293                 sec%transport(13,1) = sec%transport(13,1)+transports_2d(2,jsec,jseg) 
     1294              ELSE 
     1295                 sec%transport(14,1) = sec%transport(14,1)+transports_2d(2,jsec,jseg) 
     1296              ENDIF 
     1297 
     1298           ENDIF !end of ice case 
     1299#endif 
     1300  
     1301        ENDDO !end of loop on the segment 
     1302 
     1303     ELSE  !if sec%nb_point =0 
     1304        sec%transport(1:2,:)=0. 
     1305        IF (sec%llstrpond) sec%transport(3:10,:)=0. 
     1306        IF (sec%ll_ice_section) sec%transport( 11:14,:)=0. 
     1307     ENDIF !end of sec%nb_point =0 case 
     1308 
     1309  END SUBROUTINE dia_dct_sum 
     1310  
     1311  SUBROUTINE dia_dct_sum_h(sec,jsec) 
     1312     !!------------------------------------------------------------- 
     1313     !! Exactly as dia_dct_sum but for hourly files containing data summed at each time step 
     1314     !! 
     1315     !! Purpose: Average the transport over nn_dctwri time steps  
     1316     !! and sum over the density/salinity/temperature/depth classes 
     1317     !! 
     1318     !! Method:  
     1319     !!           Sum over relevant grid cells to obtain values 
     1320     !!           for each 
     1321     !!              There are several loops:                  
     1322     !!              loop on the segment between 2 nodes 
     1323     !!              loop on the level jk 
     1324     !!              loop on the density/temperature/salinity/level classes 
     1325     !!              test on the density/temperature/salinity/level 
     1326     !! 
     1327     !!  ** Method  :Transport through a given section is equal to the sum of transports 
     1328     !!              computed on each proc. 
     1329     !!              On each proc,transport is equal to the sum of transport computed through 
     1330     !!              segments linking each point of sec%listPoint  with the next one.    
     1331     !! 
     1332     !!------------------------------------------------------------- 
     1333     !! * arguments 
     1334     TYPE(SECTION),INTENT(INOUT) :: sec 
     1335     INTEGER      ,INTENT(IN)    :: jsec        ! numeric identifier of section 
     1336 
     1337     TYPE(POINT_SECTION) :: k 
     1338     INTEGER  :: jk,jseg,jclass                        !loop on level/segment/classes  
     1339     REAL(wp) :: ztn, zsn, zrhoi, zrhop, zsshn, zfsdep ! temperature/salinity/ssh/potential density /depth at u/v point 
     1340     !!------------------------------------------------------------- 
     1341 
     1342     !! Sum the relevant segments to obtain values for each class 
     1343     IF(sec%nb_point .NE. 0)THEN    
     1344 
     1345        !--------------------------------------! 
     1346        ! LOOP ON THE SEGMENT BETWEEN 2 NODES  ! 
     1347        !--------------------------------------! 
     1348        DO jseg=1,MAX(sec%nb_point-1,0) 
     1349            
     1350           !------------------------------------------------------------------------------------------- 
     1351           ! Select the appropriate coordinate for computing the velocity of the segment 
     1352           ! 
     1353           !                      CASE(0)                                    Case (2) 
     1354           !                      -------                                    -------- 
     1355           !  listPoint(jseg)                 listPoint(jseg+1)       listPoint(jseg)  F(i,j)       
     1356           !      F(i,j)----------V(i+1,j)-------F(i+1,j)                               | 
     1357           !                                                                            | 
     1358           !                                                                            | 
     1359           !                                                                            | 
     1360           !                      Case (3)                                            U(i,j) 
     1361           !                      --------                                              | 
     1362           !                                                                            | 
     1363           !  listPoint(jseg+1) F(i,j+1)                                                | 
     1364           !                        |                                                   | 
     1365           !                        |                                                   | 
     1366           !                        |                                 listPoint(jseg+1) F(i,j-1) 
     1367           !                        |                                             
     1368           !                        |                                             
     1369           !                     U(i,j+1)                                             
     1370           !                        |                                       Case(1)      
     1371           !                        |                                       ------       
     1372           !                        |                                             
     1373           !                        |                 listPoint(jseg+1)             listPoint(jseg)                            
     1374           !                        |                 F(i-1,j)-----------V(i,j) -------f(jseg)                            
     1375           ! listPoint(jseg)     F(i,j) 
     1376           !  
     1377           !------------------------------------------------------------------------------------------- 
     1378 
     1379           SELECT CASE( sec%direction(jseg) ) 
     1380           CASE(0)  ;   k = sec%listPoint(jseg) 
     1381           CASE(1)  ;   k = POINT_SECTION(sec%listPoint(jseg)%I+1,sec%listPoint(jseg)%J) 
     1382           CASE(2)  ;   k = sec%listPoint(jseg) 
     1383           CASE(3)  ;   k = POINT_SECTION(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J+1) 
     1384           END SELECT 
     1385 
     1386           !---------------------------| 
     1387           !     LOOP ON THE LEVEL     | 
     1388           !---------------------------| 
     1389           !Sum of the transport on the vertical  
     1390           DO jk=1,mbathy(k%I,k%J) 
     1391 
     1392              ! compute temperature, salinity, insitu & potential density, ssh and depth at U/V point 
     1393              SELECT CASE( sec%direction(jseg) ) 
     1394              CASE(0,1) 
     1395                 ztn   = interp(k%I,k%J,jk,'V',tsn(:,:,:,jp_tem) ) 
     1396                 zsn   = interp(k%I,k%J,jk,'V',tsn(:,:,:,jp_sal) ) 
     1397                 zrhop = interp(k%I,k%J,jk,'V',rhop) 
     1398                 zrhoi = interp(k%I,k%J,jk,'V',rhd*rau0+rau0) 
     1399                 zsshn =  0.5*( sshn(k%I,k%J)    + sshn(k%I,k%J+1)    ) * vmask(k%I,k%J,1) 
     1400              CASE(2,3) 
     1401                 ztn   = interp(k%I,k%J,jk,'U',tsn(:,:,:,jp_tem) ) 
     1402                 zsn   = interp(k%I,k%J,jk,'U',tsn(:,:,:,jp_sal) ) 
     1403                 zrhop = interp(k%I,k%J,jk,'U',rhop) 
     1404                 zrhoi = interp(k%I,k%J,jk,'U',rhd*rau0+rau0) 
     1405                 zsshn =  0.5*( sshn(k%I,k%J)    + sshn(k%I+1,k%J)    ) * umask(k%I,k%J,1)  
     1406              END SELECT 
     1407 
     1408              zfsdep= gdept(k%I,k%J,jk) 
     1409  
     1410              !------------------------------- 
     1411              !  LOOP ON THE DENSITY CLASSES | 
     1412              !------------------------------- 
     1413              !The computation is made for each density/heat/salt/... class 
     1414              DO jclass=1,MAX(1,sec%nb_class-1) 
     1415 
     1416                 !----------------------------------------------! 
     1417                 !TEST ON THE DENSITY/SALINITY/TEMPERATURE/LEVEL!  
     1418                 !----------------------------------------------! 
     1419  
     1420                 IF ( (                                                    & 
     1421                    ((( zrhop .GE. (sec%zsigp(jclass)+1000.  )) .AND.      & 
     1422                    (   zrhop .LE. (sec%zsigp(jclass+1)+1000. ))) .OR.     & 
     1423                    ( sec%zsigp(jclass) .EQ. 99.)) .AND.                   & 
     1424 
     1425                    ((( zrhoi .GE. (sec%zsigi(jclass) + 1000.  )) .AND.    & 
     1426                    (   zrhoi .LE. (sec%zsigi(jclass+1)+1000. ))) .OR.     & 
     1427                    ( sec%zsigi(jclass) .EQ. 99.)) .AND.                   & 
     1428 
     1429                    ((( zsn .GT. sec%zsal(jclass)) .AND.                   & 
     1430                    (   zsn .LE. sec%zsal(jclass+1))) .OR.                 & 
     1431                    ( sec%zsal(jclass) .EQ. 99.)) .AND.                    & 
     1432 
     1433                    ((( ztn .GE. sec%ztem(jclass)) .AND.                   & 
     1434                    (   ztn .LE. sec%ztem(jclass+1))) .OR.                 & 
     1435                    ( sec%ztem(jclass) .EQ.99.)) .AND.                     & 
     1436 
     1437                    ((( zfsdep .GE. sec%zlay(jclass)) .AND.                & 
     1438                    (   zfsdep .LE. sec%zlay(jclass+1))) .OR.              & 
     1439                    ( sec%zlay(jclass) .EQ. 99. ))                         & 
     1440                                                                   ))   THEN 
     1441 
     1442                    !SUM THE TRANSPORTS FOR EACH CLASSES FOR THE POSITIVE AND NEGATIVE DIRECTIONS 
     1443                    !---------------------------------------------------------------------------- 
     1444                    IF (transports_3d_h(1,jsec,jseg,jk) .GE. 0.0) THEN  
     1445                       sec%transport_h(1,jclass) = sec%transport_h(1,jclass)+transports_3d_h(1,jsec,jseg,jk) 
     1446                    ELSE 
     1447                       sec%transport_h(2,jclass) = sec%transport_h(2,jclass)+transports_3d_h(1,jsec,jseg,jk) 
     1448                    ENDIF 
     1449                    IF( sec%llstrpond )THEN 
     1450 
     1451                       IF( transports_3d_h(1,jsec,jseg,jk) .NE. 0._wp ) THEN 
     1452 
     1453                          IF (transports_3d_h(2,jsec,jseg,jk)/transports_3d_h(1,jsec,jseg,jk) .GE. 0.0 ) THEN 
     1454                             sec%transport_h(3,jclass) = sec%transport_h(3,jclass)+transports_3d_h(2,jsec,jseg,jk)/transports_3d_h(1,jsec,jseg,jk) 
     1455                          ELSE 
     1456                             sec%transport_h(4,jclass) = sec%transport_h(4,jclass)+transports_3d_h(2,jsec,jseg,jk)/transports_3d_h(1,jsec,jseg,jk) 
     1457                          ENDIF 
     1458 
     1459                          IF ( transports_3d_h(3,jsec,jseg,jk)/transports_3d_h(1,jsec,jseg,jk) .GE. 0.0 ) THEN 
     1460                             sec%transport_h(5,jclass) = sec%transport_h(5,jclass)+transports_3d_h(3,jsec,jseg,jk)/transports_3d_h(1,jsec,jseg,jk) 
     1461                          ELSE 
     1462                             sec%transport_h(6,jclass) = sec%transport_h(6,jclass)+transports_3d_h(3,jsec,jseg,jk)/transports_3d_h(1,jsec,jseg,jk) 
     1463                          ENDIF 
     1464 
     1465                       ENDIF 
     1466 
     1467                       IF ( transports_3d_h(4,jsec,jseg,jk) .GE. 0.0 ) THEN 
     1468                          sec%transport_h(7,jclass) = sec%transport_h(7,jclass)+transports_3d_h(4,jsec,jseg,jk) 
     1469                       ELSE 
     1470                          sec%transport_h(8,jclass) = sec%transport_h(8,jclass)+transports_3d_h(4,jsec,jseg,jk) 
     1471                       ENDIF 
     1472 
     1473                       IF ( transports_3d_h(5,jsec,jseg,jk) .GE. 0.0 ) THEN 
     1474                          sec%transport_h( 9,jclass) = sec%transport_h( 9,jclass)+transports_3d_h(5,jsec,jseg,jk) 
     1475                       ELSE 
     1476                          sec%transport_h(10,jclass) = sec%transport_h(10,jclass)+transports_3d_h(5,jsec,jseg,jk) 
     1477                       ENDIF 
     1478 
     1479                    ELSE 
     1480                       sec%transport_h( 3,jclass) = 0._wp 
     1481                       sec%transport_h( 4,jclass) = 0._wp 
     1482                       sec%transport_h( 5,jclass) = 0._wp 
     1483                       sec%transport_h( 6,jclass) = 0._wp 
     1484                       sec%transport_h( 7,jclass) = 0._wp 
     1485                       sec%transport_h( 8,jclass) = 0._wp 
     1486                       sec%transport_h( 9,jclass) = 0._wp 
     1487                       sec%transport_h(10,jclass) = 0._wp 
     1488                    ENDIF 
     1489 
     1490                 ENDIF ! end of test if point is in class 
     1491    
     1492              ENDDO ! end of loop on the classes 
     1493 
     1494           ENDDO ! loop over jk 
     1495 
     1496#if defined key_lim2 || defined key_lim3 
     1497 
     1498           !ICE CASE     
     1499           IF( sec%ll_ice_section )THEN 
     1500 
     1501              IF ( transports_2d_h(1,jsec,jseg) .GE. 0.0 ) THEN 
     1502                 sec%transport_h(11,1) = sec%transport_h(11,1)+transports_2d_h(1,jsec,jseg) 
     1503              ELSE 
     1504                 sec%transport_h(12,1) = sec%transport_h(12,1)+transports_2d_h(1,jsec,jseg) 
     1505              ENDIF 
     1506 
     1507              IF ( transports_2d_h(3,jsec,jseg) .GE. 0.0 ) THEN 
     1508                 sec%transport_h(13,1) = sec%transport_h(13,1)+transports_2d_h(2,jsec,jseg) 
     1509              ELSE 
     1510                 sec%transport_h(14,1) = sec%transport_h(14,1)+transports_2d_h(2,jsec,jseg) 
     1511              ENDIF 
     1512 
     1513           ENDIF !end of ice case 
     1514#endif 
     1515  
     1516        ENDDO !end of loop on the segment 
     1517 
     1518     ELSE  !if sec%nb_point =0 
     1519        sec%transport_h(1:2,:)=0. 
     1520        IF (sec%llstrpond) sec%transport_h(3:10,:)=0. 
     1521        IF (sec%ll_ice_section) sec%transport_h( 11:14,:)=0. 
     1522     ENDIF !end of sec%nb_point =0 case 
     1523 
     1524  END SUBROUTINE dia_dct_sum_h 
     1525  
     1526  SUBROUTINE dia_dct_wri_NOOS(kt,ksec,sec) 
     1527     !!------------------------------------------------------------- 
     1528     !! Write transport output in numdct using NOOS formatting  
     1529     !!  
     1530     !! Purpose: Write  transports in ascii files 
     1531     !!  
     1532     !! Method: 
     1533     !!        1. Write volume transports in "volume_transport" 
     1534     !!           Unit: Sv : area * Velocity / 1.e6  
     1535     !!  
     1536     !!        2. Write heat transports in "heat_transport" 
     1537     !!           Unit: Peta W : area * Velocity * T * rhau * Cp / 1.e15 
     1538     !!  
     1539     !!        3. Write salt transports in "salt_transport" 
     1540     !!           Unit: 10^9 g m^3 / s : area * Velocity * S / 1.e6 
     1541     !! 
     1542     !!-------------------------------------------------------------  
     1543     !!arguments 
     1544     INTEGER, INTENT(IN)          :: kt          ! time-step 
     1545     TYPE(SECTION), INTENT(INOUT) :: sec         ! section to write    
     1546     INTEGER ,INTENT(IN)          :: ksec        ! section number 
     1547 
     1548     !!local declarations 
     1549     INTEGER               :: jclass,ji             ! Dummy loop 
     1550     CHARACTER(len=2)      :: classe             ! Classname  
     1551     REAL(wp)              :: zbnd1,zbnd2        ! Class bounds 
     1552     REAL(wp)              :: zslope             ! section's slope coeff 
     1553     ! 
     1554     REAL(wp), POINTER, DIMENSION(:):: zsumclasses ! 1D workspace  
     1555     !!-------------------------------------------------------------  
     1556     CALL wrk_alloc(nb_type , zsumclasses )   
     1557 
     1558     zsumclasses(:)=0._wp 
     1559     zslope = sec%slopeSection        
     1560 
     1561     WRITE(numdct_NOOS,'(I4,a1,I2,a1,I2,a12,i3,a17,i3,a10,a25)') nyear,'.',nmonth,'.',nday,'   Transect:',ksec-1,'   No. of layers:',sec%nb_class-1,'   Name: ',sec%name 
     1562 
     1563     DO jclass=1,MAX(1,sec%nb_class-1) 
     1564        zsumclasses(1:nb_type)=zsumclasses(1:nb_type)+sec%transport(1:nb_type,jclass) 
     1565     ENDDO 
     1566  
     1567     classe   = 'total   ' 
     1568     zbnd1   = 0._wp 
     1569     zbnd2   = 0._wp 
     1570 
     1571     IF ( zslope .gt. 0._wp .and. zslope .ne. 10000._wp ) THEN 
     1572        WRITE(numdct_NOOS,'(9e12.4E2)') -(zsumclasses( 1)+zsumclasses( 2)), -zsumclasses( 2),-zsumclasses( 1),   & 
     1573                                        -(zsumclasses( 7)+zsumclasses( 8)), -zsumclasses( 8),-zsumclasses( 7),   & 
     1574                                        -(zsumclasses( 9)+zsumclasses(10)), -zsumclasses(10),-zsumclasses( 9) 
     1575     ELSE 
     1576        WRITE(numdct_NOOS,'(9e12.4E2)')   zsumclasses( 1)+zsumclasses( 2) ,  zsumclasses( 1), zsumclasses( 2),   & 
     1577                                          zsumclasses( 7)+zsumclasses( 8) ,  zsumclasses( 7), zsumclasses( 8),   & 
     1578                                          zsumclasses( 9)+zsumclasses(10) ,  zsumclasses( 9), zsumclasses(10) 
     1579     ENDIF  
     1580 
     1581     DO jclass=1,MAX(1,sec%nb_class-1) 
     1582    
     1583        classe   = 'N       ' 
     1584        zbnd1   = 0._wp 
     1585        zbnd2   = 0._wp 
     1586 
     1587        !insitu density classes transports 
     1588        IF( ( sec%zsigi(jclass)   .NE. 99._wp ) .AND. & 
     1589            ( sec%zsigi(jclass+1) .NE. 99._wp )       )THEN 
     1590           classe = 'DI       ' 
     1591           zbnd1 = sec%zsigi(jclass) 
     1592           zbnd2 = sec%zsigi(jclass+1) 
     1593        ENDIF 
     1594        !potential density classes transports 
     1595        IF( ( sec%zsigp(jclass)   .NE. 99._wp ) .AND. & 
     1596            ( sec%zsigp(jclass+1) .NE. 99._wp )       )THEN 
     1597           classe = 'DP      ' 
     1598           zbnd1 = sec%zsigp(jclass) 
     1599           zbnd2 = sec%zsigp(jclass+1) 
     1600        ENDIF 
     1601        !depth classes transports 
     1602        IF( ( sec%zlay(jclass)    .NE. 99._wp ) .AND. & 
     1603            ( sec%zlay(jclass+1)  .NE. 99._wp )       )THEN  
     1604           classe = 'Z       ' 
     1605           zbnd1 = sec%zlay(jclass) 
     1606           zbnd2 = sec%zlay(jclass+1) 
     1607        ENDIF 
     1608        !salinity classes transports 
     1609        IF( ( sec%zsal(jclass) .NE. 99._wp    ) .AND. & 
     1610            ( sec%zsal(jclass+1) .NE. 99._wp  )       )THEN 
     1611           classe = 'S       ' 
     1612           zbnd1 = sec%zsal(jclass) 
     1613           zbnd2 = sec%zsal(jclass+1)    
     1614        ENDIF 
     1615        !temperature classes transports 
     1616        IF( ( sec%ztem(jclass) .NE. 99._wp     ) .AND. & 
     1617            ( sec%ztem(jclass+1) .NE. 99._wp     )       ) THEN 
     1618           classe = 'T       ' 
     1619           zbnd1 = sec%ztem(jclass) 
     1620           zbnd2 = sec%ztem(jclass+1) 
     1621        ENDIF 
     1622                   
     1623        !write volume transport per class 
     1624        IF ( zslope .gt. 0._wp .and. zslope .ne. 10000._wp ) THEN 
     1625           WRITE(numdct_NOOS,'(9e12.4E2)') -(sec%transport( 1,jclass)+sec%transport( 2,jclass)),-sec%transport( 2,jclass),-sec%transport( 1,jclass), & 
     1626                                           -(sec%transport( 7,jclass)+sec%transport( 8,jclass)),-sec%transport( 8,jclass),-sec%transport( 7,jclass), & 
     1627                                           -(sec%transport( 9,jclass)+sec%transport(10,jclass)),-sec%transport(10,jclass),-sec%transport( 9,jclass) 
     1628        ELSE 
     1629           WRITE(numdct_NOOS,'(9e12.4E2)')   sec%transport( 1,jclass)+sec%transport( 2,jclass) , sec%transport( 1,jclass), sec%transport( 2,jclass), & 
     1630                                             sec%transport( 7,jclass)+sec%transport( 8,jclass) , sec%transport( 7,jclass), sec%transport( 8,jclass), & 
     1631                                             sec%transport( 9,jclass)+sec%transport(10,jclass) , sec%transport( 9,jclass), sec%transport(10,jclass) 
     1632        ENDIF 
     1633 
     1634     ENDDO 
     1635 
     1636     CALL wrk_dealloc(nb_type , zsumclasses )   
     1637 
     1638  END SUBROUTINE dia_dct_wri_NOOS 
     1639 
     1640  SUBROUTINE dia_dct_wri_NOOS_h(hr,ksec,sec) 
     1641     !!------------------------------------------------------------- 
     1642     !! As routine dia_dct_wri_NOOS but for hourly output files 
     1643     !! 
     1644     !! Write transport output in numdct using NOOS formatting  
     1645     !!  
     1646     !! Purpose: Write  transports in ascii files 
     1647     !!  
     1648     !! Method: 
     1649     !!        1. Write volume transports in "volume_transport" 
     1650     !!           Unit: Sv : area * Velocity / 1.e6  
     1651     !! 
     1652     !!-------------------------------------------------------------  
     1653     !!arguments 
     1654     INTEGER, INTENT(IN)          :: hr          ! hour 
     1655     TYPE(SECTION), INTENT(INOUT) :: sec         ! section to write    
     1656     INTEGER ,INTENT(IN)          :: ksec        ! section number 
     1657 
     1658     !!local declarations 
     1659     INTEGER               :: jclass,jhr            ! Dummy loop 
     1660     CHARACTER(len=2)      :: classe             ! Classname  
     1661     REAL(wp)              :: zbnd1,zbnd2        ! Class bounds 
     1662     REAL(wp)              :: zslope             ! section's slope coeff 
     1663     ! 
     1664     REAL(wp), POINTER, DIMENSION(:):: zsumclasses ! 1D workspace  
     1665     !!-------------------------------------------------------------  
     1666 
     1667     CALL wrk_alloc(nb_type , zsumclasses )  
     1668 
     1669     zsumclasses(:)=0._wp 
     1670     zslope = sec%slopeSection        
     1671 
     1672     DO jclass=1,MAX(1,sec%nb_class-1) 
     1673        zsumclasses(1:nb_type)=zsumclasses(1:nb_type)+sec%transport_h(1:nb_type,jclass) 
     1674     ENDDO 
     1675  
     1676     !write volume transport per class 
     1677     IF ( zslope .gt. 0._wp .and. zslope .ne. 10000._wp ) THEN 
     1678        z_hr_output(ksec,hr,1)=-(zsumclasses(1)+zsumclasses(2)) 
     1679     ELSE 
     1680        z_hr_output(ksec,hr,1)= (zsumclasses(1)+zsumclasses(2)) 
     1681     ENDIF 
     1682 
     1683     DO jclass=1,MAX(1,sec%nb_class-1) 
     1684        IF ( zslope .gt. 0._wp .and. zslope .ne. 10000._wp ) THEN 
     1685           z_hr_output(ksec,hr,jclass+1)=-(sec%transport_h(1,jclass)+sec%transport_h(2,jclass)) 
     1686        ELSE 
     1687           z_hr_output(ksec,hr,jclass+1)= (sec%transport_h(1,jclass)+sec%transport_h(2,jclass)) 
     1688        ENDIF 
     1689     ENDDO 
     1690 
     1691     IF ( hr .eq. 48._wp ) THEN 
     1692        WRITE(numdct_NOOS_h,'(I4,a1,I2,a1,I2,a12,i3,a17,i3)') nyear,'.',nmonth,'.',nday,'   Transect:',ksec-1,'   No. of layers:',sec%nb_class-1 
     1693        DO jhr=25,48 
     1694           WRITE(numdct_NOOS_h,'(11F12.1)')  z_hr_output(ksec,jhr,1), (z_hr_output(ksec,jhr,jclass+1), jclass=1,MAX(1,10) ) 
     1695        ENDDO 
     1696     ENDIF 
     1697 
     1698     CALL wrk_dealloc(nb_type , zsumclasses ) 
     1699 
     1700  END SUBROUTINE dia_dct_wri_NOOS_h 
     1701 
    8961702  SUBROUTINE dia_dct_wri(kt,ksec,sec) 
    8971703     !!------------------------------------------------------------- 
     
    9171723 
    9181724     !!local declarations 
    919      INTEGER               :: jcl,ji             ! Dummy loop 
     1725     INTEGER               :: jclass             ! Dummy loop 
    9201726     CHARACTER(len=2)      :: classe             ! Classname  
    9211727     REAL(wp)              :: zbnd1,zbnd2        ! Class bounds 
    9221728     REAL(wp)              :: zslope             ! section's slope coeff 
    9231729     ! 
    924      REAL(wp), POINTER, DIMENSION(:):: zsumclass ! 1D workspace  
     1730     REAL(wp), POINTER, DIMENSION(:):: zsumclasses ! 1D workspace  
    9251731     !!-------------------------------------------------------------  
    926      CALL wrk_alloc(nb_type_class , zsumclass )   
    927  
    928      zsumclass(:)=0._wp 
     1732     CALL wrk_alloc(nb_type , zsumclasses )   
     1733 
     1734     zsumclasses(:)=0._wp 
    9291735     zslope = sec%slopeSection        
    930  
    931   
    932      DO jcl=1,MAX(1,sec%nb_class-1) 
    933  
    934         ! Mean computation 
    935         sec%transport(:,jcl)=sec%transport(:,jcl)/(nn_dctwri/nn_dct) 
     1736  
     1737     DO jclass=1,MAX(1,sec%nb_class-1) 
     1738 
    9361739        classe   = 'N       ' 
    9371740        zbnd1   = 0._wp 
    9381741        zbnd2   = 0._wp 
    939         zsumclass(1:nb_type_class)=zsumclass(1:nb_type_class)+sec%transport(1:nb_type_class,jcl) 
     1742        zsumclasses(1:nb_type)=zsumclasses(1:nb_type)+sec%transport(1:nb_type,jclass) 
    9401743 
    9411744    
    9421745        !insitu density classes transports 
    943         IF( ( sec%zsigi(jcl)   .NE. 99._wp ) .AND. & 
    944             ( sec%zsigi(jcl+1) .NE. 99._wp )       )THEN 
     1746        IF( ( sec%zsigi(jclass)   .NE. 99._wp ) .AND. & 
     1747            ( sec%zsigi(jclass+1) .NE. 99._wp )       )THEN 
    9451748           classe = 'DI       ' 
    946            zbnd1 = sec%zsigi(jcl) 
    947            zbnd2 = sec%zsigi(jcl+1) 
     1749           zbnd1 = sec%zsigi(jclass) 
     1750           zbnd2 = sec%zsigi(jclass+1) 
    9481751        ENDIF 
    9491752        !potential density classes transports 
    950         IF( ( sec%zsigp(jcl)   .NE. 99._wp ) .AND. & 
    951             ( sec%zsigp(jcl+1) .NE. 99._wp )       )THEN 
     1753        IF( ( sec%zsigp(jclass)   .NE. 99._wp ) .AND. & 
     1754            ( sec%zsigp(jclass+1) .NE. 99._wp )       )THEN 
    9521755           classe = 'DP      ' 
    953            zbnd1 = sec%zsigp(jcl) 
    954            zbnd2 = sec%zsigp(jcl+1) 
     1756           zbnd1 = sec%zsigp(jclass) 
     1757           zbnd2 = sec%zsigp(jclass+1) 
    9551758        ENDIF 
    9561759        !depth classes transports 
    957         IF( ( sec%zlay(jcl)    .NE. 99._wp ) .AND. & 
    958             ( sec%zlay(jcl+1)  .NE. 99._wp )       )THEN  
     1760        IF( ( sec%zlay(jclass)    .NE. 99._wp ) .AND. & 
     1761            ( sec%zlay(jclass+1)  .NE. 99._wp )       )THEN  
    9591762           classe = 'Z       ' 
    960            zbnd1 = sec%zlay(jcl) 
    961            zbnd2 = sec%zlay(jcl+1) 
     1763           zbnd1 = sec%zlay(jclass) 
     1764           zbnd2 = sec%zlay(jclass+1) 
    9621765        ENDIF 
    9631766        !salinity classes transports 
    964         IF( ( sec%zsal(jcl) .NE. 99._wp    ) .AND. & 
    965             ( sec%zsal(jcl+1) .NE. 99._wp  )       )THEN 
     1767        IF( ( sec%zsal(jclass) .NE. 99._wp    ) .AND. & 
     1768            ( sec%zsal(jclass+1) .NE. 99._wp  )       )THEN 
    9661769           classe = 'S       ' 
    967            zbnd1 = sec%zsal(jcl) 
    968            zbnd2 = sec%zsal(jcl+1)    
     1770           zbnd1 = sec%zsal(jclass) 
     1771           zbnd2 = sec%zsal(jclass+1)    
    9691772        ENDIF 
    9701773        !temperature classes transports 
    971         IF( ( sec%ztem(jcl) .NE. 99._wp     ) .AND. & 
    972             ( sec%ztem(jcl+1) .NE. 99._wp     )       ) THEN 
     1774        IF( ( sec%ztem(jclass) .NE. 99._wp     ) .AND. & 
     1775            ( sec%ztem(jclass+1) .NE. 99._wp     )       ) THEN 
    9731776           classe = 'T       ' 
    974            zbnd1 = sec%ztem(jcl) 
    975            zbnd2 = sec%ztem(jcl+1) 
     1777           zbnd1 = sec%ztem(jclass) 
     1778           zbnd2 = sec%ztem(jclass+1) 
    9761779        ENDIF 
    9771780                   
    9781781        !write volume transport per class 
    9791782        WRITE(numdct_vol,118) ndastp,kt,ksec,sec%name,zslope, & 
    980                               jcl,classe,zbnd1,zbnd2,& 
    981                               sec%transport(1,jcl),sec%transport(2,jcl), & 
    982                               sec%transport(1,jcl)+sec%transport(2,jcl) 
     1783                              jclass,classe,zbnd1,zbnd2,& 
     1784                              sec%transport(1,jclass),sec%transport(2,jclass), & 
     1785                              sec%transport(1,jclass)+sec%transport(2,jclass) 
    9831786 
    9841787        IF( sec%llstrpond )THEN 
    9851788 
    9861789           !write heat transport per class: 
    987            WRITE(numdct_heat,119) ndastp,kt,ksec,sec%name,zslope,  & 
    988                               jcl,classe,zbnd1,zbnd2,& 
    989                               sec%transport(7,jcl)*1000._wp*rcp/1.e15,sec%transport(8,jcl)*1000._wp*rcp/1.e15, & 
    990                               ( sec%transport(7,jcl)+sec%transport(8,jcl) )*1000._wp*rcp/1.e15 
     1790           WRITE(numdct_temp,119) ndastp,kt,ksec,sec%name,zslope,  & 
     1791                              jclass,classe,zbnd1,zbnd2,& 
     1792                              sec%transport(7,jclass)*1000._wp*rcp/1.e15,sec%transport(8,jclass)*1000._wp*rcp/1.e15, & 
     1793                              ( sec%transport(7,jclass)+sec%transport(8,jclass) )*1000._wp*rcp/1.e15 
    9911794           !write salt transport per class 
    992            WRITE(numdct_salt,119) ndastp,kt,ksec,sec%name,zslope,  & 
    993                               jcl,classe,zbnd1,zbnd2,& 
    994                               sec%transport(9,jcl)*1000._wp/1.e9,sec%transport(10,jcl)*1000._wp/1.e9,& 
    995                               (sec%transport(9,jcl)+sec%transport(10,jcl))*1000._wp/1.e9 
     1795           WRITE(numdct_sal ,119) ndastp,kt,ksec,sec%name,zslope,  & 
     1796                              jclass,classe,zbnd1,zbnd2,& 
     1797                              sec%transport(9,jclass)*1000._wp/1.e9,sec%transport(10,jclass)*1000._wp/1.e9,& 
     1798                              (sec%transport(9,jclass)+sec%transport(10,jclass))*1000._wp/1.e9 
    9961799        ENDIF 
    9971800 
     
    10001803     zbnd1 = 0._wp 
    10011804     zbnd2 = 0._wp 
    1002      jcl=0 
     1805     jclass=0 
    10031806 
    10041807     !write total volume transport 
    10051808     WRITE(numdct_vol,118) ndastp,kt,ksec,sec%name,zslope, & 
    1006                            jcl,"total",zbnd1,zbnd2,& 
    1007                            zsumclass(1),zsumclass(2),zsumclass(1)+zsumclass(2) 
     1809                           jclass,"total",zbnd1,zbnd2,& 
     1810                           zsumclasses(1),zsumclasses(2),zsumclasses(1)+zsumclasses(2) 
    10081811 
    10091812     IF( sec%llstrpond )THEN 
    10101813 
    10111814        !write total heat transport 
    1012         WRITE(numdct_heat,119) ndastp,kt,ksec,sec%name,zslope, & 
    1013                            jcl,"total",zbnd1,zbnd2,& 
    1014                            zsumclass(7)* 1000._wp*rcp/1.e15,zsumclass(8)* 1000._wp*rcp/1.e15,& 
    1015                            (zsumclass(7)+zsumclass(8) )* 1000._wp*rcp/1.e15 
     1815        WRITE(numdct_temp,119) ndastp,kt,ksec,sec%name,zslope, & 
     1816                           jclass,"total",zbnd1,zbnd2,& 
     1817                           zsumclasses(7)* 1000._wp*rcp/1.e15,zsumclasses(8)* 1000._wp*rcp/1.e15,& 
     1818                           (zsumclasses(7)+zsumclasses(8) )* 1000._wp*rcp/1.e15 
    10161819        !write total salt transport 
    1017         WRITE(numdct_salt,119) ndastp,kt,ksec,sec%name,zslope, & 
    1018                            jcl,"total",zbnd1,zbnd2,& 
    1019                            zsumclass(9)*1000._wp/1.e9,zsumclass(10)*1000._wp/1.e9,& 
    1020                            (zsumclass(9)+zsumclass(10))*1000._wp/1.e9 
     1820        WRITE(numdct_sal ,119) ndastp,kt,ksec,sec%name,zslope, & 
     1821                           jclass,"total",zbnd1,zbnd2,& 
     1822                           zsumclasses(9)*1000._wp/1.e9,zsumclasses(10)*1000._wp/1.e9,& 
     1823                           (zsumclasses(9)+zsumclasses(10))*1000._wp/1.e9 
    10211824     ENDIF 
    10221825 
     
    10251828        !write total ice volume transport 
    10261829        WRITE(numdct_vol,118) ndastp,kt,ksec,sec%name,zslope,& 
    1027                               jcl,"ice_vol",zbnd1,zbnd2,& 
    1028                               sec%transport(9,1),sec%transport(10,1),& 
    1029                               sec%transport(9,1)+sec%transport(10,1) 
     1830                              jclass,"ice_vol",zbnd1,zbnd2,& 
     1831                              sec%transport(11,1),sec%transport(12,1),& 
     1832                              sec%transport(11,1)+sec%transport(12,1) 
    10301833        !write total ice surface transport 
    10311834        WRITE(numdct_vol,118) ndastp,kt,ksec,sec%name,zslope,& 
    1032                               jcl,"ice_surf",zbnd1,zbnd2,& 
    1033                               sec%transport(11,1),sec%transport(12,1), & 
    1034                               sec%transport(11,1)+sec%transport(12,1)  
     1835                              jclass,"ice_surf",zbnd1,zbnd2,& 
     1836                              sec%transport(13,1),sec%transport(14,1), & 
     1837                              sec%transport(13,1)+sec%transport(14,1)  
    10351838     ENDIF 
    10361839                                               
     
    10381841119 FORMAT(I8,1X,I8,1X,I4,1X,A30,1X,f9.2,1X,I4,3X,A8,1X,2F12.4,5X,3E15.6) 
    10391842 
    1040      CALL wrk_dealloc(nb_type_class , zsumclass )   
     1843     CALL wrk_dealloc(nb_type , zsumclasses )   
    10411844  END SUBROUTINE dia_dct_wri 
    10421845 
  • branches/UKMO/CO5_package_branch/NEMOGCM/NEMO/OPA_SRC/DIA/diaptr.F90

    r7363 r7367  
    332332      !!---------------------------------------------------------------------- 
    333333      USE oce,     vt  =>   ua   ! use ua as workspace 
    334       USE oce,     vs  =>   ua   ! use ua as workspace 
     334      USE oce,     vs  =>   va   ! use va as workspace 
    335335      IMPLICIT none 
    336336      !! 
     
    378378                     zv = ( vn(ji,jj,jk) + vn(ji,jj-1,jk) ) * 0.5_wp 
    379379#endif  
    380                      vt(:,jj,jk) = zv * tsn(:,jj,jk,jp_tem) 
    381                      vs(:,jj,jk) = zv * tsn(:,jj,jk,jp_sal) 
     380                     vt(ji,jj,jk) = zv * tsn(ji,jj,jk,jp_tem) 
     381                     vs(ji,jj,jk) = zv * tsn(ji,jj,jk,jp_sal) 
    382382                  END DO 
    383383               END DO 
  • branches/UKMO/CO5_package_branch/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90

    r7363 r7367  
    4444   USE iom 
    4545   USE ioipsl 
     46   USE diafoam, ONLY: dia_wri_foam  
     47!CEOD   USE insitu_tem, ONLY: insitu_t, theta2t 
     48   USE bartrop_uv, ONLY: un_dm, vn_dm, bartrop_vel 
    4649#if defined key_lim2 
    4750   USE limwri_2  
     51   USE ice_2           ! LIM_2 ice model variables 
     52   USE lbclnk          ! ocean lateral boundary conditions (or mpp link) 
     53#endif 
     54#if defined key_lim3 
     55   USE ice_3           ! LIM_3 ice model variables 
     56   USE lbclnk          ! ocean lateral boundary conditions (or mpp link) 
     57#endif 
     58   USE daymod          ! calendar 
     59   USE insitu_tem, ONLY: insitu_t, theta2t 
     60#if defined key_top  
     61   USE par_trc         ! biogeochemical variables  
     62   USE trc 
     63#endif 
     64#if defined key_spm 
     65   USE spm_con, ONLY:  Eps0XS  
     66#endif  
     67#if defined key_zdftke  
     68   USE zdftke, ONLY: en 
     69#endif 
     70   USE zdf_oce, ONLY: avt, avm 
     71#if defined key_zdfgls 
     72   USE zdfgls, ONLY: mxln, en 
    4873#endif 
    4974   USE lib_mpp         ! MPP library 
     
    5479   PRIVATE 
    5580 
     81   PUBLIC   dia_wri_tmb_init        ! Called by nemogcm module 
    5682   PUBLIC   dia_wri                 ! routines called by step.F90 
    5783   PUBLIC   dia_wri_state 
    5884   PUBLIC   dia_wri_alloc           ! Called by nemogcm module 
     85   PUBLIC   dia_wri_tide_init       ! Called by nemogcm module 
    5986 
    6087   INTEGER ::   nid_T, nz_T, nh_T, ndim_T, ndim_hT   ! grid_T file 
     
    6592   INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_hT, ndex_hU, ndex_hV 
    6693   INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_T, ndex_U, ndex_V 
     94 
     95   !! * variables for calculating 25-hourly means 
     96   REAL(wp),SAVE, ALLOCATABLE,   DIMENSION(:,:,:) ::   tn_25h  , sn_25h  , insitu_t_25h 
     97   REAL(wp),SAVE, ALLOCATABLE,   DIMENSION(:,:)   ::   sshn_25h, hmld_kara_25h 
     98   REAL(wp),SAVE, ALLOCATABLE,   DIMENSION(:,:,:) ::   un_25h  , vn_25h  , wn_25h 
     99   REAL(wp),SAVE, ALLOCATABLE,   DIMENSION(:,:,:) ::   avt_25h , avm_25h 
     100#if defined key_zdfgls || key_zdftke 
     101   REAL(wp),SAVE, ALLOCATABLE,   DIMENSION(:,:,:) ::   en_25h   
     102#endif 
     103#if defined key_zdfgls  
     104   REAL(wp),SAVE, ALLOCATABLE,   DIMENSION(:,:,:) ::   mxln_25h 
     105#endif 
     106   INTEGER, SAVE :: cnt_25h     ! Counter for 25 hour means 
     107 
     108 
    67109 
    68110   !! * Substitutions 
     
    125167      REAL(wp), POINTER, DIMENSION(:,:)   :: z2d       ! 2D workspace 
    126168      REAL(wp), POINTER, DIMENSION(:,:,:) :: z3d      ! 3D workspace 
     169      REAL(wp), POINTER, DIMENSION(:,:,:) :: zwtmb    ! temporary workspace for tmb 
    127170      !!---------------------------------------------------------------------- 
    128171      !  
     
    131174      CALL wrk_alloc( jpi , jpj      , z2d ) 
    132175      CALL wrk_alloc( jpi , jpj, jpk , z3d ) 
     176      CALL wrk_alloc( jpi , jpj, 3 , zwtmb ) 
    133177      ! 
    134178      ! Output the initial state and forcings 
     
    138182      ENDIF 
    139183 
     184      IF (ln_diatide) THEN 
     185         CALL dia_wri_tide(kt) 
     186      ENDIF 
     187 
    140188      CALL iom_put( "toce"   , tsn(:,:,:,jp_tem)                     )    ! temperature 
     189      CALL theta2t ! in-situ temperature conversion 
     190!CEOD      CALL iom_put( "tinsitu", insitu_t(:,:,:)                       )    ! in-situ temperature 
    141191      CALL iom_put( "soce"   , tsn(:,:,:,jp_sal)                     )    ! salinity 
    142192      CALL iom_put( "sst"    , tsn(:,:,1,jp_tem)                     )    ! sea surface temperature 
     
    146196      CALL iom_put( "uoce"   , un                                    )    ! i-current       
    147197      CALL iom_put( "voce"   , vn                                    )    ! j-current 
    148        
     198      CALL iom_put( "ssu"    , un(:,:,1)                             )    ! sea surface U velocity 
     199      CALL iom_put( "ssv"    , vn(:,:,1)                             )    ! sea surface V velocity 
     200      IF( cp_cfg == "natl" .OR. cp_cfg == "ind12" ) CALL bartrop_vel ! barotropic velocity conversion 
     201!These don't exist independently in this branch so we remove them to get a CO5 
     202!that works on the Cray  
     203!CEOD      CALL iom_put( "uocebt" , un_dm                                 )    ! barotropic i-current  
     204!CEOD      CALL iom_put( "vocebt" , vn_dm                                 )    ! barotropic j-current 
    149205      CALL iom_put( "avt"    , avt                                   )    ! T vert. eddy diff. coef. 
    150206      CALL iom_put( "avm"    , avmu                                  )    ! T vert. eddy visc. coef. 
    151207      IF( lk_zdfddm ) THEN 
    152208         CALL iom_put( "avs" , fsavs(:,:,:)                          )    ! S vert. eddy diff. coef. 
     209      ENDIF 
     210      ! 
     211      ! If we want tmb values  
     212 
     213      IF (ln_diatmb) THEN 
     214         CALL dia_wri_calctmb(  tsn(:,:,:,jp_tem),zwtmb )  
     215         !ssh already output but here we output it masked 
     216         CALL iom_put( "sshnmasked" , sshn(:,:)*tmask(:,:,1)+missing_val*(1-tmask(:,:,1 ) ) )   ! tmb Temperature 
     217         CALL iom_put( "top_temp" , zwtmb(:,:,1) )    ! tmb Temperature 
     218         CALL iom_put( "mid_temp" , zwtmb(:,:,2) )    ! tmb Temperature 
     219         CALL iom_put( "bot_temp" , zwtmb(:,:,3) )    ! tmb Temperature 
     220!         CALL iom_put( "sotrefml" , hmld_tref(:,:) )    ! "T criterion Mixed Layer Depth 
     221 
     222         CALL dia_wri_calctmb(  tsn(:,:,:,jp_sal),zwtmb )  
     223         CALL iom_put( "top_sal" , zwtmb(:,:,1) )    ! tmb Salinity  
     224         CALL iom_put( "mid_sal" , zwtmb(:,:,2) )    ! tmb Salinity 
     225         CALL iom_put( "bot_sal" , zwtmb(:,:,3) )    ! tmb Salinity 
     226 
     227         CALL dia_wri_calctmb(  un(:,:,:),zwtmb )  
     228         CALL iom_put( "top_u" , zwtmb(:,:,1) )    ! tmb  U Velocity 
     229         CALL iom_put( "mid_u" , zwtmb(:,:,2) )    ! tmb  U Velocity 
     230         CALL iom_put( "bot_u" , zwtmb(:,:,3) )    ! tmb  U Velocity 
     231!Called in  dynspg_ts.F90        CALL iom_put( "baro_u" , un_b )    ! Barotropic  U Velocity 
     232 
     233         CALL dia_wri_calctmb(  vn(:,:,:),zwtmb )  
     234         CALL iom_put( "top_v" , zwtmb(:,:,1) )    ! tmb  V Velocity 
     235         CALL iom_put( "mid_v" , zwtmb(:,:,2) )    ! tmb  V Velocity 
     236         CALL iom_put( "bot_v" , zwtmb(:,:,3) )    ! tmb  V Velocity 
     237!Called in  dynspg_ts.F90       CALL iom_put( "baro_v" , vn_b )    ! Barotropic  V Velocity 
    153238      ENDIF 
    154239 
     
    171256         z3d(:,:,jpk) = 0.e0 
    172257         DO jk = 1, jpkm1 
    173             z3d(:,:,jk) = rau0 * un(:,:,jk) * e1u(:,:) * fse3u(:,:,jk) 
     258            z3d(:,:,jk) = rau0 * un(:,:,jk) * e2u(:,:) * fse3u(:,:,jk) 
    174259         END DO 
    175260         CALL iom_put( "u_masstr", z3d )                  ! mass transport in i-direction 
     
    186271         CALL iom_put( "u_heattr", z2d )                  ! heat transport in i-direction 
    187272         DO jk = 1, jpkm1 
    188             z3d(:,:,jk) = rau0 * vn(:,:,jk) * e2v(:,:) * fse3v(:,:,jk) 
     273            z3d(:,:,jk) = rau0 * vn(:,:,jk) * e1v(:,:) * fse3v(:,:,jk) 
    189274         END DO 
    190275         CALL iom_put( "v_masstr", z3d )                  ! mass transport in j-direction 
     
    251336      ENDIF 
    252337      ! 
     338      ! -1. Alternative routine  
     339      !------------------------  
     340      IF (ln_diafoam) THEN  
     341         CALL dia_wri_foam( kt )  
     342         RETURN  
     343      END IF  
     344      !  
    253345      ! 0. Initialisation 
    254346      ! ----------------- 
     
    673765#endif 
    674766 
     767   SUBROUTINE dia_wri_calctmb( infield,outtmb ) 
     768      !!--------------------------------------------------------------------- 
     769      !!                  ***  ROUTINE dia_tmb  *** 
     770      !!                    
     771      !! ** Purpose :   Write diagnostics for Top,Mid, and Bottom of water Column 
     772      !! 
     773      !! ** Method  :    
     774      !!      use mbathy to find surface, mid and bottom of model levels 
     775      !! 
     776      !! History : 
     777      !!   3.4  !  04-13  (E. O'Dea) Routine taken from old dia_wri_foam 
     778      !!---------------------------------------------------------------------- 
     779      !! * Modules used 
     780 
     781      ! Routine to map 3d field to top, middle, bottom 
     782      IMPLICIT NONE 
     783 
     784      ! Routine arguments 
     785      REAL(wp), DIMENSION(jpi, jpj, jpk), INTENT(IN   ) :: infield    ! Input 3d field and mask 
     786      REAL(wp), DIMENSION(jpi, jpj, 3  ), INTENT(  OUT) :: outtmb     ! Output top, middle, bottom 
     787 
     788      ! Local variables 
     789      INTEGER :: ji,jj,jk  ! Dummy loop indices 
     790 
     791      ! Calculate top 
     792      outtmb(:,:,1) = infield(:,:,1)*tmask(:,:,1) + missing_val*(1-tmask(:,:,1)) 
     793 
     794      ! Calculate middle 
     795      DO ji = 1,jpi 
     796         DO jj = 1,jpj 
     797            jk              = max(1,mbathy(ji,jj)/2) 
     798            outtmb(ji,jj,2) = infield(ji,jj,jk)*tmask(ji,jj,jk) + missing_val*(1-tmask(ji,jj,jk)) 
     799         END DO 
     800      END DO 
     801 
     802      ! Calculate bottom 
     803      DO ji = 1,jpi 
     804         DO jj = 1,jpj 
     805            jk              = max(1,mbathy(ji,jj) - 1) 
     806            outtmb(ji,jj,3) = infield(ji,jj,jk)*tmask(ji,jj,jk)  + missing_val*(1-tmask(ji,jj,jk)) 
     807         END DO 
     808      END DO 
     809 
     810   END SUBROUTINE dia_wri_calctmb 
     811 
     812   SUBROUTINE dia_wri_tmb_init 
     813      !!--------------------------------------------------------------------------- 
     814      !!                  ***  ROUTINE dia_wri_tmb_init  *** 
     815      !!      
     816      !! ** Purpose: Initialization of tmb namelist  
     817      !!         
     818      !! ** Method : Read namelist 
     819      !!   History 
     820      !!   3.4  !  04-13  (E. O'Dea) Routine to initialize dia_wri_tmb 
     821      !!--------------------------------------------------------------------------- 
     822      !! 
     823      INTEGER            ::   ierror   ! local integer 
     824      !! 
     825      NAMELIST/nam_diatmb/ ln_diatmb 
     826      !! 
     827      !!---------------------------------------------------------------------- 
     828      ! 
     829      REWIND ( numnam )              ! Read Namelist nam_diatmb 
     830      READ   ( numnam, nam_diatmb ) 
     831      ! 
     832      IF(lwp) THEN                   ! Control print 
     833         WRITE(numout,*) 
     834         WRITE(numout,*) 'dia_wri_tmb_init : Output Top, Middle, Bottom Diagnostics' 
     835         WRITE(numout,*) '~~~~~~~~~~~~' 
     836         WRITE(numout,*) '   Namelist nam_diatmb : set tmb outputs ' 
     837         WRITE(numout,*) '      Switch for TMB diagnostics (T) or not (F)  ln_diatmb  = ', ln_diatmb 
     838      ENDIF 
     839 
     840    END SUBROUTINE dia_wri_tmb_init 
     841 
     842 
    675843   SUBROUTINE dia_wri_state( cdfile_name, kt ) 
    676844      !!--------------------------------------------------------------------- 
     
    798966   END SUBROUTINE dia_wri_state 
    799967   !!====================================================================== 
     968   !!====================================================================== 
     969 
     970   SUBROUTINE dia_wri_tide( kt ) 
     971      !!--------------------------------------------------------------------- 
     972      !!                  ***  ROUTINE dia_tide  *** 
     973      !!                    
     974      !! ** Purpose :   Write diagnostics with M2/S2 tide removed 
     975      !! 
     976      !! ** Method  :    
     977      !!      25hr mean outputs for shelf seas 
     978      !! 
     979      !! History : 
     980      !!   ?.0  !  07-04  (A. Hines) New routine, developed from dia_wri_foam 
     981      !!   3.4  !  02-13  (J. Siddorn) Routine taken from old dia_wri_foam 
     982      !!---------------------------------------------------------------------- 
     983      !! * Modules used 
     984 
     985      IMPLICIT NONE 
     986 
     987      !! * Arguments 
     988      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index 
     989 
     990 
     991      !! * Local declarations 
     992      INTEGER ::   ji, jj, jk 
     993 
     994      LOGICAL ::   ll_print = .FALSE.    ! =T print and flush numout 
     995      REAL(wp)                         ::   zsto, zout, zmax, zjulian, zdt, zmdi  ! temporary integers 
     996      INTEGER                          ::   i_steps                               ! no of timesteps per hour 
     997      REAL(wp), DIMENSION(jpi,jpj    ) ::   zw2d, un_dm, vn_dm                    ! temporary workspace 
     998      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zw3d                                  ! temporary workspace 
     999      REAL(wp), DIMENSION(jpi,jpj,3)   ::   zwtmb                                 ! temporary workspace 
     1000      INTEGER                          ::   nyear0, nmonth0,nday0                 ! start year,month,day 
     1001!#if defined key_top 
     1002!      CHARACTER (len=20) :: cltra, cltrau 
     1003!      CHARACTER (len=80) :: cltral 
     1004!      INTEGER            :: jn, jl 
     1005!#endif 
     1006!#if defined key_spm   
     1007!      ! variables needed to calculate visibility field from sediment fields 
     1008!      REAL(wp), DIMENSION(jpi,jpj,jpk) :: vis3d    ! derived 3D visibility field 
     1009!      REAL(wp) :: epsessX = 0.07d-03               ! attenuation coefficient applied to the sediment (as used in ERSEM) 
     1010!      REAL(wp) :: tiny = 1.0d-15                   ! to prevent division by zero in visibility calculation 
     1011!#endif    
     1012         
     1013      !!---------------------------------------------------------------------- 
     1014      
     1015      ! 0. Initialisation 
     1016      ! ----------------- 
     1017      ! Define frequency of summing to create 25 h mean 
     1018      zdt = rdt 
     1019      IF( nacc == 1 ) zdt = rdtmin 
     1020       
     1021      IF( MOD( 3600,INT(zdt) ) == 0 ) THEN 
     1022         i_steps = 3600/INT(zdt) 
     1023      ELSE 
     1024         CALL ctl_stop('STOP', 'dia_wri_tide: timestep must give MOD(3600,rdt) = 0 otherwise no hourly values are possible') 
     1025      ENDIF  
     1026           
     1027#if defined key_lim3 || defined key_lim2 
     1028      CALL ctl_stop('STOP', 'dia_wri_tide not setup yet to do tidemean ice') 
     1029#endif 
     1030#if defined key_spm || defined key_MOersem 
     1031      CALL ctl_stop('STOP', 'dia_wri_tide not setup yet to do tidemean ERSEM') 
     1032#endif 
     1033                
     1034      ! local variable for debugging 
     1035      ll_print = ll_print .AND. lwp 
     1036 
     1037      ! Sum of 25 hourly instantaneous values to give a 25h mean from 24hours 
     1038      ! every day 
     1039      IF( MOD( kt, i_steps ) == 0  .and. kt .ne. nn_it000 ) THEN 
     1040 
     1041         IF (lwp) THEN 
     1042              WRITE(numout,*) 'dia_wri_tide : Summing instantaneous hourly diagnostics at timestep ',kt 
     1043              WRITE(numout,*) '~~~~~~~~~~~~ ' 
     1044         ENDIF 
     1045 
     1046         tn_25h(:,:,:)        = tn_25h(:,:,:) + tsn(:,:,:,jp_tem) 
     1047         sn_25h(:,:,:)        = sn_25h(:,:,:) + tsn(:,:,:,jp_sal) 
     1048         CALL theta2t 
     1049         insitu_t_25h(:,:,:)  = insitu_t_25h(:,:,:) + insitu_t(:,:,:) 
     1050         sshn_25h(:,:)        = sshn_25h(:,:) + sshn (:,:) 
     1051!         hmld_kara_25h(:,:)   = hmld_kara_25h(:,:) + hmld_kara(:,:) 
     1052#if defined key_lim3 || defined key_lim2 
     1053         hsnif_25h(:,:)       = hsnif_25h(:,:) + hsnif(:,:) 
     1054         hicif_25h(:,:)       = hicif_25h(:,:) + hicif(:,:) 
     1055         frld_25h(:,:)        = frld_25h(:,:) + frld(:,:) 
     1056#endif  
     1057#if defined key_spm || defined key_MOersem 
     1058         trn_25h(:,:,:,:)     = trn_25h(:,:,:,:) + trn (:,:,:,:) 
     1059         trc3d_25h(:,:,:,:)   = trc3d_25h(:,:,:,:) + trc3d(:,:,:,:) 
     1060         trc2d_25h(:,:,:)     = trc2d_25h(:,:,:) + trc2d(:,:,:) 
     1061#endif 
     1062         un_25h(:,:,:)        = un_25h(:,:,:) + un(:,:,:) 
     1063         vn_25h(:,:,:)        = vn_25h(:,:,:) + vn(:,:,:) 
     1064         wn_25h(:,:,:)        = wn_25h(:,:,:) + wn(:,:,:) 
     1065         avt_25h(:,:,:)       = avt_25h(:,:,:) + avt(:,:,:) 
     1066         avm_25h(:,:,:)       = avm_25h(:,:,:) + avm(:,:,:) 
     1067# if defined key_zdfgls || defined key_zdftke 
     1068         en_25h(:,:,:)        = en_25h(:,:,:) + en(:,:,:) 
     1069#endif 
     1070# if defined key_zdfgls 
     1071         mxln_25h(:,:,:)      = mxln_25h(:,:,:) + mxln(:,:,:) 
     1072#endif 
     1073         cnt_25h = cnt_25h + 1 
     1074 
     1075         IF (lwp) THEN 
     1076            WRITE(numout,*) 'dia_tide : Summed the following number of hourly values so far',cnt_25h 
     1077         ENDIF 
     1078 
     1079      ENDIF ! MOD( kt, i_steps ) == 0 
     1080 
     1081         ! Write data for 25 hour mean output streams 
     1082      IF( cnt_25h .EQ. 25 .AND.  MOD( kt, i_steps*24) == 0 .AND. kt .NE. nn_it000 ) THEN 
     1083 
     1084            IF(lwp) THEN 
     1085               WRITE(numout,*) 'dia_wri_tide : Writing 25 hour mean tide diagnostics at timestep', kt 
     1086               WRITE(numout,*) '~~~~~~~~~~~~ ' 
     1087            ENDIF 
     1088 
     1089            tn_25h(:,:,:)        = tn_25h(:,:,:) / 25.0_wp 
     1090            sn_25h(:,:,:)        = sn_25h(:,:,:) / 25.0_wp 
     1091            insitu_t_25h(:,:,:)  = insitu_t_25h(:,:,:) / 25.0_wp 
     1092            sshn_25h(:,:)        = sshn_25h(:,:) / 25.0_wp 
     1093!            hmld_kara_25h(:,:)   = hmld_kara_25h(:,:) / 25.0_wp 
     1094#if defined key_lim3 || defined key_lim2 
     1095            hsnif_25h(:,:)       = hsnif_25h(:,:) / 25.0_wp 
     1096            hicif_25h(:,:)       = hicif_25h(:,:) / 25.0_wp 
     1097            frld_25h(:,:)        = frld_25h(:,:) / 25.0_wp 
     1098#endif  
     1099#if defined key_spm || defined key_MOersem 
     1100            trn_25h(:,:,:,:)     = trn_25h(:,:,:,:) / 25.0_wp 
     1101            trc3d_25h(:,:,:,:)   = trc3d_25h(:,:,:,:) / 25.0_wp 
     1102            trc2d_25h(:,:,:)     = trc2d_25h(:,:,:) / 25.0_wp 
     1103#endif 
     1104            un_25h(:,:,:)        = un_25h(:,:,:) / 25.0_wp 
     1105            vn_25h(:,:,:)        = vn_25h(:,:,:) / 25.0_wp 
     1106            wn_25h(:,:,:)        = wn_25h(:,:,:) / 25.0_wp 
     1107            avt_25h(:,:,:)       = avt_25h(:,:,:) / 25.0_wp 
     1108            avm_25h(:,:,:)       = avm_25h(:,:,:) / 25.0_wp 
     1109# if defined key_zdfgls || defined key_zdftke 
     1110            en_25h(:,:,:)        = en_25h(:,:,:) / 25.0_wp 
     1111#endif 
     1112# if defined key_zdfgls 
     1113            mxln_25h(:,:,:)       = mxln_25h(:,:,:) / 25.0_wp 
     1114#endif 
     1115 
     1116            IF (lwp)  WRITE(numout,*) 'dia_wri_tide : Mean calculated by dividing 25 hour sums and writing output'  
     1117  
     1118            zmdi=missing_val                 ! for masking 
     1119            ! write tracers (instantaneous) 
     1120            zw3d(:,:,:) = tn_25h(:,:,:)*tmask(:,:,:) + zmdi*(1.0-tmask(:,:,:)) 
     1121            CALL iom_put("temper25h", zw3d)   ! potential temperature 
     1122            CALL theta2t                                                                    ! calculate insitu temp 
     1123            zw3d(:,:,:) = insitu_t_25h(:,:,:)*tmask(:,:,:) + zmdi*(1.0-tmask(:,:,:)) 
     1124            CALL iom_put("tempis25h", zw3d)   ! in-situ temperature 
     1125            zw3d(:,:,:) = sn_25h(:,:,:)*tmask(:,:,:) + zmdi*(1.0-tmask(:,:,:)) 
     1126            CALL iom_put( "salin25h", zw3d  )   ! salinity 
     1127            zw2d(:,:) = sshn_25h(:,:)*tmask(:,:,1) + zmdi*(1.0-tmask(:,:,1)) 
     1128            CALL iom_put( "ssh25h", zw2d )   ! sea surface  
     1129!            zw2d(:,:) = hmld_kara_25h(:,:)*tmask(:,:,1) + zmdi*(1.0-tmask(:,:,1)) 
     1130!            CALL iom_put( "kara25h", zw2d )   ! mixed layer  
     1131 
     1132#if defined key_lim3 || defined key_lim2 
     1133            ! Write ice model variables (instantaneous) 
     1134            zw2d(:,:) = hsnif_25h(:,:)*tmask(:,:,1) + zmdi*(1.0-tmask(:,:,1)) 
     1135            CALL iom_put("isnowthi", zw2d )   ! ice thickness 
     1136            zw2d(:,:) = hicif_25h(:,:)*tmask(:,:,1) + zmdi*(1.0-tmask(:,:,1)) 
     1137            CALL iom_put("iicethic", zw2d )   ! ice thickness 
     1138            zw2d(:,:) = (1.0-frld_25h(:,:))*tmask(:,:,1) + zmdi*(1.0-tmask(:,:,1)) 
     1139            CALL iom_put("iiceconc", zw2d )   ! ice concetration 
     1140#endif 
     1141#if defined key_spm || key_MOersem 
     1142            ! output biogeochemical variables: 
     1143            ! output main tracers 
     1144            DO jn = 1, jptra 
     1145               cltra = ctrcnm(jn)      ! short title for tracer 
     1146               zw3d(:,:,:) = trn_25h(:,:,:,jn)*tmask(:,:,:) + zmdi*(1.0-tmask(:,:,:)) 
     1147               IF( lutsav(jn) )  CALL iom_put( cltra, zw3d  )   ! temperature 
     1148            END DO 
     1149            ! more 3D horizontal arrays from diagnostics 
     1150            DO jl = 1, jpdia3d 
     1151               cltra = ctrc3d(jl)   ! short title for 3D diagnostic 
     1152               zw3d(:,:,:) = trc3d_25h(:,:,:,jl)*tmask(:,:,:) + zmdi*(1.0-tmask(:,:,:)) 
     1153               CALL iom_put( cltra, zw3d  )    
     1154            END DO 
     1155            ! more 2D horizontal arrays from diagnostics 
     1156            DO jl = 1, jpdia2d 
     1157               cltra = ctrc2d(jl)   ! short title for 2D diagnostic 
     1158               zw2d(:,:) = trc2d_25h(:,:,jl)*tmask(:,:,1) + zmdi*(1.0-tmask(:,:,1)) 
     1159               CALL iom_put(cltra, zw2d ) 
     1160            END DO 
     1161#endif 
     1162 
     1163            ! Write velocities (instantaneous) 
     1164            zw3d(:,:,:) = un_25h(:,:,:)*umask(:,:,:) + zmdi*(1.0-umask(:,:,:)) 
     1165            CALL iom_put("vozocrtx25h", zw3d)    ! i-current 
     1166            zw3d(:,:,:) = vn_25h(:,:,:)*vmask(:,:,:) + zmdi*(1.0-vmask(:,:,:)) 
     1167            CALL iom_put("vomecrty25h", zw3d  )   ! j-current 
     1168 
     1169            zw3d(:,:,:) = wn_25h(:,:,:)*tmask(:,:,:) + zmdi*(1.0-tmask(:,:,:)) 
     1170            CALL iom_put("vomecrtz25h", zw3d )   ! k-current 
     1171            zw3d(:,:,:) = avt_25h(:,:,:)*tmask(:,:,:) + zmdi*(1.0-tmask(:,:,:)) 
     1172            CALL iom_put("avt25h", zw3d )   ! diffusivity 
     1173            zw3d(:,:,:) = avm_25h(:,:,:)*tmask(:,:,:) + zmdi*(1.0-tmask(:,:,:)) 
     1174            CALL iom_put("avm25h", zw3d)   ! viscosity 
     1175#if defined key_zdftke || defined key_zdfgls  
     1176            zw3d(:,:,:) = en_25h(:,:,:)*tmask(:,:,:) + zmdi*(1.0-tmask(:,:,:)) 
     1177            CALL iom_put("tke25h", zw3d)   ! tke 
     1178#endif 
     1179#if defined key_zdfgls  
     1180            zw3d(:,:,:) = mxln_25h(:,:,:)*tmask(:,:,:) + zmdi*(1.0-tmask(:,:,:)) 
     1181            CALL iom_put( "mxln25h",zw3d)  
     1182#endif 
     1183 
     1184            ! After the write reset the values to cnt=1 and sum values equal current value  
     1185            tn_25h(:,:,:) = tsn(:,:,:,jp_tem) 
     1186            sn_25h(:,:,:) = tsn(:,:,:,jp_sal) 
     1187            CALL theta2t 
     1188            insitu_t_25h(:,:,:) = insitu_t(:,:,:) 
     1189            sshn_25h(:,:) = sshn (:,:) 
     1190!            hmld_kara_25h(:,:) = hmld_kara(:,:) 
     1191#if defined key_lim3 || defined key_lim2 
     1192            hsnif_25h(:,:) = hsnif(:,:) 
     1193            hicif_25h(:,:) = hicif(:,:) 
     1194            frld_25h(:,:) = frld(:,:) 
     1195#endif  
     1196#if defined key_spm || defined key_MOersem 
     1197            trn_25h(:,:,:,:) = trn (:,:,:,:) 
     1198            trc3d_25h(:,:,:,:) = trc3d(:,:,:,:) 
     1199            trc2d_25h(:,:,:) = trc2d(:,:,:) 
     1200#endif 
     1201            un_25h(:,:,:) = un(:,:,:) 
     1202            vn_25h(:,:,:) = vn(:,:,:) 
     1203            wn_25h(:,:,:) = wn(:,:,:) 
     1204            avt_25h(:,:,:) = avt(:,:,:) 
     1205            avm_25h(:,:,:) = avm(:,:,:) 
     1206# if defined key_zdfgls || defined key_zdftke 
     1207            en_25h(:,:,:) = en(:,:,:) 
     1208#endif 
     1209# if defined key_zdfgls 
     1210            mxln_25h(:,:,:) = mxln(:,:,:) 
     1211#endif 
     1212            cnt_25h = 1 
     1213            IF (lwp)  WRITE(numout,*) 'dia_wri_tide : After 25hr mean write, reset sum to current value and cnt_25h to one for overlapping average',cnt_25h 
     1214 
     1215      ENDIF !  cnt_25h .EQ. 25 .AND.  MOD( kt, i_steps * 24) == 0 .AND. kt .NE. nn_it000 
     1216 
     1217   END SUBROUTINE dia_wri_tide 
     1218   !!====================================================================== 
     1219 
     1220   SUBROUTINE dia_wri_tide_init 
     1221      !!--------------------------------------------------------------------------- 
     1222      !!                  ***  ROUTINE dia_wri_tide_init  *** 
     1223      !!      
     1224      !! ** Purpose: Initialization of 25hour mean variables for detided output   
     1225      !!         
     1226      !! ** Method : Read namelist, allocate and assign initial values 
     1227      !!   History 
     1228      !!   3.4  !  03-13  (E. O'Dea) Routine to initialize dia_wri_tide   
     1229      !!--------------------------------------------------------------------------- 
     1230      !! 
     1231      INTEGER            ::   ierror   ! local integer 
     1232      !! 
     1233      NAMELIST/nam_diatide/ ln_diatide 
     1234      !! 
     1235      !!---------------------------------------------------------------------- 
     1236      ! 
     1237      REWIND ( numnam )              ! Read Namelist nam_tiatide 
     1238      READ   ( numnam, nam_diatide ) 
     1239      ! 
     1240      IF(lwp) THEN                   ! Control print 
     1241         WRITE(numout,*) 
     1242         WRITE(numout,*) 'dia_wri_tide_init : Output 25 hour Mean Diagnostics' 
     1243         WRITE(numout,*) '~~~~~~~~~~~~' 
     1244         WRITE(numout,*) '   Namelist nam_diatide : set 25hour mean outputs ' 
     1245         WRITE(numout,*) '      Switch for 25 hour mean diagnostics (T) or not (F)  ln_diatide  = ', ln_diatide 
     1246      ENDIF 
     1247      IF( .NOT. ln_diatide )   RETURN 
     1248 
     1249      ! ------------------- ! 
     1250      ! 1 - Allocate memory ! 
     1251      ! ------------------- ! 
     1252      ALLOCATE( tn_25h(jpi,jpj,jpk), STAT=ierror ) 
     1253      IF( ierror > 0 ) THEN 
     1254         CALL ctl_stop( 'dia_tide: unable to allocate tn_25h' )   ;   RETURN 
     1255      ENDIF 
     1256      ALLOCATE( sn_25h(jpi,jpj,jpk), STAT=ierror ) 
     1257      IF( ierror > 0 ) THEN 
     1258         CALL ctl_stop( 'dia_tide: unable to allocate sn_25h' )   ;   RETURN 
     1259      ENDIF 
     1260      ALLOCATE( insitu_t_25h(jpi,jpj,jpk), STAT=ierror ) 
     1261      IF( ierror > 0 ) THEN 
     1262         CALL ctl_stop( 'dia_tide: unable to allocate insitu_t_25h' )   ;   RETURN 
     1263      ENDIF 
     1264      ALLOCATE( un_25h(jpi,jpj,jpk), STAT=ierror ) 
     1265      IF( ierror > 0 ) THEN 
     1266         CALL ctl_stop( 'dia_tide: unable to allocate un_25h' )   ;   RETURN 
     1267      ENDIF 
     1268      ALLOCATE( vn_25h(jpi,jpj,jpk), STAT=ierror ) 
     1269      IF( ierror > 0 ) THEN 
     1270         CALL ctl_stop( 'dia_tide: unable to allocate vn_25h' )   ;   RETURN 
     1271      ENDIF 
     1272      ALLOCATE( wn_25h(jpi,jpj,jpk), STAT=ierror ) 
     1273      IF( ierror > 0 ) THEN 
     1274         CALL ctl_stop( 'dia_tide: unable to allocate wn_25h' )   ;   RETURN 
     1275      ENDIF 
     1276      ALLOCATE( avt_25h(jpi,jpj,jpk), STAT=ierror ) 
     1277      IF( ierror > 0 ) THEN 
     1278         CALL ctl_stop( 'dia_tide: unable to allocate avt_25h' )   ;   RETURN 
     1279      ENDIF 
     1280      ALLOCATE( avm_25h(jpi,jpj,jpk), STAT=ierror ) 
     1281      IF( ierror > 0 ) THEN 
     1282         CALL ctl_stop( 'dia_tide: unable to allocate avm_25h' )   ;   RETURN 
     1283      ENDIF 
     1284# if defined key_zdfgls || defined key_zdftke 
     1285      ALLOCATE( en_25h(jpi,jpj,jpk), STAT=ierror ) 
     1286      IF( ierror > 0 ) THEN 
     1287         CALL ctl_stop( 'dia_tide: unable to allocate en_25h' )   ;   RETURN 
     1288      ENDIF 
     1289#endif 
     1290# if defined key_zdfgls  
     1291      ALLOCATE( mxln_25h(jpi,jpj,jpk), STAT=ierror ) 
     1292      IF( ierror > 0 ) THEN 
     1293         CALL ctl_stop( 'dia_tide: unable to allocate mxln_25h' )   ;   RETURN 
     1294      ENDIF 
     1295#endif 
     1296      ALLOCATE( sshn_25h(jpi,jpj), STAT=ierror ) 
     1297      IF( ierror > 0 ) THEN 
     1298         CALL ctl_stop( 'dia_tide: unable to allocate sshn_25h' )   ;   RETURN 
     1299      ENDIF 
     1300      ALLOCATE( hmld_kara_25h(jpi,jpj), STAT=ierror ) 
     1301      IF( ierror > 0 ) THEN 
     1302         CALL ctl_stop( 'dia_tide: unable to allocate hmld_kara_25h' )   ;   RETURN 
     1303      ENDIF 
     1304      ! ------------------------- ! 
     1305      ! 2 - Assign Initial Values ! 
     1306      ! ------------------------- ! 
     1307      cnt_25h = 1  ! sets the first value of sum at timestep 1 (note - should strictly be at timestep zero so before values used where possible)  
     1308      tn_25h(:,:,:) = tsb(:,:,:,jp_tem) 
     1309      sn_25h(:,:,:) = tsb(:,:,:,jp_sal) 
     1310      CALL theta2t 
     1311      insitu_t_25h(:,:,:) = insitu_t(:,:,:) 
     1312      sshn_25h(:,:) = sshb(:,:) 
     1313!         hmld_kara_25h(:,:) = hmld_kara(:,:) 
     1314      un_25h(:,:,:) = ub(:,:,:) 
     1315      vn_25h(:,:,:) = vb(:,:,:) 
     1316      wn_25h(:,:,:) = wn(:,:,:) 
     1317      avt_25h(:,:,:) = avt(:,:,:) 
     1318      avm_25h(:,:,:) = avm(:,:,:) 
     1319# if defined key_zdfgls || defined key_zdftke 
     1320         en_25h(:,:,:) = en(:,:,:) 
     1321#endif 
     1322# if defined key_zdfgls 
     1323         mxln_25h(:,:,:) = mxln(:,:,:) 
     1324#endif 
     1325#if defined key_lim3 || defined key_lim2 
     1326         CALL ctl_stop('STOP', 'dia_wri_tide not setup yet to do tidemean ice') 
     1327#endif  
     1328 
     1329      ! -------------------------- ! 
     1330      ! 3 - Return to dia_wri_tide ! 
     1331      ! -------------------------- ! 
     1332 
     1333 
     1334    END SUBROUTINE dia_wri_tide_init 
     1335 
     1336 
    8001337END MODULE diawri 
  • branches/UKMO/CO5_package_branch/NEMOGCM/NEMO/OPA_SRC/DOM/closea.F90

    r7363 r7367  
    77   !!             8.5  !  02-06  (E. Durand, G. Madec)  F90 
    88   !!             9.0  !  06-07  (G. Madec)  add clo_rnf, clo_ups, clo_bat 
     9   !!        NEMO 3.4  !  03-12  (P.G. Fogli) sbc_clo bug fix & mpp reproducibility 
    910   !!---------------------------------------------------------------------- 
    1011 
     
    2021   USE in_out_manager  ! I/O manager 
    2122   USE sbc_oce         ! ocean surface boundary conditions 
    22    USE lib_mpp         ! distributed memory computing library 
    23    USE lbclnk          ! ??? 
     23   USE lib_fortran,    ONLY: glob_sum, DDPDD 
     24   USE lbclnk          ! lateral boundary condition - MPP exchanges 
     25   USE lib_mpp         ! MPP library 
     26   USE timing 
    2427 
    2528   IMPLICIT NONE 
     
    8588         SELECT CASE ( jp_cfg ) 
    8689         !                                           ! ======================= 
     90         CASE ( 1 )                                  ! ORCA_R1 configuration 
     91            !                                        ! ======================= 
     92            ncsnr(1)   = 1    ; ncstt(1)   = 0           ! Caspian Sea 
     93            ncsi1(1)   = 332  ; ncsj1(1)   = 203 
     94            ncsi2(1)   = 344  ; ncsj2(1)   = 235 
     95            ncsir(1,1) = 1    ; ncsjr(1,1) = 1 
     96            !                                         
     97            !                                        ! ======================= 
    8798         CASE ( 2 )                                  !  ORCA_R2 configuration 
    8899            !                                        ! ======================= 
     
    177188      INTEGER, INTENT(in) ::   kt   ! ocean model time step 
    178189      ! 
    179       INTEGER                     ::   ji, jj, jc, jn   ! dummy loop indices 
    180       REAL(wp)                    ::   zze2 
    181       REAL(wp), DIMENSION (jpncs) ::   zfwf  
    182       !!---------------------------------------------------------------------- 
    183       ! 
     190      INTEGER             ::   ji, jj, jc, jn   ! dummy loop indices 
     191      REAL(wp), PARAMETER ::   rsmall = 1.e-20_wp    ! Closed sea correction epsilon 
     192      REAL(wp)            ::   zze2, ztmp, zcorr     !  
     193      COMPLEX(wp)         ::   ctmp  
     194      REAL(wp), DIMENSION(jpncs) ::   zfwf   ! 1D workspace 
     195      !!---------------------------------------------------------------------- 
     196      ! 
     197      IF( nn_timing == 1 )  CALL timing_start('sbc_clo') 
    184198      !                                                   !------------------! 
    185199      IF( kt == nit000 ) THEN                             !  Initialisation  ! 
     
    189203         IF(lwp) WRITE(numout,*)'~~~~~~~' 
    190204 
    191          ! Total surface of ocean 
    192          surf(jpncs+1) = SUM( e1t(:,:) * e2t(:,:) * tmask_i(:,:) ) 
    193  
    194          DO jc = 1, jpncs 
    195             surf(jc) =0.e0 
    196             DO jj = ncsj1(jc), ncsj2(jc) 
    197                DO ji = ncsi1(jc), ncsi2(jc) 
    198                   surf(jc) = surf(jc) + e1t(ji,jj) * e2t(ji,jj) * tmask_i(ji,jj)      ! surface of closed seas 
     205         surf(:) = 0.e0_wp 
     206         ! 
     207         surf(jpncs+1) = glob_sum( e1e2t(:,:) )   ! surface of the global ocean 
     208         ! 
     209         !                                        ! surface of closed seas  
     210         IF( lk_mpp_rep ) THEN                         ! MPP reproductible calculation 
     211            DO jc = 1, jpncs 
     212               ctmp = CMPLX( 0.e0, 0.e0, wp ) 
     213               DO jj = ncsj1(jc), ncsj2(jc) 
     214                  DO ji = ncsi1(jc), ncsi2(jc) 
     215                     ztmp = e1e2t(ji,jj) * tmask_i(ji,jj) 
     216                     CALL DDPDD( CMPLX( ztmp, 0.e0, wp ), ctmp ) 
     217                  END DO  
    199218               END DO  
    200             END DO  
    201          END DO  
    202          IF( lk_mpp )   CALL mpp_sum ( surf, jpncs+1 )       ! mpp: sum over all the global domain 
     219               IF( lk_mpp )   CALL mpp_sum( ctmp ) 
     220               surf(jc) = REAL(ctmp,wp) 
     221            END DO 
     222         ELSE                                          ! Standard calculation            
     223            DO jc = 1, jpncs 
     224               DO jj = ncsj1(jc), ncsj2(jc) 
     225                  DO ji = ncsi1(jc), ncsi2(jc) 
     226                     surf(jc) = surf(jc) + e1e2t(ji,jj) * tmask_i(ji,jj)      ! surface of closed seas 
     227                  END DO  
     228               END DO  
     229            END DO  
     230            IF( lk_mpp )   CALL mpp_sum ( surf, jpncs )       ! mpp: sum over all the global domain 
     231         ENDIF 
    203232 
    204233         IF(lwp) WRITE(numout,*)'     Closed sea surfaces' 
     
    215244      !                                                   !--------------------! 
    216245      !                                                   !  update emp, emps  ! 
    217       zfwf = 0.e0                                         !--------------------! 
    218       DO jc = 1, jpncs 
    219          DO jj = ncsj1(jc), ncsj2(jc) 
    220             DO ji = ncsi1(jc), ncsi2(jc) 
    221                zfwf(jc) = zfwf(jc) + e1t(ji,jj) * e2t(ji,jj) * ( emp(ji,jj)-rnf(ji,jj) ) * tmask_i(ji,jj)  
    222             END DO   
    223          END DO  
    224       END DO 
    225       IF( lk_mpp )   CALL mpp_sum ( zfwf(:) , jpncs )       ! mpp: sum over all the global domain 
     246      zfwf = 0.e0_wp                                      !--------------------! 
     247      IF( lk_mpp_rep ) THEN                         ! MPP reproductible calculation 
     248         DO jc = 1, jpncs 
     249            ctmp = CMPLX( 0.e0, 0.e0, wp ) 
     250            DO jj = ncsj1(jc), ncsj2(jc) 
     251               DO ji = ncsi1(jc), ncsi2(jc) 
     252                  ztmp = e1e2t(ji,jj) * ( emp(ji,jj)-rnf(ji,jj) ) * tmask_i(ji,jj) 
     253                  CALL DDPDD( CMPLX( ztmp, 0.e0, wp ), ctmp ) 
     254               END DO   
     255            END DO  
     256            IF( lk_mpp )   CALL mpp_sum( ctmp ) 
     257            zfwf(jc) = REAL(ctmp,wp) 
     258         END DO 
     259      ELSE                                          ! Standard calculation            
     260         DO jc = 1, jpncs 
     261            DO jj = ncsj1(jc), ncsj2(jc) 
     262               DO ji = ncsi1(jc), ncsi2(jc) 
     263                  zfwf(jc) = zfwf(jc) + e1e2t(ji,jj) * ( emp(ji,jj)-rnf(ji,jj) ) * tmask_i(ji,jj)  
     264               END DO   
     265            END DO  
     266         END DO 
     267         IF( lk_mpp )   CALL mpp_sum ( zfwf(:) , jpncs )       ! mpp: sum over all the global domain 
     268      ENDIF 
    226269 
    227270      IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN      ! Black Sea case for ORCA_R2 configuration 
    228          zze2    = ( zfwf(3) + zfwf(4) ) / 2. 
     271         zze2    = ( zfwf(3) + zfwf(4) ) * 0.5_wp 
    229272         zfwf(3) = zze2 
    230273         zfwf(4) = zze2 
    231274      ENDIF 
    232275 
     276      zcorr = 0._wp 
     277 
    233278      DO jc = 1, jpncs 
    234279         ! 
    235          IF( ncstt(jc) == 0 ) THEN  
    236             ! water/evap excess is shared by all open ocean 
    237             emp (:,:) = emp (:,:) + zfwf(jc) / surf(jpncs+1) 
    238             emps(:,:) = emps(:,:) + zfwf(jc) / surf(jpncs+1) 
    239          ELSEIF( ncstt(jc) == 1 ) THEN  
    240             ! Excess water in open sea, at outflow location, excess evap shared 
    241             IF ( zfwf(jc) <= 0.e0 ) THEN  
    242                 DO jn = 1, ncsnr(jc) 
     280         ! The following if avoids the redistribution of the round off 
     281         IF ( ABS(zfwf(jc) / surf(jpncs+1) ) > rsmall) THEN 
     282            ! 
     283            IF( ncstt(jc) == 0 ) THEN           ! water/evap excess is shared by all open ocean 
     284               emp (:,:) = emp (:,:) + zfwf(jc) / surf(jpncs+1) 
     285               emps(:,:) = emps(:,:) + zfwf(jc) / surf(jpncs+1) 
     286               ! accumulate closed seas correction 
     287               zcorr     = zcorr     + zfwf(jc) / surf(jpncs+1) 
     288               ! 
     289            ELSEIF( ncstt(jc) == 1 ) THEN       ! Excess water in open sea, at outflow location, excess evap shared 
     290               IF ( zfwf(jc) <= 0.e0_wp ) THEN  
     291                   DO jn = 1, ncsnr(jc) 
     292                     ji = mi0(ncsir(jc,jn)) 
     293                     jj = mj0(ncsjr(jc,jn)) ! Location of outflow in open ocean 
     294                     IF (      ji > 1 .AND. ji < jpi   & 
     295                         .AND. jj > 1 .AND. jj < jpj ) THEN  
     296                         emp (ji,jj) = emp (ji,jj) + zfwf(jc) / ( REAL(ncsnr(jc)) * e1e2t(ji,jj) ) 
     297                         emps(ji,jj) = emps(ji,jj) + zfwf(jc) / ( REAL(ncsnr(jc)) * e1e2t(ji,jj) ) 
     298                     ENDIF  
     299                   END DO  
     300               ELSE  
     301                   emp (:,:) = emp (:,:) + zfwf(jc) / surf(jpncs+1) 
     302                   emps(:,:) = emps(:,:) + zfwf(jc) / surf(jpncs+1) 
     303                   ! accumulate closed seas correction 
     304                   zcorr     = zcorr     + zfwf(jc) / surf(jpncs+1) 
     305               ENDIF 
     306            ELSEIF( ncstt(jc) == 2 ) THEN       ! Excess e-p-r (either sign) goes to open ocean, at outflow location 
     307               DO jn = 1, ncsnr(jc) 
    243308                  ji = mi0(ncsir(jc,jn)) 
    244309                  jj = mj0(ncsjr(jc,jn)) ! Location of outflow in open ocean 
    245                   IF (      ji > 1 .AND. ji < jpi   & 
    246                       .AND. jj > 1 .AND. jj < jpj ) THEN  
    247                       emp (ji,jj) = emp (ji,jj) + zfwf(jc) /   & 
    248                          (FLOAT(ncsnr(jc)) * e1t(ji,jj) * e2t(ji,jj)) 
    249                       emps(ji,jj) = emps(ji,jj) + zfwf(jc) /   & 
    250                           (FLOAT(ncsnr(jc)) * e1t(ji,jj) * e2t(ji,jj)) 
    251                   END IF  
    252                 END DO  
    253             ELSE  
    254                 emp (:,:) = emp (:,:) + zfwf(jc) / surf(jpncs+1) 
    255                 emps(:,:) = emps(:,:) + zfwf(jc) / surf(jpncs+1) 
    256             ENDIF 
    257          ELSEIF( ncstt(jc) == 2 ) THEN  
    258             ! Excess e-p+r (either sign) goes to open ocean, at outflow location 
    259             IF(      ji > 1 .AND. ji < jpi    & 
    260                .AND. jj > 1 .AND. jj < jpj ) THEN  
    261                 DO jn = 1, ncsnr(jc) 
    262                   ji = mi0(ncsir(jc,jn)) 
    263                   jj = mj0(ncsjr(jc,jn)) ! Location of outflow in open ocean 
    264                   emp (ji,jj) = emp (ji,jj) + zfwf(jc)   & 
    265                       / (FLOAT(ncsnr(jc)) *  e1t(ji,jj) * e2t(ji,jj) ) 
    266                   emps(ji,jj) = emps(ji,jj) + zfwf(jc)   & 
    267                       / (FLOAT(ncsnr(jc)) *  e1t(ji,jj) * e2t(ji,jj) ) 
    268                 END DO  
     310                  IF(      ji > 1 .AND. ji < jpi    & 
     311                     .AND. jj > 1 .AND. jj < jpj ) THEN  
     312                     emp (ji,jj) = emp (ji,jj) + zfwf(jc) / ( REAL(ncsnr(jc)) *  e1e2t(ji,jj) ) 
     313                     emps(ji,jj) = emps(ji,jj) + zfwf(jc) / ( REAL(ncsnr(jc)) *  e1e2t(ji,jj) ) 
     314                  ENDIF  
     315               END DO  
    269316            ENDIF  
    270          ENDIF  
    271          ! 
    272          DO jj = ncsj1(jc), ncsj2(jc) 
    273             DO ji = ncsi1(jc), ncsi2(jc) 
    274                emp (ji,jj) = emp (ji,jj) - zfwf(jc) / surf(jc) 
    275                emps(ji,jj) = emps(ji,jj) - zfwf(jc) / surf(jc) 
    276             END DO   
    277          END DO  
    278          ! 
     317            ! 
     318            DO jj = ncsj1(jc), ncsj2(jc) 
     319               DO ji = ncsi1(jc), ncsi2(jc) 
     320                  emp (ji,jj) = emp (ji,jj) - zfwf(jc) / surf(jc) 
     321                  emps(ji,jj) = emps(ji,jj) - zfwf(jc) / surf(jc) 
     322               END DO   
     323            END DO  
     324            ! 
     325         END IF 
    279326      END DO  
    280       ! 
    281       CALL lbc_lnk( emp , 'T', 1. ) 
    282       CALL lbc_lnk( emps, 'T', 1. ) 
     327 
     328      IF ( ABS(zcorr) > rsmall ) THEN      ! remove the global correction from the closed seas 
     329         DO jc = 1, jpncs                  ! only if it is large enough 
     330            DO jj = ncsj1(jc), ncsj2(jc) 
     331               DO ji = ncsi1(jc), ncsi2(jc) 
     332                  emp (ji,jj) = emp (ji,jj) - zcorr 
     333                  emps(ji,jj) = emps(ji,jj) - zcorr 
     334               END DO   
     335             END DO  
     336          END DO 
     337      ENDIF 
     338      ! 
     339      emp (:,:) = emp (:,:) * tmask(:,:,1) 
     340      emps(:,:) = emps(:,:) * tmask(:,:,1) 
     341      ! 
     342      CALL lbc_lnk( emp , 'T', 1._wp ) 
     343      CALL lbc_lnk( emps, 'T', 1._wp ) 
     344      ! 
     345      IF( nn_timing == 1 )  CALL timing_stop('sbc_clo') 
    283346      ! 
    284347   END SUBROUTINE sbc_clo 
    285     
    286     
     348 
     349 
    287350   SUBROUTINE clo_rnf( p_rnfmsk ) 
    288351      !!--------------------------------------------------------------------- 
     
    308371               ii = mi0( ncsir(jc,jn) ) 
    309372               ij = mj0( ncsjr(jc,jn) ) 
    310                p_rnfmsk(ii,ij) = MAX( p_rnfmsk(ii,ij), 1.0 ) 
     373               p_rnfmsk(ii,ij) = MAX( p_rnfmsk(ii,ij), 1.0_wp ) 
    311374            END DO  
    312375         ENDIF  
     
    336399         DO jj = ncsj1(jc), ncsj2(jc) 
    337400            DO ji = ncsi1(jc), ncsi2(jc) 
    338                p_upsmsk(ji,jj) = 0.5            ! mixed upstream/centered scheme over closed seas 
     401               p_upsmsk(ji,jj) = 0.5_wp         ! mixed upstream/centered scheme over closed seas 
    339402            END DO  
    340403         END DO  
     
    374437   !!====================================================================== 
    375438END MODULE closea 
     439 
  • branches/UKMO/CO5_package_branch/NEMOGCM/NEMO/OPA_SRC/DOM/daymod.F90

    r7363 r7367  
    116116 
    117117      ! number of seconds since the beginning of current year/month/week/day at the middle of the time-step 
    118       nsec_year  = nday_year * nsecd - ndt05   ! 1 time step before the middle of the first time step 
    119       nsec_month = nday      * nsecd - ndt05   ! because day will be called at the beginning of step 
    120       nsec_week  = idweek    * nsecd - ndt05 
    121       nsec_day   =             nsecd - ndt05 
     118      nsec_year  = nday_year * nsecd - ndt   ! 1 time step before the middle of the first time step 
     119      nsec_month = nday      * nsecd - ndt   ! because day will be called at the beginning of step 
     120      nsec_week  = idweek    * nsecd - ndt 
     121      nsec_day   =             nsecd - ndt 
    122122 
    123123      ! control print 
     
    219219      IF( ABS(adatrj  - REAL(NINT(adatrj ),wp)) < zprec )   adatrj  = REAL(NINT(adatrj ),wp)   ! avoid truncation error 
    220220       
    221       IF( nsec_day > nsecd ) THEN                       ! New day 
     221      IF( nsec_day >= nsecd ) THEN                       ! New day 
    222222         ! 
    223223         nday      = nday + 1 
  • branches/UKMO/CO5_package_branch/NEMOGCM/NEMO/OPA_SRC/DOM/dom_oce.F90

    r7363 r7367  
    5252   REAL(wp), PUBLIC ::   rdtmax          !: maximum time step on tracers 
    5353   REAL(wp), PUBLIC ::   rdth            !: depth variation of tracer step 
    54    INTEGER , PUBLIC ::   nclosea         !: =0 suppress closed sea/lake from the ORCA domain or not (=1) 
    5554 
    5655   !                                                  !!! associated variables 
     
    125124   LOGICAL, PUBLIC ::   ln_zps     =  .FALSE.   !: z-coordinate - partial step 
    126125   LOGICAL, PUBLIC ::   ln_sco     =  .FALSE.   !: s-coordinate or hybrid z-s coordinate 
     126   LOGICAL, PUBLIC ::   ln_read_zenv     =  .FALSE.   !: Whether to read zenv or calculate it 
    127127 
    128128   !! All coordinates 
     
    173173   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   scosrf, scobot   !: ocean surface and bottom topographies  
    174174   !                                        !  (if deviating from coordinate surfaces in HYBRID) 
     175   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   rx1              !: maximum grid stiffness ratio 
    175176   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hifv  , hiff     !: interface depth between stretching at  V--F 
    176177   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hift  , hifu     !: and quasi-uniform spacing              T--U  points (m) 
     178   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   zenv    !:       Envelope Batymetry, calcualted or read in  
    177179 
    178180   !!---------------------------------------------------------------------- 
     
    295297         &      scosrf(jpi,jpj) , scobot(jpi,jpj) ,     & 
    296298         &      hifv  (jpi,jpj) , hiff  (jpi,jpj) ,     & 
    297          &      hift  (jpi,jpj) , hifu  (jpi,jpj) , STAT=ierr(8) ) 
     299         &      hift  (jpi,jpj) , hifu  (jpi,jpj) , rx1(jpi,jpj), STAT=ierr(8) ) 
    298300 
    299301      ALLOCATE( mbathy(jpi,jpj) , bathy(jpi,jpj) ,                     & 
    300302         &     tmask_i(jpi,jpj) , bmask(jpi,jpj) ,                     & 
    301          &     mbkt   (jpi,jpj) , mbku (jpi,jpj) , mbkv(jpi,jpj) , STAT=ierr(9) ) 
     303         &     mbkt   (jpi,jpj) , mbku (jpi,jpj) , mbkv(jpi,jpj) , zenv(jpi,jpj), STAT=ierr(9) ) 
    302304 
    303305      ALLOCATE( tmask(jpi,jpj,jpk) , umask(jpi,jpj,jpk),     &  
  • branches/UKMO/CO5_package_branch/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90

    r7363 r7367  
    3636   USE dyncor_c1d      ! Coriolis term (c1d case)         (cor_c1d routine) 
    3737   USE timing          ! Timing 
     38   USE lbclnk 
    3839 
    3940   IMPLICIT NONE 
     
    8485                             CALL dom_zgr      ! Vertical mesh and bathymetry 
    8586                             CALL dom_msk      ! Masks 
     87      IF( ln_sco )           CALL dom_stiff    ! Maximum stiffness ratio/hydrostatic consistency       
    8688      IF( lk_vvl         )   CALL dom_vvl      ! Vertical variable mesh 
    8789      ! 
     
    123125      !!---------------------------------------------------------------------- 
    124126      USE ioipsl 
     127      NAMELIST/namrun/ ln_NOOS   
    125128      NAMELIST/namrun/ nn_no   , cn_exp    , cn_ocerst_in, cn_ocerst_out, ln_rstart , nn_rstctl,   & 
     129         &             nn_stocklist,                                                               & 
    126130         &             nn_it000, nn_itend  , nn_date0    , nn_leapy     , nn_istate , nn_stock ,   & 
    127          &             nn_write, ln_dimgnnn, ln_mskland  , ln_clobber   , nn_chunksz 
     131         &             nn_write, ln_dimgnnn, ln_mskland  , ln_clobber   , nn_chunksz,              & 
     132         &             ln_diafoam, nn_diafoam, ln_depwri  
    128133      NAMELIST/namdom/ nn_bathy , rn_e3zps_min, rn_e3zps_rat, nn_msh    , rn_hmin,   & 
    129134         &             nn_acc   , rn_atfp     , rn_rdt      , rn_rdtmin ,            & 
    130135         &             rn_rdtmax, rn_rdth     , nn_baro     , nn_closea 
    131136      NAMELIST/namcla/ nn_cla 
     137      NAMELIST/namrun/ ln_rstdate 
    132138#if defined key_netcdf4 
    133139      NAMELIST/namnc4/ nn_nchunks_i, nn_nchunks_j, nn_nchunks_k, ln_nc4zip 
    134140#endif 
    135       !!---------------------------------------------------------------------- 
    136  
     141      NAMELIST/namrun/ ln_diatide 
     142      !!---------------------------------------------------------------------- 
     143      NAMELIST/namrun/ ln_diatmb 
     144 
     145 
     146      NAMELIST/namrun/ cn_rst_dir ! moved here to allow merge with CO5 branches (ln_NOOS)                                                                  
    137147      REWIND( numnam )              ! Namelist namrun : parameters of the run 
    138148      READ  ( numnam, namrun ) 
     
    152162         WRITE(numout,*) '      leap year calendar (0/1)        nn_leapy   = ', nn_leapy 
    153163         WRITE(numout,*) '      initial state output            nn_istate  = ', nn_istate 
    154          WRITE(numout,*) '      frequency of restart file       nn_stock   = ', nn_stock 
     164         IF ( ALL( nn_stocklist == 0 ) ) THEN   
     165            WRITE(numout,*) '   frequency of restart file       nn_stock   = ', nn_stock 
     166         ELSE   
     167            WRITE(numout,*) '   list of restart times           nn_stocklist = ', nn_stocklist(1:10)   
     168         ENDIF   
    155169         WRITE(numout,*) '      frequency of output file        nn_write   = ', nn_write 
    156170         WRITE(numout,*) '      multi file dimgout              ln_dimgnnn = ', ln_dimgnnn 
     171         WRITE(numout,*) '      use date in restart name        ln_rstdate = ', ln_rstdate  
    157172         WRITE(numout,*) '      mask land points                ln_mskland = ', ln_mskland 
     173         WRITE(numout,*) '      NOOS transect diagnostics       ln_NOOS    = ', ln_NOOS 
    158174         WRITE(numout,*) '      overwrite an existing file      ln_clobber = ', ln_clobber 
     175         WRITE(numout,*) '      restart directory               cn_rst_dir = ', cn_rst_dir  
    159176         WRITE(numout,*) '      NetCDF chunksize (bytes)        nn_chunksz = ', nn_chunksz 
     177         WRITE(numout,*) '      Met Office FOAM diagnostics     ln_diafoam = ', ln_diafoam   
     178         WRITE(numout,*) '      FOAM diagnostic choices         nn_diafoam = ', nn_diafoam   
     179         WRITE(numout,*) '      depths file output logical      ln_depwri  = ', ln_depwri  
    160180      ENDIF 
    161181 
     
    169189      ninist = nn_istate 
    170190      nstock = nn_stock 
     191      nstock_list = nn_stocklist 
    171192      nwrite = nn_write 
    172193 
     
    238259      rdtmax    = rn_rdtmin 
    239260      rdth      = rn_rdth 
    240       nclosea   = nn_closea 
    241261 
    242262      REWIND( numnam )              ! Namelist cross land advection 
     
    274294 
    275295 
     296   !!====================================================================== 
    276297   SUBROUTINE dom_ctl 
    277298      !!---------------------------------------------------------------------- 
     
    323344   END SUBROUTINE dom_ctl 
    324345 
     346   SUBROUTINE dom_stiff 
     347      !!---------------------------------------------------------------------- 
     348      !!                  ***  ROUTINE dom_stiff  *** 
     349      !!                      
     350      !! ** Purpose :   Diagnose maximum grid stiffness/hydrostatic consistency 
     351      !! 
     352      !! ** Method  :   Compute Haney (1991) hydrostatic condition ratio 
     353      !!                Save the maximum in the vertical direction 
     354      !!                (this number is only relevant in s-coordinates) 
     355      !! 
     356      !!                Haney, R. L., 1991: On the pressure gradient force 
     357      !!                over steep topography in sigma coordinate ocean models.  
     358      !!                J. Phys. Oceanogr., 21, 610???619. 
     359      !!---------------------------------------------------------------------- 
     360      INTEGER  ::   ji, jj, jk  
     361      REAL(wp) ::   zrxmax 
     362      REAL(wp), DIMENSION(4) :: zr1 
     363      !!---------------------------------------------------------------------- 
     364      rx1(:,:) = 0.e0 
     365      zrxmax   = 0.e0 
     366      zr1(:)   = 0.e0 
     367       
     368      DO ji = 2, jpim1 
     369         DO jj = 2, jpjm1 
     370            DO jk = 1, jpkm1 
     371               zr1(1) = umask(ji-1,jj  ,jk) *abs( (gdepw(ji  ,jj  ,jk  )-gdepw(ji-1,jj  ,jk  )  &  
     372                    &                         +gdepw(ji  ,jj  ,jk+1)-gdepw(ji-1,jj  ,jk+1)) & 
     373                    &                        /(gdepw(ji  ,jj  ,jk  )+gdepw(ji-1,jj  ,jk  )  & 
     374                    &                         -gdepw(ji  ,jj  ,jk+1)-gdepw(ji-1,jj  ,jk+1) + rsmall) ) 
     375               zr1(2) = umask(ji  ,jj  ,jk) *abs( (gdepw(ji+1,jj  ,jk  )-gdepw(ji  ,jj  ,jk  )  & 
     376                    &                         +gdepw(ji+1,jj  ,jk+1)-gdepw(ji  ,jj  ,jk+1)) & 
     377                    &                        /(gdepw(ji+1,jj  ,jk  )+gdepw(ji  ,jj  ,jk  )  & 
     378                    &                         -gdepw(ji+1,jj  ,jk+1)-gdepw(ji  ,jj  ,jk+1) + rsmall) ) 
     379               zr1(3) = vmask(ji  ,jj  ,jk) *abs( (gdepw(ji  ,jj+1,jk  )-gdepw(ji  ,jj  ,jk  )  & 
     380                    &                         +gdepw(ji  ,jj+1,jk+1)-gdepw(ji  ,jj  ,jk+1)) & 
     381                    &                        /(gdepw(ji  ,jj+1,jk  )+gdepw(ji  ,jj  ,jk  )  & 
     382                    &                         -gdepw(ji  ,jj+1,jk+1)-gdepw(ji  ,jj  ,jk+1) + rsmall) ) 
     383               zr1(4) = vmask(ji  ,jj-1,jk) *abs( (gdepw(ji  ,jj  ,jk  )-gdepw(ji  ,jj-1,jk  )  & 
     384                    &                         +gdepw(ji  ,jj  ,jk+1)-gdepw(ji  ,jj-1,jk+1)) & 
     385                    &                        /(gdepw(ji  ,jj  ,jk  )+gdepw(ji  ,jj-1,jk  )  & 
     386                    &                         -gdepw(ji,  jj  ,jk+1)-gdepw(ji  ,jj-1,jk+1) + rsmall) ) 
     387               zrxmax = MAXVAL(zr1(1:4)) 
     388               rx1(ji,jj) = MAX(rx1(ji,jj), zrxmax) 
     389            END DO 
     390         END DO 
     391      END DO 
     392 
     393      CALL lbc_lnk( rx1, 'T', 1. ) 
     394 
     395      zrxmax = MAXVAL(rx1) 
     396 
     397      IF( lk_mpp )   CALL mpp_max( zrxmax ) ! max over the global domain 
     398 
     399      IF(lwp) THEN 
     400         WRITE(numout,*) 
     401         WRITE(numout,*) 'dom_stiff : maximum grid stiffness ratio: ', zrxmax 
     402         WRITE(numout,*) '~~~~~~~~~' 
     403      ENDIF 
     404 
     405   END SUBROUTINE dom_stiff 
     406 
    325407   !!====================================================================== 
    326408END MODULE domain 
  • branches/UKMO/CO5_package_branch/NEMOGCM/NEMO/OPA_SRC/DOM/domwri.F90

    r7363 r7367  
    172172             
    173173      IF( ln_sco ) THEN                                         ! s-coordinate 
    174          CALL iom_rstput( 0, 0, inum4, 'hbatt', hbatt )         !    ! depth 
    175          CALL iom_rstput( 0, 0, inum4, 'hbatu', hbatu )  
     174         CALL iom_rstput( 0, 0, inum4, 'hbatt', hbatt ) 
     175         CALL iom_rstput( 0, 0, inum4, 'hbatu', hbatu ) 
    176176         CALL iom_rstput( 0, 0, inum4, 'hbatv', hbatv ) 
    177177         CALL iom_rstput( 0, 0, inum4, 'hbatf', hbatf ) 
     
    187187         CALL iom_rstput( 0, 0, inum4, 'e3v', e3v ) 
    188188         CALL iom_rstput( 0, 0, inum4, 'e3w', e3w ) 
    189          ! 
    190          CALL iom_rstput( 0, 0, inum4, 'gdept_0' , gdept_0 )    !    ! stretched system 
    191          CALL iom_rstput( 0, 0, inum4, 'gdepw_0' , gdepw_0 ) 
     189         CALL iom_rstput( 0, 0, inum4, 'rx1', rx1 )             !    ! Max. grid stiffness ratio 
     190         ! 
     191         CALL iom_rstput( 0, 0, inum4, 'gdept' , gdept )    !    ! stretched system 
     192         CALL iom_rstput( 0, 0, inum4, 'gdepw' , gdepw ) 
    192193      ENDIF 
    193194       
  • branches/UKMO/CO5_package_branch/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90

    r7363 r7367  
    1515   !!            3.2  ! 2009-07  (R. Benshila) Suppression of rigid-lid option 
    1616   !!            3.3  ! 2010-11  (G. Madec) add mbk. arrays associated to the deepest ocean level 
     17   !!            3.4  ! 2012-08  (J. Siddorn) added Siddorn adn Furner stretching functio 
    1718   !!---------------------------------------------------------------------- 
    1819 
     
    2829   !!       zgr_sco      : s-coordinate 
    2930   !!       fssig        : sigma coordinate non-dimensional function 
     31   !!       fgamma       : Siddorn and Furner stretching function 
    3032   !!       dfssig       : derivative of the sigma coordinate function    !!gm  (currently missing!) 
    3133   !!--------------------------------------------------------------------- 
     
    4749 
    4850   !                                       !!* Namelist namzgr_sco * 
     51   LOGICAL  ::   ln_s_sh94   = .false.      ! use hybrid s-sig Song and Haidvogel 1994 stretching function fssig1 (ln_sco=T) 
     52   LOGICAL  ::   ln_s_sf12   = .true.       ! use hybrid s-z-sig Siddorn and Furner 2012 stretching function fgamma (ln_sco=T) 
     53   LOGICAL  ::   ln_sigcrit  = .false.      ! use sigma coordinates below critical depth (T) or Z coordinates (F) for Siddorn & Furner stretch  
     54   ! 
    4955   REAL(wp) ::   rn_sbot_min =  300._wp     ! minimum depth of s-bottom surface (>0) (m) 
    5056   REAL(wp) ::   rn_sbot_max = 5250._wp     ! maximum depth of s-bottom surface (= ocean depth) (>0) (m) 
     57   REAL(wp) ::   rn_rmax     =    0.15_wp   ! maximum cut-off r-value allowed (0<rn_rmax<1) 
     58   REAL(wp) ::   rn_hc       =  150._wp     ! Critical depth for transition from sigma to stretched coordinates 
     59   ! Song and Haidvogel 1994 stretching parameters 
    5160   REAL(wp) ::   rn_theta    =    6.00_wp   ! surface control parameter (0<=rn_theta<=20) 
    5261   REAL(wp) ::   rn_thetb    =    0.75_wp   ! bottom control parameter  (0<=rn_thetb<= 1) 
    53    REAL(wp) ::   rn_rmax     =    0.15_wp   ! maximum cut-off r-value allowed (0<rn_rmax<1) 
    54    LOGICAL  ::   ln_s_sigma  = .false.      ! use hybrid s-sigma -coordinate & stretching function fssig1 (ln_sco=T) 
    55    REAL(wp) ::   rn_bb       =    0.80_wp   ! stretching parameter for song and haidvogel stretching 
     62   REAL(wp) ::   rn_bb       =    0.80_wp   ! stretching parameter  
    5663   !                                        ! ( rn_bb=0; top only, rn_bb =1; top and bottom) 
    57    REAL(wp) ::   rn_hc       =  150._wp     ! Critical depth for s-sigma coordinates 
     64   ! Siddorn and Furner stretching parameters 
     65   REAL(wp) ::   rn_alpha    =    4.4_wp    ! control parameter ( > 1 stretch towards surface, < 1 towards seabed) 
     66   REAL(wp) ::   rn_efold    =    0.0_wp    !  efold length scale for transition to stretched coord 
     67   REAL(wp) ::   rn_zs       =    1.0_wp    !  depth of surface grid box 
     68                           !  bottom cell depth (Zb) is a linear function of water depth Zb = H*a + b 
     69   REAL(wp) ::   rn_zb_a     =    0.024_wp  !  bathymetry scaling factor for calculating Zb 
     70   REAL(wp) ::   rn_zb_b     =   -0.2_wp    !  offset for calculating Zb 
    5871 
    5972  !! * Substitutions 
     
    8699      INTEGER ::   ioptio = 0   ! temporary integer 
    87100      ! 
    88       NAMELIST/namzgr/ ln_zco, ln_zps, ln_sco 
     101      NAMELIST/namzgr/ ln_zco, ln_zps, ln_sco, ln_read_zenv 
    89102      !!---------------------------------------------------------------------- 
    90103      ! 
     
    102115         WRITE(numout,*) '             z-coordinate - partial steps   ln_zps = ', ln_zps 
    103116         WRITE(numout,*) '             s- or hybrid z-s-coordinate    ln_sco = ', ln_sco 
     117         WRITE(numout,*) '             Read Zenv from Bathy T/F ln_read_zenv = ', ln_read_zenv 
    104118      ENDIF 
    105119 
     
    243257         END DO 
    244258      ELSE                                ! Madec & Imbard 1996 function 
     259# if key_levels==1  
     260         !Hard wire a deep and shallow level  
     261         !NOTE this configuration is for use with NEMOVAR,   
     262         !it is not set-up for NEMO  
     263         CALL ctl_warn("Single level model, depth of first layer set to 1cm."//&  
     264             &  "\nThis configuration is designed to be used with NEMOVAR only")  
     265         gdepw_0(1)=0  
     266         gdept_0(1)=0.01  
     267         gdepw_0(2)=7000  
     268         gdept_0(2)=14000  
     269         e3w_0(:)=7000  
     270         e3t_0(1)=6999.99  
     271         e3t_0(2)=7000  
     272# else  
    245273         IF( .NOT. ldbletanh ) THEN 
    246274            DO jk = 1, jpk 
     
    267295            END DO 
    268296         ENDIF 
     297# endif 
    269298         gdepw_0(1) = 0._wp                    ! force first w-level to be exactly at zero 
    270299      ENDIF 
     
    422451            CALL iom_close( inum ) 
    423452            mbathy(:,:) = INT( bathy(:,:) ) 
    424             !                                                ! ===================== 
     453            ! 
    425454            IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN    ! ORCA R2 configuration 
    426                !                                             ! ===================== 
     455               ! 
    427456               IF( nn_cla == 0 ) THEN 
    428457                  ii0 = 140   ;   ii1 = 140                  ! Gibraltar Strait open  
     
    453482            CALL iom_open ( 'bathy_meter.nc', inum )  
    454483            CALL iom_get  ( inum, jpdom_data, 'Bathymetry', bathy ) 
     484            IF ( ln_read_zenv ) THEN                  ! Whether we should read zenv or not 
     485                CALL iom_get  ( inum, jpdom_data, 'zenv', zenv ) 
     486            ENDIF 
    455487            CALL iom_close( inum ) 
    456             !                                                ! ===================== 
     488            !                                                 
    457489            IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN    ! ORCA R2 configuration 
    458                !                                             ! ===================== 
     490               ! 
    459491              IF( nn_cla == 0 ) THEN 
    460492                 ii0 = 140   ;   ii1 = 140                   ! Gibraltar Strait open  
     
    489521      ENDIF 
    490522      ! 
    491       !                                               ! =========================== ! 
    492       IF( nclosea == 0 ) THEN                         !   NO closed seas or lakes   ! 
    493          DO jl = 1, jpncs                             ! =========================== ! 
    494             DO jj = ncsj1(jl), ncsj2(jl) 
    495                DO ji = ncsi1(jl), ncsi2(jl) 
    496                   mbathy(ji,jj) = 0                   ! suppress closed seas and lakes from bathymetry 
    497                   bathy (ji,jj) = 0._wp                
    498                END DO 
    499             END DO 
    500          END DO 
    501       ENDIF 
    502       ! 
    503       !                                               ! =========================== ! 
    504       !                                               !     set a minimum depth     ! 
    505       !                                               ! =========================== ! 
    506       IF ( .not. ln_sco ) THEN 
     523      IF( nn_closea == 0 )   CALL clo_bat( bathy, mbathy )    !==  NO closed seas or lakes  ==! 
     524      !                        
     525      IF ( .not. ln_sco ) THEN                                !==  set a minimum depth  ==! 
    507526         IF( rn_hmin < 0._wp ) THEN    ;   ik = - INT( rn_hmin )                                      ! from a nb of level 
    508527         ELSE                          ;   ik = MINLOC( gdepw_0, mask = gdepw_0 > rn_hmin, dim = 1 )  ! from a depth 
     
    10471066   END SUBROUTINE zgr_zps 
    10481067 
    1049  
    1050    FUNCTION fssig( pk ) RESULT( pf ) 
    1051       !!---------------------------------------------------------------------- 
    1052       !!                 ***  ROUTINE eos_init  *** 
    1053       !!        
    1054       !! ** Purpose :   provide the analytical function in s-coordinate 
    1055       !!           
    1056       !! ** Method  :   the function provide the non-dimensional position of 
    1057       !!                T and W (i.e. between 0 and 1) 
    1058       !!                T-points at integer values (between 1 and jpk) 
    1059       !!                W-points at integer values - 1/2 (between 0.5 and jpk-0.5) 
    1060       !!---------------------------------------------------------------------- 
    1061       REAL(wp), INTENT(in) ::   pk   ! continuous "k" coordinate 
    1062       REAL(wp)             ::   pf   ! sigma value 
    1063       !!---------------------------------------------------------------------- 
    1064       ! 
    1065       pf =   (   TANH( rn_theta * ( -(pk-0.5_wp) / REAL(jpkm1) + rn_thetb )  )   & 
    1066          &     - TANH( rn_thetb * rn_theta                                )  )   & 
    1067          & * (   COSH( rn_theta                           )                      & 
    1068          &     + COSH( rn_theta * ( 2._wp * rn_thetb - 1._wp ) )  )              & 
    1069          & / ( 2._wp * SINH( rn_theta ) ) 
    1070       ! 
    1071    END FUNCTION fssig 
    1072  
    1073  
    1074    FUNCTION fssig1( pk1, pbb ) RESULT( pf1 ) 
    1075       !!---------------------------------------------------------------------- 
    1076       !!                 ***  ROUTINE eos_init  *** 
    1077       !! 
    1078       !! ** Purpose :   provide the Song and Haidvogel version of the analytical function in s-coordinate 
    1079       !! 
    1080       !! ** Method  :   the function provides the non-dimensional position of 
    1081       !!                T and W (i.e. between 0 and 1) 
    1082       !!                T-points at integer values (between 1 and jpk) 
    1083       !!                W-points at integer values - 1/2 (between 0.5 and jpk-0.5) 
    1084       !!---------------------------------------------------------------------- 
    1085       REAL(wp), INTENT(in) ::   pk1   ! continuous "k" coordinate 
    1086       REAL(wp), INTENT(in) ::   pbb   ! Stretching coefficient 
    1087       REAL(wp)             ::   pf1   ! sigma value 
    1088       !!---------------------------------------------------------------------- 
    1089       ! 
    1090       IF ( rn_theta == 0 ) then      ! uniform sigma 
    1091          pf1 = - ( pk1 - 0.5_wp ) / REAL( jpkm1 ) 
    1092       ELSE                        ! stretched sigma 
    1093          pf1 =   ( 1._wp - pbb ) * ( SINH( rn_theta*(-(pk1-0.5_wp)/REAL(jpkm1)) ) ) / SINH( rn_theta )              & 
    1094             &  + pbb * (  (TANH( rn_theta*( (-(pk1-0.5_wp)/REAL(jpkm1)) + 0.5_wp) ) - TANH( 0.5_wp * rn_theta )  )  & 
    1095             &        / ( 2._wp * TANH( 0.5_wp * rn_theta ) )  ) 
    1096       ENDIF 
    1097       ! 
    1098    END FUNCTION fssig1 
    1099  
    1100  
    11011068   SUBROUTINE zgr_sco 
    11021069      !!---------------------------------------------------------------------- 
     
    11231090      !!            esigt(k) = fsdsig(k    ) 
    11241091      !!            esigw(k) = fsdsig(k-0.5) 
    1125       !!      This routine is given as an example, it must be modified 
    1126       !!      following the user s desiderata. nevertheless, the output as 
     1092      !!      Three options for stretching are give, and they can be modified 
     1093      !!      following the users requirements. Nevertheless, the output as 
    11271094      !!      well as the way to compute the model levels and scale factors 
    1128       !!      must be respected in order to insure second order a!!uracy 
     1095      !!      must be respected in order to insure second order accuracy 
    11291096      !!      schemes. 
    11301097      !! 
    1131       !! Reference : Madec, Lott, Delecluse and Crepon, 1996. JPO, 26, 1393-1408. 
     1098      !!      The three methods for stretching available are: 
     1099      !!  
     1100      !!           s_sh94 (Song and Haidvogel 1994) 
     1101      !!                a sinh/tanh function that allows sigma and stretched sigma 
     1102      !! 
     1103      !!           s_sf12 (Siddorn and Furner 2012?) 
     1104      !!                allows the maintenance of fixed surface and or 
     1105      !!                bottom cell resolutions (cf. geopotential coordinates)  
     1106      !!                within an analytically derived stretched S-coordinate framework. 
     1107      !!  
     1108      !!          s_tanh  (Madec et al 1996) 
     1109      !!                a cosh/tanh function that gives stretched coordinates         
     1110      !! 
    11321111      !!---------------------------------------------------------------------- 
    11331112      ! 
    11341113      INTEGER  ::   ji, jj, jk, jl           ! dummy loop argument 
    11351114      INTEGER  ::   iip1, ijp1, iim1, ijm1   ! temporary integers 
    1136       REAL(wp) ::   zcoeft, zcoefw, zrmax, ztaper   ! temporary scalars 
    1137       ! 
    1138       REAL(wp), POINTER, DIMENSION(:,:  ) :: zenv, ztmp, zmsk, zri, zrj, zhbat 
     1115      REAL(wp) ::   zrmax, ztaper   ! temporary scalars 
     1116      ! 
     1117      REAL(wp), POINTER, DIMENSION(:,:  ) :: ztmp, zmsk, zri, zrj, zhbat 
    11391118      REAL(wp), POINTER, DIMENSION(:,:,:) :: gsigw3, gsigt3, gsi3w3 
    11401119      REAL(wp), POINTER, DIMENSION(:,:,:) :: esigt3, esigw3, esigtu3, esigtv3, esigtf3, esigwu3, esigwv3            
    11411120 
    1142       NAMELIST/namzgr_sco/ rn_sbot_max, rn_sbot_min, rn_theta, rn_thetb, rn_rmax, ln_s_sigma, rn_bb, rn_hc 
    1143       !!---------------------------------------------------------------------- 
     1121      NAMELIST/namzgr_sco/ln_s_sh94, ln_s_sf12, ln_sigcrit, rn_sbot_min, rn_sbot_max, rn_hc, rn_rmax,rn_theta, & 
     1122                           rn_thetb, rn_bb, rn_alpha, rn_efold, rn_zs, rn_zb_a, rn_zb_b 
     1123     !!---------------------------------------------------------------------- 
    11441124      ! 
    11451125      IF( nn_timing == 1 )  CALL timing_start('zgr_sco') 
    11461126      ! 
    1147       CALL wrk_alloc( jpi, jpj,      zenv, ztmp, zmsk, zri, zrj, zhbat                           ) 
     1127      CALL wrk_alloc( jpi, jpj,      ztmp, zmsk, zri, zrj, zhbat                           ) 
    11481128      CALL wrk_alloc( jpi, jpj, jpk, gsigw3, gsigt3, gsi3w3                                      ) 
    11491129      CALL wrk_alloc( jpi, jpj, jpk, esigt3, esigw3, esigtu3, esigtv3, esigtf3, esigwu3, esigwv3 ) 
     
    11571137         WRITE(numout,*) '~~~~~~~~~~~' 
    11581138         WRITE(numout,*) '   Namelist namzgr_sco' 
    1159          WRITE(numout,*) '      sigma-stretching coeffs ' 
    1160          WRITE(numout,*) '      maximum depth of s-bottom surface (>0)       rn_sbot_max   = ' ,rn_sbot_max 
    1161          WRITE(numout,*) '      minimum depth of s-bottom surface (>0)       rn_sbot_min   = ' ,rn_sbot_min 
    1162          WRITE(numout,*) '      surface control parameter (0<=rn_theta<=20)  rn_theta      = ', rn_theta 
    1163          WRITE(numout,*) '      bottom  control parameter (0<=rn_thetb<= 1)  rn_thetb      = ', rn_thetb 
    1164          WRITE(numout,*) '      maximum cut-off r-value allowed              rn_rmax       = ', rn_rmax 
    1165          WRITE(numout,*) '      Hybrid s-sigma-coordinate                    ln_s_sigma    = ', ln_s_sigma 
    1166          WRITE(numout,*) '      stretching parameter (song and haidvogel)    rn_bb         = ', rn_bb 
    1167          WRITE(numout,*) '      Critical depth                               rn_hc         = ', rn_hc 
    1168       ENDIF 
    1169  
    1170       gsigw3  = 0._wp   ;   gsigt3  = 0._wp   ;   gsi3w3  = 0._wp 
    1171       esigt3  = 0._wp   ;   esigw3  = 0._wp  
    1172       esigtu3 = 0._wp   ;   esigtv3 = 0._wp   ;   esigtf3 = 0._wp 
    1173       esigwu3 = 0._wp   ;   esigwv3 = 0._wp 
     1139         WRITE(numout,*) '     stretching coeffs ' 
     1140         WRITE(numout,*) '        maximum depth of s-bottom surface (>0)       rn_sbot_max   = ',rn_sbot_max 
     1141         WRITE(numout,*) '        minimum depth of s-bottom surface (>0)       rn_sbot_min   = ',rn_sbot_min 
     1142         WRITE(numout,*) '        Critical depth                               rn_hc         = ',rn_hc 
     1143         WRITE(numout,*) '        maximum cut-off r-value allowed              rn_rmax       = ',rn_rmax 
     1144         WRITE(numout,*) '     Song and Haidvogel 1994 stretching              ln_s_sh94     = ',ln_s_sh94 
     1145         WRITE(numout,*) '        Song and Haidvogel 1994 stretching coefficients' 
     1146         WRITE(numout,*) '        surface control parameter (0<=rn_theta<=20)  rn_theta      = ',rn_theta 
     1147         WRITE(numout,*) '        bottom  control parameter (0<=rn_thetb<= 1)  rn_thetb      = ',rn_thetb 
     1148         WRITE(numout,*) '        stretching parameter (song and haidvogel)    rn_bb         = ',rn_bb 
     1149         WRITE(numout,*) '     Siddorn and Furner 2012 stretching              ln_s_sf12     = ',ln_s_sf12 
     1150         WRITE(numout,*) '        Siddorn and Furner 2012 stretching coefficients' 
     1151         WRITE(numout,*) '        stretchin parameter ( >1 surface; <1 bottom) rn_alpha      = ',rn_alpha 
     1152         WRITE(numout,*) '        e-fold length scale for transition region    rn_efold      = ',rn_efold 
     1153         WRITE(numout,*) '        Surface cell depth (Zs) (m)                  rn_zs         = ',rn_zs 
     1154         WRITE(numout,*) '        Bathymetry multiplier for Zb                 rn_zb_a       = ',rn_zb_a 
     1155         WRITE(numout,*) '        Offset for Zb                                rn_zb_b       = ',rn_zb_b 
     1156         WRITE(numout,*) '        Bottom cell (Zb) (m) = H*rn_zb_a + rn_zb_b' 
     1157      ENDIF 
    11741158 
    11751159      hift(:,:) = rn_sbot_min                     ! set the minimum depth for the s-coordinate 
     
    11901174      !                                        ! ============================= 
    11911175      ! use r-value to create hybrid coordinates 
     1176 
     1177      ! Smooth the bathymetry (if required) 
     1178      scosrf(:,:) = 0._wp             ! ocean surface depth (here zero: no under ice-shelf sea) 
     1179      scobot(:,:) = bathy(:,:)        ! ocean bottom  depth 
     1180      IF( ln_read_zenv) THEN 
     1181         WRITE(numout,*) '      Zenv is not calculated but read from Bathy File  ln_read_zenv        = ', ln_read_zenv 
     1182      ELSE 
     1183      IF ( jpnij .gt.1)  CALL ctl_stop( ' zgr_zps : ln_read_zenv=false and jpnij > 1,  Calculating zenv on more than one Proc is not safe, calculate on one proc only ' )    
     1184 
    11921185      DO jj = 1, jpj 
    11931186         DO ji = 1, jpi 
     
    11961189      END DO 
    11971190      !  
    1198       ! Smooth the bathymetry (if required) 
    1199       scosrf(:,:) = 0._wp             ! ocean surface depth (here zero: no under ice-shelf sea) 
    1200       scobot(:,:) = bathy(:,:)        ! ocean bottom  depth 
    12011191      ! 
    12021192      jl = 0 
     
    12701260      ! 
    12711261      !                                        ! envelop bathymetry saved in hbatt 
     1262      ENDIF       ! End of IF block for reading from a file or calculating zenv 
    12721263      hbatt(:,:) = zenv(:,:)  
    12731264      IF( MINVAL( gphit(:,:) ) * MAXVAL( gphit(:,:) ) <= 0._wp ) THEN 
     
    13651356      ! non-dimensional "sigma" for model level depth at w- and t-levels 
    13661357 
    1367       IF( ln_s_sigma ) THEN        ! Song and Haidvogel style stretched sigma for depths 
    1368          !                         ! below rn_hc, with uniform sigma in shallower waters 
    1369          DO ji = 1, jpi 
    1370             DO jj = 1, jpj 
    1371  
    1372                IF( hbatt(ji,jj) > rn_hc ) THEN    !deep water, stretched sigma 
    1373                   DO jk = 1, jpk 
    1374                      gsigw3(ji,jj,jk) = -fssig1( REAL(jk,wp)-0.5_wp, rn_bb ) 
    1375                      gsigt3(ji,jj,jk) = -fssig1( REAL(jk,wp)       , rn_bb ) 
    1376                   END DO 
    1377                ELSE ! shallow water, uniform sigma 
    1378                   DO jk = 1, jpk 
    1379                      gsigw3(ji,jj,jk) =   REAL(jk-1,wp)            / REAL(jpk-1,wp) 
    1380                      gsigt3(ji,jj,jk) = ( REAL(jk-1,wp) + 0.5_wp ) / REAL(jpk-1,wp) 
    1381                   END DO 
    1382                ENDIF 
    1383                IF( nprint == 1 .AND. lwp )   WRITE(numout,*) 'gsigw3 1 jpk    ', gsigw3(ji,jj,1), gsigw3(ji,jj,jpk) 
    1384                ! 
    1385                DO jk = 1, jpkm1 
    1386                   esigt3(ji,jj,jk  ) = gsigw3(ji,jj,jk+1) - gsigw3(ji,jj,jk) 
    1387                   esigw3(ji,jj,jk+1) = gsigt3(ji,jj,jk+1) - gsigt3(ji,jj,jk) 
    1388                END DO 
    1389                esigw3(ji,jj,1  ) = 2._wp * ( gsigt3(ji,jj,1  ) - gsigw3(ji,jj,1  ) ) 
    1390                esigt3(ji,jj,jpk) = 2._wp * ( gsigt3(ji,jj,jpk) - gsigw3(ji,jj,jpk) ) 
    1391                ! 
    1392                ! Coefficients for vertical depth as the sum of e3w scale factors 
    1393                gsi3w3(ji,jj,1) = 0.5_wp * esigw3(ji,jj,1) 
    1394                DO jk = 2, jpk 
    1395                   gsi3w3(ji,jj,jk) = gsi3w3(ji,jj,jk-1) + esigw3(ji,jj,jk) 
    1396                END DO 
    1397                ! 
    1398                DO jk = 1, jpk 
    1399                   zcoeft = ( REAL(jk,wp) - 0.5_wp ) / REAL(jpkm1,wp) 
    1400                   zcoefw = ( REAL(jk,wp) - 1.0_wp ) / REAL(jpkm1,wp) 
    1401                   gdept (ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*gsigt3(ji,jj,jk)+rn_hc*zcoeft ) 
    1402                   gdepw (ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*gsigw3(ji,jj,jk)+rn_hc*zcoefw ) 
    1403                   gdep3w(ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*gsi3w3(ji,jj,jk)+rn_hc*zcoeft ) 
    1404                END DO 
    1405                ! 
    1406             END DO   ! for all jj's 
    1407          END DO    ! for all ji's 
    1408  
    1409          DO ji = 1, jpim1 
    1410             DO jj = 1, jpjm1 
    1411                DO jk = 1, jpk 
    1412                   esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*esigt3(ji,jj,jk)+hbatt(ji+1,jj)*esigt3(ji+1,jj,jk) )   & 
    1413                      &              / ( hbatt(ji,jj)+hbatt(ji+1,jj) ) 
    1414                   esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*esigt3(ji,jj,jk)+hbatt(ji,jj+1)*esigt3(ji,jj+1,jk) )   & 
    1415                      &              / ( hbatt(ji,jj)+hbatt(ji,jj+1) ) 
    1416                   esigtf3(ji,jj,jk) = ( hbatt(ji,jj)*esigt3(ji,jj,jk)+hbatt(ji+1,jj)*esigt3(ji+1,jj,jk)     & 
    1417                      &                + hbatt(ji,jj+1)*esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)*esigt3(ji+1,jj+1,jk) )   & 
    1418                      &              / ( hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) ) 
    1419                   esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*esigw3(ji,jj,jk)+hbatt(ji+1,jj)*esigw3(ji+1,jj,jk) )   & 
    1420                      &              / ( hbatt(ji,jj)+hbatt(ji+1,jj) ) 
    1421                   esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*esigw3(ji,jj,jk)+hbatt(ji,jj+1)*esigw3(ji,jj+1,jk) )   & 
    1422                      &              / ( hbatt(ji,jj)+hbatt(ji,jj+1) ) 
    1423                   ! 
    1424                   e3t(ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*esigt3 (ji,jj,jk) + rn_hc/FLOAT(jpkm1) ) 
    1425                   e3u(ji,jj,jk) = ( (hbatu(ji,jj)-rn_hc)*esigtu3(ji,jj,jk) + rn_hc/FLOAT(jpkm1) ) 
    1426                   e3v(ji,jj,jk) = ( (hbatv(ji,jj)-rn_hc)*esigtv3(ji,jj,jk) + rn_hc/FLOAT(jpkm1) ) 
    1427                   e3f(ji,jj,jk) = ( (hbatf(ji,jj)-rn_hc)*esigtf3(ji,jj,jk) + rn_hc/FLOAT(jpkm1) ) 
    1428                   ! 
    1429                   e3w (ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*esigw3 (ji,jj,jk) + rn_hc/FLOAT(jpkm1) ) 
    1430                   e3uw(ji,jj,jk) = ( (hbatu(ji,jj)-rn_hc)*esigwu3(ji,jj,jk) + rn_hc/FLOAT(jpkm1) ) 
    1431                   e3vw(ji,jj,jk) = ( (hbatv(ji,jj)-rn_hc)*esigwv3(ji,jj,jk) + rn_hc/FLOAT(jpkm1) ) 
    1432                END DO 
    1433             END DO 
    1434          END DO 
    1435  
    1436          CALL lbc_lnk( e3t , 'T', 1._wp ) 
    1437          CALL lbc_lnk( e3u , 'U', 1._wp ) 
    1438          CALL lbc_lnk( e3v , 'V', 1._wp ) 
    1439          CALL lbc_lnk( e3f , 'F', 1._wp ) 
    1440          CALL lbc_lnk( e3w , 'W', 1._wp ) 
    1441          CALL lbc_lnk( e3uw, 'U', 1._wp ) 
    1442          CALL lbc_lnk( e3vw, 'V', 1._wp ) 
    1443  
    1444          ! 
    1445       ELSE   ! not ln_s_sigma 
    1446          ! 
    1447          DO jk = 1, jpk 
    1448            gsigw(jk) = -fssig( REAL(jk,wp)-0.5_wp ) 
    1449            gsigt(jk) = -fssig( REAL(jk,wp)        ) 
    1450          END DO 
    1451          IF( nprint == 1 .AND. lwp )   WRITE(numout,*) 'gsigw 1 jpk    ', gsigw(1), gsigw(jpk) 
    1452          ! 
    1453          ! Coefficients for vertical scale factors at w-, t- levels 
    1454 !!gm bug :  define it from analytical function, not like juste bellow.... 
    1455 !!gm        or betteroffer the 2 possibilities.... 
    1456          DO jk = 1, jpkm1 
    1457             esigt(jk  ) = gsigw(jk+1) - gsigw(jk) 
    1458             esigw(jk+1) = gsigt(jk+1) - gsigt(jk) 
    1459          END DO 
    1460          esigw( 1 ) = 2._wp * ( gsigt(1  ) - gsigw(1  ) )  
    1461          esigt(jpk) = 2._wp * ( gsigt(jpk) - gsigw(jpk) ) 
    1462  
    1463 !!gm  original form 
    1464 !!org DO jk = 1, jpk 
    1465 !!org    esigt(jk)=fsdsig( FLOAT(jk)     ) 
    1466 !!org    esigw(jk)=fsdsig( FLOAT(jk)-0.5 ) 
    1467 !!org END DO 
    1468 !!gm 
    1469          ! 
    1470          ! Coefficients for vertical depth as the sum of e3w scale factors 
    1471          gsi3w(1) = 0.5_wp * esigw(1) 
    1472          DO jk = 2, jpk 
    1473             gsi3w(jk) = gsi3w(jk-1) + esigw(jk) 
    1474          END DO 
    1475 !!gm: depuw, depvw can be suppressed (modif in ldfslp) and depw=dep3w can be set (save 3 3D arrays) 
    1476          DO jk = 1, jpk 
    1477             zcoeft = ( REAL(jk,wp) - 0.5_wp ) / REAL(jpkm1,wp) 
    1478             zcoefw = ( REAL(jk,wp) - 1.0_wp ) / REAL(jpkm1,wp) 
    1479             gdept (:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*gsigt(jk) + hift(:,:)*zcoeft ) 
    1480             gdepw (:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*gsigw(jk) + hift(:,:)*zcoefw ) 
    1481             gdep3w(:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*gsi3w(jk) + hift(:,:)*zcoeft ) 
    1482          END DO 
    1483 !!gm: e3uw, e3vw can be suppressed  (modif in dynzdf, dynzdf_iso, zdfbfr) (save 2 3D arrays) 
    1484          DO jj = 1, jpj 
    1485             DO ji = 1, jpi 
    1486                DO jk = 1, jpk 
    1487                  e3t(ji,jj,jk) = ( (hbatt(ji,jj)-hift(ji,jj))*esigt(jk) + hift(ji,jj)/REAL(jpkm1,wp) ) 
    1488                  e3u(ji,jj,jk) = ( (hbatu(ji,jj)-hifu(ji,jj))*esigt(jk) + hifu(ji,jj)/REAL(jpkm1,wp) ) 
    1489                  e3v(ji,jj,jk) = ( (hbatv(ji,jj)-hifv(ji,jj))*esigt(jk) + hifv(ji,jj)/REAL(jpkm1,wp) ) 
    1490                  e3f(ji,jj,jk) = ( (hbatf(ji,jj)-hiff(ji,jj))*esigt(jk) + hiff(ji,jj)/REAL(jpkm1,wp) ) 
    1491                  ! 
    1492                  e3w (ji,jj,jk) = ( (hbatt(ji,jj)-hift(ji,jj))*esigw(jk) + hift(ji,jj)/REAL(jpkm1,wp) ) 
    1493                  e3uw(ji,jj,jk) = ( (hbatu(ji,jj)-hifu(ji,jj))*esigw(jk) + hifu(ji,jj)/REAL(jpkm1,wp) ) 
    1494                  e3vw(ji,jj,jk) = ( (hbatv(ji,jj)-hifv(ji,jj))*esigw(jk) + hifv(ji,jj)/REAL(jpkm1,wp) ) 
    1495                END DO 
    1496             END DO 
    1497          END DO 
    1498          ! 
    1499       ENDIF ! ln_s_sigma 
    1500  
    1501  
     1358 
     1359!======================================================================== 
     1360! Song and Haidvogel  1994 (ln_s_sh94=T) 
     1361! Siddorn and Furner 2012 (ln_sf12=T) 
     1362! or  tanh function       (both false)                     
     1363!======================================================================== 
     1364      IF      ( ln_s_sh94 ) THEN  
     1365                           CALL s_sh94() 
     1366      ELSE IF ( ln_s_sf12 ) THEN 
     1367                           CALL s_sf12() 
     1368      ELSE                  
     1369                           CALL s_tanh() 
     1370      ENDIF  
     1371 
     1372      CALL lbc_lnk( e3t , 'T', 1._wp ) 
     1373      CALL lbc_lnk( e3u , 'U', 1._wp ) 
     1374      CALL lbc_lnk( e3v , 'V', 1._wp ) 
     1375      CALL lbc_lnk( e3f , 'F', 1._wp ) 
     1376      CALL lbc_lnk( e3w , 'W', 1._wp ) 
     1377      CALL lbc_lnk( e3uw, 'U', 1._wp ) 
     1378      CALL lbc_lnk( e3vw, 'V', 1._wp ) 
     1379 
     1380      fsdepw(:,:,:) = gdepw (:,:,:) 
     1381      fsde3w(:,:,:) = gdep3w(:,:,:) 
    15021382      ! 
    15031383      where (e3t   (:,:,:).eq.0.0)  e3t(:,:,:) = 1.0 
     
    15571437            &                          ' w ', MAXVAL( fse3w (:,:,:) ) 
    15581438      ENDIF 
    1559       ! 
     1439      !  END DO 
    15601440      IF(lwp) THEN                                  ! selected vertical profiles 
    15611441         WRITE(numout,*) 
     
    15871467      ENDIF 
    15881468 
    1589 !!gm bug?  no more necessary?  if ! defined key_helsinki 
    1590       DO jk = 1, jpk 
     1469!================================================================================ 
     1470! check the coordinate makes sense 
     1471!================================================================================ 
     1472      DO ji = 1, jpi 
    15911473         DO jj = 1, jpj 
    1592             DO ji = 1, jpi 
    1593                IF( fse3w(ji,jj,jk) <= 0._wp .OR. fse3t(ji,jj,jk) <= 0._wp ) THEN 
    1594                   WRITE(ctmp1,*) 'zgr_sco :   e3w   or e3t   =< 0  at point (i,j,k)= ', ji, jj, jk 
    1595                   CALL ctl_stop( ctmp1 ) 
    1596                ENDIF 
    1597                IF( fsdepw(ji,jj,jk) < 0._wp .OR. fsdept(ji,jj,jk) < 0._wp ) THEN 
    1598                   WRITE(ctmp1,*) 'zgr_sco :   gdepw or gdept =< 0  at point (i,j,k)= ', ji, jj, jk 
    1599                   CALL ctl_stop( ctmp1 ) 
    1600                ENDIF 
    1601             END DO 
    1602          END DO 
    1603       END DO 
    1604 !!gm bug    #endif 
    1605       ! 
    1606       CALL wrk_dealloc( jpi, jpj,      zenv, ztmp, zmsk, zri, zrj, zhbat                           ) 
     1474 
     1475            IF( hbatt(ji,jj) > 0._wp) THEN 
     1476               DO jk = 1, mbathy(ji,jj) 
     1477                 ! check coordinate is monotonically increasing 
     1478                 IF (fse3w(ji,jj,jk) <= 0._wp .OR. fse3t(ji,jj,jk) <= 0._wp ) THEN 
     1479                    WRITE(ctmp1,*) 'ERROR zgr_sco :   e3w   or e3t   =< 0  at point (i,j,k)= ', ji, jj, jk 
     1480                    WRITE(numout,*) 'ERROR zgr_sco :   e3w   or e3t   =< 0  at point (i,j,k)= ', ji, jj, jk 
     1481                    WRITE(numout,*) 'e3w',fse3w(ji,jj,:) 
     1482                    WRITE(numout,*) 'e3t',fse3t(ji,jj,:) 
     1483                    CALL ctl_stop( ctmp1 ) 
     1484                 ENDIF 
     1485                 ! and check it has never gone negative 
     1486                 IF( fsdepw(ji,jj,jk) < 0._wp .OR. fsdept(ji,jj,jk) < 0._wp ) THEN 
     1487                    WRITE(ctmp1,*) 'ERROR zgr_sco :   gdepw or gdept =< 0  at point (i,j,k)= ', ji, jj, jk 
     1488                    WRITE(numout,*) 'ERROR zgr_sco :   gdepw   or gdept   =< 0  at point (i,j,k)= ', ji, jj, jk 
     1489                    WRITE(numout,*) 'gdepw',fsdepw(ji,jj,:) 
     1490                    WRITE(numout,*) 'gdept',fsdept(ji,jj,:) 
     1491                    CALL ctl_stop( ctmp1 ) 
     1492                 ENDIF 
     1493                 ! and check it never exceeds the total depth 
     1494                 IF( fsdepw(ji,jj,jk) > hbatt(ji,jj) ) THEN 
     1495                    WRITE(ctmp1,*) 'ERROR zgr_sco :   gdepw > hbatt  at point (i,j,k)= ', ji, jj, jk 
     1496                    WRITE(numout,*) 'ERROR zgr_sco :   gdepw > hbatt  at point (i,j,k)= ', ji, jj, jk 
     1497                    WRITE(numout,*) 'gdepw',fsdepw(ji,jj,:) 
     1498                    CALL ctl_stop( ctmp1 ) 
     1499                 ENDIF 
     1500               END DO 
     1501 
     1502               DO jk = 1, mbathy(ji,jj)-1 
     1503                 ! and check it never exceeds the total depth 
     1504                IF( fsdept(ji,jj,jk) > hbatt(ji,jj) ) THEN 
     1505                    WRITE(ctmp1,*) 'ERROR zgr_sco :   gdept > hbatt  at point (i,j,k)= ', ji, jj, jk 
     1506                    WRITE(numout,*) 'ERROR zgr_sco :   gdept > hbatt  at point (i,j,k)= ', ji, jj, jk 
     1507                    WRITE(numout,*) 'gdept',fsdept(ji,jj,:) 
     1508                    CALL ctl_stop( ctmp1 ) 
     1509                 ENDIF 
     1510               END DO 
     1511 
     1512            ENDIF 
     1513 
     1514         END DO 
     1515      END DO 
     1516      ! 
     1517      CALL wrk_dealloc( jpi, jpj,       ztmp, zmsk, zri, zrj, zhbat                           ) 
    16071518      CALL wrk_dealloc( jpi, jpj, jpk, gsigw3, gsigt3, gsi3w3                                      ) 
    16081519      CALL wrk_dealloc( jpi, jpj, jpk, esigt3, esigw3, esigtu3, esigtv3, esigtf3, esigwu3, esigwv3 ) 
     
    16121523   END SUBROUTINE zgr_sco 
    16131524 
     1525!!====================================================================== 
     1526   SUBROUTINE s_sh94() 
     1527 
     1528      !!---------------------------------------------------------------------- 
     1529      !!                  ***  ROUTINE s_sh94  *** 
     1530      !!                      
     1531      !! ** Purpose :   stretch the s-coordinate system 
     1532      !! 
     1533      !! ** Method  :   s-coordinate stretch using the Song and Haidvogel 1994 
     1534      !!                mixed S/sigma coordinate 
     1535      !! 
     1536      !! Reference : Song and Haidvogel 1994.  
     1537      !!---------------------------------------------------------------------- 
     1538      ! 
     1539      INTEGER  ::   ji, jj, jk           ! dummy loop argument 
     1540      REAL(wp) ::   zcoeft, zcoefw   ! temporary scalars 
     1541      ! 
     1542      REAL(wp), POINTER, DIMENSION(:,:,:) :: gsigw3, gsigt3, gsi3w3 
     1543      REAL(wp), POINTER, DIMENSION(:,:,:) :: esigt3, esigw3, esigtu3, esigtv3, esigtf3, esigwu3, esigwv3            
     1544 
     1545      CALL wrk_alloc( jpi, jpj, jpk, gsigw3, gsigt3, gsi3w3                                      ) 
     1546      CALL wrk_alloc( jpi, jpj, jpk, esigt3, esigw3, esigtu3, esigtv3, esigtf3, esigwu3, esigwv3 ) 
     1547 
     1548      gsigw3  = 0._wp   ;   gsigt3  = 0._wp   ;   gsi3w3  = 0._wp 
     1549      esigt3  = 0._wp   ;   esigw3  = 0._wp  
     1550      esigtu3 = 0._wp   ;   esigtv3 = 0._wp   ;   esigtf3 = 0._wp 
     1551      esigwu3 = 0._wp   ;   esigwv3 = 0._wp 
     1552 
     1553      DO ji = 1, jpi 
     1554         DO jj = 1, jpj 
     1555 
     1556            IF( hbatt(ji,jj) > rn_hc ) THEN    !deep water, stretched sigma 
     1557               DO jk = 1, jpk 
     1558                  gsigw3(ji,jj,jk) = -fssig1( REAL(jk,wp)-0.5_wp, rn_bb ) 
     1559                  gsigt3(ji,jj,jk) = -fssig1( REAL(jk,wp)       , rn_bb ) 
     1560               END DO 
     1561            ELSE ! shallow water, uniform sigma 
     1562               DO jk = 1, jpk 
     1563                  gsigw3(ji,jj,jk) =   REAL(jk-1,wp)            / REAL(jpk-1,wp) 
     1564                  gsigt3(ji,jj,jk) = ( REAL(jk-1,wp) + 0.5_wp ) / REAL(jpk-1,wp) 
     1565                  END DO 
     1566            ENDIF 
     1567            ! 
     1568            DO jk = 1, jpkm1 
     1569               esigt3(ji,jj,jk  ) = gsigw3(ji,jj,jk+1) - gsigw3(ji,jj,jk) 
     1570               esigw3(ji,jj,jk+1) = gsigt3(ji,jj,jk+1) - gsigt3(ji,jj,jk) 
     1571            END DO 
     1572            esigw3(ji,jj,1  ) = 2._wp * ( gsigt3(ji,jj,1  ) - gsigw3(ji,jj,1  ) ) 
     1573            esigt3(ji,jj,jpk) = 2._wp * ( gsigt3(ji,jj,jpk) - gsigw3(ji,jj,jpk) ) 
     1574            ! 
     1575            ! Coefficients for vertical depth as the sum of e3w scale factors 
     1576            gsi3w3(ji,jj,1) = 0.5_wp * esigw3(ji,jj,1) 
     1577            DO jk = 2, jpk 
     1578               gsi3w3(ji,jj,jk) = gsi3w3(ji,jj,jk-1) + esigw3(ji,jj,jk) 
     1579            END DO 
     1580            ! 
     1581            DO jk = 1, jpk 
     1582               zcoeft = ( REAL(jk,wp) - 0.5_wp ) / REAL(jpkm1,wp) 
     1583               zcoefw = ( REAL(jk,wp) - 1.0_wp ) / REAL(jpkm1,wp) 
     1584               gdept (ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*gsigt3(ji,jj,jk)+rn_hc*zcoeft ) 
     1585               gdepw (ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*gsigw3(ji,jj,jk)+rn_hc*zcoefw ) 
     1586               gdep3w(ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*gsi3w3(ji,jj,jk)+rn_hc*zcoeft ) 
     1587            END DO 
     1588           ! 
     1589         END DO   ! for all jj's 
     1590      END DO    ! for all ji's 
     1591 
     1592      DO ji = 1, jpim1 
     1593         DO jj = 1, jpjm1 
     1594            DO jk = 1, jpk 
     1595               esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*esigt3(ji,jj,jk)+hbatt(ji+1,jj)*esigt3(ji+1,jj,jk) )   & 
     1596                  &              / ( hbatt(ji,jj)+hbatt(ji+1,jj) ) 
     1597               esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*esigt3(ji,jj,jk)+hbatt(ji,jj+1)*esigt3(ji,jj+1,jk) )   & 
     1598                  &              / ( hbatt(ji,jj)+hbatt(ji,jj+1) ) 
     1599               esigtf3(ji,jj,jk) = ( hbatt(ji,jj)*esigt3(ji,jj,jk)+hbatt(ji+1,jj)*esigt3(ji+1,jj,jk)     & 
     1600                  &                + hbatt(ji,jj+1)*esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)*esigt3(ji+1,jj+1,jk) )   & 
     1601                  &              / ( hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) ) 
     1602               esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*esigw3(ji,jj,jk)+hbatt(ji+1,jj)*esigw3(ji+1,jj,jk) )   & 
     1603                  &              / ( hbatt(ji,jj)+hbatt(ji+1,jj) ) 
     1604               esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*esigw3(ji,jj,jk)+hbatt(ji,jj+1)*esigw3(ji,jj+1,jk) )   & 
     1605                  &              / ( hbatt(ji,jj)+hbatt(ji,jj+1) ) 
     1606               ! 
     1607               e3t(ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*esigt3 (ji,jj,jk) + rn_hc/REAL(jpkm1,wp) ) 
     1608               e3u(ji,jj,jk) = ( (hbatu(ji,jj)-rn_hc)*esigtu3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) ) 
     1609               e3v(ji,jj,jk) = ( (hbatv(ji,jj)-rn_hc)*esigtv3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) ) 
     1610               e3f(ji,jj,jk) = ( (hbatf(ji,jj)-rn_hc)*esigtf3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) ) 
     1611               ! 
     1612               e3w (ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*esigw3 (ji,jj,jk) + rn_hc/REAL(jpkm1,wp) ) 
     1613               e3uw(ji,jj,jk) = ( (hbatu(ji,jj)-rn_hc)*esigwu3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) ) 
     1614               e3vw(ji,jj,jk) = ( (hbatv(ji,jj)-rn_hc)*esigwv3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) ) 
     1615            END DO 
     1616        END DO 
     1617      END DO 
     1618 
     1619      CALL wrk_dealloc( jpi, jpj, jpk, gsigw3, gsigt3, gsi3w3                                      ) 
     1620      CALL wrk_dealloc( jpi, jpj, jpk, esigt3, esigw3, esigtu3, esigtv3, esigtf3, esigwu3, esigwv3 ) 
     1621 
     1622   END SUBROUTINE s_sh94 
     1623 
     1624   SUBROUTINE s_sf12 
     1625 
     1626      !!---------------------------------------------------------------------- 
     1627      !!                  ***  ROUTINE s_sf12 ***  
     1628      !!                      
     1629      !! ** Purpose :   stretch the s-coordinate system 
     1630      !! 
     1631      !! ** Method  :   s-coordinate stretch using the Siddorn and Furner 2012? 
     1632      !!                mixed S/sigma/Z coordinate 
     1633      !! 
     1634      !!                This method allows the maintenance of fixed surface and or 
     1635      !!                bottom cell resolutions (cf. geopotential coordinates)  
     1636      !!                within an analytically derived stretched S-coordinate framework. 
     1637      !! 
     1638      !! 
     1639      !! Reference : Siddorn and Furner 2012 (submitted Ocean modelling). 
     1640      !!---------------------------------------------------------------------- 
     1641      ! 
     1642      INTEGER  ::   ji, jj, jk           ! dummy loop argument 
     1643      REAL(wp) ::   fsmth          ! smoothing around critical depth 
     1644      REAL(wp) ::   zss, zbb       ! Surface and bottom cell thickness in sigma space 
     1645      ! 
     1646      REAL(wp), POINTER, DIMENSION(:,:,:) :: gsigw3, gsigt3, gsi3w3 
     1647      REAL(wp), POINTER, DIMENSION(:,:,:) :: esigt3, esigw3, esigtu3, esigtv3, esigtf3, esigwu3, esigwv3            
     1648 
     1649      !<