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

Changeset 7567


Ignore:
Timestamp:
2017-01-16T20:11:00+01:00 (8 years ago)
Author:
hadjt
Message:

CO6 version adapted for shelf seas climate projections, including added diagnostics

Location:
branches/UKMO/CO6_shelfclimate/NEMOGCM/NEMO
Files:
7 added
49 edited

Legend:

Unmodified
Added
Removed

  • branches/UKMO/CO6_shelfclimate/NEMOGCM/NEMO/OPA_SRC/BDY/bdy_oce.F90

    r7566 r7567  
    7171      REAL, POINTER, DIMENSION(:,:)   ::  ht_s  !: now snow thickness 
    7272#endif 
     73#if defined key_top 
     74      CHARACTER(LEN=20)                   :: cn_obc  !: type of boundary condition to apply 
     75      REAL(wp)                            :: rn_fac  !: multiplicative scaling factor 
     76      REAL(wp), POINTER, DIMENSION(:,:)   :: trc     !: now field of the tracer 
     77      LOGICAL                             :: dmp     !: obc damping term 
     78#endif 
     79 
    7380   END TYPE OBC_DATA 
    7481 
    … …  
    8390   LOGICAL                    ::   ln_mask_file             !: =T read bdymask from file 
    8491   LOGICAL                    ::   ln_vol                   !: =T volume correction              
     92   !JT 
     93   LOGICAL, DIMENSION(jp_bdy) ::   ln_sponge                !: =T use sponge layer  
     94   !JT 
    8595   ! 
    8696   INTEGER                    ::   nb_bdy                   !: number of open boundary sets 
    … …  
    101111   LOGICAL, DIMENSION(jp_bdy) ::   ln_tra_dmp               !: =T Tracer damping 
    102112   LOGICAL, DIMENSION(jp_bdy) ::   ln_dyn3d_dmp             !: =T Baroclinic velocity damping 
     113    
     114!   !JT 
     115   LOGICAL, DIMENSION(jp_bdy) ::   ln_ssh_bdy               !: =T USE SSH BDY - name list switch 
     116!   !JT 
     117    
    103118   REAL(wp),    DIMENSION(jp_bdy) ::   rn_time_dmp              !: Damping time scale in days 
    104119   REAL(wp),    DIMENSION(jp_bdy) ::   rn_time_dmp_out          !: Damping time scale in days at radiation outflow points 
     120   !JT 
     121   REAL(wp)                   ::   rn_sponge                  !: multiplier of diffusion for sponge layer 
     122   !JT 
    105123 
    106124   CHARACTER(len=20), DIMENSION(jp_bdy) ::   cn_ice_lim       ! Choice of boundary condition for sea ice variables  
    … …  
    118136   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), TARGET ::   bdyumask   !: Mask defining computational domain at U-points 
    119137   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), TARGET ::   bdyvmask   !: Mask defining computational domain at V-points 
     138   !JT 
     139   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   sponge_factor !: Multiplier for diffusion for sponge layer 
     140   !JT 
    120141 
    121142   REAL(wp)                                    ::   bdysurftot !: Lateral surface of unstructured open boundary 
    … …  
    147168      !!---------------------------------------------------------------------- 
    148169      ! 
    149       ALLOCATE( bdytmask(jpi,jpj) , bdyumask(jpi,jpj), bdyvmask(jpi,jpj),     &   
     170      !JT ALLOCATE( bdytmask(jpi,jpj) , bdyumask(jpi,jpj), bdyvmask(jpi,jpj),     &   
     171      ALLOCATE( bdytmask(jpi,jpj) , bdyumask(jpi,jpj), bdyvmask(jpi,jpj),sponge_factor(jpi,jpj),     &      
    150172         &      STAT=bdy_oce_alloc ) 
    151173      ! 
    … …  
    154176      bdyumask(:,:) = 1._wp 
    155177      bdyvmask(:,:) = 1._wp 
     178      !JT 
     179      sponge_factor(:,:) = 1._wp 
     180      !JT 
    156181      !  
    157182      IF( lk_mpp             )   CALL mpp_sum ( bdy_oce_alloc ) 

  • branches/UKMO/CO6_shelfclimate/NEMOGCM/NEMO/OPA_SRC/BDY/bdydta.F90

    r7566 r7567  
    3737#endif 
    3838   USE sbcapr 
     39#if defined key_top 
     40   USE par_trc 
     41   USE trc, ONLY: trn 
     42#endif 
    3943 
    4044   IMPLICIT NONE 
    … …  
    394398#endif 
    395399      ! end jchanut tschanges 
     400       
     401       
     402      !JT use sshn (ssh now) if ln_ssh_bdy set to false in the name list 
     403      DO ib_bdy = 1, nb_bdy    
     404        nblen => idx_bdy(ib_bdy)%nblen 
     405        dta => dta_bdy(ib_bdy) 
     406          
     407        ilen1(:) = nblen(:) 
     408        !JT IF( .NOT. dta%ll_ssh ) THEN  
     409        IF( .NOT. ln_ssh_bdy(ib_bdy) ) THEN  
     410          igrd = 1 ! t Grid 
     411          DO ib = 1, ilen1(igrd) 
     412              ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 
     413              ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 
     414              dta_bdy(ib_bdy)%ssh(ib) = sshn(ii,ij) * tmask(ii,ij,1)          
     415          END DO  
     416        END IF 
     417      END DO  
     418      !JT 
    396419 
    397420      IF ( ln_apr_obc ) THEN 
    … …  
    782805            IF( dta%ll_v2d ) ALLOCATE( dta%v2d(nblen(3)) ) 
    783806         ENDIF 
    784          IF ( nn_dyn2d_dta(ib_bdy) .eq. 1 .or. nn_dyn2d_dta(ib_bdy) .eq. 3 ) THEN 
    785             IF( dta%ll_ssh ) THEN 
    786                if(lwp) write(numout,*) '++++++ dta%ssh pointing to fnow' 
    787                jfld = jfld + 1 
    788                dta%ssh => bf(jfld)%fnow(:,1,1) 
    789             ENDIF 
     807         IF ( nn_dyn2d_dta(ib_bdy) .eq. 1 .or. nn_dyn2d_dta(ib_bdy) .eq. 3 ) THEN          
     808            !JT IF( dta%ll_ssh ) THEN 
     809            !JT    if(lwp) write(numout,*) '++++++ dta%ssh pointing to fnow' 
     810            !JT    jfld = jfld + 1 
     811            !JT    dta%ssh => bf(jfld)%fnow(:,1,1)    
     812            !JT ENDIF 
     813             
     814            !JT  
     815            !JT allocate ssh if dta%ll_ssh set too false, as may still use it 
     816            IF (dta%ll_ssh) THEN 
     817                IF( dta%ll_ssh ) THEN 
     818                  if(lwp) write(numout,*) '++++++ dta%ssh pointing to fnow' 
     819                  jfld = jfld + 1 
     820                  dta%ssh => bf(jfld)%fnow(:,1,1) 
     821                ENDIF 
     822            ELSE 
     823              if(lwp) write(numout,*) '++++++ dta%ssh allocated space' 
     824              !ALLOCATE( dta_bdy(ib_bdy)%ssh(nblen(1)) )             
     825              ALLOCATE( dta%ssh(nblen(1)) )             
     826            ENDIF 
     827            !JT if  
     828             
    790829            IF ( dta%ll_u2d ) THEN 
    791830               IF ( ln_full_vel_array(ib_bdy) ) THEN 
    … …  
    814853            IF( dta%ll_u3d ) ALLOCATE( dta_bdy(ib_bdy)%u3d(nblen(2),jpk) ) 
    815854            IF( dta%ll_v3d ) ALLOCATE( dta_bdy(ib_bdy)%v3d(nblen(3),jpk) ) 
    816          ENDIF 
     855         ENDIF         
    817856         IF ( nn_dyn3d_dta(ib_bdy) .eq. 1 .or. & 
    818857           &  ( ln_full_vel_array(ib_bdy) .and. ( nn_dyn2d_dta(ib_bdy) .eq. 1 .or. nn_dyn2d_dta(ib_bdy) .eq. 3 ) ) ) THEN 

  • branches/UKMO/CO6_shelfclimate/NEMOGCM/NEMO/OPA_SRC/BDY/bdydyn3d.F90

    r7566 r7567  
    3333   !!---------------------------------------------------------------------- 
    3434   !! NEMO/OPA 3.3 , NEMO Consortium (2010) 
    35    !! $Id$ 
     35   !! $Id$  
    3636   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 
    3737   !!---------------------------------------------------------------------- 
    … …  
    5959         CASE('specified') 
    6060            CALL bdy_dyn3d_spe( idx_bdy(ib_bdy), dta_bdy(ib_bdy), kt, ib_bdy ) 
     61         CASE('zerograd')  
     62            CALL bdy_dyn3d_zgrad( idx_bdy(ib_bdy), dta_bdy(ib_bdy), kt, ib_bdy ) 
    6163         CASE('zero') 
    6264            CALL bdy_dyn3d_zro( idx_bdy(ib_bdy), dta_bdy(ib_bdy), kt, ib_bdy ) 
     65         CASE('neumann') 
     66            CALL bdy_dyn3d_nmn( idx_bdy(ib_bdy), ib_bdy ) 
    6367         CASE('orlanski') 
    6468            CALL bdy_dyn3d_orlanski( idx_bdy(ib_bdy), dta_bdy(ib_bdy), ib_bdy, ll_npo=.false. ) 
    … …  
    117121 
    118122   END SUBROUTINE bdy_dyn3d_spe 
     123 
     124   SUBROUTINE bdy_dyn3d_zgrad( idx, dta, kt , ib_bdy ) 
     125      !!---------------------------------------------------------------------- 
     126      !!                  ***  SUBROUTINE bdy_dyn3d_zgrad  *** 
     127      !! 
     128      !! ** Purpose : - Enforce a zero gradient of normal velocity 
     129      !! 
     130      !!---------------------------------------------------------------------- 
     131      INTEGER                     ::   kt 
     132      TYPE(OBC_INDEX), INTENT(in) ::   idx  ! OBC indices 
     133      TYPE(OBC_DATA),  INTENT(in) ::   dta  ! OBC external data 
     134      INTEGER,         INTENT(in) ::   ib_bdy  ! BDY set index 
     135      !! 
     136      INTEGER  ::   jb, jk         ! dummy loop indices 
     137      INTEGER  ::   ii, ij, igrd   ! local integers 
     138      REAL(wp) ::   zwgt           ! boundary weight 
     139      INTEGER  ::   fu, fv 
     140      !!---------------------------------------------------------------------- 
     141      ! 
     142      IF( nn_timing == 1 ) CALL timing_start('bdy_dyn3d_zgrad') 
     143      ! 
     144      igrd = 2                      ! Copying tangential velocity into bdy points 
     145      DO jb = 1, idx%nblenrim(igrd) 
     146         DO jk = 1, jpkm1 
     147            ii   = idx%nbi(jb,igrd) 
     148            ij   = idx%nbj(jb,igrd) 
     149            fu   = ABS( ABS (NINT( idx%flagu(jb,igrd) ) ) - 1 ) 
     150            ua(ii,ij,jk) = ua(ii,ij,jk) * REAL( 1 - fu) + ( ua(ii,ij+fu,jk) * umask(ii,ij+fu,jk) & 
     151                        &+ ua(ii,ij-fu,jk) * umask(ii,ij-fu,jk) ) * umask(ii,ij,jk) * REAL( fu ) 
     152         END DO 
     153      END DO 
     154      ! 
     155      igrd = 3                      ! Copying tangential velocity into bdy points 
     156      DO jb = 1, idx%nblenrim(igrd) 
     157         DO jk = 1, jpkm1 
     158            ii   = idx%nbi(jb,igrd) 
     159            ij   = idx%nbj(jb,igrd) 
     160            fv   = ABS( ABS (NINT( idx%flagv(jb,igrd) ) ) - 1 ) 
     161            va(ii,ij,jk) = va(ii,ij,jk) * REAL( 1 - fv ) + ( va(ii+fv,ij,jk) * vmask(ii+fv,ij,jk) & 
     162                        &+ va(ii-fv,ij,jk) * vmask(ii-fv,ij,jk) ) * vmask(ii,ij,jk) * REAL( fv ) 
     163         END DO 
     164      END DO 
     165      CALL lbc_bdy_lnk( ua, 'U', -1., ib_bdy )   ! Boundary points should be updated   
     166      CALL lbc_bdy_lnk( va, 'V', -1., ib_bdy )    
     167      ! 
     168      IF( kt .eq. nit000 ) CLOSE( unit = 102 ) 
     169 
     170      IF( nn_timing == 1 ) CALL timing_stop('bdy_dyn3d_zgrad') 
     171 
     172   END SUBROUTINE bdy_dyn3d_zgrad 
    119173 
    120174   SUBROUTINE bdy_dyn3d_zro( idx, dta, kt, ib_bdy ) 
    … …  
    303357   END SUBROUTINE bdy_dyn3d_dmp 
    304358 
     359   SUBROUTINE bdy_dyn3d_nmn( idx, ib_bdy ) 
     360      !!---------------------------------------------------------------------- 
     361      !!                 ***  SUBROUTINE bdy_dyn3d_nmn  *** 
     362      !!              
     363      !!              - Apply Neumann condition to baroclinic velocities.  
     364      !!              - Wrapper routine for bdy_nmn 
     365      !!  
     366      !! 
     367      !!---------------------------------------------------------------------- 
     368      TYPE(OBC_INDEX),              INTENT(in) ::   idx  ! OBC indices 
     369      INTEGER,                      INTENT(in) ::   ib_bdy  ! BDY set index 
     370 
     371      INTEGER  ::   jb, igrd                               ! dummy loop indices 
     372      !!---------------------------------------------------------------------- 
     373 
     374      IF( nn_timing == 1 ) CALL timing_start('bdy_dyn3d_nmn') 
     375      ! 
     376      !! Note that at this stage the ub and ua arrays contain the baroclinic velocities.  
     377      ! 
     378      igrd = 2      ! Neumann bc on u-velocity;  
     379      !             
     380      CALL bdy_nmn( idx, igrd, ua ) 
     381 
     382      igrd = 3      ! Neumann bc on v-velocity 
     383      !   
     384      CALL bdy_nmn( idx, igrd, va ) 
     385      ! 
     386      CALL lbc_bdy_lnk( ua, 'U', -1., ib_bdy )    ! Boundary points should be updated 
     387      CALL lbc_bdy_lnk( va, 'V', -1., ib_bdy ) 
     388      ! 
     389      IF( nn_timing == 1 ) CALL timing_stop('bdy_dyn3d_nmn') 
     390      ! 
     391   END SUBROUTINE bdy_dyn3d_nmn 
    305392#else 
    306393   !!---------------------------------------------------------------------- 

  • branches/UKMO/CO6_shelfclimate/NEMOGCM/NEMO/OPA_SRC/BDY/bdyini.F90

    r7566 r7567  
    4949   !!---------------------------------------------------------------------- 
    5050   !! NEMO/OPA 4.0 , NEMO Consortium (2011) 
    51    !! $Id$ 
     51   !! $Id$  
    5252   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
    5353   !!---------------------------------------------------------------------- 
    … …  
    102102         &             cn_ice_lim, nn_ice_lim_dta,                           & 
    103103         &             rn_ice_tem, rn_ice_sal, rn_ice_age,                 & 
    104          &             ln_vol, nn_volctl, nn_rimwidth 
     104         &             ln_vol, nn_volctl, ln_sponge, rn_sponge, nn_rimwidth 
    105105      !! 
    106106      NAMELIST/nambdy_index/ ctypebdy, nbdyind, nbdybeg, nbdyend 
     107       
     108       
     109!      ! JT 
     110      NAMELIST/nambdy_ssh/ ln_ssh_bdy 
     111!      ! JT 
    107112      INTEGER  ::   ios                 ! Local integer output status for namelist read 
    108113      !!---------------------------------------------------------------------- 
    … …  
    132137902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nambdy in configuration namelist', lwp ) 
    133138      IF(lwm) WRITE ( numond, nambdy ) 
     139       
     140       
     141       
     142       
     143       
     144       
     145       
     146       
     147       
     148      !JT Read nambdy_ssh namelist  
     149      REWIND( numnam_ref )              ! Namelist nambdy in reference namelist :Unstructured open boundaries   
     150      READ  ( numnam_ref, nambdy_ssh, IOSTAT = ios, ERR = 905) 
     151905   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nambdy_ssh in reference namelist', lwp ) 
     152 
     153      REWIND( numnam_cfg )              ! Namelist nambdy in configuration namelist :Unstructured open boundaries 
     154      READ  ( numnam_cfg, nambdy_ssh, IOSTAT = ios, ERR = 906) 
     155906   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nambdy_ssh in configuration namelist', lwp ) 
     156      IF(lwm) WRITE ( numond, nambdy_ssh ) 
     157       
     158      IF(lwp) WRITE(numout,*) 
     159      IF(lwp) WRITE(numout,*) 'nambdy_ssh : use of ssh boundaries' 
     160      IF(lwp) WRITE(numout,*) '~~~~~~~~' 
     161      IF(lwp) WRITE(numout,*) '      ln_ssh_bdy: ' 
     162      DO ib_bdy = 1,nb_bdy 
     163        IF(lwp) WRITE(numout,*) '      ln_ssh_bdy(',ib_bdy,'): ',ln_ssh_bdy(ib_bdy) 
     164      ENDDO 
     165      IF(lwp) WRITE(numout,*) '~~~~~~~~' 
     166      IF(lwp) WRITE(numout,*)  
     167      !JT 
     168 
     169 
     170 
     171 
     172 
     173 
    134174 
    135175      ! ----------------------------------------- 
    … …  
    185225          CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for cn_dyn2d' ) 
    186226        END SELECT 
     227         
     228        !JT override dta_bdy(ib_bdy)%ll_ssh with namelist value (ln_ssh_bdy) 
     229        IF(lwp) WRITE(numout,*) 'nambdy_ssh : use of ssh boundaries' 
     230        IF(lwp) WRITE(numout,*) '~~~~~~~~' 
     231        IF(lwp) WRITE(numout,*) '      ib_bdy: ',ib_bdy 
     232        IF(lwp) WRITE(numout,*) '      Prior to Implementation of nambdy_ssh' 
     233        IF(lwp) WRITE(numout,*) '      dta_bdy(ib_bdy)%ll_ssh: ',dta_bdy(ib_bdy)%ll_ssh 
     234         
     235        dta_bdy(ib_bdy)%ll_ssh = ln_ssh_bdy(ib_bdy) 
     236         
     237        IF(lwp) WRITE(numout,*) '      After to Implementation of nambdy_ssh' 
     238        IF(lwp) WRITE(numout,*) '      dta_bdy(ib_bdy)%ll_ssh: ',dta_bdy(ib_bdy)%ll_ssh 
     239        IF(lwp) WRITE(numout,*) '~~~~~~~~' 
     240         
     241        !JT          
     242         
    187243        IF( cn_dyn2d(ib_bdy) /= 'none' ) THEN 
    188244           SELECT CASE( nn_dyn2d_dta(ib_bdy) )                   !  
    … …  
    213269             dta_bdy(ib_bdy)%ll_u3d = .true. 
    214270             dta_bdy(ib_bdy)%ll_v3d = .true. 
     271          CASE('neumann') 
     272             IF(lwp) WRITE(numout,*) '      Neumann conditions' 
     273             dta_bdy(ib_bdy)%ll_u3d = .false. 
     274             dta_bdy(ib_bdy)%ll_v3d = .false. 
     275          CASE('zerograd') 
     276             IF(lwp) WRITE(numout,*) '      Zero gradient for baroclinic velocities' 
     277             dta_bdy(ib_bdy)%ll_u3d = .false. 
     278             dta_bdy(ib_bdy)%ll_v3d = .false. 
    215279          CASE('zero') 
    216280             IF(lwp) WRITE(numout,*) '      Zero baroclinic velocities (runoff case)' 
    … …  
    365429        IF(lwp) WRITE(numout,*) '      Width of relaxation zone = ', nn_rimwidth(ib_bdy) 
    366430        IF(lwp) WRITE(numout,*) 
     431         
     432        IF( ln_sponge(ib_bdy) ) THEN                     ! check sponge layer choice  
     433          IF(lwp) WRITE(numout,*) '      Sponge layer applied at open boundaries'  
     434          IF(lwp) WRITE(numout,*) '      Multiplier for diffusion in sponge layer : ', rn_sponge  
     435          IF(lwp) WRITE(numout,*)  
     436        ELSE  
     437          IF(lwp) WRITE(numout,*) '      No Sponge layer applied at open boundaries'  
     438          IF(lwp) WRITE(numout,*)  
     439        ENDIF  
     440  
     441         
     442         
    367443 
    368444      ENDDO 
    … …  
    10921168            END DO 
    10931169         END DO  
     1170           
     1171          
     1172          !JT 
     1173         ! Compute multiplier for diffusion for sponge layer  
     1174         ! -------------------------------------------------  
     1175         IF( ln_sponge(ib_bdy) ) THEN  
     1176            igrd=1 
     1177            DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd)  
     1178               nbi => idx_bdy(ib_bdy)%nbi(ib,igrd)  
     1179               nbj => idx_bdy(ib_bdy)%nbj(ib,igrd)  
     1180               nbr => idx_bdy(ib_bdy)%nbr(ib,igrd)  
     1181               sponge_factor(nbi,nbj) = 1.0 + (rn_sponge-1.0) * ( 1.- TANH( FLOAT( nbr - 1 ) *0.5 ) )  
     1182            END DO  
     1183         ENDIF  
     1184          !JT 
     1185 
    10941186 
    10951187         ! Compute damping coefficients 

  • branches/UKMO/CO6_shelfclimate/NEMOGCM/NEMO/OPA_SRC/BDY/bdylib.F90

    r7566 r7567  
    2626   PUBLIC   bdy_orlanski_2d     ! routine called where? 
    2727   PUBLIC   bdy_orlanski_3d     ! routine called where? 
     28   PUBLIC   bdy_nmn     ! routine called where? 
    2829 
    2930   !!---------------------------------------------------------------------- 
    3031   !! NEMO/OPA 3.3 , NEMO Consortium (2010) 
    31    !! $Id$ 
     32   !! $Id$  
    3233   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 
    3334   !!---------------------------------------------------------------------- 
    … …  
    354355   END SUBROUTINE bdy_orlanski_3d 
    355356 
     357   SUBROUTINE bdy_nmn( idx, igrd, phia ) 
     358      !!---------------------------------------------------------------------- 
     359      !!                 ***  SUBROUTINE bdy_nmn  *** 
     360      !!                     
     361      !! ** Purpose : Duplicate the value at open boundaries, zero gradient. 
     362      !!  
     363      !!---------------------------------------------------------------------- 
     364      INTEGER,                    INTENT(in)     ::   igrd     ! grid index 
     365      REAL(wp), DIMENSION(:,:,:), INTENT(inout)  ::   phia     ! model after 3D field (to be updated) 
     366      TYPE(OBC_INDEX), INTENT(in) ::   idx  ! OBC indices 
     367      !!  
     368      REAL(wp) ::   zcoef, zcoef1, zcoef2 
     369      REAL(wp), POINTER, DIMENSION(:,:,:)        :: pmask      ! land/sea mask for field 
     370      REAL(wp), POINTER, DIMENSION(:,:)        :: bdypmask      ! land/sea mask for field 
     371      INTEGER  ::   ib, ik   ! dummy loop indices 
     372      INTEGER  ::   ii, ij, ip, jp   ! 2D addresses 
     373      !!---------------------------------------------------------------------- 
     374      ! 
     375      IF( nn_timing == 1 ) CALL timing_start('bdy_nmn') 
     376      ! 
     377      SELECT CASE(igrd) 
     378         CASE(1) 
     379            pmask => tmask(:,:,:) 
     380            bdypmask => bdytmask(:,:) 
     381         CASE(2) 
     382            pmask => umask(:,:,:) 
     383            bdypmask => bdyumask(:,:) 
     384         CASE(3) 
     385            pmask => vmask(:,:,:) 
     386            bdypmask => bdyvmask(:,:) 
     387         CASE DEFAULT ;   CALL ctl_stop( 'unrecognised value for igrd in bdy_nmn' ) 
     388      END SELECT 
     389      DO ib = 1, idx%nblenrim(igrd) 
     390         ii = idx%nbi(ib,igrd) 
     391         ij = idx%nbj(ib,igrd) 
     392         DO ik = 1, jpkm1 
     393            ! search the sense of the gradient 
     394            zcoef1 = bdypmask(ii-1,ij  )*pmask(ii-1,ij,ik) +  bdypmask(ii+1,ij  )*pmask(ii+1,ij,ik) 
     395            zcoef2 = bdypmask(ii  ,ij-1)*pmask(ii,ij-1,ik) +  bdypmask(ii  ,ij+1)*pmask(ii,ij+1,ik) 
     396            IF ( nint(zcoef1+zcoef2) == 0) THEN 
     397               ! corner **** we probably only want to set the tangentail component for the dynamics here 
     398               zcoef = pmask(ii-1,ij,ik) + pmask(ii+1,ij,ik) +  pmask(ii,ij-1,ik) +  pmask(ii,ij+1,ik) 
     399               IF (zcoef > .5_wp) THEN ! Only set none isolated points. 
     400                 phia(ii,ij,ik) = phia(ii-1,ij  ,ik) * pmask(ii-1,ij  ,ik) + & 
     401                   &              phia(ii+1,ij  ,ik) * pmask(ii+1,ij  ,ik) + & 
     402                   &              phia(ii  ,ij-1,ik) * pmask(ii  ,ij-1,ik) + & 
     403                   &              phia(ii  ,ij+1,ik) * pmask(ii  ,ij+1,ik) 
     404                 phia(ii,ij,ik) = ( phia(ii,ij,ik) / zcoef ) * pmask(ii,ij,ik) 
     405               ELSE 
     406                 phia(ii,ij,ik) = phia(ii,ij  ,ik) * pmask(ii,ij  ,ik) 
     407               ENDIF 
     408            ELSEIF ( nint(zcoef1+zcoef2) == 2) THEN 
     409               ! oblique corner **** we probably only want to set the normal component for the dynamics here 
     410               zcoef = pmask(ii-1,ij,ik)*bdypmask(ii-1,ij  ) + pmask(ii+1,ij,ik)*bdypmask(ii+1,ij  ) + & 
     411                   &   pmask(ii,ij-1,ik)*bdypmask(ii,ij -1 ) +  pmask(ii,ij+1,ik)*bdypmask(ii,ij+1  ) 
     412               phia(ii,ij,ik) = phia(ii-1,ij  ,ik) * pmask(ii-1,ij  ,ik)*bdypmask(ii-1,ij  ) + & 
     413                   &            phia(ii+1,ij  ,ik) * pmask(ii+1,ij  ,ik)*bdypmask(ii+1,ij  )  + & 
     414                   &            phia(ii  ,ij-1,ik) * pmask(ii  ,ij-1,ik)*bdypmask(ii,ij -1 ) + & 
     415                   &            phia(ii  ,ij+1,ik) * pmask(ii  ,ij+1,ik)*bdypmask(ii,ij+1  ) 
     416  
     417               phia(ii,ij,ik) = ( phia(ii,ij,ik) / MAX(1._wp, zcoef) ) * pmask(ii,ij,ik) 
     418            ELSE 
     419               ip = nint(bdypmask(ii+1,ij  )*pmask(ii+1,ij,ik) - bdypmask(ii-1,ij  )*pmask(ii-1,ij,ik)) 
     420               jp = nint(bdypmask(ii  ,ij+1)*pmask(ii,ij+1,ik) - bdypmask(ii  ,ij-1)*pmask(ii,ij-1,ik)) 
     421               phia(ii,ij,ik) = phia(ii+ip,ij+jp,ik) * pmask(ii+ip,ij+jp,ik) * pmask(ii,ij,ik) 
     422            ENDIF 
     423         END DO 
     424      END DO 
     425      ! 
     426      IF( nn_timing == 1 ) CALL timing_stop('bdy_nmn') 
     427      ! 
     428   END SUBROUTINE bdy_nmn 
    356429 
    357430#else 
    … …  
    366439      WRITE(*,*) 'bdy_orlanski_3d: You should not have seen this print! error?', kt 
    367440   END SUBROUTINE bdy_orlanski_3d 
     441   SUBROUTINE bdy_nmn( idx, igrd, phia )      ! Empty routine 
     442      WRITE(*,*) 'bdy_nmn: You should not have seen this print! error?', kt 
     443   END SUBROUTINE bdy_nmn 
    368444#endif 
    369445 

  • branches/UKMO/CO6_shelfclimate/NEMOGCM/NEMO/OPA_SRC/BDY/bdytides.F90

    r7566 r7567  
    102102 
    103103      REWIND(numnam_cfg) 
     104      REWIND(numnam_ref)   ! slwa 
    104105 
    105106      DO ib_bdy = 1, nb_bdy 

  • branches/UKMO/CO6_shelfclimate/NEMOGCM/NEMO/OPA_SRC/DIA/diaar5.F90

    r7566 r7567  
    3939   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:  ) ::   thick0       ! ocean thickness (interior domain) 
    4040   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   sn0          ! initial salinity 
     41   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   tn0          ! initial temperature 
    4142       
    4243   !! * Substitutions 
    … …  
    5657      !!---------------------------------------------------------------------- 
    5758      ! 
    58       ALLOCATE( area(jpi,jpj), thick0(jpi,jpj) , sn0(jpi,jpj,jpk) , STAT=dia_ar5_alloc ) 
     59      ALLOCATE( area(jpi,jpj), thick0(jpi,jpj) , sn0(jpi,jpj,jpk), tn0(jpi,jpj,jpk) , STAT=dia_ar5_alloc ) 
    5960      ! 
    6061      IF( lk_mpp             )   CALL mpp_sum ( dia_ar5_alloc ) 
    … …  
    121122      zssh_steric = - zarho / area_tot 
    122123      CALL iom_put( 'sshthster', zssh_steric ) 
     124      CALL iom_put( 'sshthster_mat', -zbotpres ) 
     125 
     126      !                      
     127      ztsn(:,:,:,jp_tem) = tn0(:,:,:)                    ! thermohaline ssh 
     128      ztsn(:,:,:,jp_sal) = tsn(:,:,:,jp_sal)  
     129      CALL eos( ztsn, zrhd, fsdept_n(:,:,:) )                       ! now in situ density using initial temperature 
     130      ! 
     131      zbotpres(:,:) = 0._wp                        ! no atmospheric surface pressure, levitating sea-ice 
     132      DO jk = 1, jpkm1 
     133         zbotpres(:,:) = zbotpres(:,:) + fse3t(:,:,jk) * zrhd(:,:,jk) 
     134      END DO 
     135      IF( .NOT.lk_vvl ) THEN 
     136         IF ( ln_isfcav ) THEN 
     137            DO ji=1,jpi 
     138               DO jj=1,jpj 
     139                  zbotpres(ji,jj) = zbotpres(ji,jj) + sshn(ji,jj) * zrhd(ji,jj,mikt(ji,jj)) + riceload(ji,jj) 
     140               END DO 
     141            END DO 
     142         ELSE 
     143            zbotpres(:,:) = zbotpres(:,:) + sshn(:,:) * zrhd(:,:,1) 
     144         END IF 
     145      END IF 
     146      !                                          
     147      zarho = SUM( area(:,:) * zbotpres(:,:) )  
     148      IF( lk_mpp )   CALL mpp_sum( zarho ) 
     149      zssh_steric = - zarho / area_tot 
     150      CALL iom_put( 'sshhlster', zssh_steric ) 
     151      !JT 
     152      CALL iom_put( 'sshhlster_mat', -zbotpres ) 
     153      !JT 
     154       
     155 
     156       
     157      !JT 
    123158       
    124159      !                                         ! steric sea surface height 
    … …  
    147182      zssh_steric = - zarho / area_tot 
    148183      CALL iom_put( 'sshsteric', zssh_steric ) 
     184      !JT 
     185      CALL iom_put( 'sshsteric_mat', -zbotpres ) 
     186      !JT 
    149187       
    150188      !                                         ! ocean bottom pressure 
    … …  
    211249      REAL(wp) ::   zztmp   
    212250      REAL(wp), POINTER, DIMENSION(:,:,:,:) ::   zsaldta   ! Jan/Dec levitus salinity 
     251      REAL(wp), POINTER, DIMENSION(:,:,:,:) ::   ztemdta   ! Jan/Dec levitus salinity 
    213252      ! reading initial file 
    214253      LOGICAL  ::   ln_tsd_init      !: T & S data flag 
    … …  
    234273      ! 
    235274      CALL wrk_alloc( jpi , jpj , jpk, jpts, zsaldta ) 
     275      CALL wrk_alloc( jpi , jpj , jpk, jpts, ztemdta ) 
    236276      !                                      ! allocate dia_ar5 arrays 
    237277      IF( dia_ar5_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'dia_ar5_init : unable to allocate arrays' ) 
    … …  
    253293      CALL iom_get  ( inum, jpdom_data, TRIM(sn_sal%clvar), zsaldta(:,:,:,2), 12 ) 
    254294      CALL iom_close( inum ) 
     295 
     296      CALL iom_open ( TRIM( cn_dir )//TRIM(sn_tem%clname), inum ) 
     297      CALL iom_get  ( inum, jpdom_data, TRIM(sn_tem%clvar), ztemdta(:,:,:,1), 1  ) 
     298      CALL iom_get  ( inum, jpdom_data, TRIM(sn_tem%clvar), ztemdta(:,:,:,2), 12 ) 
     299      CALL iom_close( inum ) 
    255300      sn0(:,:,:) = 0.5_wp * ( zsaldta(:,:,:,1) + zsaldta(:,:,:,2) )         
    256       sn0(:,:,:) = sn0(:,:,:) * tmask(:,:,:) 
     301      tn0(:,:,:) = 0.5_wp * ( ztemdta(:,:,:,1) + ztemdta(:,:,:,2) )         
     302      sn0(:,:,:) = sn0(:,:,:) * tmask(:,:,:)   
     303      tn0(:,:,:) = tn0(:,:,:) * tmask(:,:,:) 
    257304      IF( ln_zps ) THEN               ! z-coord. partial steps 
    258305         DO jj = 1, jpj               ! interpolation of salinity at the last ocean level (i.e. the partial step) 
    … …  
    262309                  zztmp = ( gdept_1d(ik) - gdept_0(ji,jj,ik) ) / ( gdept_1d(ik) - gdept_1d(ik-1) ) 
    263310                  sn0(ji,jj,ik) = ( 1._wp - zztmp ) * sn0(ji,jj,ik) + zztmp * sn0(ji,jj,ik-1) 
     311                  tn0(ji,jj,ik) = ( 1._wp - zztmp ) * tn0(ji,jj,ik) + zztmp * tn0(ji,jj,ik-1) 
    264312               ENDIF 
    265313            END DO 
    … …  
    268316      ! 
    269317      CALL wrk_dealloc( jpi , jpj , jpk, jpts, zsaldta ) 
     318      CALL wrk_dealloc( jpi , jpj , jpk, jpts, ztemdta ) 
    270319      ! 
    271320      IF( nn_timing == 1 )   CALL timing_stop('dia_ar5_init') 

  • branches/UKMO/CO6_shelfclimate/NEMOGCM/NEMO/OPA_SRC/DIA/diadct.F90

    r7566 r7567  
    5454  PUBLIC   dia_dct      ! routine called by step.F90 
    5555  PUBLIC   dia_dct_init ! routine called by opa.F90 
    56   PUBLIC   diadct_alloc ! routine called by nemo_init in nemogcm.F90  
     56  PUBLIC   diadct_alloc ! routine called by nemo_init in nemogcm.F90 
    5757  PRIVATE  readsec 
    5858  PRIVATE  removepoints 
    5959  PRIVATE  transport 
    6060  PRIVATE  dia_dct_wri 
     61  PRIVATE  dia_dct_wri_NOOS 
    6162 
    6263#include "domzgr_substitute.h90" 
    … …  
    6970  INTEGER :: nn_dctwri     ! Frequency of output 
    7071  INTEGER :: nn_secdebug   ! 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 
    7174    
    72   INTEGER, PARAMETER :: nb_class_max  = 10 
    73   INTEGER, PARAMETER :: nb_sec_max    = 150 
    74   INTEGER, PARAMETER :: nb_point_max  = 2000 
    75   INTEGER, PARAMETER :: nb_type_class = 10 
    76   INTEGER, PARAMETER :: nb_3d_vars    = 3  
    77   INTEGER, PARAMETER :: nb_2d_vars    = 2  
     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 
    7881  INTEGER            :: nb_sec  
    7982 
    … …  
    101104                                                     ztem            ,&! temperature classes(99 if you don't want) 
    102105                                                     zlay              ! level classes      (99 if you don't want) 
    103      REAL(wp), DIMENSION(nb_type_class,nb_class_max)  :: transport     ! transport output 
     106     REAL(wp), DIMENSION(nb_type,nb_class_max)        :: transport     ! transport output 
     107     REAL(wp), DIMENSION(nb_type,nb_class_max)        :: transport_h   ! transport output 
    104108     REAL(wp)                                         :: slopeSection  ! slope of the section 
    105109     INTEGER                                          :: nb_point      ! number of points in the section 
    … …  
    108112 
    109113  TYPE(SECTION),DIMENSION(nb_sec_max) :: secs ! Array of sections 
    110   
    111   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) ::  transports_3d  
    112   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:)   ::  transports_2d   
     114 
     115  REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) ::  transports_3d 
     116  REAL(wp), ALLOCATABLE, DIMENSION(:,:,:)   ::  transports_2d 
     117  REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) ::  transports_3d_h 
     118  REAL(wp), ALLOCATABLE, DIMENSION(:,:,:)   ::  transports_2d_h 
     119  REAL(wp), ALLOCATABLE, DIMENSION(:,:,:)   ::  z_hr_output 
    113120 
    114121   !! $Id$ 
    115122CONTAINS 
    116123 
    117   
    118   INTEGER FUNCTION diadct_alloc()  
    119      !!----------------------------------------------------------------------  
    120      !!                   ***  FUNCTION diadct_alloc  ***  
    121      !!----------------------------------------------------------------------  
    122      INTEGER :: ierr(2)  
    123      !!----------------------------------------------------------------------  
    124  
    125      ALLOCATE(transports_3d(nb_3d_vars,nb_sec_max,nb_point_max,jpk), STAT=ierr(1) )  
    126      ALLOCATE(transports_2d(nb_2d_vars,nb_sec_max,nb_point_max)    , STAT=ierr(2) )  
    127  
    128      diadct_alloc = MAXVAL( ierr )  
    129      IF( diadct_alloc /= 0 )   CALL ctl_warn('diadct_alloc: failed to allocate arrays')  
    130   
    131   END FUNCTION diadct_alloc  
     124  INTEGER FUNCTION diadct_alloc() 
     125     !!--------------------------------------------------------------------nb_sec-- 
     126     !!                   ***  FUNCTION diadct_alloc  *** 
     127     !!---------------------------------------------------------------------- 
     128     INTEGER :: ierr(2) 
     129     !!---------------------------------------------------------------------- 
     130     ! 
     131     ALLOCATE(transports_3d(nb_3d_vars,nb_sec_max,nb_point_max,jpk),   STAT=ierr(1) ) 
     132     ALLOCATE(transports_2d(nb_2d_vars,nb_sec_max,nb_point_max)    ,   STAT=ierr(2) ) 
     133     ALLOCATE(transports_3d_h(nb_3d_vars,nb_sec_max,nb_point_max,jpk), STAT=ierr(3) ) 
     134     ALLOCATE(transports_2d_h(nb_2d_vars,nb_sec_max,nb_point_max)    , STAT=ierr(4) ) 
     135     ALLOCATE(z_hr_output(nb_sec_max,24,nb_class_max)                , STAT=ierr(5) ) 
     136        ! 
     137     diadct_alloc = MAXVAL( ierr ) 
     138     IF( diadct_alloc /= 0 )   CALL ctl_warn('diadct_alloc: failed to allocate arrays') 
     139     ! 
     140  END FUNCTION diadct_alloc 
     141 
    132142 
    133143  SUBROUTINE dia_dct_init 
    … …  
    152162902  IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdct in configuration namelist', lwp ) 
    153163     IF(lwm) WRITE ( numond, namdct ) 
     164      
     165      
     166      
     167      
     168      
     169      
     170      
     171      
     172     !NAMELIST/namdct/nn_dct,nn_dctwri,nn_secdebug 
     173 
     174     !IF( nn_timing == 1 )   CALL timing_start('dia_dct_init') 
     175 
     176     !read namelist 
     177     !REWIND( numnam ) 
     178     !READ  ( numnam, namdct ) 
     179 
     180     IF( ln_NOOS ) THEN 
     181        nn_dct=3600./rdt         ! hard coded for NOOS transects, to give 25 hour means  
     182        nn_dctwri=86400./rdt 
     183        nn_dct_h=1       ! hard coded for NOOS transects, to give hourly data 
     184        nn_dctwri_h=3600./rdt 
     185     ENDIF 
    154186 
    155187     IF( lwp ) THEN 
    … …  
    157189        WRITE(numout,*) "diadct_init: compute transports through sections " 
    158190        WRITE(numout,*) "~~~~~~~~~~~~~~~~~~~~~" 
    159         WRITE(numout,*) "       Frequency of computation: nn_dct    = ",nn_dct 
    160         WRITE(numout,*) "       Frequency of write:       nn_dctwri = ",nn_dctwri 
     191        IF( ln_NOOS ) THEN 
     192           WRITE(numout,*) "       Frequency of computation hard coded to be every hour: nn_dct    = ",nn_dct 
     193           WRITE(numout,*) "       Frequency of write hard coded to average 25 instantaneous hour values: nn_dctwri = ",nn_dctwri 
     194           WRITE(numout,*) "       Frequency of hourly computation hard coded to be every timestep: nn_dct_h  = ",nn_dct_h 
     195           WRITE(numout,*) "       Frequency of hourly write hard coded to every hour: nn_dctwri_h = ",nn_dctwri_h 
     196        ELSE 
     197           WRITE(numout,*) "       Frequency of computation: nn_dct    = ",nn_dct 
     198           WRITE(numout,*) "       Frequency of write:       nn_dctwri = ",nn_dctwri 
     199        ENDIF 
    161200 
    162201        IF      ( nn_secdebug .GE. 1 .AND. nn_secdebug .LE. nb_sec_max )THEN 
    … …  
    174213     !Read section_ijglobal.diadct 
    175214     CALL readsec 
     215      
     216     !IF (lwp) write(numout,*) 'dct after readsec' 
     217      
    176218 
    177219     !open output file 
    178      IF( lwm ) THEN 
    179         CALL ctl_opn( numdct_vol,  'volume_transport', 'NEW', 'FORMATTED', 'SEQUENTIAL', -1, numout,  .FALSE. ) 
    180         CALL ctl_opn( numdct_heat, 'heat_transport'  , 'NEW', 'FORMATTED', 'SEQUENTIAL', -1, numout,  .FALSE. ) 
    181         CALL ctl_opn( numdct_salt, 'salt_transport'  , 'NEW', 'FORMATTED', 'SEQUENTIAL', -1, numout,  .FALSE. ) 
     220     IF( lwp ) THEN 
     221        IF( ln_NOOS ) THEN 
     222           CALL ctl_opn( numdct_NOOS  ,'NOOS_transport'  , 'NEW', 'FORMATTED', 'SEQUENTIAL', -1, numout,  .FALSE. ) 
     223           CALL ctl_opn( numdct_NOOS_h,'NOOS_transport_h', 'NEW', 'FORMATTED', 'SEQUENTIAL', -1, numout,  .FALSE. ) 
     224        ELSE 
     225           CALL ctl_opn( numdct_vol , 'volume_transport', 'NEW', 'FORMATTED', 'SEQUENTIAL', -1, numout,  .FALSE. ) 
     226           CALL ctl_opn( numdct_temp, 'heat_transport'  , 'NEW', 'FORMATTED', 'SEQUENTIAL', -1, numout,  .FALSE. ) 
     227           CALL ctl_opn( numdct_sal , 'salt_transport'  , 'NEW', 'FORMATTED', 'SEQUENTIAL', -1, numout,  .FALSE. ) 
     228        ENDIF 
    182229     ENDIF 
    183230 
    184      ! Initialise arrays to zero  
    185      transports_3d(:,:,:,:)=0.0  
    186      transports_2d(:,:,:)  =0.0  
     231     ! Initialise arrays to zero 
     232     transports_3d(:,:,:,:)  =0._wp 
     233     transports_2d(:,:,:)    =0._wp 
     234     transports_3d_h(:,:,:,:)=0._wp 
     235     transports_2d_h(:,:,:)  =0._wp 
     236     z_hr_output(:,:,:)      =0._wp 
    187237 
    188238     IF( nn_timing == 1 )   CALL timing_stop('dia_dct_init') 
    … …  
    195245     !!               ***  ROUTINE diadct  ***   
    196246     !! 
    197      !!  Purpose :: Compute section transports and write it in numdct files  
    198      !!    
    199      !!  Method  :: All arrays initialised to zero in dct_init  
    200      !!             Each nn_dct time step call subroutine 'transports' for  
    201      !!               each section to sum the transports over each grid cell.  
    202      !!             Each nn_dctwri time step:  
    203      !!               Divide the arrays by the number of summations to gain  
    204      !!               an average value  
    205      !!               Call dia_dct_sum to sum relevant grid boxes to obtain  
    206      !!               totals for each class (density, depth, temp or sal)  
    207      !!               Call dia_dct_wri to write the transports into file  
    208      !!               Reinitialise all relevant arrays to zero  
     247     !!  Purpose :: Compute section transports and write it in numdct files 
     248     !!   
     249     !!  Method  :: All arrays initialised to zero in dct_init 
     250     !!             Each nn_dct time step call subroutine 'transports' for 
     251     !!               each section to sum the transports. 
     252     !!             Each nn_dctwri time step: 
     253     !!               Divide the arrays by the number of summations to gain 
     254     !!               an average value 
     255     !!               Call dia_dct_sum to sum relevant grid boxes to obtain 
     256     !!               totals for each class (density, depth, temp or sal) 
     257     !!               Call dia_dct_wri to write the transports into file 
     258     !!               Reinitialise all relevant arrays to zero 
    209259     !!--------------------------------------------------------------------- 
    210260     !! * Arguments 
    … …  
    213263     !! * Local variables 
    214264     INTEGER             :: jsec,            &! loop on sections 
    215                             itotal            ! nb_sec_max*nb_type_class*nb_class_max 
     265                            itotal            ! nb_sec_max*nb_type*nb_class_max 
    216266     LOGICAL             :: lldebug =.FALSE.  ! debug a section   
    217       
    218      INTEGER , DIMENSION(1)             :: ish   ! tmp array for mpp_sum 
    219      INTEGER , DIMENSION(3)             :: ish2  !   " 
    220      REAL(wp), POINTER, DIMENSION(:)    :: zwork !   "   
    221      REAL(wp), POINTER, DIMENSION(:,:,:):: zsum  !   " 
     267 
     268      
     269     INTEGER , DIMENSION(1)             :: ish      ! tmp array for mpp_sum 
     270     INTEGER , DIMENSION(3)             :: ish2     !   " 
     271     REAL(wp), POINTER, DIMENSION(:)    :: zwork    !   "  
     272     REAL(wp), POINTER, DIMENSION(:,:,:):: zsum     !   " 
     273 
     274 
    222275 
    223276     !!---------------------------------------------------------------------     
     277      
     278      
     279     !WRITE(numout,*) "diadct ind ii: ",kt,nproc, narea, mig(1),mig(jpi),nlci 
     280     !WRITE(numout,*) "diadct ind jj: ",kt,nproc, narea, mjg(1),mig(jpj),nlcj 
     281 
     282      
    224283     IF( nn_timing == 1 )   CALL timing_start('dia_dct') 
    225284 
    226285     IF( lk_mpp )THEN 
    227         itotal = nb_sec_max*nb_type_class*nb_class_max 
    228         CALL wrk_alloc( itotal                                , zwork )  
    229         CALL wrk_alloc( nb_sec_max,nb_type_class,nb_class_max , zsum  ) 
     286        itotal = nb_sec_max*nb_type*nb_class_max 
     287        CALL wrk_alloc( itotal                          , zwork )  
     288        CALL wrk_alloc( nb_sec_max,nb_type,nb_class_max , zsum  ) 
    230289     ENDIF     
    231290  
    … …  
    239298         WRITE(numout,*) "~~~~~~~~~~~~~~~~~~~~~~~~~" 
    240299         WRITE(numout,*) "nb_sec = ",nb_sec 
     300         WRITE(numout,*) "nn_dct = ",nn_dct 
     301         WRITE(numout,*) "ln_NOOS = ",ln_NOOS 
     302         WRITE(numout,*) "nb_sec = ",nb_sec 
     303         WRITE(numout,*) "nb_sec_max = ",nb_sec_max 
     304         WRITE(numout,*) "nb_type = ",nb_type 
     305         WRITE(numout,*) "nb_class_max = ",nb_class_max 
    241306     ENDIF 
    242307 
    243308  
    244309     ! Compute transport and write only at nn_dctwri 
    245      IF( MOD(kt,nn_dct)==0 ) THEN  
     310     IF ( MOD(kt,nn_dct)==0 .or. &               ! compute transport every nn_dct time steps 
     311         (ln_NOOS .and. kt==nn_it000 ) )  THEN   ! also include first time step when calculating NOOS 25 hour averages 
    246312 
    247313        DO jsec=1,nb_sec 
    248314 
    249            !debug this section computing ? 
    250315           lldebug=.FALSE. 
    251316           IF( (jsec==nn_secdebug .OR. nn_secdebug==-1) .AND.  kt==nit000+nn_dct-1 .AND. lwp ) lldebug=.TRUE.  
    … …  
    253318           !Compute transport through section   
    254319           CALL transport(secs(jsec),lldebug,jsec)  
     320            
     321           !IF( lwp ) WRITE(numout,*) "diadct: call transport subroutine (kt, jsec) : ",kt,jsec 
    255322 
    256323        ENDDO 
     324        !IF( lwp ) WRITE(numout,*) "diadct: called transport subroutine (kt, nb_sec) : ",kt, nb_sec 
    257325              
    258326        IF( MOD(kt,nn_dctwri)==0 )THEN 
    259327 
    260            IF( lwp .AND. kt==nit000+nn_dctwri-1 )WRITE(numout,*)"      diadct: average transports and write at kt = ",kt          
     328           IF( lwp .AND. kt==nit000+nn_dctwri-1 )WRITE(numout,*)"      diadct: average and write at kt = ",kt          
    261329   
    262            !! divide arrays by nn_dctwri/nn_dct to obtain average  
    263            transports_3d(:,:,:,:)=transports_3d(:,:,:,:)/(nn_dctwri/nn_dct)  
    264            transports_2d(:,:,:)  =transports_2d(:,:,:)  /(nn_dctwri/nn_dct)  
    265   
    266            ! Sum over each class  
    267            DO jsec=1,nb_sec  
    268               CALL dia_dct_sum(secs(jsec),jsec)  
    269            ENDDO  
    270  
     330           !! divide arrays by nn_dctwri/nn_dct to obtain average 
     331           transports_3d(:,:,:,:)=transports_3d(:,:,:,:)/(nn_dctwri/nn_dct) 
     332           transports_2d(:,:,:)  =transports_2d(:,:,:)  /(nn_dctwri/nn_dct) 
     333 
     334           ! Sum over each class 
     335           DO jsec=1,nb_sec 
     336              CALL dia_dct_sum(secs(jsec),jsec) 
     337           ENDDO 
     338  
    271339           !Sum on all procs  
    272340           IF( lk_mpp )THEN 
    273               ish(1)  =  nb_sec_max*nb_type_class*nb_class_max  
    274               ish2    = (/nb_sec_max,nb_type_class,nb_class_max/) 
     341              zsum(:,:,:)=0.0_wp 
     342              ish(1)  =  nb_sec_max*nb_type*nb_class_max  
     343              ish2    = (/nb_sec_max,nb_type,nb_class_max/) 
    275344              DO jsec=1,nb_sec ; zsum(jsec,:,:) = secs(jsec)%transport(:,:) ; ENDDO 
    276345              zwork(:)= RESHAPE(zsum(:,:,:), ish ) 
    … …  
    283352           DO jsec=1,nb_sec 
    284353 
    285               IF( lwm )CALL dia_dct_wri(kt,jsec,secs(jsec)) 
     354              IF( lwp .and. .not. ln_NOOS )CALL dia_dct_wri(kt,jsec,secs(jsec)) 
     355              IF( lwp .and.       ln_NOOS )CALL dia_dct_wri_NOOS(kt,jsec,secs(jsec))   ! use NOOS specific formatting 
     356             
    286357             
    287358              !nullify transports values after writing 
    288               transports_3d(:,jsec,:,:)=0. 
    289               transports_2d(:,jsec,:  )=0. 
     359              transports_3d(:,jsec,:,:)=0.0 
     360              transports_2d(:,jsec,:  )=0.0 
    290361              secs(jsec)%transport(:,:)=0.   
     362              IF ( ln_NOOS ) CALL transport(secs(jsec),lldebug,jsec)  ! reinitialise for next 25 hour instantaneous average (overlapping values) 
    291363 
    292364           ENDDO 
    … …  
    295367 
    296368     ENDIF 
     369     IF ( MOD(kt,nn_dct_h)==0 ) THEN            ! compute transport every nn_dct_h time steps 
     370     !IF ( MOD(kt,nn_dct)==0 .or. &               ! compute transport every nn_dct time steps 
     371     !    (ln_NOOS .and. kt==nn_it000 ) )  THEN   ! also include first time step when calculating NOOS 25 hour averages 
     372 
     373        DO jsec=1,nb_sec 
     374 
     375           lldebug=.FALSE. 
     376           IF( (jsec==nn_secdebug .OR. nn_secdebug==-1) .AND.  kt==nit000+nn_dct_h-1 .AND. lwp ) lldebug=.TRUE.  
     377 
     378           !Compute transport through section   
     379           CALL transport_h(secs(jsec),lldebug,jsec)  
     380 
     381        ENDDO 
     382              
     383        IF( MOD(kt,nn_dctwri_h)==0 )THEN 
     384 
     385           IF( lwp .AND. kt==nit000+nn_dctwri_h-1 )WRITE(numout,*)"      diadct: average and write hourly files at kt = ",kt          
     386   
     387           !! divide arrays by nn_dctwri/nn_dct to obtain average 
     388           transports_3d_h(:,:,:,:)=transports_3d_h(:,:,:,:)/(nn_dctwri_h/nn_dct_h) 
     389           transports_2d_h(:,:,:)  =transports_2d_h(:,:,:)  /(nn_dctwri_h/nn_dct_h) 
     390 
     391           ! Sum over each class 
     392           DO jsec=1,nb_sec 
     393              CALL dia_dct_sum_h(secs(jsec),jsec) 
     394           ENDDO 
     395  
     396           !Sum on all procs  
     397          IF( lk_mpp )THEN 
     398              ish(1)  =  nb_sec_max*nb_type*nb_class_max  
     399              ish2    = (/nb_sec_max,nb_type,nb_class_max/) 
     400              DO jsec=1,nb_sec ; zsum(jsec,:,:) = secs(jsec)%transport_h(:,:) ; ENDDO 
     401              zwork(:)= RESHAPE(zsum(:,:,:), ish ) 
     402              CALL mpp_sum(zwork, ish(1)) 
     403              zsum(:,:,:)= RESHAPE(zwork,ish2) 
     404              DO jsec=1,nb_sec ; secs(jsec)%transport_h(:,:) = zsum(jsec,:,:) ; ENDDO 
     405           ENDIF 
     406 
     407           !Write the transport 
     408           DO jsec=1,nb_sec 
     409 
     410              IF( lwp .and. ln_NOOS ) THEN 
     411                !WRITE(numout,*) "diadct: call dia_dct_wri_NOOS_h (kt,nn_dctwri_h, kt/nn_dctwri_h,jsec) : ",kt,nn_dctwri_h, kt/nn_dctwri_h,jsec 
     412                CALL dia_dct_wri_NOOS_h(kt/nn_dctwri_h,jsec,secs(jsec))   ! use NOOS specific formatting 
     413                !WRITE(numout,*) "diadct: called dia_dct_wri_NOOS_h (kt,jsec) : ",kt,jsec 
     414              endif 
     415              !nullify transports values after writing 
     416              transports_3d_h(:,jsec,:,:)=0.0 
     417              transports_2d_h(:,jsec,:)=0.0 
     418              secs(jsec)%transport_h(:,:)=0.   
     419              IF ( ln_NOOS ) CALL transport_h(secs(jsec),lldebug,jsec)  ! reinitialise for next 25 hour instantaneous average (overlapping values) 
     420 
     421           ENDDO 
     422 
     423        ENDIF  
     424 
     425     ENDIF     
    297426 
    298427     IF( lk_mpp )THEN 
    299         itotal = nb_sec_max*nb_type_class*nb_class_max 
    300         CALL wrk_dealloc( itotal                                , zwork )  
    301         CALL wrk_dealloc( nb_sec_max,nb_type_class,nb_class_max , zsum  ) 
     428        itotal = nb_sec_max*nb_type*nb_class_max 
     429        CALL wrk_dealloc( itotal                          , zwork )  
     430        CALL wrk_dealloc( nb_sec_max,nb_type,nb_class_max , zsum  ) 
    302431     ENDIF     
    303432 
    … …  
    320449     INTEGER :: isec, iiglo, ijglo, iiloc, ijloc,iost,i1 ,i2  ! temporary  integer 
    321450     INTEGER :: jsec, jpt                                     ! dummy loop indices 
     451                                                              ! heat/salt tranport is actived 
    322452 
    323453     INTEGER, DIMENSION(2) :: icoord  
    … …  
    336466     !open input file 
    337467     !--------------- 
    338      CALL ctl_opn( numdct_in, 'section_ijglobal.diadct', 'OLD', 'UNFORMATTED', 'SEQUENTIAL', -1, numout, .FALSE. ) 
    339   
     468     !IF (lwp) THEN 
     469       !write(numout,*) 'dct low-level pre open: little endian ' 
     470       !OPEN(UNIT=107,FILE='section_ijglobal.diadct', FORM='UNFORMATTED', ACCESS='SEQUENTIAL', STATUS='OLD',convert='LITTLE_ENDIAN') 
     471        
     472       write(numout,*) 'dct low-level pre open: big endian :',nproc,narea 
     473       OPEN(UNIT=107,FILE='section_ijglobal.diadct', FORM='UNFORMATTED', ACCESS='SEQUENTIAL', STATUS='OLD',convert='BIG_ENDIAN') 
     474        
     475       !write(numout,*) 'dct low-level pre open: SWAP ' 
     476       !OPEN(UNIT=107,FILE='section_ijglobal.diadct', FORM='UNFORMATTED', ACCESS='SEQUENTIAL', STATUS='OLD',convert='SWAP') 
     477        
     478       !write(numout,*) 'dct low-level pre open: NATIVE ' 
     479       !OPEN(UNIT=107,FILE='section_ijglobal.diadct', FORM='UNFORMATTED', ACCESS='SEQUENTIAL', STATUS='OLD',convert='NATIVE') 
     480        
     481       write(numout,*) 'dct low-level opened :',nproc,narea 
     482       READ(107) isec 
     483       write(numout,*) 'dct low-level isec ', isec,nproc,narea 
     484       CLOSE(107) 
     485     !ENDIF 
     486      
     487      
     488     CALL ctl_opn( numdct_in, 'section_ijglobal.diadct', 'OLD', 'UNFORMATTED', 'SEQUENTIAL', -1, numout, .TRUE. ) 
     489  
     490      
    340491     !--------------- 
    341492     !Read input file 
    … …  
    350501        !--------------- 
    351502        secs(jsec)%name='' 
    352         secs(jsec)%llstrpond    = .FALSE.  ; secs(jsec)%ll_ice_section = .FALSE. 
    353         secs(jsec)%ll_date_line = .FALSE.  ; secs(jsec)%nb_class       = 0 
    354         secs(jsec)%zsigi        = 99._wp   ; secs(jsec)%zsigp          = 99._wp 
    355         secs(jsec)%zsal         = 99._wp   ; secs(jsec)%ztem           = 99._wp 
    356         secs(jsec)%zlay         = 99._wp          
    357         secs(jsec)%transport    =  0._wp   ; secs(jsec)%nb_point       = 0 
     503        !IF( lwp )  write(numout,*)  secs(jsec)%name 
     504        secs(jsec)%llstrpond      = .FALSE. 
     505        !IF( lwp )  write(numout,*) secs(jsec)%llstrpond  
     506        secs(jsec)%ll_ice_section = .FALSE. 
     507        !IF( lwp )  write(numout,*) secs 
     508        secs(jsec)%ll_date_line   = .FALSE. 
     509        !IF( lwp )  write(numout,*) secs(jsec)%ll_date_line 
     510        secs(jsec)%nb_class       = 0 
     511        !IF( lwp )  write(numout,*) secs(jsec)%nb_class 
     512        secs(jsec)%zsigi          = 99._wp 
     513        !IF( lwp )  write(numout,*) secs(jsec)%zsigi 
     514        secs(jsec)%zsigp          = 99._wp 
     515        !IF( lwp )  write(numout,*)  secs(jsec)%zsigp 
     516        secs(jsec)%zsal           = 99._wp 
     517        !IF( lwp )  write(numout,*) secs(jsec)%zsal 
     518        secs(jsec)%ztem           = 99._wp 
     519        !IF( lwp )  write(numout,*) secs(jsec)%ztem 
     520        secs(jsec)%zlay           = 99._wp 
     521        !IF( lwp )  write(numout,*) secs(jsec)%zlay 
     522        secs(jsec)%transport      =  0._wp 
     523        !IF( lwp )  write(numout,*) secs(jsec)%transport 
     524        secs(jsec)%transport_h    =  0._wp 
     525        !IF( lwp )  write(numout,*) secs(jsec)%transport_h 
     526        secs(jsec)%nb_point       = 0 
     527        !IF( lwp )  write(numout,*) secs(jsec)%nb_point  
    358528 
    359529        !read section's number / name / computing choices / classes / slopeSection / points number 
    360530        !----------------------------------------------------------------------------------------- 
    361         READ(numdct_in,iostat=iost)isec 
    362         IF (iost .NE. 0 )EXIT !end of file  
     531         
     532        !write(numout,*) 'dct isec ', isec 
     533        !write(numout,*) 'dct numdct_in ', numdct_in 
     534        READ(numdct_in,iostat=iost) isec 
     535        !write(numout,*) 'dct iost ', iost 
     536        !IF (iost .NE. 0 ) EXIT 
     537        IF (iost .NE. 0 ) then 
     538          write(numout,*) 'unable to read section_ijglobal.diadct. iost = ',iost 
     539          !nb_sec = 2 
     540          EXIT !end of file  
     541        ENDIF 
     542         
     543        !write(numout,*) 'dct isec', isec 
     544         
    363545        WRITE(cltmp,'(a,i4.4,a,i4.4)')'diadct: read sections : Problem of section number: isec= ',isec,' and jsec= ',jsec 
     546        !WRITE(numout,'(a,i4.4,a,i4.4)')'diadct: read sections : Problem of section number: isec= ',isec,' and jsec= ',jsec 
     547         
     548        !write(numout,*) 'dct isec, jsec', isec,jsec 
     549         
     550         
    364551        IF( jsec .NE. isec )  CALL ctl_stop( cltmp ) 
    365552 
    366         IF( lwp .AND. ( jsec==nn_secdebug .OR. nn_secdebug==-1 ) )WRITE(numout,*)"isec ",isec  
     553        !IF( lwp .AND. ( jsec==nn_secdebug .OR. nn_secdebug==-1 ) )WRITE(numout,*)"isec ",isec  
    367554 
    368555        READ(numdct_in)secs(jsec)%name 
    … …  
    388575 
    389576            WRITE(numout,*)       "   Section name :                       ",TRIM(secs(jsec)%name) 
    390             WRITE(numout,*)       "      Compute heat and salt transport ? ",secs(jsec)%llstrpond 
     577            WRITE(numout,*)       "      Compute temperature and salinity transports ? ",secs(jsec)%llstrpond 
    391578            WRITE(numout,*)       "      Compute ice transport ?           ",secs(jsec)%ll_ice_section 
    392579            WRITE(numout,*)       "      Section crosses date-line ?       ",secs(jsec)%ll_date_line 
    … …  
    399586            WRITE(numout,clformat)"      Temperature classes :             ",secs(jsec)%ztem 
    400587            WRITE(numout,clformat)"      Depth classes :                   ",secs(jsec)%zlay 
    401         ENDIF                
     588        ENDIF      
     589         
    402590 
    403591        IF( iptglo .NE. 0 )THEN 
    … …  
    411599              coordtemp(jpt)%I = i1  
    412600              coordtemp(jpt)%J = i2 
     601              !WRITE(numout,*)'diadct: coordtemp:', jpt,i1,i2 
    413602           ENDDO 
    414603           READ(numdct_in)directemp(1:iptglo) 
    … …  
    432621 
    433622              IF( iiglo==jpidta .AND. nimpp==1 ) iiglo = 2 
    434  
     623               
     624              !WRITE(numout,*)"diadct readsec: reg/glo domain:",narea,nlci,nimpp,mig(1),nlcj,njmpp,mjg(1) 
     625 
     626      !diadct init: reg/glo domain: 3*1,  24,  22,  24,  22,  2*1 
    435627              iiloc=iiglo-jpizoom+1-nimpp+1   ! local coordinates of the point 
    436628              ijloc=ijglo-jpjzoom+1-njmpp+1   !  " 
    437  
     629              !WRITE(numout,*)"      List of limits of each processor:",jpt,nimpp,iiloc,iiglo,jpizoom+1,nimpp+1,ijloc,ijglo,jpjzoom+1,njmpp+1,nlei,nlej 
     630 
     631                
     632              !WRITE(*,*)"diadct readsec: ii: ",narea, jpt,iiglo,mig(1),mig(jpi),nlci 
     633              !WRITE(*,*)"diadct readsec: jj: ",narea, jpt,ijglo,mjg(1),mjg(jpj),nlcj 
     634              !WRITE(*,*)"diadct readsec: log: ",narea, iiglo .GE. mig(1),iiglo .LE. mig(1)+nlci-1 ,ijglo .GE. mjg(1),ijglo .LE. mjg(1)+nlcj-1   
     635              !WRITE(*,*)"diadct readsec: log: ",narea, iiglo .GE. mig(1),iiglo .LE. mig(jpi),ijglo .GE. mjg(1),ijglo .LE. mjg(jpj) 
    438636              !verify if the point is on the local domain:(1,nlei)*(1,nlej) 
    439637              IF( iiloc .GE. 1 .AND. iiloc .LE. nlei .AND. & 
    440638                  ijloc .GE. 1 .AND. ijloc .LE. nlej       )THEN 
     639              !IF( iiglo .GE. mig(1) .AND. iiglo .LE. mig(1)+nlci-1 .AND. & 
     640              !    ijglo .GE. mjg(1) .AND. ijglo .LE. mjg(1)+nlcj-1       )THEN 
     641                   
     642                   
     643              !IF( iiglo .GE. mig(1) .AND. iiglo .LE. mig(jpi) .AND. & 
     644              !    ijglo .GE. mjg(1) .AND. ijglo .LE. mjg(jpj)       )THEN 
     645                  WRITE(*,*)"diadct readsec: assigned proc!",narea,nproc,jpt 
     646                   
     647                   
    441648                 iptloc = iptloc + 1                                                 ! count local points 
    442649                 secs(jsec)%listPoint(iptloc) = POINT_SECTION(mi0(iiglo),mj0(ijglo)) ! store local coordinates 
    … …  
    459666           ENDIF 
    460667 
    461               IF(jsec==nn_secdebug .AND. secs(jsec)%nb_point .NE. 0)THEN 
     668           IF(jsec==nn_secdebug .AND. secs(jsec)%nb_point .NE. 0)THEN 
    462669              DO jpt = 1,iptloc 
    463670                 iiglo = secs(jsec)%listPoint(jpt)%I + jpizoom - 1 + nimpp - 1 
    464671                 ijglo = secs(jsec)%listPoint(jpt)%J + jpjzoom - 1 + njmpp - 1 
    465672              ENDDO 
    466               ENDIF 
     673           ENDIF 
    467674 
    468675           !remove redundant points between processors 
    … …  
    501708 
    502709     ENDDO !end of the loop on jsec 
     710      
     711     !IF (lwp) write(numout,*) 'dct end of readsec loop, after exit' 
    503712  
    504713     nb_sec = jsec-1   !number of section read in the file 
    505714 
    506715     CALL wrk_dealloc( nb_point_max, directemp ) 
     716      
     717     !IF (lwp) write(numout,*) 'dct after dealloc' 
     718      
    507719     ! 
    508720  END SUBROUTINE readsec 
    … …  
    586798     CALL wrk_dealloc(    nb_point_max, idirec ) 
    587799     CALL wrk_dealloc( 2, nb_point_max, icoord ) 
     800 
    588801  END SUBROUTINE removepoints 
    589802 
    … …  
    592805     !!                     ***  ROUTINE transport  *** 
    593806     !! 
    594      !!  Purpose ::  Compute the transport for each point in a section  
     807     !!  Purpose ::  Compute the transport for each point in a section 
     808     !! 
     809     !!  Method  ::  Loop over each segment, and each vertical level and add the transport 
     810     !!              Be aware :           
     811     !!              One section is a sum of segments 
     812     !!              One segment is defined by 2 consecutive points in sec%listPoint 
     813     !!              All points of sec%listPoint are positioned on the F-point of the cell 
    595814     !!  
    596      !!  Method  ::  Loop over each segment, and each vertical level and add the transport  
    597      !!              Be aware :            
    598      !!              One section is a sum of segments  
    599      !!              One segment is defined by 2 consecutive points in sec%listPoint  
    600      !!              All points of sec%listPoint are positioned on the F-point of the cell  
    601      !!  
    602      !!              There are two loops:                   
    603      !!              loop on the segment between 2 nodes  
    604      !!              loop on the level jk !! 
    605      !!  
    606      !!  Output  ::  Arrays containing the volume,density,heat,salt transports for each i 
    607      !!              point in a section, summed over each nn_dct.  
     815     !!              There are two loops:                  
     816     !!              loop on the segment between 2 nodes 
     817     !!              loop on the level jk 
     818     !! 
     819     !! ** Output: Arrays containing the volume,density,salinity,temperature etc 
     820     !!            transports for each point in a section, summed over each nn_dct. 
    608821     !! 
    609822     !!------------------------------------------------------------------------------------------- 
    … …  
    614827     
    615828     !! * Local variables 
    616      INTEGER             :: jk, jseg, jclass,jl,                 &!loop on level/segment/classes/ice categories 
    617                             isgnu, isgnv                          !  
    618      REAL(wp)            :: zumid, zvmid,                        &!U/V velocity on a cell segment  
    619                             zumid_ice, zvmid_ice,                &!U/V ice velocity  
    620                             zTnorm                                !transport of velocity through one cell's sides  
    621      REAL(wp)            :: ztn, zsn, zrhoi, zrhop, zsshn, zfsdep !temperature/salinity/potential density/ssh/depth at u/v point 
     829     INTEGER             :: jk,jseg,jclass,    &!loop on level/segment/classes  
     830                            isgnu  , isgnv      ! 
     831     REAL(wp):: zumid        , zvmid        ,  &!U/V velocity on a cell segment 
     832                zumid_ice    , zvmid_ice    ,  &!U/V ice velocity 
     833                zTnorm                          !transport of velocity through one cell's sides 
     834     REAL(wp):: ztn, zsn, zrhoi, zrhop, zsshn, zfsdep ! temperature/salinity/ssh/potential density /depth at u/v point 
    622835 
    623836     TYPE(POINT_SECTION) :: k 
    624837     !!-------------------------------------------------------- 
    625838 
    626      IF( ld_debug )WRITE(numout,*)'      Compute transport' 
     839     IF( ld_debug )WRITE(numout,*)'      Compute transport (jsec,sec%nb_point,sec%slopeSection) : ', jsec,sec%nb_point,sec%slopeSection 
    627840 
    628841     !---------------------------! 
    … …  
    701914           END SELECT 
    702915 
    703            !---------------------------|  
    704            !     LOOP ON THE LEVEL     |  
    705            !---------------------------|  
    706            !Sum of the transport on the vertical   
    707            DO jk=1,mbathy(k%I,k%J)  
     916           !---------------------------| 
     917           !     LOOP ON THE LEVEL     | 
     918           !---------------------------| 
     919           !Sum of the transport on the vertical  
     920           DO jk=1,mbathy(k%I,k%J) 
    708921  
    709922              ! compute temperature, salinity, insitu & potential density, ssh and depth at U/V point  
    710               SELECT CASE( sec%direction(jseg) )  
    711               CASE(0,1)  
    712                  ztn   = interp(k%I,k%J,jk,'V',tsn(:,:,:,jp_tem) )  
    713                  zsn   = interp(k%I,k%J,jk,'V',tsn(:,:,:,jp_sal) )  
    714                  zrhop = interp(k%I,k%J,jk,'V',rhop)  
    715                  zrhoi = interp(k%I,k%J,jk,'V',rhd*rau0+rau0)  
    716                  zsshn =  0.5*( sshn(k%I,k%J)    + sshn(k%I,k%J+1)    ) * vmask(k%I,k%J,1)  
    717               CASE(2,3)  
    718                  ztn   = interp(k%I,k%J,jk,'U',tsn(:,:,:,jp_tem) )  
    719                  zsn   = interp(k%I,k%J,jk,'U',tsn(:,:,:,jp_sal) )  
    720                  zrhop = interp(k%I,k%J,jk,'U',rhop)  
    721                  zrhoi = interp(k%I,k%J,jk,'U',rhd*rau0+rau0)  
    722                  zsshn =  0.5*( sshn(k%I,k%J)    + sshn(k%I+1,k%J)    ) * umask(k%I,k%J,1)   
    723               END SELECT  
    724   
    725               zfsdep= fsdept(k%I,k%J,jk)  
    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               !zTnorm=transport through one cell;  
    738               !velocity* cell's length * cell's thickness  
    739               zTnorm=zumid*e2u(k%I,k%J)*  fse3u(k%I,k%J,jk)+     &  
    740                      zvmid*e1v(k%I,k%J)*  fse3v(k%I,k%J,jk)  
     923              SELECT CASE( sec%direction(jseg) ) 
     924              CASE(0,1) 
     925                 ztn   = interp(k%I,k%J,jk,'V',tsn(:,:,:,jp_tem) ) 
     926                 zsn   = interp(k%I,k%J,jk,'V',tsn(:,:,:,jp_sal) ) 
     927                 zrhop = interp(k%I,k%J,jk,'V',rhop) 
     928                 zrhoi = interp(k%I,k%J,jk,'V',rhd*rau0+rau0) 
     929                 zsshn =  0.5*( sshn(k%I,k%J)    + sshn(k%I,k%J+1)    ) * vmask(k%I,k%J,1) 
     930              CASE(2,3) 
     931                 ztn   = interp(k%I,k%J,jk,'U',tsn(:,:,:,jp_tem) ) 
     932                 zsn   = interp(k%I,k%J,jk,'U',tsn(:,:,:,jp_sal) ) 
     933                 zrhop = interp(k%I,k%J,jk,'U',rhop) 
     934                 zrhoi = interp(k%I,k%J,jk,'U',rhd*rau0+rau0) 
     935                 zsshn =  0.5*( sshn(k%I,k%J)    + sshn(k%I+1,k%J)    ) * umask(k%I,k%J,1)  
     936              END SELECT 
     937 
     938              zfsdep= fsdept(k%I,k%J,jk) 
     939  
     940              !compute velocity with the correct direction 
     941              SELECT CASE( sec%direction(jseg) ) 
     942              CASE(0,1)   
     943                 zumid=0. 
     944                 zvmid=isgnv*vn(k%I,k%J,jk)*vmask(k%I,k%J,jk) 
     945              CASE(2,3) 
     946                 zumid=isgnu*un(k%I,k%J,jk)*umask(k%I,k%J,jk) 
     947                 zvmid=0. 
     948              END SELECT 
     949 
     950              !zTnorm=transport through one cell; 
     951              !velocity* cell's length * cell's thickness 
     952              zTnorm=zumid*e2u(k%I,k%J)*  fse3u(k%I,k%J,jk)+     & 
     953                     zvmid*e1v(k%I,k%J)*  fse3v(k%I,k%J,jk) 
    741954 
    742955#if ! defined key_vvl 
    743               !add transport due to free surface  
    744               IF( jk==1 )THEN  
    745                  zTnorm = zTnorm + zumid* e2u(k%I,k%J) * zsshn * umask(k%I,k%J,jk) + &  
    746                                    zvmid* e1v(k%I,k%J) * zsshn * vmask(k%I,k%J,jk)  
    747               ENDIF  
     956              !add transport due to free surface 
     957              IF( jk==1 )THEN 
     958                 zTnorm = zTnorm + zumid* e2u(k%I,k%J) * zsshn * umask(k%I,k%J,jk) + & 
     959                                   zvmid* e1v(k%I,k%J) * zsshn * vmask(k%I,k%J,jk) 
     960              ENDIF 
    748961#endif 
    749               !COMPUTE TRANSPORT   
    750   
    751               transports_3d(1,jsec,jseg,jk) = transports_3d(1,jsec,jseg,jk) + zTnorm  
    752    
    753               IF ( sec%llstrpond ) THEN  
    754                  transports_3d(2,jsec,jseg,jk) = transports_3d(2,jsec,jseg,jk)  + zTnorm * ztn * zrhop * rcp 
    755                  transports_3d(3,jsec,jseg,jk) = transports_3d(3,jsec,jseg,jk)  + zTnorm * zsn * zrhop * 0.001 
     962              !COMPUTE TRANSPORT  
     963 
     964              transports_3d(1,jsec,jseg,jk) = transports_3d(1,jsec,jseg,jk) + zTnorm 
     965  
     966              IF ( sec%llstrpond ) THEN 
     967                 transports_3d(2,jsec,jseg,jk) = transports_3d(2,jsec,jseg,jk)  + zTnorm * zrhoi 
     968                 transports_3d(3,jsec,jseg,jk) = transports_3d(3,jsec,jseg,jk)  + zTnorm * zrhop 
     969                 transports_3d(4,jsec,jseg,jk) = transports_3d(4,jsec,jseg,jk)  + zTnorm * 4.e+3_wp * (ztn+273.15) * 1026._wp 
     970                 transports_3d(5,jsec,jseg,jk) = transports_3d(5,jsec,jseg,jk)  + zTnorm * 0.001 * zsn * 1026._wp 
    756971              ENDIF 
    757     
     972 
    758973           ENDDO !end of loop on the level 
    759974 
    … …  
    779994    
    780995              zTnorm=zumid_ice*e2u(k%I,k%J)+zvmid_ice*e1v(k%I,k%J) 
    781  
    782 #if defined key_lim2    
    783               transports_2d(1,jsec,jseg) = transports_2d(1,jsec,jseg) + (zTnorm)*   &  
    784                                    (1.0 - frld(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J))  &  
    785                                   *(hsnif(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J) +  &  
    786                                     hicif(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J)) 
    787               transports_2d(2,jsec,jseg) = transports_2d(2,jsec,jseg) + (zTnorm)*   &  
    788                                     (1.0 -  frld(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J)) 
    789 #endif 
    790 #if defined key_lim3 
    791               DO jl=1,jpl 
    792                  transports_2d(1,jsec,jseg) = transports_2d(1,jsec,jseg) + (zTnorm)*     & 
    793                                    a_i(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J,jl) * & 
    794                                   ( ht_i(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J,jl) +  & 
    795                                     ht_s(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J,jl) ) 
    796                                     
    797                  transports_2d(2,jsec,jseg) = transports_2d(2,jsec,jseg) + (zTnorm)*   & 
    798                                    a_i(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J,jl) 
    799               ENDDO 
    800 #endif 
     996    
     997              transports_2d(1,jsec,jseg) = transports_2d(1,jsec,jseg) + (zTnorm)*   & 
     998                                   (1.0 - frld(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J))  & 
     999                                  *(hsnif(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J) +  & 
     1000                                    hicif(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J)) & 
     1001                                   +zice_vol_pos 
     1002              transports_2d(2,jsec,jseg) = transports_2d(2,jsec,jseg) + (zTnorm)*   & 
     1003                                    (1.0 -  frld(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J))  & 
     1004                                   +zice_surf_pos 
    8011005    
    8021006           ENDIF !end of ice case 
    … …  
    8051009        ENDDO !end of loop on the segment 
    8061010 
    807      ENDIF !end of sec%nb_point =0 case 
     1011     ENDIF   !end of sec%nb_point =0 case 
    8081012     ! 
    8091013  END SUBROUTINE transport 
    810    
    811   SUBROUTINE dia_dct_sum(sec,jsec)  
     1014 
     1015  SUBROUTINE transport_h(sec,ld_debug,jsec) 
     1016     !!------------------------------------------------------------------------------------------- 
     1017     !!                     ***  ROUTINE hourly_transport  *** 
     1018     !! 
     1019     !!              Exactly as routine transport but for data summed at 
     1020     !!              each time step and averaged each hour 
     1021     !! 
     1022     !!  Purpose ::  Compute the transport for each point in a section 
     1023     !! 
     1024     !!  Method  ::  Loop over each segment, and each vertical level and add the transport 
     1025     !!              Be aware :           
     1026     !!              One section is a sum of segments 
     1027     !!              One segment is defined by 2 consecutive points in sec%listPoint 
     1028     !!              All points of sec%listPoint are positioned on the F-point of the cell 
     1029     !!  
     1030     !!              There are two loops:                  
     1031     !!              loop on the segment between 2 nodes 
     1032     !!              loop on the level jk 
     1033     !! 
     1034     !! ** Output: Arrays containing the volume,density,salinity,temperature etc 
     1035     !!            transports for each point in a section, summed over each nn_dct. 
     1036     !! 
     1037     !!------------------------------------------------------------------------------------------- 
     1038     !! * Arguments 
     1039     TYPE(SECTION),INTENT(INOUT) :: sec 
     1040     LOGICAL      ,INTENT(IN)    :: ld_debug 
     1041     INTEGER      ,INTENT(IN)    :: jsec        ! numeric identifier of section 
     1042     
     1043     !! * Local variables 
     1044     INTEGER             :: jk,jseg,jclass,    &!loop on level/segment/classes  
     1045                            isgnu  , isgnv      ! 
     1046     REAL(wp):: zumid        , zvmid        ,  &!U/V velocity on a cell segment 
     1047                zumid_ice    , zvmid_ice    ,  &!U/V ice velocity 
     1048                zTnorm                          !transport of velocity through one cell's sides 
     1049     REAL(wp):: ztn, zsn, zrhoi, zrhop, zsshn, zfsdep ! temperature/salinity/ssh/potential density /depth at u/v point 
     1050 
     1051     TYPE(POINT_SECTION) :: k 
     1052     !!-------------------------------------------------------- 
     1053 
     1054     IF( ld_debug )WRITE(numout,*)'      Compute transport' 
     1055 
     1056     !---------------------------! 
     1057     !  COMPUTE TRANSPORT        ! 
     1058     !---------------------------! 
     1059     IF(sec%nb_point .NE. 0)THEN    
     1060 
     1061        !---------------------------------------------------------------------------------------------------- 
     1062        !Compute sign for velocities: 
     1063        ! 
     1064        !convention: 
     1065        !   non horizontal section: direction + is toward left hand of section 
     1066        !       horizontal section: direction + is toward north of section 
     1067        ! 
     1068        ! 
     1069        !       slopeSection < 0     slopeSection > 0       slopeSection=inf            slopeSection=0 
     1070        !       ----------------      -----------------     ---------------             -------------- 
     1071        ! 
     1072        !   isgnv=1         direction +       
     1073        !  ______         _____             ______                                                    
     1074        !        |           //|            |                  |                         direction +    
     1075        !        | isgnu=1  // |            |isgnu=1           |isgnu=1                     /|\ 
     1076        !        |_______  //         ______|    \\            | ---\                        | 
     1077        !               |             | isgnv=-1  \\ |         | ---/ direction +       ____________ 
     1078        !               |             |          __\\|         |                     
     1079        !               |             |     direction +        |                      isgnv=1                                  
     1080        !                                                       
     1081        !---------------------------------------------------------------------------------------------------- 
     1082        isgnu = 1 
     1083        IF( sec%slopeSection .GT. 0 ) THEN  ; isgnv = -1  
     1084        ELSE                                ; isgnv =  1 
     1085        ENDIF 
     1086 
     1087        IF( ld_debug )write(numout,*)"isgnu isgnv ",isgnu,isgnv 
     1088 
     1089        !--------------------------------------! 
     1090        ! LOOP ON THE SEGMENT BETWEEN 2 NODES  ! 
     1091        !--------------------------------------! 
     1092        DO jseg=1,MAX(sec%nb_point-1,0) 
     1093               
     1094           !------------------------------------------------------------------------------------------- 
     1095           ! Select the appropriate coordinate for computing the velocity of the segment 
     1096           ! 
     1097           !                      CASE(0)                                    Case (2) 
     1098           !                      -------                                    -------- 
     1099           !  listPoint(jseg)                 listPoint(jseg+1)       listPoint(jseg)  F(i,j)       
     1100           !      F(i,j)----------V(i+1,j)-------F(i+1,j)                               | 
     1101           !                                                                            | 
     1102           !                                                                            | 
     1103           !                                                                            | 
     1104           !                      Case (3)                                            U(i,j) 
     1105           !                      --------                                              | 
     1106           !                                                                            | 
     1107           !  listPoint(jseg+1) F(i,j+1)                                                | 
     1108           !                        |                                                   | 
     1109           !                        |                                                   | 
     1110           !                        |                                 listPoint(jseg+1) F(i,j-1) 
     1111           !                        |                                             
     1112           !                        |                                             
     1113           !                     U(i,j+1)                                             
     1114           !                        |                                       Case(1)      
     1115           !                        |                                       ------       
     1116           !                        |                                             
     1117           !                        |                 listPoint(jseg+1)             listPoint(jseg)                            
     1118           !                        |                 F(i-1,j)-----------V(i,j) -------f(jseg)                            
     1119           ! listPoint(jseg)     F(i,j) 
     1120           !  
     1121           !------------------------------------------------------------------------------------------- 
     1122 
     1123           SELECT CASE( sec%direction(jseg) ) 
     1124           CASE(0)  ;   k = sec%listPoint(jseg) 
     1125           CASE(1)  ;   k = POINT_SECTION(sec%listPoint(jseg)%I+1,sec%listPoint(jseg)%J) 
     1126           CASE(2)  ;   k = sec%listPoint(jseg) 
     1127           CASE(3)  ;   k = POINT_SECTION(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J+1) 
     1128           END SELECT 
     1129 
     1130           !---------------------------| 
     1131           !     LOOP ON THE LEVEL     | 
     1132           !---------------------------| 
     1133           !Sum of the transport on the vertical  
     1134           DO jk=1,mbathy(k%I,k%J) 
     1135 
     1136              ! compute temparature, salinity, insitu & potential density, ssh and depth at U/V point 
     1137              SELECT CASE( sec%direction(jseg) ) 
     1138              CASE(0,1) 
     1139                 ztn   = interp(k%I,k%J,jk,'V',tsn(:,:,:,jp_tem) ) 
     1140                 zsn   = interp(k%I,k%J,jk,'V',tsn(:,:,:,jp_sal) ) 
     1141                 zrhop = interp(k%I,k%J,jk,'V',rhop) 
     1142                 zrhoi = interp(k%I,k%J,jk,'V',rhd*rau0+rau0) 
     1143                 zsshn =  0.5*( sshn(k%I,k%J)    + sshn(k%I,k%J+1)    ) * vmask(k%I,k%J,1) 
     1144              CASE(2,3) 
     1145                 ztn   = interp(k%I,k%J,jk,'U',tsn(:,:,:,jp_tem) ) 
     1146                 zsn   = interp(k%I,k%J,jk,'U',tsn(:,:,:,jp_sal) ) 
     1147                 zrhop = interp(k%I,k%J,jk,'U',rhop) 
     1148                 zrhoi = interp(k%I,k%J,jk,'U',rhd*rau0+rau0) 
     1149                 zsshn =  0.5*( sshn(k%I,k%J)    + sshn(k%I+1,k%J)    ) * umask(k%I,k%J,1)  
     1150              END SELECT 
     1151 
     1152              !JT zfsdep= gdept(k%I,k%J,jk) 
     1153              !JT zfsdep= gdept_0(k%I,k%J,jk) 
     1154              zfsdep= fsdept(k%I,k%J,jk) 
     1155  
     1156              !compute velocity with the correct direction 
     1157              SELECT CASE( sec%direction(jseg) ) 
     1158              CASE(0,1)   
     1159                 zumid=0. 
     1160                 zvmid=isgnv*vn(k%I,k%J,jk)*vmask(k%I,k%J,jk) 
     1161              CASE(2,3) 
     1162                 zumid=isgnu*un(k%I,k%J,jk)*umask(k%I,k%J,jk) 
     1163                 zvmid=0. 
     1164              END SELECT 
     1165 
     1166              !zTnorm=transport through one cell; 
     1167              !velocity* cell's length * cell's thickness 
     1168              zTnorm=zumid*e2u(k%I,k%J)*  fse3u(k%I,k%J,jk)+     & 
     1169                     zvmid*e1v(k%I,k%J)*  fse3v(k%I,k%J,jk) 
     1170 
     1171#if ! defined key_vvl 
     1172              !add transport due to free surface 
     1173              IF( jk==1 )THEN 
     1174                 zTnorm = zTnorm + zumid* e2u(k%I,k%J) * zsshn * umask(k%I,k%J,jk) + & 
     1175                                   zvmid* e1v(k%I,k%J) * zsshn * vmask(k%I,k%J,jk) 
     1176              ENDIF 
     1177#endif 
     1178              !COMPUTE TRANSPORT  
     1179 
     1180              transports_3d_h(1,jsec,jseg,jk) = transports_3d_h(1,jsec,jseg,jk) + zTnorm 
     1181  
     1182              IF ( sec%llstrpond ) THEN 
     1183                 transports_3d_h(2,jsec,jseg,jk) = transports_3d_h(2,jsec,jseg,jk)  + zTnorm * zrhoi 
     1184                 transports_3d_h(3,jsec,jseg,jk) = transports_3d_h(3,jsec,jseg,jk)  + zTnorm * zrhop 
     1185                 transports_3d_h(4,jsec,jseg,jk) = transports_3d_h(4,jsec,jseg,jk)  + zTnorm * 4.e+3_wp * (ztn+273.15) * 1026._wp 
     1186                 transports_3d_h(5,jsec,jseg,jk) = transports_3d_h(5,jsec,jseg,jk)  + zTnorm * 0.001 * zsn * 1026._wp 
     1187              ENDIF 
     1188 
     1189           ENDDO !end of loop on the level 
     1190 
     1191#if defined key_lim2 || defined key_lim3 
     1192 
     1193           !ICE CASE     
     1194           !------------ 
     1195           IF( sec%ll_ice_section )THEN 
     1196              SELECT CASE (sec%direction(jseg)) 
     1197              CASE(0) 
     1198                 zumid_ice = 0 
     1199                 zvmid_ice =  isgnv*0.5*(v_ice(k%I,k%J+1)+v_ice(k%I+1,k%J+1)) 
     1200              CASE(1) 
     1201                 zumid_ice = 0 
     1202                 zvmid_ice =  isgnv*0.5*(v_ice(k%I,k%J+1)+v_ice(k%I+1,k%J+1)) 
     1203              CASE(2) 
     1204                 zvmid_ice = 0 
     1205                 zumid_ice =  isgnu*0.5*(u_ice(k%I+1,k%J)+u_ice(k%I+1,k%J+1)) 
     1206              CASE(3) 
     1207                 zvmid_ice = 0 
     1208                 zumid_ice =  isgnu*0.5*(u_ice(k%I+1,k%J)+u_ice(k%I+1,k%J+1)) 
     1209              END SELECT 
     1210    
     1211              zTnorm=zumid_ice*e2u(k%I,k%J)+zvmid_ice*e1v(k%I,k%J) 
     1212    
     1213              transports_2d_h(1,jsec,jseg) = transports_2d_h(1,jsec,jseg) + (zTnorm)*   & 
     1214                                   (1.0 - frld(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J))  & 
     1215                                  *(hsnif(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J) +  & 
     1216                                    hicif(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J)) & 
     1217                                   +zice_vol_pos 
     1218              transports_2d_h(2,jsec,jseg) = transports_2d_h(2,jsec,jseg) + (zTnorm)*   & 
     1219                                    (1.0 -  frld(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J))  & 
     1220                                   +zice_surf_pos 
     1221    
     1222           ENDIF !end of ice case 
     1223#endif 
     1224  
     1225        ENDDO !end of loop on the segment 
     1226 
     1227     ENDIF   !end of sec%nb_point =0 case 
     1228     ! 
     1229  END SUBROUTINE transport_h 
     1230  
     1231  SUBROUTINE dia_dct_sum(sec,jsec) 
     1232     !!------------------------------------------------------------- 
     1233     !! Purpose: Average the transport over nn_dctwri time steps  
     1234     !! and sum over the density/salinity/temperature/depth classes 
     1235     !! 
     1236     !! Method:  
     1237     !!           Sum over relevant grid cells to obtain values 
     1238     !!           for each 
     1239     !!              There are several loops:                  
     1240     !!              loop on the segment between 2 nodes 
     1241     !!              loop on the level jk 
     1242     !!              loop on the density/temperature/salinity/level classes 
     1243     !!              test on the density/temperature/salinity/level 
     1244     !! 
     1245     !!  ** Method  :Transport through a given section is equal to the sum of transports 
     1246     !!              computed on each proc. 
     1247     !!              On each proc,transport is equal to the sum of transport computed through 
     1248     !!              segments linking each point of sec%listPoint  with the next one.    
     1249     !! 
     1250     !!------------------------------------------------------------- 
     1251     !! * arguments 
     1252     TYPE(SECTION),INTENT(INOUT) :: sec 
     1253     INTEGER      ,INTENT(IN)    :: jsec        ! numeric identifier of section 
     1254 
     1255     TYPE(POINT_SECTION) :: k 
     1256     INTEGER  :: jk,jseg,jclass                        !loop on level/segment/classes  
     1257     REAL(wp) :: ztn, zsn, zrhoi, zrhop, zsshn, zfsdep ! temperature/salinity/ssh/potential density /depth at u/v point 
     1258     !!------------------------------------------------------------- 
     1259 
     1260     !! Sum the relevant segments to obtain values for each class 
     1261     IF(sec%nb_point .NE. 0)THEN    
     1262 
     1263        !--------------------------------------! 
     1264        ! LOOP ON THE SEGMENT BETWEEN 2 NODES  ! 
     1265        !--------------------------------------! 
     1266        DO jseg=1,MAX(sec%nb_point-1,0) 
     1267            
     1268           !------------------------------------------------------------------------------------------- 
     1269           ! Select the appropriate coordinate for computing the velocity of the segment 
     1270           ! 
     1271           !                      CASE(0)                                    Case (2) 
     1272           !                      -------                                    -------- 
     1273           !  listPoint(jseg)                 listPoint(jseg+1)       listPoint(jseg)  F(i,j)       
     1274           !      F(i,j)----------V(i+1,j)-------F(i+1,j)                               | 
     1275           !                                                                            | 
     1276           !                                                                            | 
     1277           !                                                                            | 
     1278           !                      Case (3)                                            U(i,j) 
     1279           !                      --------                                              | 
     1280           !                                                                            | 
     1281           !  listPoint(jseg+1) F(i,j+1)                                                | 
     1282           !                        |                                                   | 
     1283           !                        |                                                   | 
     1284           !                        |                                 listPoint(jseg+1) F(i,j-1) 
     1285           !                        |                                             
     1286           !                        |                                             
     1287           !                     U(i,j+1)                                             
     1288           !                        |                                       Case(1)      
     1289           !                        |                                       ------       
     1290           !                        |                                             
     1291           !                        |                 listPoint(jseg+1)             listPoint(jseg)                            
     1292           !                        |                 F(i-1,j)-----------V(i,j) -------f(jseg)                            
     1293           ! listPoint(jseg)     F(i,j) 
     1294           !  
     1295           !------------------------------------------------------------------------------------------- 
     1296 
     1297           SELECT CASE( sec%direction(jseg) ) 
     1298           CASE(0)  ;   k = sec%listPoint(jseg) 
     1299           CASE(1)  ;   k = POINT_SECTION(sec%listPoint(jseg)%I+1,sec%listPoint(jseg)%J) 
     1300           CASE(2)  ;   k = sec%listPoint(jseg) 
     1301           CASE(3)  ;   k = POINT_SECTION(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J+1) 
     1302           END SELECT 
     1303 
     1304           !---------------------------| 
     1305           !     LOOP ON THE LEVEL     | 
     1306           !---------------------------| 
     1307           !Sum of the transport on the vertical  
     1308           DO jk=1,mbathy(k%I,k%J) 
     1309 
     1310              ! compute temperature, salinity, insitu & potential density, ssh and depth at U/V point 
     1311              SELECT CASE( sec%direction(jseg) ) 
     1312              CASE(0,1) 
     1313                 ztn   = interp(k%I,k%J,jk,'V',tsn(:,:,:,jp_tem) ) 
     1314                 zsn   = interp(k%I,k%J,jk,'V',tsn(:,:,:,jp_sal) ) 
     1315                 zrhop = interp(k%I,k%J,jk,'V',rhop) 
     1316                 zrhoi = interp(k%I,k%J,jk,'V',rhd*rau0+rau0) 
     1317                 zsshn =  0.5*( sshn(k%I,k%J)    + sshn(k%I,k%J+1)    ) * vmask(k%I,k%J,1) 
     1318              CASE(2,3) 
     1319                 ztn   = interp(k%I,k%J,jk,'U',tsn(:,:,:,jp_tem) ) 
     1320                 zsn   = interp(k%I,k%J,jk,'U',tsn(:,:,:,jp_sal) ) 
     1321                 zrhop = interp(k%I,k%J,jk,'U',rhop) 
     1322                 zrhoi = interp(k%I,k%J,jk,'U',rhd*rau0+rau0) 
     1323                 zsshn =  0.5*( sshn(k%I,k%J)    + sshn(k%I+1,k%J)    ) * umask(k%I,k%J,1)  
     1324              END SELECT 
     1325 
     1326              !JT zfsdep= gdept(k%I,k%J,jk) 
     1327              !JT zfsdep= gdept_0(k%I,k%J,jk) 
     1328              zfsdep= fsdept(k%I,k%J,jk)  
     1329  
     1330              !------------------------------- 
     1331              !  LOOP ON THE DENSITY CLASSES | 
     1332              !------------------------------- 
     1333              !The computation is made for each density/temperature/salinity/depth class 
     1334              DO jclass=1,MAX(1,sec%nb_class-1) 
     1335  
     1336                 !----------------------------------------------! 
     1337                 !TEST ON THE DENSITY/SALINITY/TEMPERATURE/LEVEL!  
     1338                 !----------------------------------------------! 
     1339  
     1340                 IF ( (                                                    & 
     1341                    ((( zrhop .GE. (sec%zsigp(jclass)+1000.  )) .AND.      & 
     1342                    (   zrhop .LE. (sec%zsigp(jclass+1)+1000. ))) .OR.     & 
     1343                    ( sec%zsigp(jclass) .EQ. 99.)) .AND.                   & 
     1344 
     1345                    ((( zrhoi .GE. (sec%zsigi(jclass) + 1000.  )) .AND.    & 
     1346                    (   zrhoi .LE. (sec%zsigi(jclass+1)+1000. ))) .OR.     & 
     1347                    ( sec%zsigi(jclass) .EQ. 99.)) .AND.                   & 
     1348 
     1349                    ((( zsn .GT. sec%zsal(jclass)) .AND.                   & 
     1350                    (   zsn .LE. sec%zsal(jclass+1))) .OR.                 & 
     1351                    ( sec%zsal(jclass) .EQ. 99.)) .AND.                    & 
     1352 
     1353                    ((( ztn .GE. sec%ztem(jclass)) .AND.                   & 
     1354                    (   ztn .LE. sec%ztem(jclass+1))) .OR.                 & 
     1355                    ( sec%ztem(jclass) .EQ.99.)) .AND.                     & 
     1356 
     1357                    ((( zfsdep .GE. sec%zlay(jclass)) .AND.                & 
     1358                    (   zfsdep .LE. sec%zlay(jclass+1))) .OR.              & 
     1359                    ( sec%zlay(jclass) .EQ. 99. ))                         & 
     1360                                                                   ))   THEN 
     1361 
     1362                    !SUM THE TRANSPORTS FOR EACH CLASSES FOR THE POSITIVE AND NEGATIVE DIRECTIONS 
     1363                    !---------------------------------------------------------------------------- 
     1364                    IF (transports_3d(1,jsec,jseg,jk) .GE. 0.0) THEN  
     1365                       sec%transport(1,jclass) = sec%transport(1,jclass)+transports_3d(1,jsec,jseg,jk) 
     1366                    ELSE 
     1367                       sec%transport(2,jclass) = sec%transport(2,jclass)+transports_3d(1,jsec,jseg,jk) 
     1368                    ENDIF 
     1369                    IF( sec%llstrpond )THEN 
     1370 
     1371                       IF( transports_3d(1,jsec,jseg,jk) .NE. 0._wp ) THEN 
     1372 
     1373                          IF (transports_3d(2,jsec,jseg,jk)/transports_3d(1,jsec,jseg,jk) .GE. 0.0 ) THEN 
     1374                             sec%transport(3,jclass) = sec%transport(3,jclass)+transports_3d(2,jsec,jseg,jk)/transports_3d(1,jsec,jseg,jk) 
     1375                          ELSE 
     1376                             sec%transport(4,jclass) = sec%transport(4,jclass)+transports_3d(2,jsec,jseg,jk)/transports_3d(1,jsec,jseg,jk) 
     1377                          ENDIF 
     1378 
     1379                          IF ( transports_3d(3,jsec,jseg,jk)/transports_3d(1,jsec,jseg,jk) .GE. 0.0 ) THEN 
     1380                             sec%transport(5,jclass) = sec%transport(5,jclass)+transports_3d(3,jsec,jseg,jk)/transports_3d(1,jsec,jseg,jk) 
     1381                          ELSE 
     1382                             sec%transport(6,jclass) = sec%transport(6,jclass)+transports_3d(3,jsec,jseg,jk)/transports_3d(1,jsec,jseg,jk) 
     1383                          ENDIF 
     1384 
     1385                       ENDIF 
     1386 
     1387                       IF ( transports_3d(4,jsec,jseg,jk) .GE. 0.0 ) THEN 
     1388                          sec%transport(7,jclass) = sec%transport(7,jclass)+transports_3d(4,jsec,jseg,jk) 
     1389                       ELSE 
     1390                          sec%transport(8,jclass) = sec%transport(8,jclass)+transports_3d(4,jsec,jseg,jk) 
     1391                       ENDIF 
     1392 
     1393                       IF ( transports_3d(5,jsec,jseg,jk) .GE. 0.0 ) THEN 
     1394                          sec%transport( 9,jclass) = sec%transport( 9,jclass)+transports_3d(5,jsec,jseg,jk) 
     1395                       ELSE 
     1396                          sec%transport(10,jclass) = sec%transport(10,jclass)+transports_3d(5,jsec,jseg,jk) 
     1397                       ENDIF 
     1398 
     1399                    ELSE 
     1400                       sec%transport( 3,jclass) = 0._wp 
     1401                       sec%transport( 4,jclass) = 0._wp 
     1402                       sec%transport( 5,jclass) = 0._wp 
     1403                       sec%transport( 6,jclass) = 0._wp 
     1404                       sec%transport( 7,jclass) = 0._wp 
     1405                       sec%transport( 8,jclass) = 0._wp 
     1406                       sec%transport( 9,jclass) = 0._wp 
     1407                       sec%transport(10,jclass) = 0._wp 
     1408                    ENDIF 
     1409 
     1410                 ENDIF ! end of test if point is in class 
     1411    
     1412              ENDDO ! end of loop on the classes 
     1413 
     1414           ENDDO ! loop over jk 
     1415 
     1416#if defined key_lim2 || defined key_lim3 
     1417 
     1418           !ICE CASE     
     1419           IF( sec%ll_ice_section )THEN 
     1420 
     1421              IF ( transports_2d(1,jsec,jseg) .GE. 0.0 ) THEN 
     1422                 sec%transport(11,1) = sec%transport(11,1)+transports_2d(1,jsec,jseg) 
     1423              ELSE 
     1424                 sec%transport(12,1) = sec%transport(12,1)+transports_2d(1,jsec,jseg) 
     1425              ENDIF 
     1426 
     1427              IF ( transports_2d(3,jsec,jseg) .GE. 0.0 ) THEN 
     1428                 sec%transport(13,1) = sec%transport(13,1)+transports_2d(2,jsec,jseg) 
     1429              ELSE 
     1430                 sec%transport(14,1) = sec%transport(14,1)+transports_2d(2,jsec,jseg) 
     1431              ENDIF 
     1432 
     1433           ENDIF !end of ice case 
     1434#endif 
     1435  
     1436        ENDDO !end of loop on the segment 
     1437 
     1438     ELSE  !if sec%nb_point =0 
     1439        sec%transport(1:2,:)=0. 
     1440        IF (sec%llstrpond) sec%transport(3:10,:)=0. 
     1441        IF (sec%ll_ice_section) sec%transport( 11:14,:)=0. 
     1442     ENDIF !end of sec%nb_point =0 case 
     1443 
     1444  END SUBROUTINE dia_dct_sum 
     1445  
     1446  SUBROUTINE dia_dct_sum_h(sec,jsec) 
     1447     !!------------------------------------------------------------- 
     1448     !! Exactly as dia_dct_sum but for hourly files containing data summed at each time step 
     1449     !! 
     1450     !! Purpose: Average the transport over nn_dctwri time steps  
     1451     !! and sum over the density/salinity/temperature/depth classes 
     1452     !! 
     1453     !! Method:  
     1454     !!           Sum over relevant grid cells to obtain values 
     1455     !!           for each 
     1456     !!              There are several loops:                  
     1457     !!              loop on the segment between 2 nodes 
     1458     !!              loop on the level jk 
     1459     !!              loop on the density/temperature/salinity/level classes 
     1460     !!              test on the density/temperature/salinity/level 
     1461     !! 
     1462     !!  ** Method  :Transport through a given section is equal to the sum of transports 
     1463     !!              computed on each proc. 
     1464     !!              On each proc,transport is equal to the sum of transport computed through 
     1465     !!              segments linking each point of sec%listPoint  with the next one.    
     1466     !! 
     1467     !!------------------------------------------------------------- 
     1468     !! * arguments 
     1469     TYPE(SECTION),INTENT(INOUT) :: sec 
     1470     INTEGER      ,INTENT(IN)    :: jsec        ! numeric identifier of section 
     1471 
     1472     TYPE(POINT_SECTION) :: k 
     1473     INTEGER  :: jk,jseg,jclass                        !loop on level/segment/classes  
     1474     REAL(wp) :: ztn, zsn, zrhoi, zrhop, zsshn, zfsdep ! temperature/salinity/ssh/potential density /depth at u/v point 
     1475     !!------------------------------------------------------------- 
     1476 
     1477     !! Sum the relevant segments to obtain values for each class 
     1478     IF(sec%nb_point .NE. 0)THEN    
     1479 
     1480        !--------------------------------------! 
     1481        ! LOOP ON THE SEGMENT BETWEEN 2 NODES  ! 
     1482        !--------------------------------------! 
     1483        DO jseg=1,MAX(sec%nb_point-1,0) 
     1484            
     1485           !------------------------------------------------------------------------------------------- 
     1486           ! Select the appropriate coordinate for computing the velocity of the segment 
     1487           ! 
     1488           !                      CASE(0)                                    Case (2) 
     1489           !                      -------                                    -------- 
     1490           !  listPoint(jseg)                 listPoint(jseg+1)       listPoint(jseg)  F(i,j)       
     1491           !      F(i,j)----------V(i+1,j)-------F(i+1,j)                               | 
     1492           !                                                                            | 
     1493           !                                                                            | 
     1494           !                                                                            | 
     1495           !                      Case (3)                                            U(i,j) 
     1496           !                      --------                                              | 
     1497           !                                                                            | 
     1498           !  listPoint(jseg+1) F(i,j+1)                                                | 
     1499           !                        |                                                   | 
     1500           !                        |                                                   | 
     1501           !                        |                                 listPoint(jseg+1) F(i,j-1) 
     1502           !                        |                                             
     1503           !                        |                                             
     1504           !                     U(i,j+1)                                             
     1505           !                        |                                       Case(1)      
     1506           !                        |                                       ------       
     1507           !                        |                                             
     1508           !                        |                 listPoint(jseg+1)             listPoint(jseg)                            
     1509           !                        |                 F(i-1,j)-----------V(i,j) -------f(jseg)                            
     1510           ! listPoint(jseg)     F(i,j) 
     1511           !  
     1512           !------------------------------------------------------------------------------------------- 
     1513 
     1514           SELECT CASE( sec%direction(jseg) ) 
     1515           CASE(0)  ;   k = sec%listPoint(jseg) 
     1516           CASE(1)  ;   k = POINT_SECTION(sec%listPoint(jseg)%I+1,sec%listPoint(jseg)%J) 
     1517           CASE(2)  ;   k = sec%listPoint(jseg) 
     1518           CASE(3)  ;   k = POINT_SECTION(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J+1) 
     1519           END SELECT 
     1520 
     1521           !---------------------------| 
     1522           !     LOOP ON THE LEVEL     | 
     1523           !---------------------------| 
     1524           !Sum of the transport on the vertical  
     1525           DO jk=1,mbathy(k%I,k%J) 
     1526 
     1527              ! compute temperature, salinity, insitu & potential density, ssh and depth at U/V point 
     1528              SELECT CASE( sec%direction(jseg) ) 
     1529              CASE(0,1) 
     1530                 ztn   = interp(k%I,k%J,jk,'V',tsn(:,:,:,jp_tem) ) 
     1531                 zsn   = interp(k%I,k%J,jk,'V',tsn(:,:,:,jp_sal) ) 
     1532                 zrhop = interp(k%I,k%J,jk,'V',rhop) 
     1533                 zrhoi = interp(k%I,k%J,jk,'V',rhd*rau0+rau0) 
     1534                 zsshn =  0.5*( sshn(k%I,k%J)    + sshn(k%I,k%J+1)    ) * vmask(k%I,k%J,1) 
     1535              CASE(2,3) 
     1536                 ztn   = interp(k%I,k%J,jk,'U',tsn(:,:,:,jp_tem) ) 
     1537                 zsn   = interp(k%I,k%J,jk,'U',tsn(:,:,:,jp_sal) ) 
     1538                 zrhop = interp(k%I,k%J,jk,'U',rhop) 
     1539                 zrhoi = interp(k%I,k%J,jk,'U',rhd*rau0+rau0) 
     1540                 zsshn =  0.5*( sshn(k%I,k%J)    + sshn(k%I+1,k%J)    ) * umask(k%I,k%J,1)  
     1541              END SELECT 
     1542 
     1543              !JT zfsdep= gdept(k%I,k%J,jk) 
     1544              !JT zfsdep= gdept_0(k%I,k%J,jk) 
     1545              zfsdep= fsdept(k%I,k%J,jk) 
     1546  
     1547              !------------------------------- 
     1548              !  LOOP ON THE DENSITY CLASSES | 
     1549              !------------------------------- 
     1550              !The computation is made for each density/heat/salt/... class 
     1551              DO jclass=1,MAX(1,sec%nb_class-1) 
     1552 
     1553                 !----------------------------------------------! 
     1554                 !TEST ON THE DENSITY/SALINITY/TEMPERATURE/LEVEL!  
     1555                 !----------------------------------------------! 
     1556  
     1557                 IF ( (                                                    & 
     1558                    ((( zrhop .GE. (sec%zsigp(jclass)+1000.  )) .AND.      & 
     1559                    (   zrhop .LE. (sec%zsigp(jclass+1)+1000. ))) .OR.     & 
     1560                    ( sec%zsigp(jclass) .EQ. 99.)) .AND.                   & 
     1561 
     1562                    ((( zrhoi .GE. (sec%zsigi(jclass) + 1000.  )) .AND.    & 
     1563                    (   zrhoi .LE. (sec%zsigi(jclass+1)+1000. ))) .OR.     & 
     1564                    ( sec%zsigi(jclass) .EQ. 99.)) .AND.                   & 
     1565 
     1566                    ((( zsn .GT. sec%zsal(jclass)) .AND.                   & 
     1567                    (   zsn .LE. sec%zsal(jclass+1))) .OR.                 & 
     1568                    ( sec%zsal(jclass) .EQ. 99.)) .AND.                    & 
     1569 
     1570                    ((( ztn .GE. sec%ztem(jclass)) .AND.                   & 
     1571                    (   ztn .LE. sec%ztem(jclass+1))) .OR.                 & 
     1572                    ( sec%ztem(jclass) .EQ.99.)) .AND.                     & 
     1573 
     1574                    ((( zfsdep .GE. sec%zlay(jclass)) .AND.                & 
     1575                    (   zfsdep .LE. sec%zlay(jclass+1))) .OR.              & 
     1576                    ( sec%zlay(jclass) .EQ. 99. ))                         & 
     1577                                                                   ))   THEN 
     1578 
     1579                    !SUM THE TRANSPORTS FOR EACH CLASSES FOR THE POSITIVE AND NEGATIVE DIRECTIONS 
     1580                    !---------------------------------------------------------------------------- 
     1581                    IF (transports_3d_h(1,jsec,jseg,jk) .GE. 0.0) THEN  
     1582                       sec%transport_h(1,jclass) = sec%transport_h(1,jclass)+transports_3d_h(1,jsec,jseg,jk) 
     1583                    ELSE 
     1584                       sec%transport_h(2,jclass) = sec%transport_h(2,jclass)+transports_3d_h(1,jsec,jseg,jk) 
     1585                    ENDIF 
     1586                    IF( sec%llstrpond )THEN 
     1587 
     1588                       IF( transports_3d_h(1,jsec,jseg,jk) .NE. 0._wp ) THEN 
     1589 
     1590                          IF (transports_3d_h(2,jsec,jseg,jk)/transports_3d_h(1,jsec,jseg,jk) .GE. 0.0 ) THEN 
     1591                             sec%transport_h(3,jclass) = sec%transport_h(3,jclass)+transports_3d_h(2,jsec,jseg,jk)/transports_3d_h(1,jsec,jseg,jk) 
     1592                          ELSE 
     1593                             sec%transport_h(4,jclass) = sec%transport_h(4,jclass)+transports_3d_h(2,jsec,jseg,jk)/transports_3d_h(1,jsec,jseg,jk) 
     1594                          ENDIF 
     1595 
     1596                          IF ( transports_3d_h(3,jsec,jseg,jk)/transports_3d_h(1,jsec,jseg,jk) .GE. 0.0 ) THEN 
     1597                             sec%transport_h(5,jclass) = sec%transport_h(5,jclass)+transports_3d_h(3,jsec,jseg,jk)/transports_3d_h(1,jsec,jseg,jk) 
     1598                          ELSE 
     1599                             sec%transport_h(6,jclass) = sec%transport_h(6,jclass)+transports_3d_h(3,jsec,jseg,jk)/transports_3d_h(1,jsec,jseg,jk) 
     1600                          ENDIF 
     1601 
     1602                       ENDIF 
     1603 
     1604                       IF ( transports_3d_h(4,jsec,jseg,jk) .GE. 0.0 ) THEN 
     1605                          sec%transport_h(7,jclass) = sec%transport_h(7,jclass)+transports_3d_h(4,jsec,jseg,jk) 
     1606                       ELSE 
     1607                          sec%transport_h(8,jclass) = sec%transport_h(8,jclass)+transports_3d_h(4,jsec,jseg,jk) 
     1608                       ENDIF 
     1609 
     1610                       IF ( transports_3d_h(5,jsec,jseg,jk) .GE. 0.0 ) THEN 
     1611                          sec%transport_h( 9,jclass) = sec%transport_h( 9,jclass)+transports_3d_h(5,jsec,jseg,jk) 
     1612                       ELSE 
     1613                          sec%transport_h(10,jclass) = sec%transport_h(10,jclass)+transports_3d_h(5,jsec,jseg,jk) 
     1614                       ENDIF 
     1615 
     1616                    ELSE 
     1617                       sec%transport_h( 3,jclass) = 0._wp 
     1618                       sec%transport_h( 4,jclass) = 0._wp 
     1619                       sec%transport_h( 5,jclass) = 0._wp 
     1620                       sec%transport_h( 6,jclass) = 0._wp 
     1621                       sec%transport_h( 7,jclass) = 0._wp 
     1622                       sec%transport_h( 8,jclass) = 0._wp 
     1623                       sec%transport_h( 9,jclass) = 0._wp 
     1624                       sec%transport_h(10,jclass) = 0._wp 
     1625                    ENDIF 
     1626 
     1627                 ENDIF ! end of test if point is in class 
     1628    
     1629              ENDDO ! end of loop on the classes 
     1630 
     1631           ENDDO ! loop over jk 
     1632 
     1633#if defined key_lim2 || defined key_lim3 
     1634 
     1635           !ICE CASE     
     1636           IF( sec%ll_ice_section )THEN 
     1637 
     1638              IF ( transports_2d_h(1,jsec,jseg) .GE. 0.0 ) THEN 
     1639                 sec%transport_h(11,1) = sec%transport_h(11,1)+transports_2d_h(1,jsec,jseg) 
     1640              ELSE 
     1641                 sec%transport_h(12,1) = sec%transport_h(12,1)+transports_2d_h(1,jsec,jseg) 
     1642              ENDIF 
     1643 
     1644              IF ( transports_2d_h(3,jsec,jseg) .GE. 0.0 ) THEN 
     1645                 sec%transport_h(13,1) = sec%transport_h(13,1)+transports_2d_h(2,jsec,jseg) 
     1646              ELSE 
     1647                 sec%transport_h(14,1) = sec%transport_h(14,1)+transports_2d_h(2,jsec,jseg) 
     1648              ENDIF 
     1649 
     1650           ENDIF !end of ice case 
     1651#endif 
     1652  
     1653        ENDDO !end of loop on the segment 
     1654 
     1655     ELSE  !if sec%nb_point =0 
     1656        sec%transport_h(1:2,:)=0. 
     1657        IF (sec%llstrpond) sec%transport_h(3:10,:)=0. 
     1658        IF (sec%ll_ice_section) sec%transport_h( 11:14,:)=0. 
     1659     ENDIF !end of sec%nb_point =0 case 
     1660 
     1661  END SUBROUTINE dia_dct_sum_h 
     1662  
     1663  SUBROUTINE dia_dct_wri_NOOS(kt,ksec,sec) 
     1664     !!------------------------------------------------------------- 
     1665     !! Write transport output in numdct using NOOS formatting  
     1666     !!  
     1667     !! Purpose: Write  transports in ascii files 
     1668     !!  
     1669     !! Method: 
     1670     !!        1. Write volume transports in "volume_transport" 
     1671     !!           Unit: Sv : area * Velocity / 1.e6  
     1672     !!  
     1673     !!        2. Write heat transports in "heat_transport" 
     1674     !!           Unit: Peta W : area * Velocity * T * rhau * Cp / 1.e15 
     1675     !!  
     1676     !!        3. Write salt transports in "salt_transport" 
     1677     !!           Unit: 10^9 g m^3 / s : area * Velocity * S / 1.e6 
     1678     !! 
    8121679     !!-------------------------------------------------------------  
    813      !! Purpose: Average the transport over nn_dctwri time steps   
    814      !! and sum over the density/salinity/temperature/depth classes  
     1680     !!arguments 
     1681     INTEGER, INTENT(IN)          :: kt          ! time-step 
     1682     TYPE(SECTION), INTENT(INOUT) :: sec         ! section to write    
     1683     INTEGER ,INTENT(IN)          :: ksec        ! section number 
     1684 
     1685     !!local declarations 
     1686     INTEGER               :: jclass,ji             ! Dummy loop 
     1687     CHARACTER(len=2)      :: classe             ! Classname  
     1688     REAL(wp)              :: zbnd1,zbnd2        ! Class bounds 
     1689     REAL(wp)              :: zslope             ! section's slope coeff 
     1690     ! 
     1691     REAL(wp), POINTER, DIMENSION(:):: zsumclasses ! 1D workspace  
     1692     !!-------------------------------------------------------------  
     1693      
     1694      
     1695     IF( lwp ) THEN 
     1696        WRITE(numout,*) " " 
     1697        WRITE(numout,*) "dia_dct_wri_NOOS: write transports through sections at timestep: ", kt 
     1698        WRITE(numout,*) "~~~~~~~~~~~~~~~~~~~~~" 
     1699     ENDIF    
     1700         
     1701     CALL wrk_alloc(nb_type , zsumclasses )   
     1702 
     1703     zsumclasses(:)=0._wp 
     1704     zslope = sec%slopeSection        
     1705 
     1706     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 
     1707 
     1708     DO jclass=1,MAX(1,sec%nb_class-1) 
     1709        zsumclasses(1:nb_type)=zsumclasses(1:nb_type)+sec%transport(1:nb_type,jclass) 
     1710     ENDDO 
     1711  
     1712     classe   = 'total   ' 
     1713     zbnd1   = 0._wp 
     1714     zbnd2   = 0._wp 
     1715 
     1716     IF ( zslope .gt. 0._wp .and. zslope .ne. 10000._wp ) THEN 
     1717        WRITE(numdct_NOOS,'(9e12.4E2)') -(zsumclasses( 1)+zsumclasses( 2)), -zsumclasses( 2),-zsumclasses( 1),   & 
     1718                                        -(zsumclasses( 7)+zsumclasses( 8)), -zsumclasses( 8),-zsumclasses( 7),   & 
     1719                                        -(zsumclasses( 9)+zsumclasses(10)), -zsumclasses(10),-zsumclasses( 9) 
     1720     ELSE 
     1721        WRITE(numdct_NOOS,'(9e12.4E2)')   zsumclasses( 1)+zsumclasses( 2) ,  zsumclasses( 1), zsumclasses( 2),   & 
     1722                                          zsumclasses( 7)+zsumclasses( 8) ,  zsumclasses( 7), zsumclasses( 8),   & 
     1723                                          zsumclasses( 9)+zsumclasses(10) ,  zsumclasses( 9), zsumclasses(10) 
     1724     ENDIF  
     1725 
     1726     DO jclass=1,MAX(1,sec%nb_class-1) 
     1727    
     1728        classe   = 'N       ' 
     1729        zbnd1   = 0._wp 
     1730        zbnd2   = 0._wp 
     1731 
     1732        !insitu density classes transports 
     1733        IF( ( sec%zsigi(jclass)   .NE. 99._wp ) .AND. & 
     1734            ( sec%zsigi(jclass+1) .NE. 99._wp )       )THEN 
     1735           classe = 'DI       ' 
     1736           zbnd1 = sec%zsigi(jclass) 
     1737           zbnd2 = sec%zsigi(jclass+1) 
     1738        ENDIF 
     1739        !potential density classes transports 
     1740        IF( ( sec%zsigp(jclass)   .NE. 99._wp ) .AND. & 
     1741            ( sec%zsigp(jclass+1) .NE. 99._wp )       )THEN 
     1742           classe = 'DP      ' 
     1743           zbnd1 = sec%zsigp(jclass) 
     1744           zbnd2 = sec%zsigp(jclass+1) 
     1745        ENDIF 
     1746        !depth classes transports 
     1747        IF( ( sec%zlay(jclass)    .NE. 99._wp ) .AND. & 
     1748            ( sec%zlay(jclass+1)  .NE. 99._wp )       )THEN  
     1749           classe = 'Z       ' 
     1750           zbnd1 = sec%zlay(jclass) 
     1751           zbnd2 = sec%zlay(jclass+1) 
     1752        ENDIF 
     1753        !salinity classes transports 
     1754        IF( ( sec%zsal(jclass) .NE. 99._wp    ) .AND. & 
     1755            ( sec%zsal(jclass+1) .NE. 99._wp  )       )THEN 
     1756           classe = 'S       ' 
     1757           zbnd1 = sec%zsal(jclass) 
     1758           zbnd2 = sec%zsal(jclass+1)    
     1759        ENDIF 
     1760        !temperature classes transports 
     1761        IF( ( sec%ztem(jclass) .NE. 99._wp     ) .AND. & 
     1762            ( sec%ztem(jclass+1) .NE. 99._wp     )       ) THEN 
     1763           classe = 'T       ' 
     1764           zbnd1 = sec%ztem(jclass) 
     1765           zbnd2 = sec%ztem(jclass+1) 
     1766        ENDIF 
     1767                   
     1768        !write volume transport per class 
     1769        IF ( zslope .gt. 0._wp .and. zslope .ne. 10000._wp ) THEN 
     1770           WRITE(numdct_NOOS,'(9e12.4E2)') -(sec%transport( 1,jclass)+sec%transport( 2,jclass)),-sec%transport( 2,jclass),-sec%transport( 1,jclass), & 
     1771                                           -(sec%transport( 7,jclass)+sec%transport( 8,jclass)),-sec%transport( 8,jclass),-sec%transport( 7,jclass), & 
     1772                                           -(sec%transport( 9,jclass)+sec%transport(10,jclass)),-sec%transport(10,jclass),-sec%transport( 9,jclass) 
     1773        ELSE 
     1774           WRITE(numdct_NOOS,'(9e12.4E2)')   sec%transport( 1,jclass)+sec%transport( 2,jclass) , sec%transport( 1,jclass), sec%transport( 2,jclass), & 
     1775                                             sec%transport( 7,jclass)+sec%transport( 8,jclass) , sec%transport( 7,jclass), sec%transport( 8,jclass), & 
     1776                                             sec%transport( 9,jclass)+sec%transport(10,jclass) , sec%transport( 9,jclass), sec%transport(10,jclass) 
     1777        ENDIF 
     1778 
     1779     ENDDO 
     1780      
     1781     !CALL FLUSH(numdct_NOOS) 
     1782 
     1783     CALL wrk_dealloc(nb_type , zsumclasses )   
     1784 
     1785  END SUBROUTINE dia_dct_wri_NOOS 
     1786 
     1787  SUBROUTINE dia_dct_wri_NOOS_h(hr,ksec,sec) 
     1788     !!------------------------------------------------------------- 
     1789     !! As routine dia_dct_wri_NOOS but for hourly output files 
     1790     !! 
     1791     !! Write transport output in numdct using NOOS formatting  
    8151792     !!  
    816      !! Method:   Sum over relevant grid cells to obtain values   
    817      !!           for each class 
    818      !!              There are several loops:                   
    819      !!              loop on the segment between 2 nodes  
    820      !!              loop on the level jk  
    821      !!              loop on the density/temperature/salinity/level classes  
    822      !!              test on the density/temperature/salinity/level  
     1793     !! Purpose: Write  transports in ascii files 
    8231794     !!  
    824      !!  Note:    Transport through a given section is equal to the sum of transports  
    825      !!           computed on each proc.  
    826      !!           On each proc,transport is equal to the sum of transport computed through  
    827      !!           segments linking each point of sec%listPoint  with the next one.     
    828      !!  
     1795     !! Method: 
     1796     !!        1. Write volume transports in "volume_transport" 
     1797     !!           Unit: Sv : area * Velocity / 1.e6  
     1798     !! 
    8291799     !!-------------------------------------------------------------  
    830      !! * arguments  
    831      TYPE(SECTION),INTENT(INOUT) :: sec  
    832      INTEGER      ,INTENT(IN)    :: jsec        ! numeric identifier of section  
    833   
    834      TYPE(POINT_SECTION) :: k  
    835      INTEGER  :: jk,jseg,jclass                        ! dummy variables for looping on level/segment/classes   
    836      REAL(wp) :: ztn, zsn, zrhoi, zrhop, zsshn, zfsdep ! temperature/salinity/ssh/potential density /depth at u/v point  
     1800     !!arguments 
     1801     INTEGER, INTENT(IN)          :: hr          ! hour 
     1802     TYPE(SECTION), INTENT(INOUT) :: sec         ! section to write    
     1803     INTEGER ,INTENT(IN)          :: ksec        ! section number 
     1804 
     1805     !!local declarations 
     1806     INTEGER               :: jclass,jhr            ! Dummy loop 
     1807     CHARACTER(len=2)      :: classe             ! Classname  
     1808     REAL(wp)              :: zbnd1,zbnd2        ! Class bounds 
     1809     REAL(wp)              :: zslope             ! section's slope coeff 
     1810     ! 
     1811     REAL(wp), POINTER, DIMENSION(:):: zsumclasses ! 1D workspace  
    8371812     !!-------------------------------------------------------------  
    838   
    839      !! Sum the relevant segments to obtain values for each class  
    840      IF(sec%nb_point .NE. 0)THEN     
    841   
    842         !--------------------------------------!  
    843         ! LOOP ON THE SEGMENT BETWEEN 2 NODES  !  
    844         !--------------------------------------!  
    845         DO jseg=1,MAX(sec%nb_point-1,0)  
    846              
    847            !-------------------------------------------------------------------------------------------  
    848            ! Select the appropriate coordinate for computing the velocity of the segment  
    849            !  
    850            !                      CASE(0)                                    Case (2)  
    851            !                      -------                                    --------  
    852            !  listPoint(jseg)                 listPoint(jseg+1)       listPoint(jseg)  F(i,j)        
    853            !      F(i,j)----------V(i+1,j)-------F(i+1,j)                               |  
    854            !                                                                            |  
    855            !                                                                            |  
    856            !                                                                            |  
    857            !                      Case (3)                                            U(i,j)  
    858            !                      --------                                              |  
    859            !                                                                            |  
    860            !  listPoint(jseg+1) F(i,j+1)                                                |  
    861            !                        |                                                   |  
    862            !                        |                                                   |  
    863            !                        |                                 listPoint(jseg+1) F(i,j-1)  
    864            !                        |                                              
    865            !                        |                                              
    866            !                     U(i,j+1)                                              
    867            !                        |                                       Case(1)       
    868            !                        |                                       ------        
    869            !                        |                                              
    870            !                        |                 listPoint(jseg+1)             listPoint(jseg)                             
    871            !                        |                 F(i-1,j)-----------V(i,j) -------f(jseg)                             
    872            ! listPoint(jseg)     F(i,j)  
    873            !   
    874            !-------------------------------------------------------------------------------------------  
    875   
    876            SELECT CASE( sec%direction(jseg) )  
    877            CASE(0)  ;   k = sec%listPoint(jseg)  
    878            CASE(1)  ;   k = POINT_SECTION(sec%listPoint(jseg)%I+1,sec%listPoint(jseg)%J)  
    879            CASE(2)  ;   k = sec%listPoint(jseg)  
    880            CASE(3)  ;   k = POINT_SECTION(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J+1)  
    881            END SELECT  
    882   
    883            !---------------------------|  
    884            !     LOOP ON THE LEVEL     |  
    885            !---------------------------|  
    886            !Sum of the transport on the vertical   
    887            DO jk=1,mbathy(k%I,k%J)  
    888   
    889               ! compute temperature, salinity, insitu & potential density, ssh and depth at U/V point  
    890               SELECT CASE( sec%direction(jseg) )  
    891               CASE(0,1)  
    892                  ztn   = interp(k%I,k%J,jk,'V',tsn(:,:,:,jp_tem) )  
    893                  zsn   = interp(k%I,k%J,jk,'V',tsn(:,:,:,jp_sal) )  
    894                  zrhop = interp(k%I,k%J,jk,'V',rhop)  
    895                  zrhoi = interp(k%I,k%J,jk,'V',rhd*rau0+rau0)  
    896  
    897               CASE(2,3)  
    898                  ztn   = interp(k%I,k%J,jk,'U',tsn(:,:,:,jp_tem) )  
    899                  zsn   = interp(k%I,k%J,jk,'U',tsn(:,:,:,jp_sal) )  
    900                  zrhop = interp(k%I,k%J,jk,'U',rhop)  
    901                  zrhoi = interp(k%I,k%J,jk,'U',rhd*rau0+rau0)  
    902                  zsshn =  0.5*( sshn(k%I,k%J)    + sshn(k%I+1,k%J)    ) * umask(k%I,k%J,1)   
    903               END SELECT  
    904   
    905               zfsdep= fsdept(k%I,k%J,jk)  
    906    
    907               !-------------------------------  
    908               !  LOOP ON THE DENSITY CLASSES |  
    909               !-------------------------------  
    910               !The computation is made for each density/temperature/salinity/depth class  
    911               DO jclass=1,MAX(1,sec%nb_class-1)  
    912   
    913                  !----------------------------------------------!  
    914                  !TEST ON THE DENSITY/SALINITY/TEMPERATURE/LEVEL!   
    915                  !----------------------------------------------!  
    916  
    917                  IF ( (                                                    &  
    918                     ((( zrhop .GE. (sec%zsigp(jclass)+1000.  )) .AND.      &  
    919                     (   zrhop .LE. (sec%zsigp(jclass+1)+1000. ))) .OR.     &  
    920                     ( sec%zsigp(jclass) .EQ. 99.)) .AND.                   &  
    921   
    922                     ((( zrhoi .GE. (sec%zsigi(jclass) + 1000.  )) .AND.    &  
    923                     (   zrhoi .LE. (sec%zsigi(jclass+1)+1000. ))) .OR.     &  
    924                     ( sec%zsigi(jclass) .EQ. 99.)) .AND.                   &  
    925   
    926                     ((( zsn .GT. sec%zsal(jclass)) .AND.                   &  
    927                     (   zsn .LE. sec%zsal(jclass+1))) .OR.                 &  
    928                     ( sec%zsal(jclass) .EQ. 99.)) .AND.                    &  
    929   
    930                     ((( ztn .GE. sec%ztem(jclass)) .AND.                   &  
    931                     (   ztn .LE. sec%ztem(jclass+1))) .OR.                 &  
    932                     ( sec%ztem(jclass) .EQ.99.)) .AND.                     &  
    933   
    934                     ((( zfsdep .GE. sec%zlay(jclass)) .AND.                &  
    935                     (   zfsdep .LE. sec%zlay(jclass+1))) .OR.              &  
    936                     ( sec%zlay(jclass) .EQ. 99. ))                         &  
    937                                                                    ))   THEN  
    938   
    939                     !SUM THE TRANSPORTS FOR EACH CLASSES FOR THE POSITIVE AND NEGATIVE DIRECTIONS  
    940                     !----------------------------------------------------------------------------  
    941                     IF (transports_3d(1,jsec,jseg,jk) .GE. 0.0) THEN   
    942                        sec%transport(1,jclass) = sec%transport(1,jclass)+transports_3d(1,jsec,jseg,jk)*1.E-6  
    943                     ELSE  
    944                        sec%transport(2,jclass) = sec%transport(2,jclass)+transports_3d(1,jsec,jseg,jk)*1.E-6  
    945                     ENDIF  
    946                     IF( sec%llstrpond )THEN  
    947   
    948                        IF ( transports_3d(2,jsec,jseg,jk) .GE. 0.0 ) THEN  
    949                           sec%transport(3,jclass) = sec%transport(3,jclass)+transports_3d(2,jsec,jseg,jk)  
    950                        ELSE  
    951                           sec%transport(4,jclass) = sec%transport(4,jclass)+transports_3d(2,jsec,jseg,jk)  
    952                        ENDIF  
    953   
    954                        IF ( transports_3d(3,jsec,jseg,jk) .GE. 0.0 ) THEN  
    955                           sec%transport(5,jclass) = sec%transport(5,jclass)+transports_3d(3,jsec,jseg,jk)  
    956                        ELSE  
    957                           sec%transport(6,jclass) = sec%transport(6,jclass)+transports_3d(3,jsec,jseg,jk)  
    958                        ENDIF  
    959   
    960                     ELSE  
    961                        sec%transport( 3,jclass) = 0._wp  
    962                        sec%transport( 4,jclass) = 0._wp  
    963                        sec%transport( 5,jclass) = 0._wp  
    964                        sec%transport( 6,jclass) = 0._wp  
    965                     ENDIF  
    966   
    967                  ENDIF ! end of test if point is in class  
    968      
    969               ENDDO ! end of loop on the classes  
    970   
    971            ENDDO ! loop over jk  
    972   
    973 #if defined key_lim2 || defined key_lim3  
    974   
    975            !ICE CASE      
    976            IF( sec%ll_ice_section )THEN  
    977   
    978               IF ( transports_2d(1,jsec,jseg) .GE. 0.0 ) THEN  
    979                  sec%transport( 7,1) = sec%transport( 7,1)+transports_2d(1,jsec,jseg)*1.E-6  
    980               ELSE  
    981                  sec%transport( 8,1) = sec%transport( 8,1)+transports_2d(1,jsec,jseg)*1.E-6  
    982               ENDIF  
    983   
    984               IF ( transports_2d(3,jsec,jseg) .GE. 0.0 ) THEN  
    985                  sec%transport( 9,1) = sec%transport( 9,1)+transports_2d(2,jsec,jseg)*1.E-6  
    986               ELSE  
    987                  sec%transport(10,1) = sec%transport(10,1)+transports_2d(2,jsec,jseg)*1.E-6  
    988               ENDIF  
    989   
    990            ENDIF !end of ice case  
    991 #endif  
    992    
    993         ENDDO !end of loop on the segment  
    994   
    995      ELSE  !if sec%nb_point =0  
    996         sec%transport(1:2,:)=0.  
    997         IF (sec%llstrpond) sec%transport(3:6,:)=0.  
    998         IF (sec%ll_ice_section) sec%transport(7:10,:)=0.  
    999      ENDIF !end of sec%nb_point =0 case  
    1000   
    1001   END SUBROUTINE dia_dct_sum  
    1002    
     1813 
     1814     IF( lwp ) THEN 
     1815        WRITE(numout,*) " " 
     1816        WRITE(numout,*) "dia_dct_wri_NOOS_h: write transports through sections at timestep: ", hr 
     1817        WRITE(numout,*) "~~~~~~~~~~~~~~~~~~~~~" 
     1818     ENDIF    
     1819      
     1820     CALL wrk_alloc(nb_type , zsumclasses )  
     1821 
     1822     zsumclasses(:)=0._wp 
     1823     zslope = sec%slopeSection        
     1824 
     1825     DO jclass=1,MAX(1,sec%nb_class-1) 
     1826        zsumclasses(1:nb_type)=zsumclasses(1:nb_type)+sec%transport_h(1:nb_type,jclass) 
     1827     ENDDO 
     1828  
     1829     !write volume transport per class 
     1830     IF ( zslope .gt. 0._wp .and. zslope .ne. 10000._wp ) THEN 
     1831        z_hr_output(ksec,hr,1)=-(zsumclasses(1)+zsumclasses(2)) 
     1832     ELSE 
     1833        z_hr_output(ksec,hr,1)= (zsumclasses(1)+zsumclasses(2)) 
     1834     ENDIF 
     1835 
     1836     DO jclass=1,MAX(1,sec%nb_class-1) 
     1837        IF ( zslope .gt. 0._wp .and. zslope .ne. 10000._wp ) THEN 
     1838           z_hr_output(ksec,hr,jclass+1)=-(sec%transport_h(1,jclass)+sec%transport_h(2,jclass)) 
     1839        ELSE 
     1840           z_hr_output(ksec,hr,jclass+1)= (sec%transport_h(1,jclass)+sec%transport_h(2,jclass)) 
     1841        ENDIF 
     1842     ENDDO 
     1843 
     1844     IF ( hr .eq. 48._wp ) THEN 
     1845        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 
     1846        DO jhr=25,48 
     1847           WRITE(numdct_NOOS_h,'(11F12.1)')  z_hr_output(ksec,jhr,1), (z_hr_output(ksec,jhr,jclass+1), jclass=1,MAX(1,10) ) 
     1848        ENDDO 
     1849     ENDIF 
     1850 
     1851     CALL wrk_dealloc(nb_type , zsumclasses ) 
     1852 
     1853  END SUBROUTINE dia_dct_wri_NOOS_h 
     1854 
    10031855  SUBROUTINE dia_dct_wri(kt,ksec,sec) 
    10041856     !!------------------------------------------------------------- 
    … …  
    10121864     !!  
    10131865     !!        2. Write heat transports in "heat_transport" 
    1014      !!           Unit: Peta W : area * Velocity * T * rhop * Cp * 1.e-15 
     1866     !!           Unit: Peta W : area * Velocity * T * rhau * Cp / 1.e15 
    10151867     !!  
    10161868     !!        3. Write salt transports in "salt_transport" 
    1017      !!           Unit: 10^9 Kg/m^2/s : area * Velocity * S * rhop * 1.e-9  
     1869     !!           Unit: 10^9 g m^3 / s : area * Velocity * S / 1.e6 
    10181870     !! 
    10191871     !!-------------------------------------------------------------  
    … …  
    10311883     REAL(wp), POINTER, DIMENSION(:):: zsumclasses ! 1D workspace  
    10321884     !!-------------------------------------------------------------  
    1033      CALL wrk_alloc(nb_type_class , zsumclasses )   
     1885     CALL wrk_alloc(nb_type , zsumclasses )   
    10341886 
    10351887     zsumclasses(:)=0._wp 
    … …  
    10421894        zbnd1   = 0._wp 
    10431895        zbnd2   = 0._wp 
    1044         zsumclasses(1:nb_type_class)=zsumclasses(1:nb_type_class)+sec%transport(1:nb_type_class,jclass) 
     1896        zsumclasses(1:nb_type)=zsumclasses(1:nb_type)+sec%transport(1:nb_type,jclass) 
    10451897 
    10461898    
    … …  
    10901942 
    10911943           !write heat transport per class: 
    1092            WRITE(numdct_heat,119) ndastp,kt,ksec,sec%name,zslope,  & 
     1944           WRITE(numdct_temp,119) ndastp,kt,ksec,sec%name,zslope,  & 
    10931945                              jclass,classe,zbnd1,zbnd2,& 
    1094                               sec%transport(3,jclass)*1.e-15,sec%transport(4,jclass)*1.e-15, & 
    1095                               ( sec%transport(3,jclass)+sec%transport(4,jclass) )*1.e-15 
     1946                              sec%transport(7,jclass)*1000._wp*rcp/1.e15,sec%transport(8,jclass)*1000._wp*rcp/1.e15, & 
     1947                              ( sec%transport(7,jclass)+sec%transport(8,jclass) )*1000._wp*rcp/1.e15 
    10961948           !write salt transport per class 
    1097            WRITE(numdct_salt,119) ndastp,kt,ksec,sec%name,zslope,  & 
     1949           WRITE(numdct_sal ,119) ndastp,kt,ksec,sec%name,zslope,  & 
    10981950                              jclass,classe,zbnd1,zbnd2,& 
    1099                               sec%transport(5,jclass)*1.e-9,sec%transport(6,jclass)*1.e-9,& 
    1100                               (sec%transport(5,jclass)+sec%transport(6,jclass))*1.e-9 
     1951                              sec%transport(9,jclass)*1000._wp/1.e9,sec%transport(10,jclass)*1000._wp/1.e9,& 
     1952                              (sec%transport(9,jclass)+sec%transport(10,jclass))*1000._wp/1.e9 
    11011953        ENDIF 
    11021954 
    … …  
    11151967 
    11161968        !write total heat transport 
    1117         WRITE(numdct_heat,119) ndastp,kt,ksec,sec%name,zslope, & 
     1969        WRITE(numdct_temp,119) ndastp,kt,ksec,sec%name,zslope, & 
    11181970                           jclass,"total",zbnd1,zbnd2,& 
    1119                            zsumclasses(3)*1.e-15,zsumclasses(4)*1.e-15,& 
    1120                            (zsumclasses(3)+zsumclasses(4) )*1.e-15 
     1971                           zsumclasses(7)* 1000._wp*rcp/1.e15,zsumclasses(8)* 1000._wp*rcp/1.e15,& 
     1972                           (zsumclasses(7)+zsumclasses(8) )* 1000._wp*rcp/1.e15 
    11211973        !write total salt transport 
    1122         WRITE(numdct_salt,119) ndastp,kt,ksec,sec%name,zslope, & 
     1974        WRITE(numdct_sal ,119) ndastp,kt,ksec,sec%name,zslope, & 
    11231975                           jclass,"total",zbnd1,zbnd2,& 
    1124                            zsumclasses(5)*1.e-9,zsumclasses(6)*1.e-9,& 
    1125                            (zsumclasses(5)+zsumclasses(6))*1.e-9 
     1976                           zsumclasses(9)*1000._wp/1.e9,zsumclasses(10)*1000._wp/1.e9,& 
     1977                           (zsumclasses(9)+zsumclasses(10))*1000._wp/1.e9 
    11261978     ENDIF 
    11271979 
    … …  
    11311983        WRITE(numdct_vol,118) ndastp,kt,ksec,sec%name,zslope,& 
    11321984                              jclass,"ice_vol",zbnd1,zbnd2,& 
    1133                               sec%transport(7,1),sec%transport(8,1),& 
    1134                               sec%transport(7,1)+sec%transport(8,1) 
     1985                              sec%transport(11,1),sec%transport(12,1),& 
     1986                              sec%transport(11,1)+sec%transport(12,1) 
    11351987        !write total ice surface transport 
    11361988        WRITE(numdct_vol,118) ndastp,kt,ksec,sec%name,zslope,& 
    11371989                              jclass,"ice_surf",zbnd1,zbnd2,& 
    1138                               sec%transport(9,1),sec%transport(10,1), & 
    1139                               sec%transport(9,1)+sec%transport(10,1)  
     1990                              sec%transport(13,1),sec%transport(14,1), & 
     1991                              sec%transport(13,1)+sec%transport(14,1)  
    11401992     ENDIF 
    11411993                                               
    … …  
    11431995119 FORMAT(I8,1X,I8,1X,I4,1X,A30,1X,f9.2,1X,I4,3X,A8,1X,2F12.4,5X,3E15.6) 
    11441996 
    1145      CALL wrk_dealloc(nb_type_class , zsumclasses )   
     1997     CALL wrk_dealloc(nb_type , zsumclasses )   
    11461998  END SUBROUTINE dia_dct_wri 
    11471999 
    … …  
    11492001  !!---------------------------------------------------------------------- 
    11502002  !! 
    1151   !!   Purpose: compute temperature/salinity/density at U-point or V-point 
     2003  !!   Purpose: compute Temperature/Salinity/density at U-point or V-point 
    11522004  !!   -------- 
    11532005  !! 
    … …  
    11582010  !!  
    11592011  !! 
    1160   !!    |    I          |    I+1           |    Z=temperature/salinity/density at U-poinT 
     2012  !!    |    I          |    I+1           |    Z=Temperature/Salinity/density at U-poinT 
    11612013  !!    |               |                  | 
    1162   !!  ----------------------------------------  1. Veritcal interpolation: compute zbis 
     2014  !!  ----------------------------------------  1. Veritcale interpolation: compute zbis 
    11632015  !!    |               |                  |       interpolation between ptab(I,J,K) and ptab(I,J,K+1) 
    11642016  !!    |               |                  |       zbis =  
    … …  
    12452097     zdep2 = fsdept(ii2,ij2,kk) - zdepu 
    12462098 
    1247      ! weights 
     2099     !weights 
    12482100     zwgt1 = SQRT( ( 0.5 * zet1 ) * ( 0.5 * zet1 ) + ( zdep1 * zdep1 ) ) 
    12492101     zwgt2 = SQRT( ( 0.5 * zet2 ) * ( 0.5 * zet2 ) + ( zdep2 * zdep2 ) ) 
    12502102   
    12512103     ! result 
    1252      interp = zmsk * ( zwgt2 *  ptab(ii1,ij1,kk) + zwgt1 *  ptab(ii1,ij1,kk) ) / ( zwgt2 + zwgt1 )    
     2104     !JT interp = zmsk * ( zwgt2 *  ptab(ii1,ij1,kk) + zwgt1 *  ptab(ii1,ij1,kk) ) / ( zwgt2 + zwgt1 )    
     2105     interp = umask(ii1,ij1,kk) * ( zwgt2 *  ptab(ii1,ij1,kk) + zwgt1 *  ptab(ii1,ij1,kk) ) / ( zwgt2 + zwgt1 )    
    12532106 
    12542107 
    … …  
    12752128           zbis = ptab(ii2,ij2,kk) + zwgt1 * ( ptab(ii2,ij2,kk-1) - ptab(ii2,ij2,kk) )  
    12762129           ! result 
    1277             interp = zmsk * ( zet2 * ptab(ii1,ij1,kk) + zet1 * zbis )/( zet1 + zet2 ) 
     2130            !JT interp = zmsk * ( zet2 * ptab(ii1,ij1,kk) + zet1 * zbis )/( zet1 + zet2 ) 
     2131            interp = umask(ii1,ij1,kk) * ( zet2 * ptab(ii1,ij1,kk) + zet1 * zbis )/( zet1 + zet2 ) 
    12782132        ELSE 
    12792133           ! zbis 
    12802134           zbis = ptab(ii1,ij1,kk) + zwgt2 * ( ptab(ii1,ij1,kk-1) - ptab(ii1,ij2,kk) ) 
    12812135           ! result 
    1282            interp = zmsk * ( zet2 * zbis + zet1 * ptab(ii2,ij2,kk) )/( zet1 + zet2 ) 
     2136           !JT interp = zmsk * ( zet2 * zbis + zet1 * ptab(ii2,ij2,kk) )/( zet1 + zet2 ) 
     2137           interp = umask(ii1,ij1,kk) * ( zet2 * zbis + zet1 * ptab(ii2,ij2,kk) )/( zet1 + zet2 ) 
    12832138        ENDIF     
    12842139 
    12852140     ELSE 
    1286         interp = zmsk * (  zet2 * ptab(ii1,ij1,kk) + zet1 * ptab(ii2,ij2,kk) )/( zet1 + zet2 ) 
     2141        !JT interp = zmsk * (  zet2 * ptab(ii1,ij1,kk) + zet1 * ptab(ii2,ij2,kk) )/( zet1 + zet2 ) 
     2142        interp = umask(ii1,ij1,kk) * (  zet2 * ptab(ii1,ij1,kk) + zet1 * ptab(ii2,ij2,kk) )/( zet1 + zet2 ) 
    12872143     ENDIF 
    12882144 

  • branches/UKMO/CO6_shelfclimate/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90

    r7566 r7567  
    4444   USE in_out_manager  ! I/O manager 
    4545   USE diadimg         ! dimg direct access file format output 
     46   USE diatmb          ! Top,middle,bottom output 
     47   USE dia25h          ! 25h Mean output 
     48   USE diaregmean      ! regionalmean 
     49   USE diapea          ! pea 
    4650   USE iom 
    4751   USE ioipsl 
    … …  
    379383      CALL wrk_dealloc( jpi , jpj      , z2d ) 
    380384      CALL wrk_dealloc( jpi , jpj, jpk , z3d ) 
     385      ! 
     386      ! If we want tmb values  
     387 
     388      IF (ln_diatmb) THEN 
     389         CALL dia_tmb 
     390      ENDIF 
     391      IF (ln_dia25h) THEN 
     392         CALL dia_25h( kt ) 
     393      ENDIF 
     394       
     395      CALL dia_pea( kt ) 
     396       
     397      IF (ln_diaregmean) THEN 
     398           
     399         !write(*,*) kt,narea,'diawri before dia_regmean' 
     400         CALL dia_regmean( kt ) 
     401         !write(*,*) kt,narea,'diawri after dia_regmean' 
     402      ENDIF 
    381403      ! 
    382404      IF( nn_timing == 1 )   CALL timing_stop('dia_wri') 

  • branches/UKMO/CO6_shelfclimate/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90

    r7566 r7567  
    135135      !!---------------------------------------------------------------------- 
    136136      USE ioipsl 
     137      NAMELIST/namrun/ ln_NOOS   
    137138      NAMELIST/namrun/ cn_ocerst_indir, cn_ocerst_outdir, nn_stocklist, ln_rst_list,               & 
    138          &             nn_no   , cn_exp    , cn_ocerst_in, cn_ocerst_out, ln_rstart , nn_rstctl,   & 
     139         &             nn_no   , cn_exp    , cn_ocerst_in, cn_ocerst_out, ln_rstdate, ln_rstart , nn_rstctl,   & 
    139140         &             nn_it000, nn_itend  , nn_date0    , nn_leapy     , nn_istate , nn_stock ,   & 
    140141         &             nn_write, ln_dimgnnn, ln_mskland  , ln_cfmeta    , ln_clobber, nn_chunksz, nn_euler 
    … …  
    174175         WRITE(numout,*) '      restart output directory        cn_ocerst_outdir= ', cn_ocerst_outdir 
    175176         WRITE(numout,*) '      restart logical                 ln_rstart  = ', ln_rstart 
     177         WRITE(numout,*) '      datestamping of restarts        ln_rstdate  = ', ln_rstdate 
    176178         WRITE(numout,*) '      start with forward time step    nn_euler   = ', nn_euler 
    177179         WRITE(numout,*) '      control of time step            nn_rstctl  = ', nn_rstctl 
    … …  
    189191         WRITE(numout,*) '      multi file dimgout              ln_dimgnnn = ', ln_dimgnnn 
    190192         WRITE(numout,*) '      mask land points                ln_mskland = ', ln_mskland 
     193         ! JT 
     194         WRITE(numout,*) '      NOOS transect diagnostics       ln_NOOS    = ', ln_NOOS 
     195         ! JT 
    191196         WRITE(numout,*) '      additional CF standard metadata ln_cfmeta  = ', ln_cfmeta 
    192197         WRITE(numout,*) '      overwrite an existing file      ln_clobber = ', ln_clobber 

  • branches/UKMO/CO6_shelfclimate/NEMOGCM/NEMO/OPA_SRC/DOM/dommsk.F90

    r7566 r7567  
    3131   USE wrk_nemo        ! Memory allocation 
    3232   USE timing          ! Timing 
     33   USE iom    ! slwa 
    3334 
    3435   IMPLICIT NONE 

  • branches/UKMO/CO6_shelfclimate/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf_bilap.F90

    r7566 r7567  
    2424   USE wrk_nemo        ! Memory Allocation 
    2525   USE timing          ! Timing 
     26#if defined key_bdy   
     27   USE bdy_oce         ! needed for extra diffusion in rim   
     28#endif 
     29 
    2630 
    2731   IMPLICIT NONE 
    … …  
    7983      REAL(wp), POINTER, DIMENSION(:,:  ) :: zcu, zcv 
    8084      REAL(wp), POINTER, DIMENSION(:,:,:) :: zuf, zut, zlu, zlv 
     85      ! 
     86      REAL(wp), DIMENSION(jpi,jpj) :: zfactor  ! multiplier for diffusion 
    8187      !!---------------------------------------------------------------------- 
    8288      ! 
    … …  
    9096         WRITE(numout,*) 'dyn_ldf_bilap : iso-level bilaplacian operator' 
    9197         WRITE(numout,*) '~~~~~~~~~~~~~' 
     98#if defined key_bdy  
     99         IF(lwp) WRITE(numout,*) '~~~~~ using sponge_factor' 
     100#endif 
    92101      ENDIF 
     102 
     103#if defined key_bdy  
     104      zfactor(:,:) = sponge_factor(:,:)  
     105#else  
     106      zfactor(:,:) = 1.0  
     107#endif  
    93108 
    94109!!bug gm this should be enough 

  • branches/UKMO/CO6_shelfclimate/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf_iso.F90

    r7566 r7567  
    3030   USE wrk_nemo        ! Memory Allocation 
    3131   USE timing          ! Timing 
     32#if defined key_bdy   
     33   USE bdy_oce         ! needed for extra diffusion in rim   
     34#endif 
    3235 
    3336   IMPLICIT NONE 
    … …  
    115118      REAL(wp) ::   zuav, zvav, zuwslpi, zuwslpj, zvwslpi, zvwslpj   !   -      - 
    116119      ! 
     120      REAL(wp), DIMENSION(jpi,jpj) :: zfactor  ! multiplier for diffusion 
     121      ! 
    117122      REAL(wp), POINTER, DIMENSION(:,:) :: ziut, zjuf, zjvt, zivf, zdku, zdk1u, zdkv, zdk1v 
    118123      !!---------------------------------------------------------------------- 
    … …  
    126131         IF(lwp) WRITE(numout,*) 'dyn_ldf_iso : iso-neutral laplacian diffusive operator or ' 
    127132         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   s-coordinate horizontal diffusive operator' 
     133#if defined key_bdy  
     134         IF(lwp) WRITE(numout,*) '~~~~~ using sponge_factor' 
     135#endif 
    128136         !                                      ! allocate dyn_ldf_bilap arrays 
    129137         IF( dyn_ldf_iso_alloc() /= 0 )   CALL ctl_stop('STOP', 'dyn_ldf_iso: failed to allocate arrays') 
    … …  
    155163      ENDIF 
    156164 
     165#if defined key_bdy  
     166      zfactor(:,:) = sponge_factor(:,:)  
     167#else  
     168      zfactor(:,:) = 1.0  
     169#endif  
    157170      !                                                ! =============== 
    158171      DO jk = 1, jpkm1                                 ! Horizontal slab 
    … …  
    199212            DO jj = 2, jpjm1 
    200213               DO ji = fs_2, jpi   ! vector opt. 
    201                   zabe1 = (fsahmt(ji,jj,jk)+ahmb0) * e2t(ji,jj) * fse3t(ji,jj,jk) / e1t(ji,jj) 
     214                  zabe1 = zfactor(ji,jj) * (fsahmt(ji,jj,jk)+ahmb0) * e2t(ji,jj) * fse3t(ji,jj,jk) / e1t(ji,jj) 
    202215 
    203216                  zmskt = 1./MAX(  umask(ji-1,jj,jk  )+umask(ji,jj,jk+1)   & 
    204217                     &           + umask(ji-1,jj,jk+1)+umask(ji,jj,jk  ), 1. ) 
    205218 
    206                   zcof1 = - aht0 * e2t(ji,jj) * zmskt * 0.5  * ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) ) 
     219                  zcof1 = - zfactor(ji,jj) * aht0 * e2t(ji,jj) * zmskt * 0.5  * ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) ) 
    207220 
    208221                  ziut(ji,jj) = (  zabe1 * ( ub(ji,jj,jk) - ub(ji-1,jj,jk) )   & 
    … …  
    216229         DO jj = 1, jpjm1 
    217230            DO ji = 1, fs_jpim1   ! vector opt. 
    218                zabe2 = ( fsahmf(ji,jj,jk) + ahmb0 ) * e1f(ji,jj) * fse3f(ji,jj,jk) / e2f(ji,jj) 
     231               zabe2 = zfactor(ji,jj) * ( fsahmf(ji,jj,jk) + ahmb0 ) * e1f(ji,jj) * fse3f(ji,jj,jk) / e2f(ji,jj) 
    219232 
    220233               zmskf = 1./MAX(  umask(ji,jj+1,jk  )+umask(ji,jj,jk+1)   & 
    221234                  &           + umask(ji,jj+1,jk+1)+umask(ji,jj,jk  ), 1. ) 
    222235 
    223                zcof2 = - aht0 * e1f(ji,jj) * zmskf * 0.5  * ( vslp(ji+1,jj,jk) + vslp(ji,jj,jk) ) 
     236               zcof2 = - zfactor(ji,jj) * aht0 * e1f(ji,jj) * zmskf * 0.5  * ( vslp(ji+1,jj,jk) + vslp(ji,jj,jk) ) 
    224237 
    225238               zjuf(ji,jj) = (  zabe2 * ( ub(ji,jj+1,jk) - ub(ji,jj,jk) )   & 
    … …  
    237250         DO jj = 2, jpjm1 
    238251            DO ji = 1, fs_jpim1   ! vector opt. 
    239                zabe1 = ( fsahmf(ji,jj,jk) + ahmb0 ) * e2f(ji,jj) * fse3f(ji,jj,jk) / e1f(ji,jj) 
     252               zabe1 = zfactor(ji,jj) * ( fsahmf(ji,jj,jk) + ahmb0 ) * e2f(ji,jj) * fse3f(ji,jj,jk) / e1f(ji,jj) 
    240253 
    241254               zmskf = 1./MAX(  vmask(ji+1,jj,jk  )+vmask(ji,jj,jk+1)   & 
    242255                  &           + vmask(ji+1,jj,jk+1)+vmask(ji,jj,jk  ), 1. ) 
    243256 
    244                zcof1 = - aht0 * e2f(ji,jj) * zmskf * 0.5 * ( uslp(ji,jj+1,jk) + uslp(ji,jj,jk) ) 
     257               zcof1 = - zfactor(ji,jj) * aht0 * e2f(ji,jj) * zmskf * 0.5 * ( uslp(ji,jj+1,jk) + uslp(ji,jj,jk) ) 
    245258 
    246259               zivf(ji,jj) = (  zabe1 * ( vb(ji+1,jj,jk) - vb(ji,jj,jk) )   & 
    … …  
    269282            DO jj = 2, jpj 
    270283               DO ji = 1, fs_jpim1   ! vector opt. 
    271                   zabe2 = (fsahmt(ji,jj,jk)+ahmb0) * e1t(ji,jj) * fse3t(ji,jj,jk) / e2t(ji,jj) 
     284                  zabe2 = zfactor(ji,jj) * (fsahmt(ji,jj,jk)+ahmb0) * e1t(ji,jj) * fse3t(ji,jj,jk) / e2t(ji,jj) 
    272285 
    273286                  zmskt = 1./MAX(  vmask(ji,jj-1,jk  )+vmask(ji,jj,jk+1)   & 
    274287                     &           + vmask(ji,jj-1,jk+1)+vmask(ji,jj,jk  ), 1. ) 
    275288 
    276                   zcof2 = - aht0 * e1t(ji,jj) * zmskt * 0.5 * ( vslp(ji,jj-1,jk) + vslp(ji,jj,jk) ) 
     289                  zcof2 = - zfactor(ji,jj) * aht0 * e1t(ji,jj) * zmskt * 0.5 * ( vslp(ji,jj-1,jk) + vslp(ji,jj,jk) ) 
    277290 
    278291                  zjvt(ji,jj) = (  zabe2 * ( vb(ji,jj,jk) - vb(ji,jj-1,jk) )   & 

  • branches/UKMO/CO6_shelfclimate/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90

    r7566 r7567  
    4141   USE timing          ! Timing     
    4242   USE sbcapr          ! surface boundary condition: atmospheric pressure 
     43   USE diatmb          ! Top,middle,bottom output 
    4344   USE dynadv, ONLY: ln_dynadv_vec 
    4445#if defined key_agrif 
    … …  
    144145      INTEGER  ::   ji, jj, jk, jn        ! dummy loop indices 
    145146      INTEGER  ::   ikbu, ikbv, noffset      ! local integers 
     147      REAL(wp) ::   zmdi 
    146148      REAL(wp) ::   zraur, z1_2dt_b, z2dt_bf    ! local scalars 
    147149      REAL(wp) ::   zx1, zy1, zx2, zy2         !   -      - 
    … …  
    169171      CALL wrk_alloc( jpi, jpj, zhf ) 
    170172      ! 
     173      zmdi=1.e+20                               !  missing data indicator for masking 
    171174      !                                         !* Local constant initialization 
    172175      z1_12 = 1._wp / 12._wp  
    … …  
    926929      CALL wrk_dealloc( jpi, jpj, zhf ) 
    927930      ! 
     931      IF ( ln_diatmb ) THEN 
     932         CALL iom_put( "baro_u" , un_b*umask(:,:,1)+zmdi*(1-umask(:,:,1 ) ) )  ! Barotropic  U Velocity 
     933         CALL iom_put( "baro_v" , vn_b*vmask(:,:,1)+zmdi*(1-vmask(:,:,1 ) ) )  ! Barotropic  V Velocity 
     934      ENDIF 
    928935      IF( nn_timing == 1 )  CALL timing_stop('dyn_spg_ts') 
    929936      ! 

  • branches/UKMO/CO6_shelfclimate/NEMOGCM/NEMO/OPA_SRC/ICB/icbrst.F90

    r7566 r7567  
    1818   !!---------------------------------------------------------------------- 
    1919   USE par_oce        ! NEMO parameters 
     20   USE phycst         ! for rday 
    2021   USE dom_oce        ! NEMO domain 
    2122   USE in_out_manager ! NEMO IO routines 
     23   USE ioipsl, ONLY : ju2ymds    ! for calendar 
    2224   USE lib_mpp        ! NEMO MPI library, lk_mpp in particular 
    2325   USE netcdf         ! netcdf routines for IO 
    … …  
    231233      INTEGER ::   jn   ! dummy loop index 
    232234      INTEGER ::   ix_dim, iy_dim, ik_dim, in_dim 
    233       CHARACTER(len=256)     :: cl_path 
    234       CHARACTER(len=256)     :: cl_filename 
     235      INTEGER             ::   iyear, imonth, iday 
     236      REAL (wp)           ::   zsec 
     237      REAL (wp)           ::   zfjulday 
     238      CHARACTER(len=256)  :: cl_path 
     239      CHARACTER(len=256)  :: cl_filename 
     240      CHARACTER(LEN=20)   ::   clkt     ! ocean time-step deine as a character 
    235241      TYPE(iceberg), POINTER :: this 
    236242      TYPE(point)  , POINTER :: pt 
    … …  
    240246      cl_path = TRIM(cn_ocerst_outdir) 
    241247      IF( cl_path(LEN_TRIM(cl_path):) /= '/' ) cl_path = TRIM(cl_path) // '/' 
     248      IF ( ln_rstdate ) THEN 
     249         zfjulday = fjulday + rdttra(1) / rday 
     250         IF( ABS(zfjulday - REAL(NINT(zfjulday),wp)) < 0.1 / rday )   zfjulday = REAL(NINT(zfjulday),wp)   ! avoid truncation error 
     251         CALL ju2ymds( zfjulday, iyear, imonth, iday, zsec )            
     252         WRITE(clkt, '(i4.4,2i2.2)') iyear, imonth, iday 
     253      ELSE 
     254         IF( kt > 999999999 ) THEN   ;   WRITE(clkt, *       ) kt 
     255         ELSE                        ;   WRITE(clkt, '(i8.8)') kt 
     256         ENDIF 
     257      ENDIF 
    242258      IF( lk_mpp ) THEN 
    243          WRITE(cl_filename,'(A,"_icebergs_",I8.8,"_restart_",I4.4,".nc")') TRIM(cexper), kt, narea-1 
     259         WRITE(cl_filename,'(A,"_icebergs_",A,"_restart_",I4.4,".nc")') TRIM(cexper), TRIM(ADJUSTL(clkt)), narea-1 
    244260      ELSE 
    245          WRITE(cl_filename,'(A,"_icebergs_",I8.8,"_restart.nc")') TRIM(cexper), kt 
     261         WRITE(cl_filename,'(A,"_icebergs_",A,"_restart.nc")') TRIM(cexper), TRIM(ADJUSTL(clkt)) 
    246262      ENDIF 
    247263      IF (nn_verbose_level >= 0) WRITE(numout,'(2a)') 'icebergs, write_restart: creating ',TRIM(cl_path)//TRIM(cl_filename) 

  • branches/UKMO/CO6_shelfclimate/NEMOGCM/NEMO/OPA_SRC/IOM/in_out_manager.F90

    r7566 r7567  
    3030   CHARACTER(lc) ::   cn_ocerst_outdir !: restart output directory 
    3131   LOGICAL       ::   ln_rstart        !: start from (F) rest or (T) a restart file 
     32   LOGICAL       ::   ln_rstdate       !: datestamping of restarts  
    3233   LOGICAL       ::   ln_rst_list      !: output restarts at list of times (T) or by frequency (F) 
    3334   INTEGER       ::   nn_no            !: job number 
    … …  
    4142   INTEGER       ::   nn_write         !: model standard output frequency 
    4243   INTEGER       ::   nn_stock         !: restart file frequency 
    43    INTEGER, DIMENSION(10) :: nn_stocklist  !: restart dump times 
     44   INTEGER, DIMENSION(100) :: nn_stocklist  !: restart dump times 
    4445   LOGICAL       ::   ln_dimgnnn       !: type of dimgout. (F): 1 file for all proc 
    4546                                                       !:                  (T): 1 file per proc 
    … …  
    4748   LOGICAL       ::   ln_cfmeta        !: output additional data to netCDF files required for compliance with the CF metadata standard 
    4849   LOGICAL       ::   ln_clobber       !: clobber (overwrite) an existing file 
     50   !JT 
     51   LOGICAL       ::   ln_NOOS          !: NOOS transects  diagnostics 
     52   !JT 
    4953   INTEGER       ::   nn_chunksz       !: chunksize (bytes) for NetCDF file (works only with iom_nf90 routines) 
    5054#if defined key_netcdf4 
    … …  
    8387   INTEGER       ::   nwrite                      !: model standard output frequency 
    8488   INTEGER       ::   nstock                      !: restart file frequency 
    85    INTEGER, DIMENSION(10) :: nstocklist           !: restart dump times 
     89   INTEGER, DIMENSION(100) :: nstocklist           !: restart dump times 
    8690 
    8791   !!---------------------------------------------------------------------- 
    … …  
    131135   INTEGER ::   numdct_in       =   -1      !: logical unit for transports computing 
    132136   INTEGER ::   numdct_vol      =   -1      !: logical unit for voulume transports output 
    133    INTEGER ::   numdct_heat     =   -1      !: logical unit for heat    transports output 
    134    INTEGER ::   numdct_salt     =   -1      !: logical unit for salt    transports output 
     137   !JT INTEGER ::   numdct_heat     =   -1      !: logical unit for heat    transports output 
     138   !JT INTEGER ::   numdct_salt     =   -1      !: logical unit for salt    transports output 
     139   !JT 
     140   INTEGER ::   numdct_temp =   -1      !: logical unit for heat    transports output 
     141   INTEGER ::   numdct_sal  =   -1      !: logical unit for salt    transports output 
     142 
     143   INTEGER ::   numdct_NOOS     =   -1      !: logical unit for NOOS    transports output 
     144   INTEGER ::   numdct_NOOS_h   =   -1      !: logical unit for NOOS hourly transports output 
     145   !JT 
    135146   INTEGER ::   numfl           =   -1      !: logical unit for floats ascii output 
    136147   INTEGER ::   numflo          =   -1      !: logical unit for floats ascii output 

  • branches/UKMO/CO6_shelfclimate/NEMOGCM/NEMO/OPA_SRC/IOM/iom.F90

    r7566 r7567  
    4444   USE ioipsl, ONLY :  ju2ymds    ! for calendar 
    4545   USE crs             ! Grid coarsening 
    46  
     46     
    4747   IMPLICIT NONE 
    4848   PUBLIC   !   must be public to be able to access iom_def through iom 
    … …  
    5555   PUBLIC iom_init, iom_swap, iom_open, iom_close, iom_setkt, iom_varid, iom_get, iom_gettime, iom_rstput, iom_put 
    5656   PUBLIC iom_getatt, iom_use, iom_context_finalize 
    57  
     57    
     58    
     59   !JT REGION MEANS 
     60   !INTEGER , PUBLIC ::   n_regions_output = 100 
     61    
     62   INTEGER , PUBLIC ::   n_regions_output 
     63   !JT REGION MEANS 
     64    
    5865   PRIVATE iom_rp0d, iom_rp1d, iom_rp2d, iom_rp3d 
    5966   PRIVATE iom_g0d, iom_g1d, iom_g2d, iom_g3d, iom_get_123d 
    … …  
    106113      REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: z_bnds 
    107114      !!---------------------------------------------------------------------- 
     115       
     116       
     117       
     118      !JT REGION MEANS 
     119      !! read namelist to see if the region mask code is called, if so read the region mask, and count the regions.  
     120       
     121       
     122      REAL(wp),  ALLOCATABLE,   DIMENSION(:,:) ::   tmpregion !: temporary region_mask 
     123      INTEGER, DIMENSION(3) ::   zdimsz   ! number of elements in each of the 3 dimensions (i.e., lon, lat, no of masks, 297,  375,  4) for an array 
     124      INTEGER               ::   zndims   ! number of dimensions in an array (i.e. 3, ) 
     125      INTEGER :: inum, nmasks,ierr,maskno,idmaskvar,tmpint 
     126      REAL(wp), ALLOCATABLE,   DIMENSION(:,:,:)  ::   tmp_region_mask_real   ! tempory region_mask of reals 
     127       
     128      LOGICAL ::   ln_diaregmean  ! region mean calculation 
     129    
     130     
     131      INTEGER ::   ios                  ! Local integer output status for namelist read 
     132      LOGICAL :: ln_diaregmean_ascii  ! region mean calculation ascii output 
     133      LOGICAL :: ln_diaregmean_bin  ! region mean calculation binary output 
     134      LOGICAL :: ln_diaregmean_nc  ! region mean calculation netcdf output 
     135     
     136       
     137       
     138      NAMELIST/nam_diaregmean/ ln_diaregmean,ln_diaregmean_ascii,ln_diaregmean_bin,ln_diaregmean_nc 
     139       
     140      ! read in Namelist.  
     141      !!---------------------------------------------------------------------- 
     142      ! 
     143      REWIND ( numnam_ref )              ! Read Namelist nam_diatmb in referdiatmbence namelist : TMB diagnostics 
     144      READ   ( numnam_ref, nam_diaregmean, IOSTAT=ios, ERR= 901 ) 
     145901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_diaregmean in reference namelist', lwp ) 
     146 
     147      REWIND( numnam_cfg )              ! Namelist nam_diatmb in configuration namelist  TMB diagnostics 
     148      READ  ( numnam_cfg, nam_diaregmean, IOSTAT = ios, ERR = 902 ) 
     149902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_diaregmean in configuration namelist', lwp ) 
     150      IF(lwm) WRITE ( numond, nam_diaregmean ) 
     151 
     152      IF (ln_diaregmean) THEN 
     153       
     154        ! Open region mask for region means, and retrieve the size of the mask (number of levels)           
     155          CALL iom_open ( 'region_mask.nc', inum ) 
     156          idmaskvar = iom_varid( inum, 'mask', kdimsz=zdimsz, kndims=zndims, ldstop = .FALSE.)           
     157          nmasks = zdimsz(3) 
     158           
     159          ! read in the region mask (which contains floating point numbers) into a temporary array of reals. 
     160          ALLOCATE( tmp_region_mask_real(jpi,jpj,nmasks),  STAT= ierr ) 
     161          IF( ierr /= 0 )   CALL ctl_stop( 'dia_regmean_init: failed to allocate tmp_region_mask_real array' ) 
     162           
     163          ! Use jpdom_unknown to read in a n layer mask. 
     164          tmp_region_mask_real(:,:,:) = 0 
     165          CALL iom_get( inum, jpdom_unknown, 'mask', tmp_region_mask_real(1:nlci,1:nlcj,1:nmasks),   & 
     166              &          kstart = (/ mig(1),mjg(1),1 /), kcount = (/ nlci,nlcj,nmasks /) ) 
     167           
     168          CALL iom_close( inum ) 
     169          !Convert the region mask of reals into one of integers.  
     170           
     171           
     172          n_regions_output = 0 
     173          DO maskno = 1,nmasks 
     174              tmpint = maxval(int(tmp_region_mask_real(:,:,maskno))) 
     175              CALL mpp_max( tmpint ) 
     176              n_regions_output = n_regions_output + (tmpint + 1) 
     177          END DO 
     178       
     179           
     180         
     181        WRITE(numout,*)  'IOM: n_regions_output: ',n_regions_output 
     182         
     183      ELSE 
     184        n_regions_output = 1 
     185      ENDIF 
     186       
     187       
     188       
     189      !JT REGION MEANS 
     190       
     191       
     192       
     193       
    108194#if ! defined key_xios2 
    109195      ALLOCATE( z_bnds(jpk,2) ) 
    … …  
    227313      CALL iom_set_axis_attr( "iax_20C", (/ REAL(20,wp) /) ) 
    228314      CALL iom_set_axis_attr( "iax_28C", (/ REAL(28,wp) /) ) 
     315       
     316       
     317       
     318      ! JT Region means.  
     319      CALL iom_set_axis_attr( "region", (/ (REAL(ji,wp), ji=1,n_regions_output) /) ) 
     320 
     321      !CALL iom_set_axis_attr( "region", (/ (REAL(ji,wp), ji=1,100) /) ) 
     322 
    229323       
    230324      ! automatic definitions of some of the xml attributs 
    … …  
    12461340      REAL(wp), DIMENSION(:)  , OPTIONAL, INTENT(in) ::   paxis 
    12471341      REAL(wp), DIMENSION(:,:), OPTIONAL, INTENT(in) ::   bounds 
     1342       
     1343      !INTEGER :: iind_JT 
     1344       
     1345       
     1346      !write(numout,*) 'IOM/iom.F90:iom_set_axis_attr: ',cdid 
     1347       
    12481348      IF ( PRESENT(paxis) ) THEN 
     1349       
     1350        !write(numout,*) 'IOM/iom.F90:iom_set_axis_attr paxis size for: ',cdid,SIZE(paxis) 
     1351        !write(numout,*) 'IOM/iom.F90:iom_set_axis_attr paxis values for: ',cdid,paxis 
     1352        !do iind_JT = 1,SIZE(paxis)         
     1353        !  write(numout,*) 'IOM/iom.F90:iom_set_axis_attr paxis individual values for: ',cdid,iind_JT,paxis(iind_JT) 
     1354        !end do 
     1355         
    12491356#if ! defined key_xios2 
    12501357         IF ( xios_is_valid_axis     (cdid) )   CALL xios_set_axis_attr     ( cdid, size=SIZE(paxis), value=paxis ) 

  • branches/UKMO/CO6_shelfclimate/NEMOGCM/NEMO/OPA_SRC/IOM/restart.F90

    r7566 r7567  
    2121   USE in_out_manager  ! I/O manager 
    2222   USE iom             ! I/O module 
     23   USE ioipsl, ONLY : ju2ymds    ! for calendar 
    2324   USE eosbn2          ! equation of state            (eos bn2 routine) 
    2425   USE trdmxl_oce      ! ocean active mixed layer tracers trends variables 
    … …  
    5455      !!---------------------------------------------------------------------- 
    5556      INTEGER, INTENT(in) ::   kt     ! ocean time-step 
     57      INTEGER             ::   iyear, imonth, iday 
     58      REAL (wp)           ::   zsec 
     59      REAL (wp)           ::   zfjulday 
    5660      !! 
    5761      CHARACTER(LEN=20)   ::   clkt     ! ocean time-step deine as a character 
    5862      CHARACTER(LEN=50)   ::   clname   ! ocean output restart file name 
    59       CHARACTER(lc)       ::   clpath   ! full path to ocean output restart file 
     63      CHARACTER(LEN=150)  ::   clpath   ! full path to ocean output restart file 
    6064      !!---------------------------------------------------------------------- 
    6165      ! 
    … …  
    8185      IF( kt == nitrst - 1 .OR. nstock == 1 .OR. ( kt == nitend .AND. .NOT. lrst_oce ) ) THEN 
    8286         IF( nitrst <= nitend .AND. nitrst > 0 ) THEN  
    83             ! beware of the format used to write kt (default is i8.8, that should be large enough...) 
    84             IF( nitrst > 999999999 ) THEN   ;   WRITE(clkt, *       ) nitrst 
    85             ELSE                            ;   WRITE(clkt, '(i8.8)') nitrst 
     87            IF ( ln_rstdate ) THEN 
     88               zfjulday = fjulday + rdttra(1) / rday 
     89               IF( ABS(zfjulday - REAL(NINT(zfjulday),wp)) < 0.1 / rday )   zfjulday = REAL(NINT(zfjulday),wp)   ! avoid truncation error 
     90               CALL ju2ymds( zfjulday, iyear, imonth, iday, zsec )            
     91               WRITE(clkt, '(i4.4,2i2.2)') iyear, imonth, iday             
     92               IF ( ln_rst_list .AND. ( kt .NE. nitend) ) THEN 
     93                 ! JT IF ( nstock_list_in_use_JT .AND. ( kt .NE. nitend - 1) ) THEN 
     94                 WRITE(clkt, '(i4.4,2i2.2,a1,i10.10)') iyear, imonth, iday,'_',kt !JT 
     95               ENDIF 
     96            ELSE 
     97               ! beware of the format used to write kt (default is i8.8, that should be large enough...) 
     98               IF( nitrst > 999999999 ) THEN   ;   WRITE(clkt, *       ) nitrst 
     99               ELSE                            ;   WRITE(clkt, '(i8.8)') nitrst 
     100               ENDIF 
    86101            ENDIF 
    87102            ! create the file 

  • branches/UKMO/CO6_shelfclimate/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_oce.F90

    r7566 r7567  
    121121   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   sprecip           !: solid precipitation                          [Kg/m2/s] 
    122122   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   fr_i              !: ice fraction = 1 - lead fraction      (between 0 to 1) 
     123   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   pressnow          !: UKMO SHELF pressure 
     124   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   apgu              !: UKMO SHELF pressure forcing 
     125   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   apgv              !: UKMO SHELF pressure forcing 
    123126#if defined key_cpl_carbon_cycle 
    124127   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   atm_co2           !: atmospheric pCO2                             [ppm] 
    … …  
    171174#endif 
    172175         &      ssu_m  (jpi,jpj) , sst_m(jpi,jpj) , frq_m(jpi,jpj) ,      & 
     176         &      pressnow(jpi,jpj), apgu(jpi,jpj)    , apgv(jpi,jpj) ,     & 
    173177         &      ssv_m  (jpi,jpj) , sss_m(jpi,jpj) , ssh_m(jpi,jpj) , STAT=ierr(4) ) 
    174178         ! 

  • branches/UKMO/CO6_shelfclimate/NEMOGCM/NEMO/OPA_SRC/SBC/sbcflx.F90

    r7566 r7567  
    2828   PUBLIC sbc_flx       ! routine called by step.F90 
    2929 
    30    INTEGER , PARAMETER ::   jpfld   = 5   ! maximum number of files to read  
     30   INTEGER , PARAMETER ::   jpfld   = 6   ! maximum number of files to read  
    3131   INTEGER , PARAMETER ::   jp_utau = 1   ! index of wind stress (i-component) file 
    3232   INTEGER , PARAMETER ::   jp_vtau = 2   ! index of wind stress (j-component) file 
    … …  
    3434   INTEGER , PARAMETER ::   jp_qsr  = 4   ! index of solar heat file 
    3535   INTEGER , PARAMETER ::   jp_emp  = 5   ! index of evaporation-precipation file 
     36   INTEGER , PARAMETER ::   jp_press = 6  ! index of pressure for UKMO shelf fluxes 
    3637   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf    ! structure of input fields (file informations, fields read) 
     38   LOGICAL , PUBLIC    ::   ln_shelf_flx = .FALSE. ! UKMO SHELF specific flux flag 
     39   INTEGER             ::   jpfld_local   ! maximum number of files to read (locally modified depending on ln_shelf_flx)  
    3740 
    3841   !! * Substitutions 
    … …  
    8285      REAL(wp) ::   zcdrag = 1.5e-3       ! drag coefficient 
    8386      REAL(wp) ::   ztx, zty, zmod, zcoef ! temporary variables 
     87      REAL     ::   cs                    ! UKMO SHELF: Friction co-efficient at surface 
     88      REAL     ::   totwindspd            ! UKMO SHELF: Magnitude of wind speed vector 
     89 
     90      REAL(wp) ::   rhoa  = 1.22         ! Air density kg/m3 
     91      REAL(wp) ::   cdrag = 1.5e-3       ! drag coefficient  
    8492      !! 
    8593      CHARACTER(len=100) ::  cn_dir                               ! Root directory for location of flx files 
    8694      TYPE(FLD_N), DIMENSION(jpfld) ::   slf_i                    ! array of namelist information structures 
    87       TYPE(FLD_N) ::   sn_utau, sn_vtau, sn_qtot, sn_qsr, sn_emp  ! informations about the fields to be read 
    88       NAMELIST/namsbc_flx/ cn_dir, sn_utau, sn_vtau, sn_qtot, sn_qsr, sn_emp 
     95      TYPE(FLD_N) ::   sn_utau, sn_vtau, sn_qtot, sn_qsr, sn_emp, sn_press  !  informations about the fields to be read 
     96      LOGICAL     ::   ln_foam_flx  = .FALSE.                     ! UKMO FOAM specific flux flag 
     97      NAMELIST/namsbc_flx/ cn_dir, sn_utau, sn_vtau, sn_qtot, sn_qsr, sn_emp,   & 
     98      &                    ln_foam_flx, sn_press, ln_shelf_flx 
    8999      !!--------------------------------------------------------------------- 
    90100      ! 
    … …  
    109119         slf_i(jp_emp ) = sn_emp 
    110120         ! 
    111          ALLOCATE( sf(jpfld), STAT=ierror )        ! set sf structure 
     121            ALLOCATE( sf(jpfld), STAT=ierror )        ! set sf structure 
     122            IF( ln_shelf_flx ) slf_i(jp_press) = sn_press 
     123    
     124            ! define local jpfld depending on shelf_flx logical 
     125            IF( ln_shelf_flx ) THEN 
     126               jpfld_local = jpfld 
     127            ELSE 
     128               jpfld_local = jpfld-1 
     129            ENDIF 
     130            ! 
    112131         IF( ierror > 0 ) THEN    
    113132            CALL ctl_stop( 'sbc_flx: unable to allocate sf structure' )   ;   RETURN   
    114133         ENDIF 
    115          DO ji= 1, jpfld 
     134         DO ji= 1, jpfld_local 
    116135            ALLOCATE( sf(ji)%fnow(jpi,jpj,1) ) 
    117136            IF( slf_i(ji)%ln_tint ) ALLOCATE( sf(ji)%fdta(jpi,jpj,1,2) ) 
    … …  
    132151         ENDIF 
    133152!CDIR COLLAPSE 
     153            !!UKMO SHELF effect of atmospheric pressure on SSH 
     154            ! If using ln_apr_dyn, this is done there so don't repeat here. 
     155            IF( ln_shelf_flx .AND. .NOT. ln_apr_dyn) THEN 
     156               DO jj = 1, jpjm1 
     157                  DO ji = 1, jpim1 
     158                     apgu(ji,jj) = (-1.0/rau0)*(sf(jp_press)%fnow(ji+1,jj,1)-sf(jp_press)%fnow(ji,jj,1))/e1u(ji,jj) 
     159                     apgv(ji,jj) = (-1.0/rau0)*(sf(jp_press)%fnow(ji,jj+1,1)-sf(jp_press)%fnow(ji,jj,1))/e2v(ji,jj) 
     160                  END DO 
     161               END DO 
     162            ENDIF ! ln_shelf_flx 
     163       
    134164         DO jj = 1, jpj                                           ! set the ocean fluxes from read fields 
    135165            DO ji = 1, jpi 
    136                utau(ji,jj) = sf(jp_utau)%fnow(ji,jj,1) 
    137                vtau(ji,jj) = sf(jp_vtau)%fnow(ji,jj,1) 
    138                qns (ji,jj) = sf(jp_qtot)%fnow(ji,jj,1) - sf(jp_qsr)%fnow(ji,jj,1) 
    139                emp (ji,jj) = sf(jp_emp )%fnow(ji,jj,1) 
     166                   IF( ln_shelf_flx ) THEN 
     167                      !! UKMO SHELF - need atmospheric pressure to calculate Haney forcing 
     168                      pressnow(ji,jj) = sf(jp_press)%fnow(ji,jj,1) 
     169                      !! UKMO SHELF flux files contain wind speed not wind stress 
     170                      totwindspd = sqrt((sf(jp_utau)%fnow(ji,jj,1))**2.0 + (sf(jp_vtau)%fnow(ji,jj,1))**2.0) 
     171                      cs = 0.63 + (0.066 * totwindspd) 
     172                      utau(ji,jj) = cs * (rhoa/rau0) * sf(jp_utau)%fnow(ji,jj,1) * totwindspd 
     173                      vtau(ji,jj) = cs * (rhoa/rau0) * sf(jp_vtau)%fnow(ji,jj,1) * totwindspd 
     174                   ELSE 
     175                      utau(ji,jj) = sf(jp_utau)%fnow(ji,jj,1) 
     176                      vtau(ji,jj) = sf(jp_vtau)%fnow(ji,jj,1) 
     177                   ENDIF 
     178                   qsr (ji,jj) = sf(jp_qsr )%fnow(ji,jj,1) 
     179                   IF( ln_foam_flx .OR. ln_shelf_flx ) THEN 
     180                      !! UKMO FOAM flux files contain non-solar heat flux (qns) rather than total heat flux (qtot) 
     181                      qns (ji,jj) = sf(jp_qtot)%fnow(ji,jj,1) 
     182                      !! UKMO FOAM flux files contain the net DOWNWARD freshwater flux P-E rather then E-P 
     183                      emp (ji,jj) = -1. * sf(jp_emp )%fnow(ji,jj,1) 
     184                   ELSE 
     185                      qns (ji,jj) = sf(jp_qtot)%fnow(ji,jj,1) - sf(jp_qsr)%fnow(ji,jj,1) 
     186                      emp (ji,jj) = sf(jp_emp )%fnow(ji,jj,1) 
     187                   ENDIF 
    140188            END DO 
    141189         END DO 
    … …  
    143191         qns(:,:) = qns(:,:) - emp(:,:) * sst_m(:,:) * rcp        ! mass flux is at SST 
    144192         ! 
     193    
     194            !! UKMO FOAM wind fluxes need lbc_lnk calls owing to a bug in interp.exe 
     195            IF( ln_foam_flx ) THEN 
     196               CALL lbc_lnk( utau(:,:), 'U', -1. ) 
     197               CALL lbc_lnk( vtau(:,:), 'V', -1. ) 
     198            ENDIF 
     199     
    145200         !                                                        ! module of wind stress and wind speed at T-point 
    146201         zcoef = 1. / ( zrhoa * zcdrag ) 
    … …  
    162217            WRITE(numout,*)  
    163218            WRITE(numout,*) '        read daily momentum, heat and freshwater fluxes OK' 
    164             DO jf = 1, jpfld 
     219            DO jf = 1, jpfld_local 
    165220               IF( jf == jp_utau .OR. jf == jp_vtau )   zfact =     1. 
    166221               IF( jf == jp_qtot .OR. jf == jp_qsr  )   zfact =     0.1 

  • branches/UKMO/CO6_shelfclimate/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90

    r7566 r7567  
    326326      !!---------------------------------------------------------------------- 
    327327      INTEGER, INTENT(in) ::   kt       ! ocean time step 
     328      REAL(wp)                         ::   zmdi  ! temporary reals 
     329       
     330       
     331      zmdi=1.e+20 !missing data indicator for masking 
     332       
    328333      !!--------------------------------------------------------------------- 
    329334      ! 
    … …  
    443448            &                    'at it= ', kt,' date= ', ndastp 
    444449         IF(lwp) WRITE(numout,*) '~~~~' 
    445          CALL iom_rstput( kt, nitrst, numrow, 'utau_b' , utau ) 
    446          CALL iom_rstput( kt, nitrst, numrow, 'vtau_b' , vtau ) 
    447          CALL iom_rstput( kt, nitrst, numrow, 'qns_b'  , qns  ) 
     450         !CALL iom_rstput( kt, nitrst, numrow, 'utau_b' , utau*tmask(:,:,1)+zmdi*(1-tmask(:,:,1 ) ) ) ! (land masked) # JT 
     451         !CALL iom_rstput( kt, nitrst, numrow, 'vtau_b' , vtau*tmask(:,:,1)+zmdi*(1-tmask(:,:,1 ) ) ) ! (land masked) # JT 
     452         !CALL iom_rstput( kt, nitrst, numrow, 'qns_b'  ,  qns*tmask(:,:,1)+zmdi*(1-tmask(:,:,1 ) )  ) ! (land masked) # JT 
     453         CALL iom_rstput( kt, nitrst, numrow, 'utau_b' , utau*tmask(:,:,1) ) ! (land masked) # JT 
     454         CALL iom_rstput( kt, nitrst, numrow, 'vtau_b' , vtau*tmask(:,:,1) ) ! (land masked) # JT 
     455         CALL iom_rstput( kt, nitrst, numrow, 'qns_b'  ,  qns*tmask(:,:,1) ) ! (land masked) # JT 
    448456         ! The 3D heat content due to qsr forcing is treated in traqsr 
    449457         ! CALL iom_rstput( kt, nitrst, numrow, 'qsr_b'  , qsr  ) 
    450          CALL iom_rstput( kt, nitrst, numrow, 'emp_b'  , emp  ) 
    451          CALL iom_rstput( kt, nitrst, numrow, 'sfx_b'  , sfx  ) 
     458         !CALL iom_rstput( kt, nitrst, numrow, 'emp_b'  , emp*tmask(:,:,1)+zmdi*(1-tmask(:,:,1 ) )  ) ! (land masked) # JT 
     459         !CALL iom_rstput( kt, nitrst, numrow, 'sfx_b'  , sfx*tmask(:,:,1)+zmdi*(1-tmask(:,:,1 ) )  ) ! (land masked) # JT 
     460         CALL iom_rstput( kt, nitrst, numrow, 'emp_b'  , emp*tmask(:,:,1)  ) ! (land masked) # JT 
     461         CALL iom_rstput( kt, nitrst, numrow, 'sfx_b'  , sfx*tmask(:,:,1)  ) ! (land masked) # JT 
    452462      ENDIF 
    453463 
    … …  
    456466      !                                                ! ---------------------------------------- ! 
    457467      IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN 
    458          CALL iom_put( "empmr" , emp  - rnf )                   ! upward water flux 
     468         CALL iom_put( "empmr" , (emp  - rnf)*tmask(:,:,1)+zmdi*(1-tmask(:,:,1 ) ) )    ! upward water flux (land masked) # JT 
     469         !CALL iom_put( "empmr" , emp  - rnf )                   ! upward water flux 
    459470         CALL iom_put( "saltflx", sfx  )                        ! downward salt flux   
    460471                                                                ! (includes virtual salt flux beneath ice  
    461472                                                                ! in linear free surface case) 
    462473         CALL iom_put( "fmmflx", fmmflx  )                      ! Freezing-melting water flux 
    463          CALL iom_put( "qt"    , qns  + qsr )                   ! total heat flux  
    464          CALL iom_put( "qns"   , qns        )                   ! solar heat flux 
    465          CALL iom_put( "qsr"   ,       qsr  )                   ! solar heat flux 
     474         !CALL iom_put( "qt"    , qns  + qsr )                   ! total heat flux  
     475         CALL iom_put( "qt"    , (qns  + qsr)*tmask(:,:,1)+zmdi*(1-tmask(:,:,1 ) ) )    ! total heat flux (land masked) # JT 
     476         CALL iom_put( "qns"   , qns*tmask(:,:,1)+zmdi*(1-tmask(:,:,1 ) )        )                   ! solar heat flux (land masked) # JT 
     477         CALL iom_put( "qsr"   ,       qsr*tmask(:,:,1)+zmdi*(1-tmask(:,:,1 ) )  )                   ! solar heat flux (land masked) # JT 
    466478         IF( nn_ice > 0 .OR. nn_components == jp_iam_opa )   CALL iom_put( "ice_cover", fr_i )   ! ice fraction  
    467          CALL iom_put( "taum"  , taum       )                   ! wind stress module  
    468          CALL iom_put( "wspd"  , wndm       )                   ! wind speed  module over free ocean or leads in presence of sea-ice 
    469       ENDIF 
    470       ! 
    471       CALL iom_put( "utau", utau )   ! i-wind stress   (stress can be updated at  
    472       CALL iom_put( "vtau", vtau )   ! j-wind stress    each time step in sea-ice) 
     479         !CALL iom_put( "taum"  , taum       )                   ! wind stress module   
     480         CALL iom_put( "taum", taum*tmask(:,:,1)+zmdi*(1-tmask(:,:,1 ) ) )   ! wind stress module (land masked) # JT 
     481         CALL iom_put( "wspd"  , wndm*tmask(:,:,1)+zmdi*(1-tmask(:,:,1 ) )       )                   ! wind speed  module over free ocean or leads in presence of sea-ice (land masked) # JT 
     482      ENDIF 
     483      ! 
     484      CALL iom_put( "utau", utau*tmask(:,:,1)+zmdi*(1-tmask(:,:,1 ) ) )   ! i-wind stress   (stress can be updated at  (land masked) # JT 
     485      CALL iom_put( "vtau", vtau*tmask(:,:,1)+zmdi*(1-tmask(:,:,1 ) ) )   ! j-wind stress    each time step in sea-ice) (land masked) # JT 
    473486      ! 
    474487      IF(ln_ctl) THEN         ! print mean trends (used for debugging) 

  • branches/UKMO/CO6_shelfclimate/NEMOGCM/NEMO/OPA_SRC/SBC/sbcssr.F90

    r7566 r7567  
    4242   LOGICAL         ::   ln_sssr_bnd     ! flag to bound erp term  
    4343   REAL(wp)        ::   rn_sssr_bnd     ! ABS(Max./Min.) value of erp term [mm/day] 
     44   LOGICAL         ::   ln_UKMO_haney   ! UKMO specific flag to calculate Haney forcing   
    4445 
    4546   REAL(wp) , ALLOCATABLE, DIMENSION(:) ::   buffer   ! Temporary buffer for exchange 
    … …  
    7980      INTEGER  ::   ierror   ! return error code 
    8081      !! 
     82      REAL(wp) ::   sst1,sst2                      ! sea surface temperature 
     83      REAL(wp) ::   e_sst1, e_sst2                 ! saturation vapour pressure 
     84      REAL(wp) ::   qs1,qs2                        ! specific humidity 
     85      REAL(wp) ::   pr_tmp                         ! temporary variable for pressure 
     86  
     87      REAL(wp), DIMENSION(jpi,jpj) ::  hny_frc1    ! Haney forcing for sensible heat, correction for latent heat    
     88      REAL(wp), DIMENSION(jpi,jpj) ::  hny_frc2    ! Haney forcing for sensible heat, correction for latent heat    
     89      !! 
    8190      CHARACTER(len=100) ::  cn_dir          ! Root directory for location of ssr files 
    8291      TYPE(FLD_N) ::   sn_sst, sn_sss        ! informations about the fields to be read 
    … …  
    95104            ! 
    96105            IF( nn_sstr == 1 ) THEN                                   !* Temperature restoring term 
    97                DO jj = 1, jpj 
    98                   DO ji = 1, jpi 
    99                      zqrp = rn_dqdt * ( sst_m(ji,jj) - sf_sst(1)%fnow(ji,jj,1) ) 
    100                      qns(ji,jj) = qns(ji,jj) + zqrp 
    101                      qrp(ji,jj) = zqrp 
    102                   END DO 
    103                END DO 
     106                  IF( ln_UKMO_haney ) THEN 
     107                     DO jj = 1, jpj 
     108                        DO ji = 1, jpi 
     109                           sst1   =  sst_m(ji,jj) 
     110                           sst2   =  sf_sst(1)%fnow(ji,jj,1)    
     111                           e_sst1 = 10**((0.7859+0.03477*sst1)/(1.+0.00412*sst1)) 
     112                           e_sst2 = 10**((0.7859+0.03477*sst2)/(1.+0.00412*sst2))          
     113                           pr_tmp = 0.01*pressnow(ji,jj)  !pr_tmp = 1012.0 
     114                           qs1    = (0.62197*e_sst1)/(pr_tmp-0.378*e_sst1) 
     115                           qs2    = (0.62197*e_sst2)/(pr_tmp-0.378*e_sst2) 
     116                           hny_frc1(ji,jj) = sst1-sst2                    
     117                           hny_frc2(ji,jj) = qs1-qs2                      
     118                          !Might need to mask off land points. 
     119                           hny_frc1(ji,jj)=-hny_frc1(ji,jj)*wndm(ji,jj)*1.42 
     120                           hny_frc2(ji,jj)=-hny_frc2(ji,jj)*wndm(ji,jj)*4688.0 
     121                           ! JT Have masked out the land points  
     122                           qns(ji,jj)=qns(ji,jj)+(hny_frc1(ji,jj)+hny_frc2(ji,jj))*tmask(ji,jj,1) 
     123                           qrp(ji,jj) = 0.e0 
     124                        END DO 
     125                     END DO 
     126                  ELSE 
     127                     DO jj = 1, jpj 
     128                        DO ji = 1, jpi 
     129                           zqrp = rn_dqdt * ( sst_m(ji,jj) - sf_sst(1)%fnow(ji,jj,1) ) 
     130                           qns(ji,jj) = qns(ji,jj) + zqrp 
     131                           qrp(ji,jj) = zqrp 
     132                        END DO 
     133                     END DO 
     134                  ENDIF 
    104135               CALL iom_put( "qrp", qrp )                             ! heat flux damping 
    105136            ENDIF 
    … …  
    163194      CHARACTER(len=100) ::  cn_dir          ! Root directory for location of ssr files 
    164195      TYPE(FLD_N) ::   sn_sst, sn_sss        ! informations about the fields to be read 
    165       NAMELIST/namsbc_ssr/ cn_dir, nn_sstr, nn_sssr, rn_dqdt, rn_deds, sn_sst, sn_sss, ln_sssr_bnd, rn_sssr_bnd 
     196      NAMELIST/namsbc_ssr/ cn_dir, nn_sstr, nn_sssr, rn_dqdt, rn_deds, sn_sst, sn_sss, ln_sssr_bnd, rn_sssr_bnd, ln_UKMO_haney 
    166197      INTEGER     ::  ios 
    167198      !!---------------------------------------------------------------------- 
    … …  
    189220         WRITE(numout,*) '      flag to bound erp term                 ln_sssr_bnd = ', ln_sssr_bnd 
    190221         WRITE(numout,*) '      ABS(Max./Min.) erp threshold           rn_sssr_bnd = ', rn_sssr_bnd, ' mm/day' 
     222         WRITE(numout,*) '      Haney forcing                          ln_UKMO_haney = ', ln_UKMO_haney 
    191223      ENDIF 
    192224      ! 

  • branches/UKMO/CO6_shelfclimate/NEMOGCM/NEMO/OPA_SRC/TRA/eosbn2.F90

    r7566 r7567  
    12411241      IF(lwm) WRITE( numond, nameos ) 
    12421242      ! 
    1243       rau0        = 1026._wp                 !: volumic mass of reference     [kg/m3] 
     1243      rau0        = 1020._wp                 !: volumic mass of reference     [kg/m3] 
     1244!     rau0        = 1026._wp                 !: volumic mass of reference     [kg/m3] 
    12441245      rcp         = 3991.86795711963_wp      !: heat capacity     [J/K] 
    12451246      ! 

  • branches/UKMO/CO6_shelfclimate/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_tvd.F90

    r7566 r7567  
    100100         IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) ) l_trd = .TRUE. 
    101101      ENDIF 
     102! slwa unless you use l_trdtra too, the above switches off trend calculations for l_trdtrc 
     103         l_trd = .FALSE. 
     104         IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) ) l_trd = .TRUE. 
     105!slwa 
    102106      ! 
    103107      IF( l_trd )  THEN 

  • branches/UKMO/CO6_shelfclimate/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_iso.F90

    r7566 r7567  
    3232   USE wrk_nemo        ! Memory Allocation 
    3333   USE timing          ! Timing 
     34#if defined key_bdy   
     35   USE bdy_oce   
     36#endif   
    3437 
    3538   IMPLICIT NONE 
    … …  
    102105      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta        ! tracer trend  
    103106      REAL(wp)                             , INTENT(in   ) ::   pahtb0     ! background diffusion coef 
     107      !JT 
     108      REAL(wp), DIMENSION(jpi,jpj)     ::   zfactor  ! multiplier for diffusion 
     109      !JT 
    104110      ! 
    105111      INTEGER  ::  ji, jj, jk, jn   ! dummy loop indices 
    … …  
    118124      ! 
    119125 
     126      ! 
     127#if defined key_bdy  
     128      zfactor(:,:) = sponge_factor(:,:)  
     129#else  
     130      zfactor(:,:) = 1.0  
     131#endif 
    120132      IF( kt == kit000 )  THEN 
    121133         IF(lwp) WRITE(numout,*) 
    122134         IF(lwp) WRITE(numout,*) 'tra_ldf_iso : rotated laplacian diffusion operator on ', cdtype 
    123135         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 
     136#if defined key_bdy  
     137         IF(lwp) WRITE(numout,*) '~~~~~ using sponge_factor' 
     138#endif 
    124139      ENDIF 
    125       ! 
    126140      !                                                          ! =========== 
    127141      DO jn = 1, kjpt                                            ! tracer loop 
    … …  
    200214            DO jj = 1 , jpjm1 
    201215               DO ji = 1, fs_jpim1   ! vector opt. 
    202                   zabe1 = ( fsahtu(ji,jj,jk) + pahtb0 ) * re2u_e1u(ji,jj) * fse3u_n(ji,jj,jk) 
    203                   zabe2 = ( fsahtv(ji,jj,jk) + pahtb0 ) * re1v_e2v(ji,jj) * fse3v_n(ji,jj,jk) 
     216                  zabe1 =  zfactor(ji,jj) * ( fsahtu(ji,jj,jk) + pahtb0 ) * re2u_e1u(ji,jj) * fse3u_n(ji,jj,jk) 
     217                  zabe2 =  zfactor(ji,jj) * ( fsahtv(ji,jj,jk) + pahtb0 ) * re1v_e2v(ji,jj) * fse3v_n(ji,jj,jk) 
    204218                  ! 
    205219                  zmsku = 1. / MAX(  tmask(ji+1,jj,jk  ) + tmask(ji,jj,jk+1)   & 

  • branches/UKMO/CO6_shelfclimate/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_lap.F90

    r7566 r7567  
    2626   USE lib_mpp         ! MPP library 
    2727   USE timing          ! Timing 
     28#if defined key_bdy   
     29   USE bdy_oce   
     30#endif   
    2831 
    2932   IMPLICIT NONE 
    … …  
    7376      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb        ! before and now tracer fields 
    7477      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta        ! tracer trend  
     78      !JT 
     79      REAL(wp), DIMENSION(jpi,jpj)     ::   zfactor  ! multiplier for diffusion 
     80      !JT 
    7581      ! 
    7682      INTEGER  ::   ji, jj, jk, jn       ! dummy loop indices 

  • branches/UKMO/CO6_shelfclimate/NEMOGCM/NEMO/OPA_SRC/TRA/traqsr.F90

    r7566 r7567  
    4646   LOGICAL , PUBLIC ::   ln_qsr_ice   !: light penetration for ice-model LIM3 (clem) 
    4747   INTEGER , PUBLIC ::   nn_chldta    !: use Chlorophyll data (=1) or not (=0) 
     48   INTEGER , PUBLIC ::   nn_kd490dta  !: use kd490dta data (=1) or not (=0) 
    4849   REAL(wp), PUBLIC ::   rn_abs       !: fraction absorbed in the very near surface (RGB & 2 bands) 
    4950   REAL(wp), PUBLIC ::   rn_si0       !: very near surface depth of extinction      (RGB & 2 bands) 
    … …  
    5455   REAL(wp) ::   xsi1r                           !: inverse of rn_si1 
    5556   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_chl   ! structure of input Chl (file informations, fields read) 
     57   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_kd490 ! structure of input kd490 (file informations, fields read) 
    5658   INTEGER, PUBLIC ::   nksr              ! levels below which the light cannot penetrate ( depth larger than 391 m) 
    5759   REAL(wp), DIMENSION(3,61) ::   rkrgb   !: tabulated attenuation coefficients for RGB absorption 
    … …  
    306308            ! 
    307309         ENDIF 
     310! slwa 
     311         IF( nn_kd490dta == 1 ) THEN                      !  use KD490 data read in   ! 
     312            !                                             ! ------------------------- ! 
     313               nksr = jpk - 1 
     314               ! 
     315               CALL fld_read( kt, 1, sf_kd490 )     ! Read kd490 data and provide it at the current time step 
     316               ! 
     317               zcoef  = ( 1. - rn_abs ) 
     318               ze0(:,:,1) = rn_abs  * qsr(:,:) 
     319               ze1(:,:,1) = zcoef * qsr(:,:) 
     320               zea(:,:,1) =         qsr(:,:) 
     321               ! 
     322               DO jk = 2, nksr+1 
     323!CDIR NOVERRCHK 
     324                  DO jj = 1, jpj 
     325!CDIR NOVERRCHK    
     326                     DO ji = 1, jpi 
     327                        zc0 = ze0(ji,jj,jk-1) * EXP( - fse3t(ji,jj,jk-1) * xsi0r     ) 
     328                        zc1 = ze1(ji,jj,jk-1) * EXP( - fse3t(ji,jj,jk-1) * sf_kd490(1)%fnow(ji,jj,1) ) 
     329                        ze0(ji,jj,jk) = zc0 
     330                        ze1(ji,jj,jk) = zc1 
     331                        zea(ji,jj,jk) = ( zc0 + zc1 ) * tmask(ji,jj,jk) 
     332                     END DO 
     333                  END DO 
     334               END DO 
     335               ! clem: store attenuation coefficient of the first ocean level 
     336               IF ( ln_qsr_ice ) THEN 
     337                  DO jj = 1, jpj 
     338                     DO ji = 1, jpi 
     339                        zzc0 = rn_abs * EXP( - fse3t(ji,jj,1) * xsi0r     ) 
     340                        zzc1 = zcoef  * EXP( - fse3t(ji,jj,1) * sf_kd490(1)%fnow(ji,jj,1) ) 
     341                        fraqsr_1lev(ji,jj) = 1.0 - ( zzc0 + zzc1 ) * tmask(ji,jj,2)  
     342                     END DO 
     343                  END DO 
     344               ENDIF 
     345               ! 
     346               DO jk = 1, nksr                                        ! compute and add qsr trend to ta 
     347                  qsr_hc(:,:,jk) = r1_rau0_rcp * ( zea(:,:,jk) - zea(:,:,jk+1) ) 
     348               END DO 
     349               zea(:,:,nksr+1:jpk) = 0.e0     !  
     350               CALL iom_put( 'qsr3d', zea )   ! Shortwave Radiation 3D distribution 
     351               ! 
     352        ENDIF   ! use KD490 data 
     353!slwa 
    308354         ! 
    309355         !                                        Add to the general trend 
    … …  
    374420      CHARACTER(len=100) ::   cn_dir   ! Root directory for location of ssr files 
    375421      TYPE(FLD_N)        ::   sn_chl   ! informations about the chlorofyl field to be read 
    376       !! 
    377       NAMELIST/namtra_qsr/  sn_chl, cn_dir, ln_traqsr, ln_qsr_rgb, ln_qsr_2bd, ln_qsr_bio, ln_qsr_ice,  & 
    378          &                  nn_chldta, rn_abs, rn_si0, rn_si1 
     422      TYPE(FLD_N)        ::   sn_kd490 ! informations about the kd490 field to be read 
     423      !! 
     424      NAMELIST/namtra_qsr/  sn_chl, sn_kd490, cn_dir, ln_traqsr, ln_qsr_rgb, ln_qsr_2bd, ln_qsr_bio, ln_qsr_ice,  & 
     425         &                  nn_chldta, rn_abs, rn_si0, rn_si1, nn_kd490dta 
    379426      !!---------------------------------------------------------------------- 
    380427 
    … …  
    409456         WRITE(numout,*) '      RGB & 2 bands: shortess depth of extinction  rn_si0 = ', rn_si0 
    410457         WRITE(numout,*) '      2 bands: longest depth of extinction         rn_si1 = ', rn_si1 
     458         WRITE(numout,*) '      read in KD490 data                       nn_kd490dta  = ', nn_kd490dta 
    411459      ENDIF 
    412460 
    … …  
    422470         IF( ln_qsr_2bd  )   ioptio = ioptio + 1 
    423471         IF( ln_qsr_bio  )   ioptio = ioptio + 1 
     472         IF( nn_kd490dta == 1 )   ioptio = ioptio + 1 
    424473         ! 
    425474         IF( ioptio /= 1 ) & 
    … …  
    431480         IF( ln_qsr_2bd                      )   nqsr =  3 
    432481         IF( ln_qsr_bio                      )   nqsr =  4 
     482         IF( nn_kd490dta == 1                )   nqsr =  5 
    433483         ! 
    434484         IF(lwp) THEN                   ! Print the choice 
    … …  
    438488            IF( nqsr ==  3 )   WRITE(numout,*) '         2 bands light penetration' 
    439489            IF( nqsr ==  4 )   WRITE(numout,*) '         bio-model light penetration' 
     490            IF( nqsr ==  5 )   WRITE(numout,*) '         KD490 light penetration' 
    440491         ENDIF 
    441492         ! 
    … …  
    447498         xsi0r = 1.e0 / rn_si0 
    448499         xsi1r = 1.e0 / rn_si1 
     500         IF( nn_kd490dta == 1 ) THEN           !* KD490 data : set sf_kd490 structure 
     501            IF(lwp) WRITE(numout,*) 
     502            IF(lwp) WRITE(numout,*) '        KD490 read in a file' 
     503            ALLOCATE( sf_kd490(1), STAT=ierror ) 
     504            IF( ierror > 0 ) THEN 
     505               CALL ctl_stop( 'tra_qsr_init: unable to allocate sf_kd490 structure' )   ;   RETURN 
     506            ENDIF 
     507            ALLOCATE( sf_kd490(1)%fnow(jpi,jpj,1)   ) 
     508            IF( sn_kd490%ln_tint )ALLOCATE( sf_kd490(1)%fdta(jpi,jpj,1,2) ) 
     509            !                                        ! fill sf_kd490 with sn_kd490 and control print 
     510            CALL fld_fill( sf_kd490, (/ sn_kd490 /), cn_dir, 'tra_qsr_init',   & 
     511               &                                         'Solar penetration function of read KD490', 'namtra_qsr' ) 
    449512         !                                ! ---------------------------------- ! 
    450          IF( ln_qsr_rgb ) THEN            !  Red-Green-Blue light penetration  ! 
     513         ELSEIF( ln_qsr_rgb ) THEN            !  Red-Green-Blue light penetration  ! 
    451514            !                             ! ---------------------------------- ! 
    452515            ! 

  • branches/UKMO/CO6_shelfclimate/NEMOGCM/NEMO/OPA_SRC/TRA/trasbc.F90

    r7566 r7567  
    2525   USE trd_oce         ! trends: ocean variables 
    2626   USE trdtra          ! trends manager: tracers  
     27   USE tradwl          ! solar radiation penetration (downwell method) 
    2728   ! 
    2829   USE in_out_manager  ! I/O manager 
    … …  
    138139 
    139140!!gm      IF( .NOT.ln_traqsr )   qsr(:,:) = 0.e0   ! no solar radiation penetration 
    140       IF( .NOT.ln_traqsr ) THEN     ! no solar radiation penetration 
     141      IF( .NOT.ln_traqsr .and. .NOT.ln_tradwl ) THEN     ! no solar radiation penetration 
    141142         qns(:,:) = qns(:,:) + qsr(:,:)      ! total heat flux in qns 
    142143         qsr(:,:) = 0.e0                     ! qsr set to zero 

  • branches/UKMO/CO6_shelfclimate/NEMOGCM/NEMO/OPA_SRC/TRD/trdtra.F90

    r7566 r7567  
    203203         DO jj = 2, jpjm1 
    204204            DO ji = fs_2, fs_jpim1   ! vector opt. 
     205#if defined key_tracer_budget 
     206!              ptrd(ji,jj,jk) = - (     pf (ji,jj,jk) - pf (ji-ii,jj-ij,jk-ik)  )  * tmask(ji,jj,jk) 
     207               ptrd(ji,jj,jk) = -      pf (ji,jj,jk) * tmask(ji,jj,jk) 
     208#else 
    205209               ptrd(ji,jj,jk) = - (     pf (ji,jj,jk) - pf (ji-ii,jj-ij,jk-ik)                        & 
    206210                 &                  - ( pun(ji,jj,jk) - pun(ji-ii,jj-ij,jk-ik) ) * ptn(ji,jj,jk)  )   & 
    207211                 &              / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )  * tmask(ji,jj,jk) 
     212#endif 
    208213            END DO 
    209214         END DO 

  • branches/UKMO/CO6_shelfclimate/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfmxl.F90

    r7566 r7567  
    2222   USE timing          ! Timing 
    2323   USE trc_oce, ONLY : lk_offline ! offline flag 
     24   USE lbclnk  
     25   USE eosbn2          ! Equation of state 
    2426 
    2527   IMPLICIT NONE 
    … …  
    2729 
    2830   PUBLIC   zdf_mxl       ! called by step.F90 
     31   PUBLIC   zdf_mxl_tref  ! called by asminc.F90 
    2932   PUBLIC   zdf_mxl_alloc ! Used in zdf_tke_init 
    3033 
    … …  
    3336   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hmlp    !: mixed layer depth  (rho=rho0+zdcrit) [m] 
    3437   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hmlpt   !: mixed layer depth at t-points        [m] 
     38   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hmld_tref  !: mixed layer depth at t-points - temperature criterion [m] 
     39   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hmld_kara  !: mixed layer depth of Kara et al   [m] 
    3540 
    3641   REAL(wp), PUBLIC ::   rho_c = 0.01_wp    !: density criterion for mixed layer depth 
    3742   REAL(wp)         ::   avt_c = 5.e-4_wp   ! Kz criterion for the turbocline depth 
     43 
     44   ! Namelist variables for  namzdf_karaml  
     45  
     46   LOGICAL   :: ln_kara            ! Logical switch for calculating kara mixed 
     47                                     ! layer 
     48   LOGICAL   :: ln_kara_write25h   ! Logical switch for 25 hour outputs 
     49   INTEGER   :: jpmld_type         ! mixed layer type              
     50   REAL(wp)  :: ppz_ref            ! depth of initial T_ref  
     51   REAL(wp)  :: ppdT_crit          ! Critical temp diff   
     52   REAL(wp)  :: ppiso_frac         ! Fraction of ppdT_crit used  
     53    
     54   !Used for 25h mean 
     55   LOGICAL, PRIVATE :: kara_25h_init = .TRUE.   !Logical used to initalise 25h  
     56                                                !outputs. Necissary, because we need to 
     57                                                !initalise the kara_25h on the zeroth 
     58                                                !timestep (i.e in the nemogcm_init call) 
     59   REAL, PRIVATE, ALLOCATABLE, DIMENSION(:,:) :: hmld_kara_25h 
    3860 
    3961   !! * Substitutions 
    … …  
    5274      zdf_mxl_alloc = 0      ! set to zero if no array to be allocated 
    5375      IF( .NOT. ALLOCATED( nmln ) ) THEN 
    54          ALLOCATE( nmln(jpi,jpj), hmld(jpi,jpj), hmlp(jpi,jpj), hmlpt(jpi,jpj), STAT= zdf_mxl_alloc ) 
     76         ALLOCATE( nmln(jpi,jpj), hmld(jpi,jpj), hmlp(jpi,jpj), hmlpt(jpi,jpj), & 
     77         &                           hmld_tref(jpi,jpj), STAT= zdf_mxl_alloc ) 
    5578         ! 
    5679         IF( lk_mpp             )   CALL mpp_sum ( zdf_mxl_alloc ) 
    … …  
    123146      END DO 
    124147      ! depth of the mixing and mixed layers 
     148 
     149      CALL zdf_mxl_kara( kt ) ! kara MLD 
     150 
    125151      DO jj = 1, jpj 
    126152         DO ji = 1, jpi 
    … …  
    146172   END SUBROUTINE zdf_mxl 
    147173 
     174 
     175   SUBROUTINE zdf_mxl_tref() 
     176      !!---------------------------------------------------------------------- 
     177      !!                  ***  ROUTINE zdf_mxl_tref  *** 
     178      !!                    
     179      !! ** Purpose :   Compute the mixed layer depth with temperature criteria. 
     180      !! 
     181      !! ** Method  :   The temperature-defined mixed layer depth is required 
     182      !!                   when assimilating SST in a 2D analysis.  
     183      !! 
     184      !! ** Action  :   hmld_tref 
     185      !!---------------------------------------------------------------------- 
     186      ! 
     187      INTEGER  ::   ji, jj, jk   ! dummy loop indices 
     188      REAL(wp) ::   t_ref               ! Reference temperature   
     189      REAL(wp) ::   temp_c = 0.2        ! temperature criterion for mixed layer depth   
     190      !!---------------------------------------------------------------------- 
     191      ! 
     192      ! Initialise array 
     193      IF( zdf_mxl_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'zdf_mxl_tref : unable to allocate arrays' ) 
     194       
     195      !For the AMM model assimiation uses a temperature based mixed layer depth   
     196      !This is defined here   
     197      DO jj = 1, jpj   
     198         DO ji = 1, jpi   
     199           hmld_tref(ji,jj)=fsdept(ji,jj,1  )    
     200           IF(ssmask(ji,jj) > 0.)THEN   
     201             t_ref=tsn(ji,jj,1,jp_tem)  
     202             DO jk=2,jpk   
     203               IF(ssmask(ji,jj)==0.)THEN   
     204                  hmld_tref(ji,jj)=fsdept(ji,jj,jk )   
     205                  EXIT   
     206               ELSEIF( ABS(tsn(ji,jj,jk,jp_tem)-t_ref) < temp_c)THEN   
     207                  hmld_tref(ji,jj)=fsdept(ji,jj,jk )   
     208               ELSE   
     209                  EXIT   
     210               ENDIF   
     211             ENDDO   
     212           ENDIF   
     213         ENDDO   
     214      ENDDO 
     215 
     216   END SUBROUTINE zdf_mxl_tref 
     217 
     218   SUBROUTINE zdf_mxl_kara( kt )  
     219      !!----------------------------------------------------------------------------------  
     220      !!                    ***  ROUTINE zdf_mxl_kara  ***  
     221      !                                                                         
     222      !   Calculate mixed layer depth according to the definition of           
     223      !   Kara et al, 2000, JGR, 105, 16803.   
     224      !  
     225      !   If mld_type=1 the mixed layer depth is calculated as the depth at which the   
     226      !   density has increased by an amount equivalent to a temperature difference of   
     227      !   0.8C at the surface.  
     228      !  
     229      !   For other values of mld_type the mixed layer is calculated as the depth at   
     230      !   which the temperature differs by 0.8C from the surface temperature.   
     231      !                                                                         
     232      !   Original version: David Acreman                                       
     233      !  
     234      !!----------------------------------------------------------------------------------- 
     235      
     236      INTEGER, INTENT( in ) ::   kt   ! ocean time-step index  
     237  
     238      NAMELIST/namzdf_karaml/ ln_kara,jpmld_type, ppz_ref, ppdT_crit, & 
     239      &                       ppiso_frac, ln_kara_write25h  
     240  
     241      ! Local variables                                                         
     242      REAL, DIMENSION(jpi,jpk) :: ppzdep      ! depth for use in calculating d(rho)  
     243      REAL(wp), DIMENSION(jpi,jpj,jpts) :: ztsn_2d  !Local version of tsn  
     244      LOGICAL :: ll_found(jpi,jpj)              ! Is T_b to be found by interpolation ?  
     245      LOGICAL :: ll_belowml(jpi,jpj,jpk)        ! Flag points below mixed layer when ll_found=F  
     246      INTEGER :: ji, jj, jk                     ! loop counter  
     247      INTEGER :: ik_ref(jpi,jpj)                ! index of reference level  
     248      INTEGER :: ik_iso(jpi,jpj)                ! index of last uniform temp level  
     249      REAL    :: zT(jpi,jpj,jpk)                ! Temperature or denisty  
     250      REAL    :: zT_ref(jpi,jpj)                ! reference temperature  
     251      REAL    :: zT_b                           ! base temperature  
     252      REAL    :: zdTdz(jpi,jpj,jpk-2)           ! gradient of zT  
     253      REAL    :: zmoddT(jpi,jpj,jpk-2)          ! Absolute temperature difference  
     254      REAL    :: zdz                            ! depth difference  
     255      REAL    :: zdT                            ! temperature difference  
     256      REAL    :: zdelta_T(jpi,jpj)              ! difference critereon  
     257      REAL    :: zRHO1(jpi,jpj), zRHO2(jpi,jpj) ! Densities 
     258      INTEGER, SAVE :: idt, i_steps 
     259      INTEGER, SAVE :: i_cnt_25h     ! Counter for 25 hour means 
     260      INTEGER :: ios                 ! Local integer output status for namelist read 
     261 
     262      
     263  
     264      !!-------------------------------------------------------------------------------------  
     265  
     266      IF( kt == nit000 ) THEN  
     267         ! Default values  
     268         ln_kara      = .FALSE. 
     269         ln_kara_write25h = .FALSE.  
     270         jpmld_type   = 1      
     271         ppz_ref      = 10.0  
     272         ppdT_crit    = 0.2   
     273         ppiso_frac   = 0.1    
     274         ! Read namelist 
     275         REWIND ( numnam_ref ) 
     276         READ   ( numnam_ref, namzdf_karaml, IOSTAT=ios, ERR= 901 ) 
     277901      IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_karaml in reference namelist', lwp ) 
     278 
     279         REWIND( numnam_cfg )              ! Namelist nam_diadiaop in configuration namelist  3D hourly diagnostics 
     280         READ  ( numnam_cfg,  namzdf_karaml, IOSTAT = ios, ERR = 902 ) 
     281902      IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_karaml in configuration namelist', lwp ) 
     282         IF(lwm) WRITE ( numond, namzdf_karaml ) 
     283  
     284 
     285         WRITE(numout,*) '===== Kara mixed layer ====='  
     286         WRITE(numout,*) 'ln_kara = ',    ln_kara 
     287         WRITE(numout,*) 'jpmld_type = ', jpmld_type  
     288         WRITE(numout,*) 'ppz_ref = ',    ppz_ref  
     289         WRITE(numout,*) 'ppdT_crit = ',  ppdT_crit  
     290         WRITE(numout,*) 'ppiso_frac = ', ppiso_frac 
     291         WRITE(numout,*) 'ln_kara_write25h = ', ln_kara_write25h 
     292         WRITE(numout,*) '============================'  
     293       
     294         IF ( .NOT.ln_kara ) THEN 
     295            WRITE(numout,*) "ln_kara not set; Kara mixed layer not calculated"  
     296         ELSE 
     297            IF (.NOT. ALLOCATED(hmld_kara) ) ALLOCATE( hmld_kara(jpi,jpj) ) 
     298            IF ( ln_kara_write25h .AND. kara_25h_init ) THEN 
     299               i_cnt_25h = 0 
     300               IF (.NOT. ALLOCATED(hmld_kara_25h) ) & 
     301               &   ALLOCATE( hmld_kara_25h(jpi,jpj) ) 
     302               hmld_kara_25h = 0._wp 
     303               IF( nacc == 1 ) THEN 
     304                  idt = INT(rdtmin) 
     305               ELSE 
     306                  idt = INT(rdt) 
     307               ENDIF 
     308               IF( MOD( 3600,idt) == 0 ) THEN  
     309                  i_steps = 3600 / idt   
     310               ELSE  
     311                  CALL ctl_stop('STOP', & 
     312                  & 'zdf_mxl_kara: timestep must give MOD(3600,rdt)'// & 
     313                  & ' = 0 otherwise no hourly values are possible')  
     314               ENDIF   
     315            ENDIF 
     316         ENDIF 
     317      ENDIF 
     318       
     319      IF ( ln_kara ) THEN 
     320          
     321         !set critical ln_kara 
     322         ztsn_2d = tsn(:,:,1,:) 
     323         ztsn_2d(:,:,jp_tem) = ztsn_2d(:,:,jp_tem) + ppdT_crit 
     324      
     325         ! Set the mixed layer depth criterion at each grid point  
     326         ppzdep = 0._wp 
     327         IF( jpmld_type == 1 ) THEN                                          
     328            CALL eos ( tsn(:,:,1,:), & 
     329            &          ppzdep(:,:), zRHO1(:,:) )  
     330            CALL eos ( ztsn_2d(:,:,:), & 
     331            &           ppzdep(:,:), zRHO2(:,:) )  
     332            zdelta_T(:,:) = abs( zRHO1(:,:) - zRHO2(:,:) ) * rau0  
     333            ! RHO from eos (2d version) doesn't calculate north or east halo:  
     334            CALL lbc_lnk( zdelta_T, 'T', 1. )  
     335            zT(:,:,:) = rhop(:,:,:)  
     336         ELSE  
     337            zdelta_T(:,:) = ppdT_crit                       
     338            zT(:,:,:) = tsn(:,:,:,jp_tem)                            
     339         ENDIF  
     340      
     341         ! Calculate the gradient of zT and absolute difference for use later  
     342         DO jk = 1 ,jpk-2  
     343            zdTdz(:,:,jk)  =    ( zT(:,:,jk+1) - zT(:,:,jk) ) / fse3w(:,:,jk+1)  
     344            zmoddT(:,:,jk) = abs( zT(:,:,jk+1) - zT(:,:,jk) )  
     345         END DO  
     346      
     347         ! Find density/temperature at the reference level (Kara et al use 10m).           
     348         ! ik_ref is the index of the box centre immediately above or at the reference level  
     349         ! Find ppz_ref in the array of model level depths and find the ref     
     350         ! density/temperature by linear interpolation.                                    
     351         ik_ref = -1 
     352         DO jk = jpkm1, 2, -1  
     353            WHERE( fsdept(:,:,jk) > ppz_ref )  
     354               ik_ref(:,:) = jk - 1  
     355               zT_ref(:,:) = zT(:,:,jk-1) + & 
     356               &             zdTdz(:,:,jk-1) * ( ppz_ref - fsdept(:,:,jk-1) )  
     357            ENDWHERE  
     358         END DO 
     359         IF ( ANY( ik_ref  < 0 ) .OR. ANY( ik_ref  > jpkm1 ) ) THEN 
     360            CALL ctl_stop( "STOP", & 
     361            & "zdf_mxl_kara: unable to find reference level for kara ML" )  
     362         ELSE 
     363            ! If the first grid box centre is below the reference level then use the  
     364            ! top model level to get zT_ref  
     365            WHERE( fsdept(:,:,1) > ppz_ref )   
     366               zT_ref = zT(:,:,1)  
     367               ik_ref = 1  
     368            ENDWHERE  
     369      
     370            ! Search for a uniform density/temperature region where adjacent levels           
     371            ! differ by less than ppiso_frac * deltaT.                                       
     372            ! ik_iso is the index of the last level in the uniform layer   
     373            ! ll_found indicates whether the mixed layer depth can be found by interpolation  
     374            ik_iso(:,:)   = ik_ref(:,:)  
     375            ll_found(:,:) = .false.  
     376            DO jj = 1, nlcj  
     377               DO ji = 1, nlci  
     378                 !CDIR NOVECTOR  
     379                  DO jk = ik_ref(ji,jj), mbathy(ji,jj)-1  
     380                     IF( zmoddT(ji,jj,jk) > ( ppiso_frac * zdelta_T(ji,jj) ) ) THEN  
     381                        ik_iso(ji,jj)   = jk  
     382                        ll_found(ji,jj) = ( zmoddT(ji,jj,jk) > zdelta_T(ji,jj) )  
     383                        EXIT  
     384                     ENDIF  
     385                  END DO  
     386               END DO  
     387            END DO  
     388      
     389            ! Use linear interpolation to find depth of mixed layer base where possible  
     390            hmld_kara(:,:) = ppz_ref  
     391            DO jj = 1, jpj  
     392               DO ji = 1, jpi  
     393                  IF( ll_found(ji,jj) .and. tmask(ji,jj,1) == 1.0 ) THEN  
     394                     zdz =  abs( zdelta_T(ji,jj) / zdTdz(ji,jj,ik_iso(ji,jj)) )  
     395                     hmld_kara(ji,jj) = fsdept(ji,jj,ik_iso(ji,jj)) + zdz  
     396                  ENDIF  
     397               END DO  
     398            END DO  
     399      
     400            ! If ll_found = .false. then calculate MLD using difference of zdelta_T     
     401            ! from the reference density/temperature  
     402      
     403            ! Prevent this section from working on land points  
     404            WHERE( tmask(:,:,1) /= 1.0 )  
     405               ll_found = .true.  
     406            ENDWHERE  
     407      
     408            DO jk = 1, jpk  
     409               ll_belowml(:,:,jk) = abs( zT(:,:,jk) - zT_ref(:,:) ) >= & 
     410               & zdelta_T(:,:) 
     411            END DO  
     412      
     413            ! Set default value where interpolation cannot be used (ll_found=false)   
     414            DO jj = 1, jpj  
     415               DO ji = 1, jpi  
     416                  IF( .NOT. ll_found(ji,jj) )  & 
     417                  &   hmld_kara(ji,jj) = fsdept(ji,jj,mbathy(ji,jj))  
     418               END DO  
     419            END DO  
     420      
     421            DO jj = 1, jpj  
     422               DO ji = 1, jpi  
     423                  !CDIR NOVECTOR  
     424                  DO jk = ik_ref(ji,jj)+1, mbathy(ji,jj)  
     425                     IF( ll_found(ji,jj) ) EXIT  
     426                     IF( ll_belowml(ji,jj,jk) ) THEN                 
     427                        zT_b = zT_ref(ji,jj) + zdelta_T(ji,jj) * & 
     428                        &      SIGN(1.0, zdTdz(ji,jj,jk-1) )  
     429                        zdT  = zT_b - zT(ji,jj,jk-1)                                       
     430                        zdz  = zdT / zdTdz(ji,jj,jk-1)                                        
     431                        hmld_kara(ji,jj) = fsdept(ji,jj,jk-1) + zdz  
     432                        EXIT                                                    
     433                     ENDIF  
     434                  END DO  
     435               END DO  
     436            END DO  
     437      
     438            hmld_kara(:,:) = hmld_kara(:,:) * tmask(:,:,1)  
     439  
     440            IF(  ln_kara_write25h  ) THEN 
     441               !Double IF required as i_steps not defined if ln_kara_write25h = 
     442               ! FALSE 
     443               IF ( ( MOD( kt, i_steps ) == 0 ) .OR.  kara_25h_init ) THEN 
     444                  hmld_kara_25h = hmld_kara_25h + hmld_kara 
     445                  i_cnt_25h = i_cnt_25h + 1 
     446                  IF ( kara_25h_init ) kara_25h_init = .FALSE. 
     447               ENDIF 
     448            ENDIF 
     449  
     450#if defined key_iomput  
     451            IF( kt /= nit000 ) THEN  
     452               CALL iom_put( "mldkara"  , hmld_kara )    
     453               IF( ( MOD( i_cnt_25h, 25) == 0 ) .AND.  ln_kara_write25h ) & 
     454                  CALL iom_put( "kara25h"  , ( hmld_kara_25h / 25._wp ) ) 
     455            ENDIF 
     456#endif 
     457  
     458         ENDIF 
     459      ENDIF 
     460        
     461   END SUBROUTINE zdf_mxl_kara  
     462 
    148463   !!====================================================================== 
    149464END MODULE zdfmxl 

  • branches/UKMO/CO6_shelfclimate/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90

    r7566 r7567  
    8585   USE stopar 
    8686   USE stopts 
     87   USE diatmb          ! Top,middle,bottom output 
     88   USE diaregmean      ! Top,middle,bottom output 
     89   USE diapea          ! Top,middle,bottom output 
     90   USE dia25h          ! 25h mean output 
    8791 
    8892   IMPLICIT NONE 
    … …  
    475479      IF( lk_asminc     )   CALL asm_inc_init   ! Initialize assimilation increments 
    476480      IF(lwp) WRITE(numout,*) 'Euler time step switch is ', neuler 
     481                            CALL dia_tmb_init  ! TMB outputs 
     482                            CALL mpp_max( nstop ) 
     483                            !write(*,*) nstop, narea, 'nstop nemogcm before dia_regmean_init' 
     484                              
     485                            CALL dia_regmean_init  ! TMB outputs 
     486                            CALL dia_pea_init  ! TMB outputs 
     487                            CALL mpp_max( nstop ) 
     488                            !write(*,*)  nstop, narea, 'nstop nemogcm finished dia_regmean_init' 
     489                            CALL dia_25h_init  ! 25h mean  outputs 
    477490      ! 
    478491   END SUBROUTINE nemo_init 
    … …  
    608621      IF( numout          /=  6 )   CLOSE( numout          )   ! standard model output file 
    609622      IF( numdct_vol      /= -1 )   CLOSE( numdct_vol      )   ! volume transports 
    610       IF( numdct_heat     /= -1 )   CLOSE( numdct_heat     )   ! heat transports 
    611       IF( numdct_salt     /= -1 )   CLOSE( numdct_salt     )   ! salt transports 
     623      !JT IF( numdct_heat     /= -1 )   CLOSE( numdct_heat     )   ! heat transports 
     624      !JT IF( numdct_salt     /= -1 )   CLOSE( numdct_salt     )   ! salt transports 
     625      IF( numdct_vol  /= -1 )   CLOSE( numdct_vol  )   ! volume transports 
     626      IF( numdct_temp /= -1 )   CLOSE( numdct_temp )   ! heat transports 
     627      IF( numdct_sal  /= -1 )   CLOSE( numdct_sal  )   ! salt transports 
     628      IF( numdct_NOOS /= -1 )   CLOSE( numdct_NOOS )   ! NOOS transports 
     629 
    612630 
    613631      ! 
    … …  
    630648      USE ldftra_oce, ONLY: ldftra_oce_alloc 
    631649      USE trc_oce   , ONLY: trc_oce_alloc 
     650      USE diainsitutem, ONLY: insitu_tem_alloc 
    632651#if defined key_diadct  
    633652      USE diadct    , ONLY: diadct_alloc  
    … …  
    646665      ierr = ierr + ldftra_oce_alloc()          ! ocean lateral  physics : tracers 
    647666      ierr = ierr + zdf_oce_alloc   ()          ! ocean vertical physics 
     667      ierr = ierr + insitu_tem_alloc() 
    648668      ! 
    649669      ierr = ierr + trc_oce_alloc   ()          ! shared TRC / TRA arrays 

  • branches/UKMO/CO6_shelfclimate/NEMOGCM/NEMO/OPA_SRC/step.F90

    r7566 r7567  
    255255                             CALL tra_sbc    ( kstp )       ! surface boundary condition 
    256256      IF( ln_traqsr      )   CALL tra_qsr    ( kstp )       ! penetrative solar radiation qsr 
     257      IF( ln_tradwl      )   CALL tra_dwl    ( kstp )       ! Polcoms Style Short Wave Radiation  
    257258      IF( ln_trabbc      )   CALL tra_bbc    ( kstp )       ! bottom heat flux 
    258259      IF( lk_trabbl      )   CALL tra_bbl    ( kstp )       ! advective (and/or diffusive) bottom boundary layer scheme 

  • branches/UKMO/CO6_shelfclimate/NEMOGCM/NEMO/OPA_SRC/step_oce.F90

    r7566 r7567  
    2525   USE sbcrnf           ! surface boundary condition: runoff variables 
    2626   USE sbccpl           ! surface boundary condition: coupled formulation (call send at end of step) 
     27   USE sbcflx           ! surface boundary condition: Fluxes 
    2728   USE sbc_oce          ! surface boundary condition: ocean 
    2829   USE sbctide          ! Tide initialisation 
    … …  
    3031 
    3132   USE traqsr           ! solar radiation penetration      (tra_qsr routine) 
     33   USE tradwl           ! POLCOMS style solar radiation    (tra_dwl routine)  
    3234   USE trasbc           ! surface boundary condition       (tra_sbc routine) 
    3335   USE trabbc           ! bottom boundary condition        (tra_bbc routine) 

  • branches/UKMO/CO6_shelfclimate/NEMOGCM/NEMO/OPA_SRC/trc_oce.F90

    r7566 r7567  
    2727   REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::   etot3         !: light absortion coefficient 
    2828   REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::   facvol        !: volume for degraded regions 
     29   REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:)   ::   rlambda2      !: Lambda2 for downwell version of Short wave Radiation 
     30   REAL(wp), PUBLIC                                      ::   rlambda       !: Lambda  for downwell version of Short wave Radiation 
    2931 
    3032#if defined key_top  
    … …  
    7880      !!                  ***  trc_oce_alloc  *** 
    7981      !!---------------------------------------------------------------------- 
    80       INTEGER ::   ierr(2)        ! Local variables 
     82      INTEGER ::   ierr(3)        ! Local variables 
    8183      !!---------------------------------------------------------------------- 
    8284      ierr(:) = 0 
    8385                     ALLOCATE( etot3 (jpi,jpj,jpk), STAT=ierr(1) ) 
    8486      IF( lk_degrad) ALLOCATE( facvol(jpi,jpj,jpk), STAT=ierr(2) ) 
     87                    ALLOCATE( rlambda2(jpi,jpj),   STAT=ierr(3) ) 
    8588      trc_oce_alloc  = MAXVAL( ierr ) 
    8689      ! 
    8790      IF( trc_oce_alloc /= 0 )   CALL ctl_warn('trc_oce_alloc: failed to allocate etot3 array') 
     91      IF( trc_oce_alloc /= 0 )   CALL ctl_warn('trc_oce_alloc: failed to allocate etot3, facvol or rlambda2 array') 
    8892   END FUNCTION trc_oce_alloc 
    8993 

  • branches/UKMO/CO6_shelfclimate/NEMOGCM/NEMO/TOP_SRC/MY_TRC/trcsms_my_trc.F90

    r7566 r7567  
    1818   USE trd_oce 
    1919   USE trdtrc 
     20   USE trcbc, only : trc_bc_read 
    2021 
    2122   IMPLICIT NONE 
    … …  
    5556 
    5657      IF( l_trdtrc )  CALL wrk_alloc( jpi, jpj, jpk, ztrmyt ) 
     58 
     59      CALL trc_bc_read  ( kt )       ! tracers: surface and lateral Boundary Conditions 
    5760 
    5861      IF( l_trdtrc ) THEN      ! Save the trends in the ixed layer 

  • branches/UKMO/CO6_shelfclimate/NEMOGCM/NEMO/TOP_SRC/MY_TRC/trcwri_my_trc.F90

    r7566 r7567  
    1919 
    2020   PUBLIC trc_wri_my_trc  
     21#if defined key_tracer_budget 
     22   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:), SAVE :: trb_temp ! slwa 
     23#endif 
     24 
    2125 
    2226#  include "top_substitute.h90" 
    2327CONTAINS 
    2428 
     29#if defined key_tracer_budget 
     30   SUBROUTINE trc_wri_my_trc (kt, fl) ! slwa 
     31#else 
    2532   SUBROUTINE trc_wri_my_trc 
     33#endif 
    2634      !!--------------------------------------------------------------------- 
    2735      !!                     ***  ROUTINE trc_wri_trc  *** 
    … …  
    2937      !! ** Purpose :   output passive tracers fields  
    3038      !!--------------------------------------------------------------------- 
     39#if defined key_tracer_budget 
     40      INTEGER, INTENT( in ), OPTIONAL     :: fl  
     41      INTEGER, INTENT( in )               :: kt 
     42      REAL(wp), DIMENSION(jpi,jpj,jpk)    :: trpool !tracer pool temporary output 
     43#endif 
    3144      CHARACTER (len=20)   :: cltra 
    32       INTEGER              :: jn 
     45      INTEGER              :: jn,jk 
    3346      !!--------------------------------------------------------------------- 
    3447  
    3548      ! write the tracer concentrations in the file 
    3649      ! --------------------------------------- 
     50 
     51 
     52#if defined key_tracer_budget 
     53      IF( PRESENT(fl)) THEN 
     54! depth integrated 
     55! for strict budgetting write this out at end of timestep as an average between 'now' and 'after' at kt 
     56         DO jn = jp_myt0, jp_myt1  
     57            trpool(:,:,:) = 0.5 * ( trn(:,:,:,jn) * fse3t_a(:,:,:) +  & 
     58                                        trb_temp(:,:,:,jn) * fse3t(:,:,:) ) 
     59! 
     60            cltra = TRIM( ctrcnm(jn) )                  ! output of tracer density  
     61            CALL iom_put( cltra, trpool(:,:,:) / (0.5* (fse3t(:,:,:) + fse3t_a(:,:,:) ) ) ) 
     62! 
     63            cltra = TRIM( ctrcnm(jn) )//"_pool"     ! volume integrated output 
     64            DO jk = 1, jpk 
     65               trpool(:,:,jk) = trpool(:,:,jk) * e1t(:,:) * e2t(:,:) 
     66            END DO 
     67            CALL iom_put( cltra, trpool) 
     68 
     69!           cltra = TRIM( ctrcnm(jn) )//"_pool"     ! volume integrated output 
     70!           DO jk = 1, jpk 
     71!              trpool(:,:,jk) = 0.5 * ( trn(:,:,jk,jn) * fse3t_a(:,:,jk) +  &  
     72!                                       trb_temp(:,:,jk,jn) * fse3t(:,:,jk) ) * &  
     73!                                       e1t(:,:) * e2t(:,:) 
     74!           END DO 
     75!           CALL iom_put( cltra, trpool) 
     76!           cltra = TRIM( ctrcnm(jn) )                  ! output of tracer density  
     77!           CALL iom_put( cltra, trpool(:,:,:) / (0.5* (fse3t(:,:,:) + fse3t_a(:,:,:) ) ) ) 
     78         END DO 
     79         CALL iom_put( "DEPTH" , 0.5* (fse3t(:,:,:) + fse3t_a(:,:,:) ) )  !  equivalent 'depth' at same time as tracer pool output 
     80      ELSE 
     81 
     82         IF( kt == nittrc000 ) THEN 
     83           ALLOCATE(trb_temp(jpi,jpj,jpk,jptra))  ! slwa 
     84         ENDIF 
     85         trb_temp(:,:,:,:)=trn(:,:,:,:) ! slwa save for tracer budget (unfiltered trn) 
     86 
     87!        DO jn = jp_myt0, jp_myt1 
     88!           cltra = TRIM( ctrcnm(jn) )                  ! short title for tracer 
     89!           CALL iom_put( cltra, trn(:,:,:,jn) )  
     90!        END DO 
     91! write out depths and areas in double precision for tracer budget calculations 
     92         CALL iom_put( "AREA" , e1t(:,:) * e2t(:,:)) 
     93!        CALL iom_put( "DEPTH" , fse3t(:,:,:) )  ! need depth at same time as tracer output 
     94 
     95      END IF 
     96#else 
    3797      DO jn = jp_myt0, jp_myt1 
    3898         cltra = TRIM( ctrcnm(jn) )                  ! short title for tracer 
    3999         CALL iom_put( cltra, trn(:,:,:,jn) ) 
    40100      END DO 
     101#endif 
    41102      ! 
    42103   END SUBROUTINE trc_wri_my_trc 

  • branches/UKMO/CO6_shelfclimate/NEMOGCM/NEMO/TOP_SRC/TRP/trcldf.F90

    r7566 r7567  
    5656      INTEGER, INTENT( in ) ::   kt   ! ocean time-step index 
    5757      !! 
    58       INTEGER            :: jn 
     58      INTEGER            :: jn, jk 
    5959      CHARACTER (len=22) :: charout 
    6060      REAL(wp), POINTER, DIMENSION(:,:,:,:) ::   ztrtrd 
    … …  
    105105        DO jn = 1, jptra 
    106106           ztrtrd(:,:,:,jn) = tra(:,:,:,jn) - ztrtrd(:,:,:,jn) 
     107#if defined key_tracer_budget 
     108           DO jk = 1, jpkm1 
     109             ztrtrd(:,:,jk,jn) = ztrtrd(:,:,jk,jn) * e1t(:,:) * e2t(:,:) * fse3t(:,:,jk)  ! slwa 
     110           END DO 
     111#endif 
    107112           CALL trd_tra( kt, 'TRC', jn, jptra_ldf, ztrtrd(:,:,:,jn) ) 
    108113        END DO 

  • branches/UKMO/CO6_shelfclimate/NEMOGCM/NEMO/TOP_SRC/TRP/trcnxt.F90

    r7566 r7567  
    3333   USE trdtra 
    3434   USE tranxt 
     35   USE trcbdy          ! BDY open boundaries 
     36   USE bdy_par, only: lk_bdy 
     37   USE iom 
    3538# if defined key_agrif 
    3639   USE agrif_top_interp 
    … …  
    9396      CHARACTER (len=22) :: charout 
    9497      REAL(wp), POINTER, DIMENSION(:,:,:,:) ::  ztrdt  
     98#if defined key_tracer_budget 
     99      REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::  ztrdt_m1 ! slwa 
     100#endif 
    95101      !!---------------------------------------------------------------------- 
    96102      ! 
    … …  
    101107         WRITE(numout,*) 'trc_nxt : time stepping on passive tracers' 
    102108      ENDIF 
     109#if defined key_tracer_budget 
     110      IF( kt == nittrc000 .AND. l_trdtrc ) THEN 
     111         ALLOCATE( ztrdt_m1(jpi,jpj,jpk,jptra) )  ! slwa 
     112         IF( ln_rsttr .AND.    &                     ! Restart: read in restart  file 
     113            iom_varid( numrtr, 'atf_trend_'//TRIM(ctrcnm(1)), ldstop = .FALSE. ) > 0 ) THEN 
     114            IF(lwp) WRITE(numout,*) '          nittrc000-nn_dttrc ATF tracer trend read in the restart file' 
     115            DO jn = 1, jptra 
     116               CALL iom_get( numrtr, jpdom_autoglo, 'atf_trend_'//TRIM(ctrcnm(jn)), ztrdt_m1(:,:,:,jn) )   ! before tracer trend for atf 
     117            END DO 
     118         ELSE           
     119           ztrdt_m1=0.0 
     120         ENDIF 
     121      ENDIF 
     122#endif 
    103123 
    104124#if defined key_agrif 
    … …  
    111131 
    112132 
    113 #if defined key_bdy 
    114 !!      CALL bdy_trc( kt )               ! BDY open boundaries 
    115 #endif 
     133      IF( lk_bdy )  CALL trc_bdy( kt )               ! BDY open boundaries 
    116134 
    117135 
    … …  
    149167               zfact = 1.e0 / r2dt(jk)   
    150168               ztrdt(:,:,jk,jn) = ( trb(:,:,jk,jn) - ztrdt(:,:,jk,jn) ) * zfact  
    151                CALL trd_tra( kt, 'TRC', jn, jptra_atf, ztrdt ) 
     169!slwa          CALL trd_tra( kt, 'TRC', jn, jptra_atf, ztrdt ) 
     170#if defined key_tracer_budget 
     171               ztrdt(:,:,jk,jn) = ztrdt(:,:,jk,jn) * e1t(:,:) * e2t(:,:) * e3t_n(:,:,jk)  ! slwa vvl 
     172               !ztrdt(:,:,jk,jn) = ztrdt(:,:,jk,jn) * e1t(:,:) * e2t(:,:) * e3t_0(:,:,jk)  ! slwa CHANGE for vvl 
     173#endif 
    152174            END DO 
     175#if defined key_tracer_budget 
     176! slwa budget code 
     177              CALL trd_tra( kt, 'TRC', jn, jptra_atf, ztrdt_m1(:,:,:,jn) ) 
     178#else 
     179              CALL trd_tra( kt, 'TRC', jn, jptra_atf, ztrdt(:,:,:,jn) ) 
     180#endif 
    153181         END DO 
     182#if defined key_tracer_budget 
     183        ztrdt_m1(:,:,:,:) = ztrdt(:,:,:,:)    ! need previous time step for budget slwa 
     184#endif 
    154185         CALL wrk_dealloc( jpi, jpj, jpk, jptra, ztrdt )  
    155186      END IF 
     187 
     188#if defined key_tracer_budget 
     189      !                                           Write in the tracer restart file 
     190      !                                          ******************************* 
     191      IF( lrst_trc ) THEN 
     192         IF(lwp) WRITE(numout,*) 
     193         IF(lwp) WRITE(numout,*) 'trc : ATF trend at last time step for tracer budget written in tracer restart file ',   & 
     194            &                    'at it= ', kt,' date= ', ndastp 
     195         IF(lwp) WRITE(numout,*) '~~~~' 
     196         DO jn = 1, jptra 
     197            CALL iom_rstput( kt, nitrst, numrtw, 'atf_trend_'//TRIM(ctrcnm(jn)), ztrdt_m1(:,:,:,jn) ) 
     198         END DO 
     199      ENDIF 
     200#endif 
     201 
    156202      ! 
    157203      IF(ln_ctl)   THEN  ! print mean trends (used for debugging) 

  • branches/UKMO/CO6_shelfclimate/NEMOGCM/NEMO/TOP_SRC/TRP/trcrad.F90

    r7566 r7567  
    1818   USE trdtra 
    1919   USE prtctl_trc          ! Print control for debbuging 
     20#if defined key_tracer_budget 
     21   USE iom 
     22#endif 
    2023 
    2124   IMPLICIT NONE 
    … …  
    110113      REAL(wp) :: zcoef, ztrcorn, ztrmasn   !    "         " 
    111114      REAL(wp), POINTER, DIMENSION(:,:,:) ::   ztrtrdb, ztrtrdn   ! workspace arrays 
     115#if defined key_tracer_budget 
     116      REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::  ztrtrdb_m1 ! slwa  
     117#endif 
    112118      REAL(wp) :: zs2rdt 
    113119      LOGICAL ::   lldebug = .FALSE. 
    … …  
    116122  
    117123      IF( l_trdtrc )  CALL wrk_alloc( jpi, jpj, jpk, ztrtrdb, ztrtrdn ) 
     124#if defined key_tracer_budget 
     125      IF( kt == nittrc000 .AND. l_trdtrc) THEN 
     126         ALLOCATE( ztrtrdb_m1(jpi,jpj,jpk,jptra) )  ! slwa 
     127         IF( ln_rsttr .AND.    &                     ! Restart: read in restart  file 
     128            iom_varid( numrtr, 'rdb_trend_'//TRIM(ctrcnm(1)), ldstop = .FALSE. ) > 0 ) THEN 
     129            IF(lwp) WRITE(numout,*) '          nittrc000-nn_dttrc RDB tracer trend read in the restart file' 
     130            DO jn = 1, jptra 
     131               CALL iom_get( numrtr, jpdom_autoglo, 'rdb_trend_'//TRIM(ctrcnm(jn)), ztrtrdb_m1(:,:,:,jn) )   ! before tracer trend for rdb 
     132            END DO 
     133         ELSE 
     134           ztrtrdb_m1=0.0 
     135         ENDIF 
     136      ENDIF 
     137#endif 
    118138       
    119139      IF( PRESENT( cpreserv )  ) THEN   !  total tracer concentration is preserved  
    … …  
    156176               ztrtrdb(:,:,:) = ( ptrb(:,:,:,jn) - ztrtrdb(:,:,:) ) * zs2rdt 
    157177               ztrtrdn(:,:,:) = ( ptrn(:,:,:,jn) - ztrtrdn(:,:,:) ) * zs2rdt  
     178#if defined key_tracer_budget 
     179! slwa budget code 
     180               DO jk = 1, jpkm1 
     181                  ztrtrdb(:,:,jk) = ztrtrdb(:,:,jk) * e1t(:,:) * e2t(:,:) * fse3t(:,:,jk) 
     182                  ztrtrdn(:,:,jk) = ztrtrdn(:,:,jk) * e1t(:,:) * e2t(:,:) * fse3t(:,:,jk) 
     183               END DO 
     184               CALL trd_tra( kt, 'TRC', jn, jptra_radb, ztrtrdb_m1(:,:,:,jn) ) 
     185               ztrtrdb_m1(:,:,:,jn)=ztrtrdb(:,:,:) 
     186#else 
    158187               CALL trd_tra( kt, 'TRC', jn, jptra_radb, ztrtrdb )       ! Asselin-like trend handling 
     188#endif 
    159189               CALL trd_tra( kt, 'TRC', jn, jptra_radn, ztrtrdn )       ! standard     trend handling 
    160190              ! 
    … …  
    187217               ztrtrdb(:,:,:) = ( ptrb(:,:,:,jn) - ztrtrdb(:,:,:) ) * zs2rdt 
    188218               ztrtrdn(:,:,:) = ( ptrn(:,:,:,jn) - ztrtrdn(:,:,:) ) * zs2rdt  
     219#if defined key_tracer_budget 
     220! slwa budget code 
     221               DO jk = 1, jpkm1 
     222                  ztrtrdb(:,:,jk) = ztrtrdb(:,:,jk) * e1t(:,:) * e2t(:,:) * fse3t(:,:,jk) 
     223                  ztrtrdn(:,:,jk) = ztrtrdn(:,:,jk) * e1t(:,:) * e2t(:,:) * fse3t(:,:,jk) 
     224               END DO 
     225               CALL trd_tra( kt, 'TRC', jn, jptra_radb, ztrtrdb_m1(:,:,:,jn) ) 
     226               ztrtrdb_m1(:,:,:,jn)=ztrtrdb(:,:,:) 
     227#else 
    189228               CALL trd_tra( kt, 'TRC', jn, jptra_radb, ztrtrdb )       ! Asselin-like trend handling 
     229#endif 
    190230               CALL trd_tra( kt, 'TRC', jn, jptra_radn, ztrtrdn )       ! standard     trend handling 
    191231              ! 
    … …  
    195235 
    196236      ENDIF 
     237 
     238#if defined key_tracer_budget 
     239      !                                           Write in the tracer restart file 
     240      !                                          ******************************* 
     241      IF( lrst_trc ) THEN 
     242         IF(lwp) WRITE(numout,*) 
     243         IF(lwp) WRITE(numout,*) 'trc : RDB trend at last time step for tracer budget written in tracer restart file ',   & 
     244            &                    'at it= ', kt,' date= ', ndastp 
     245         IF(lwp) WRITE(numout,*) '~~~~' 
     246         DO jn = 1, jptra 
     247            CALL iom_rstput( kt, nitrst, numrtw, 'rdb_trend_'//TRIM(ctrcnm(jn)), ztrtrdb_m1(:,:,:,jn) ) 
     248         END DO 
     249      ENDIF 
     250#endif 
    197251 
    198252      IF( l_trdtrc )  CALL wrk_dealloc( jpi, jpj, jpk, ztrtrdb, ztrtrdn ) 

  • branches/UKMO/CO6_shelfclimate/NEMOGCM/NEMO/TOP_SRC/TRP/trcsbc.F90

    r7566 r7567  
    113113           sbc_trc_b(:,:,:) = 0._wp 
    114114         ENDIF 
     115         sbc_trc(:,:,:) = 0._wp    !slwa initialise for vvl 
    115116      ELSE                                         ! Swap of forcing fields 
    116117         IF( ln_top_euler ) THEN 

  • branches/UKMO/CO6_shelfclimate/NEMOGCM/NEMO/TOP_SRC/TRP/trctrp.F90

    r7566 r7567  
    2727   USE trcsbc          ! surface boundary condition          (trc_sbc routine) 
    2828   USE zpshde          ! partial step: hor. derivative       (zps_hde routine) 
     29   USE trcbdy          ! BDY open boundaries 
     30   USE bdy_par, only: lk_bdy 
    2931 
    3032#if defined key_agrif 
    … …  
    6870         IF( ln_trcdmp )        CALL trc_dmp( kstp )            ! internal damping trends 
    6971         IF( ln_trcdmp_clo )    CALL trc_dmp_clo( kstp )        ! internal damping trends on closed seas only 
     72         IF( lk_bdy )           CALL trc_bdy_dmp( kstp )        ! BDY damping trends 
    7073                                CALL trc_adv( kstp )            ! horizontal & vertical advection  
    7174                                CALL trc_ldf( kstp )            ! lateral mixing 

  • branches/UKMO/CO6_shelfclimate/NEMOGCM/NEMO/TOP_SRC/trc.F90

    r7566 r7567  
    1414   USE par_oce 
    1515   USE par_trc 
     16#if defined key_bdy 
     17   USE bdy_oce, only: nb_bdy, OBC_DATA 
     18#endif 
    1619    
    1720   IMPLICIT NONE 
    … …  
    9194       CHARACTER(len = 20)  :: clunit   !: unit 
    9295       LOGICAL              :: llinit   !: read in a file or not 
     96#if defined  key_my_trc 
     97       LOGICAL              :: llsbc   !: read in a file or not 
     98       LOGICAL              :: llcbc   !: read in a file or not 
     99       LOGICAL              :: llobc   !: read in a file or not 
     100#endif 
    93101       LOGICAL              :: llsave   !: save the tracer or not 
    94102   END TYPE PTRACER 
    … …  
    191199# endif 
    192200   ! 
     201#if defined key_bdy 
     202   CHARACTER(len=20), PUBLIC, ALLOCATABLE,  SAVE,  DIMENSION(:)   ::  cn_trc_dflt          ! Default OBC condition for all tracers 
     203   CHARACTER(len=20), PUBLIC, ALLOCATABLE,  SAVE,  DIMENSION(:)   ::  cn_trc               ! Choice of boundary condition for tracers 
     204   INTEGER,           PUBLIC, ALLOCATABLE,  SAVE,  DIMENSION(:)   ::  nn_trcdmp_bdy        !: =T Tracer damping 
     205   ! External data structure of BDY for TOP. Available elements: cn_obc, ll_trc, trcnow, dmp 
     206   TYPE(OBC_DATA),    PUBLIC, ALLOCATABLE, DIMENSION(:,:), TARGET ::  trcdta_bdy           !: bdy external data (local process) 
     207#endif 
    193208 
    194209   !!---------------------------------------------------------------------- 
    … …  
    213228         &      cvol(jpi,jpj,jpk)     , rdttrc(jpk)           , trai(jptra)           ,       & 
    214229         &      ctrcnm(jptra)         , ctrcln(jptra)         , ctrcun(jptra)         ,       &  
     230#if defined key_my_trc 
     231         &      ln_trc_sbc(jptra)     , ln_trc_cbc(jptra)     , ln_trc_obc(jptra)     ,       & 
     232#endif 
     233#if defined key_bdy 
     234         &      cn_trc_dflt(nb_bdy)   , cn_trc(nb_bdy)        , nn_trcdmp_bdy(nb_bdy) ,       & 
     235         &      trcdta_bdy(jptra,nb_bdy)                                              ,       & 
     236#endif 
    215237         &      ln_trc_ini(jptra)     , ln_trc_wri(jptra)     , qsr_mean(jpi,jpj)     ,  STAT = trc_alloc  )   
    216238 

  • branches/UKMO/CO6_shelfclimate/NEMOGCM/NEMO/TOP_SRC/trcbc.F90

    r7566 r7567  
    44   !! TOP :  module for passive tracer boundary conditions 
    55   !!===================================================================== 
     6   !! History :  3.5 !  2014-04  (M. Vichi, T. Lovato)  Original 
     7   !!            3.6 !  2015-03  (T . Lovato) Revision and BDY support 
    68   !!---------------------------------------------------------------------- 
    79#if  defined key_top  
    … …  
    911   !!   'key_top'                                                TOP model  
    1012   !!---------------------------------------------------------------------- 
    11    !!   trc_dta    : read and time interpolated passive tracer data 
     13   !!   trc_bc       : read and time interpolated tracer Boundary Conditions 
    1214   !!---------------------------------------------------------------------- 
    1315   USE par_trc       !  passive tracers parameters 
    … …  
    1719   USE lib_mpp       !  MPP library 
    1820   USE fldread       !  read input fields 
     21#if defined key_bdy 
     22   USE bdy_oce, only: nb_bdy , idx_bdy, ln_coords_file, rn_time_dmp, rn_time_dmp_out 
     23#endif 
    1924 
    2025   IMPLICIT NONE 
    … …  
    3035   INTEGER  , SAVE, PUBLIC, ALLOCATABLE, DIMENSION(:)  :: n_trc_indsbc ! index of tracer with SBC data 
    3136   INTEGER  , SAVE, PUBLIC, ALLOCATABLE, DIMENSION(:)  :: n_trc_indcbc ! index of tracer with CBC data 
    32    INTEGER  , SAVE, PUBLIC                             :: ntra_obc     ! MAX( 1, nb_trcxxx ) to avoid compilation error with bounds checking 
    33    INTEGER  , SAVE, PUBLIC                             :: ntra_sbc     ! MAX( 1, nb_trcxxx ) to avoid compilation error with bounds checking 
    34    INTEGER  , SAVE, PUBLIC                             :: ntra_cbc     ! MAX( 1, nb_trcxxx ) to avoid compilation error with bounds checking 
    35    REAL(wp) , SAVE, PUBLIC, ALLOCATABLE, DIMENSION(:)  :: rf_trofac   ! multiplicative factor for OBCtracer values 
    36    TYPE(FLD), SAVE, PUBLIC, ALLOCATABLE, DIMENSION(:)  :: sf_trcobc   ! structure of data input OBC (file informations, fields read) 
    3737   REAL(wp) , SAVE, PUBLIC, ALLOCATABLE, DIMENSION(:)  :: rf_trsfac   ! multiplicative factor for SBC tracer values 
    3838   TYPE(FLD), SAVE, PUBLIC, ALLOCATABLE, DIMENSION(:)  :: sf_trcsbc   ! structure of data input SBC (file informations, fields read) 
    3939   REAL(wp) , SAVE, PUBLIC, ALLOCATABLE, DIMENSION(:)  :: rf_trcfac   ! multiplicative factor for CBC tracer values 
    4040   TYPE(FLD), SAVE, PUBLIC, ALLOCATABLE, DIMENSION(:)  :: sf_trccbc   ! structure of data input CBC (file informations, fields read) 
     41   REAL(wp) , SAVE, PUBLIC, ALLOCATABLE, DIMENSION(:,:)  :: rf_trofac    ! multiplicative factor for OBCtracer values 
     42   TYPE(FLD), SAVE, PUBLIC, ALLOCATABLE, DIMENSION(:,:), TARGET  :: sf_trcobc    ! structure of data input OBC (file informations, fields read) 
     43   TYPE(MAP_POINTER), ALLOCATABLE, DIMENSION(:,:) :: nbmap_ptr   ! array of pointers to nbmap 
    4144 
    4245   !! * Substitutions 
    4346#  include "domzgr_substitute.h90" 
    4447   !!---------------------------------------------------------------------- 
    45    !! NEMO/OPA 3.3 , NEMO Consortium (2010) 
     48   !! NEMO/OPA 3.6 , NEMO Consortium (2015) 
    4649   !! $Id$ 
    4750   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
    … …  
    6063      ! 
    6164      INTEGER,INTENT(IN) :: ntrc                           ! number of tracers 
    62       INTEGER            :: jl, jn                         ! dummy loop indices 
     65      INTEGER            :: jl, jn , ib, ibd, ii, ij, ik   ! dummy loop indices 
    6366      INTEGER            :: ierr0, ierr1, ierr2, ierr3     ! temporary integers 
    6467      INTEGER            ::  ios                           ! Local integer output status for namelist read 
     68      INTEGER            :: nblen, igrd                    ! support arrays for BDY 
    6569      CHARACTER(len=100) :: clndta, clntrc 
    6670      ! 
    67       CHARACTER(len=100) :: cn_dir 
     71      CHARACTER(len=100) :: cn_dir_sbc, cn_dir_cbc, cn_dir_obc 
     72 
    6873      TYPE(FLD_N), ALLOCATABLE, DIMENSION(:) :: slf_i  ! local array of namelist informations on the fields to read 
    69       TYPE(FLD_N), DIMENSION(jpmaxtrc) :: sn_trcobc    ! open 
     74      TYPE(FLD_N), DIMENSION(jpmaxtrc,2) :: sn_trcobc  ! open 
     75      TYPE(FLD_N), DIMENSION(jpmaxtrc) :: sn_trcobc2   ! to read in multiple (2) open bdy 
    7076      TYPE(FLD_N), DIMENSION(jpmaxtrc) :: sn_trcsbc    ! surface 
    7177      TYPE(FLD_N), DIMENSION(jpmaxtrc) :: sn_trccbc    ! coastal 
    … …  
    7480      REAL(wp)   , DIMENSION(jpmaxtrc) :: rn_trcfac    ! multiplicative factor for tracer values 
    7581      !! 
    76       NAMELIST/namtrc_bc/ cn_dir, sn_trcobc, rn_trofac, sn_trcsbc, rn_trsfac, sn_trccbc, rn_trcfac  
     82      NAMELIST/namtrc_bc/ cn_dir_sbc, cn_dir_cbc, cn_dir_obc, sn_trcobc2, rn_trofac, sn_trcsbc, rn_trsfac, sn_trccbc, rn_trcfac 
     83#if defined key_bdy 
     84      NAMELIST/namtrc_bdy/ cn_trc_dflt, cn_trc, nn_trcdmp_bdy 
     85#endif 
    7786      !!---------------------------------------------------------------------- 
    7887      IF( nn_timing == 1 )  CALL timing_start('trc_bc_init') 
    7988      ! 
     89      IF( lwp ) THEN 
     90         WRITE(numout,*) ' ' 
     91         WRITE(numout,*) 'trc_bc_init : Tracers Boundary Conditions (BC)' 
     92         WRITE(numout,*) '~~~~~~~~~~~ ' 
     93      ENDIF 
    8094      !  Initialisation and local array allocation 
    8195      ierr0 = 0  ;  ierr1 = 0  ;  ierr2 = 0  ;  ierr3 = 0   
    … …  
    107121      n_trc_indcbc(:) = 0 
    108122      ! 
    109       DO jn = 1, ntrc 
    110          IF( ln_trc_obc(jn) ) THEN 
    111              nb_trcobc       = nb_trcobc + 1  
    112              n_trc_indobc(jn) = nb_trcobc  
    113          ENDIF 
    114          IF( ln_trc_sbc(jn) ) THEN 
    115              nb_trcsbc       = nb_trcsbc + 1 
    116              n_trc_indsbc(jn) = nb_trcsbc 
    117          ENDIF 
    118          IF( ln_trc_cbc(jn) ) THEN 
    119              nb_trccbc       = nb_trccbc + 1 
    120              n_trc_indcbc(jn) = nb_trccbc 
    121          ENDIF 
    122       ENDDO 
    123       ntra_obc = MAX( 1, nb_trcobc )   ! To avoid compilation error with bounds checking 
    124       IF( lwp ) WRITE(numout,*) ' ' 
    125       IF( lwp ) WRITE(numout,*) ' Number of passive tracers to be initialized with open boundary data :', nb_trcobc 
    126       IF( lwp ) WRITE(numout,*) ' ' 
    127       ntra_sbc = MAX( 1, nb_trcsbc )   ! To avoid compilation error with bounds checking 
    128       IF( lwp ) WRITE(numout,*) ' ' 
    129       IF( lwp ) WRITE(numout,*) ' Number of passive tracers to be initialized with surface boundary data :', nb_trcsbc 
    130       IF( lwp ) WRITE(numout,*) ' ' 
    131       ntra_cbc = MAX( 1, nb_trccbc )   ! To avoid compilation error with bounds checking 
    132       IF( lwp ) WRITE(numout,*) ' ' 
    133       IF( lwp ) WRITE(numout,*) ' Number of passive tracers to be initialized with coastal boundary data :', nb_trccbc 
    134       IF( lwp ) WRITE(numout,*) ' ' 
    135  
     123      ! Read Boundary Conditions Namelists 
    136124      REWIND( numnat_ref )              ! Namelist namtrc_bc in reference namelist : Passive tracer data structure 
    137125      READ  ( numnat_ref, namtrc_bc, IOSTAT = ios, ERR = 901) 
    … …  
    139127 
    140128      REWIND( numnat_cfg )              ! Namelist namtrc_bc in configuration namelist : Passive tracer data structure 
     129#if defined key_bdy 
     130      DO ib = 1, nb_bdy 
     131#endif 
    141132      READ  ( numnat_cfg, namtrc_bc, IOSTAT = ios, ERR = 902 ) 
    142133902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namtrc_bc in configuration namelist', lwp ) 
    143134      IF(lwm) WRITE ( numont, namtrc_bc ) 
    144  
    145       ! print some information for each  
     135#if defined key_bdy 
     136        sn_trcobc(:,ib)=sn_trcobc2(:) 
     137<