Changeset 12377 for NEMO/trunk/src/OCE


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

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

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

svn merge —ignore-ancestry \

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

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

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

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

Location:
NEMO/trunk
Files:
17 deleted
164 edited
19 copied

Legend:

Unmodified
Added
Removed
  • NEMO/trunk

    • Property svn:externals
      •  

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

    r10425 r12377  
    5252CONTAINS 
    5353 
    54    SUBROUTINE asm_bkg_wri( kt ) 
     54   SUBROUTINE asm_bkg_wri( kt, Kmm ) 
    5555      !!----------------------------------------------------------------------- 
    5656      !!                  ***  ROUTINE asm_bkg_wri *** 
     
    6565      !!----------------------------------------------------------------------- 
    6666      INTEGER, INTENT( IN ) :: kt               ! Current time-step 
     67      INTEGER, INTENT( IN ) :: Kmm              ! time level index 
    6768      ! 
    6869      CHARACTER (LEN=50) :: cl_asmbkg 
     
    9899            ! 
    99100            !                                      ! Write the information 
    100             CALL iom_rstput( kt, nitbkg_r, inum, 'rdastp' , zdate             ) 
    101             CALL iom_rstput( kt, nitbkg_r, inum, 'un'     , un                ) 
    102             CALL iom_rstput( kt, nitbkg_r, inum, 'vn'     , vn                ) 
    103             CALL iom_rstput( kt, nitbkg_r, inum, 'tn'     , tsn(:,:,:,jp_tem) ) 
    104             CALL iom_rstput( kt, nitbkg_r, inum, 'sn'     , tsn(:,:,:,jp_sal) ) 
    105             CALL iom_rstput( kt, nitbkg_r, inum, 'sshn'   , sshn              ) 
    106             IF( ln_zdftke )   CALL iom_rstput( kt, nitbkg_r, inum, 'en'     , en                ) 
     101            CALL iom_rstput( kt, nitbkg_r, inum, 'rdastp' , zdate                ) 
     102            CALL iom_rstput( kt, nitbkg_r, inum, 'un'     , uu(:,:,:,Kmm)        ) 
     103            CALL iom_rstput( kt, nitbkg_r, inum, 'vn'     , vv(:,:,:,Kmm)        ) 
     104            CALL iom_rstput( kt, nitbkg_r, inum, 'tn'     , ts(:,:,:,jp_tem,Kmm) ) 
     105            CALL iom_rstput( kt, nitbkg_r, inum, 'sn'     , ts(:,:,:,jp_sal,Kmm) ) 
     106            CALL iom_rstput( kt, nitbkg_r, inum, 'sshn'   , ssh(:,:,Kmm)         ) 
     107            IF( ln_zdftke )   CALL iom_rstput( kt, nitbkg_r, inum, 'en'     , en ) 
    107108            ! 
    108109            CALL iom_close( inum ) 
     
    133134            ! 
    134135            !                                      ! Write the information 
    135             CALL iom_rstput( kt, nitdin_r, inum, 'rdastp' , zdate             ) 
    136             CALL iom_rstput( kt, nitdin_r, inum, 'un'     , un                ) 
    137             CALL iom_rstput( kt, nitdin_r, inum, 'vn'     , vn                ) 
    138             CALL iom_rstput( kt, nitdin_r, inum, 'tn'     , tsn(:,:,:,jp_tem) ) 
    139             CALL iom_rstput( kt, nitdin_r, inum, 'sn'     , tsn(:,:,:,jp_sal) ) 
    140             CALL iom_rstput( kt, nitdin_r, inum, 'sshn'   , sshn              ) 
     136            CALL iom_rstput( kt, nitdin_r, inum, 'rdastp' , zdate                ) 
     137            CALL iom_rstput( kt, nitdin_r, inum, 'un'     , uu(:,:,:,Kmm)        ) 
     138            CALL iom_rstput( kt, nitdin_r, inum, 'vn'     , vv(:,:,:,Kmm)        ) 
     139            CALL iom_rstput( kt, nitdin_r, inum, 'tn'     , ts(:,:,:,jp_tem,Kmm) ) 
     140            CALL iom_rstput( kt, nitdin_r, inum, 'sn'     , ts(:,:,:,jp_sal,Kmm) ) 
     141            CALL iom_rstput( kt, nitdin_r, inum, 'sshn'   , ssh(:,:,Kmm)         ) 
    141142#if defined key_si3 
    142143            IF( nn_ice == 2 ) THEN 
  • NEMO/trunk/src/OCE/ASM/asminc.F90

    r11536 r12377  
    9494 
    9595   !! * Substitutions 
    96 #  include "vectopt_loop_substitute.h90" 
     96#  include "do_loop_substitute.h90" 
    9797   !!---------------------------------------------------------------------- 
    9898   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
     
    102102CONTAINS 
    103103 
    104    SUBROUTINE asm_inc_init 
     104   SUBROUTINE asm_inc_init( Kbb, Kmm, Krhs ) 
    105105      !!---------------------------------------------------------------------- 
    106106      !!                    ***  ROUTINE asm_inc_init  *** 
     
    112112      !! ** Action  :  
    113113      !!---------------------------------------------------------------------- 
     114      INTEGER, INTENT(in) ::  Kbb, Kmm, Krhs  ! time level indices 
     115      ! 
    114116      INTEGER :: ji, jj, jk, jt  ! dummy loop indices 
    115117      INTEGER :: imid, inum      ! local integers 
     
    145147      ln_temnofreeze = .FALSE. 
    146148 
    147       REWIND( numnam_ref )              ! Namelist nam_asminc in reference namelist : Assimilation increment 
    148149      READ  ( numnam_ref, nam_asminc, IOSTAT = ios, ERR = 901) 
    149150901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nam_asminc in reference namelist' ) 
    150       REWIND( numnam_cfg )              ! Namelist nam_asminc in configuration namelist : Assimilation increment 
    151151      READ  ( numnam_cfg, nam_asminc, IOSTAT = ios, ERR = 902 ) 
    152152902   IF( ios >  0 )   CALL ctl_nam ( ios , 'nam_asminc in configuration namelist' ) 
     
    413413            DO jk = 1, jpkm1           ! zhdiv = e1e1 * div 
    414414               zhdiv(:,:) = 0._wp 
    415                DO jj = 2, jpjm1 
    416                   DO ji = fs_2, fs_jpim1   ! vector opt. 
    417                      zhdiv(ji,jj) = (  e2u(ji  ,jj) * e3u_n(ji  ,jj,jk) * u_bkginc(ji  ,jj,jk)    & 
    418                         &            - e2u(ji-1,jj) * e3u_n(ji-1,jj,jk) * u_bkginc(ji-1,jj,jk)    & 
    419                         &            + e1v(ji,jj  ) * e3v_n(ji,jj  ,jk) * v_bkginc(ji,jj  ,jk)    & 
    420                         &            - e1v(ji,jj-1) * e3v_n(ji,jj-1,jk) * v_bkginc(ji,jj-1,jk)  ) / e3t_n(ji,jj,jk) 
    421                   END DO 
    422                END DO 
     415               DO_2D_00_00 
     416                  zhdiv(ji,jj) = (  e2u(ji  ,jj) * e3u(ji  ,jj,jk,Kmm) * u_bkginc(ji  ,jj,jk)    & 
     417                     &            - e2u(ji-1,jj) * e3u(ji-1,jj,jk,Kmm) * u_bkginc(ji-1,jj,jk)    & 
     418                     &            + e1v(ji,jj  ) * e3v(ji,jj  ,jk,Kmm) * v_bkginc(ji,jj  ,jk)    & 
     419                     &            - e1v(ji,jj-1) * e3v(ji,jj-1,jk,Kmm) * v_bkginc(ji,jj-1,jk)  ) / e3t(ji,jj,jk,Kmm) 
     420               END_2D 
    423421               CALL lbc_lnk( 'asminc', zhdiv, 'T', 1. )   ! lateral boundary cond. (no sign change) 
    424422               ! 
    425                DO jj = 2, jpjm1 
    426                   DO ji = fs_2, fs_jpim1   ! vector opt. 
    427                      u_bkginc(ji,jj,jk) = u_bkginc(ji,jj,jk)                         & 
    428                         &               + 0.2_wp * ( zhdiv(ji+1,jj) - zhdiv(ji  ,jj) ) * r1_e1u(ji,jj) * umask(ji,jj,jk) 
    429                      v_bkginc(ji,jj,jk) = v_bkginc(ji,jj,jk)                         & 
    430                         &               + 0.2_wp * ( zhdiv(ji,jj+1) - zhdiv(ji,jj  ) ) * r1_e2v(ji,jj) * vmask(ji,jj,jk)  
    431                   END DO 
    432                END DO 
     423               DO_2D_00_00 
     424                  u_bkginc(ji,jj,jk) = u_bkginc(ji,jj,jk)                         & 
     425                     &               + 0.2_wp * ( zhdiv(ji+1,jj) - zhdiv(ji  ,jj) ) * r1_e1u(ji,jj) * umask(ji,jj,jk) 
     426                  v_bkginc(ji,jj,jk) = v_bkginc(ji,jj,jk)                         & 
     427                     &               + 0.2_wp * ( zhdiv(ji,jj+1) - zhdiv(ji,jj  ) ) * r1_e2v(ji,jj) * vmask(ji,jj,jk)  
     428               END_2D 
    433429            END DO 
    434430            ! 
     
    494490      ! 
    495491      IF( lk_asminc ) THEN                            !==  data assimilation  ==! 
    496          IF( ln_bkgwri )   CALL asm_bkg_wri( nit000 - 1 )      ! Output background fields 
     492         IF( ln_bkgwri )   CALL asm_bkg_wri( nit000 - 1, Kmm )      ! Output background fields 
    497493         IF( ln_asmdin ) THEN                                  ! Direct initialization 
    498             IF( ln_trainc )   CALL tra_asm_inc( nit000 - 1 )      ! Tracers 
    499             IF( ln_dyninc )   CALL dyn_asm_inc( nit000 - 1 )      ! Dynamics 
    500             IF( ln_sshinc )   CALL ssh_asm_inc( nit000 - 1 )      ! SSH 
     494            IF( ln_trainc )   CALL tra_asm_inc( nit000 - 1, Kbb, Kmm, ts    , Krhs )      ! Tracers 
     495            IF( ln_dyninc )   CALL dyn_asm_inc( nit000 - 1, Kbb, Kmm, uu, vv, Krhs )      ! Dynamics 
     496            IF( ln_sshinc )   CALL ssh_asm_inc( nit000 - 1, Kbb, Kmm )                    ! SSH 
    501497         ENDIF 
    502498      ENDIF 
     
    505501    
    506502    
    507    SUBROUTINE tra_asm_inc( kt ) 
     503   SUBROUTINE tra_asm_inc( kt, Kbb, Kmm, pts, Krhs ) 
    508504      !!---------------------------------------------------------------------- 
    509505      !!                    ***  ROUTINE tra_asm_inc  *** 
     
    515511      !! ** Action  :  
    516512      !!---------------------------------------------------------------------- 
    517       INTEGER, INTENT(IN) ::   kt   ! Current time step 
     513      INTEGER                                  , INTENT(in   ) :: kt             ! Current time step 
     514      INTEGER                                  , INTENT(in   ) :: Kbb, Kmm, Krhs ! Time level indices 
     515      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts,jpt), INTENT(inout) :: pts            ! active tracers and RHS of tracer equation 
    518516      ! 
    519517      INTEGER  :: ji, jj, jk 
     
    526524      ! used to prevent the applied increments taking the temperature below the local freezing point  
    527525      DO jk = 1, jpkm1 
    528         CALL eos_fzp( tsn(:,:,jk,jp_sal), fzptnz(:,:,jk), gdept_n(:,:,jk) ) 
     526        CALL eos_fzp( pts(:,:,jk,jp_sal,Kmm), fzptnz(:,:,jk), gdept(:,:,jk,Kmm) ) 
    529527      END DO 
    530528         ! 
     
    549547                  ! Do not apply negative increments if the temperature will fall below freezing 
    550548                  WHERE(t_bkginc(:,:,jk) > 0.0_wp .OR. & 
    551                      &   tsn(:,:,jk,jp_tem) + tsa(:,:,jk,jp_tem) + t_bkginc(:,:,jk) * wgtiau(it) > fzptnz(:,:,jk) )  
    552                      tsa(:,:,jk,jp_tem) = tsa(:,:,jk,jp_tem) + t_bkginc(:,:,jk) * zincwgt   
     549                     &   pts(:,:,jk,jp_tem,Kmm) + pts(:,:,jk,jp_tem,Krhs) + t_bkginc(:,:,jk) * wgtiau(it) > fzptnz(:,:,jk) )  
     550                     pts(:,:,jk,jp_tem,Krhs) = pts(:,:,jk,jp_tem,Krhs) + t_bkginc(:,:,jk) * zincwgt   
    553551                  END WHERE 
    554552               ELSE 
    555                   tsa(:,:,jk,jp_tem) = tsa(:,:,jk,jp_tem) + t_bkginc(:,:,jk) * zincwgt   
     553                  pts(:,:,jk,jp_tem,Krhs) = pts(:,:,jk,jp_tem,Krhs) + t_bkginc(:,:,jk) * zincwgt   
    556554               ENDIF 
    557555               IF (ln_salfix) THEN 
     
    559557                  ! minimum value salfixmin 
    560558                  WHERE(s_bkginc(:,:,jk) > 0.0_wp .OR. & 
    561                      &   tsn(:,:,jk,jp_sal) + tsa(:,:,jk,jp_sal) + s_bkginc(:,:,jk) * wgtiau(it) > salfixmin )  
    562                      tsa(:,:,jk,jp_sal) = tsa(:,:,jk,jp_sal) + s_bkginc(:,:,jk) * zincwgt 
     559                     &   pts(:,:,jk,jp_sal,Kmm) + pts(:,:,jk,jp_sal,Krhs) + s_bkginc(:,:,jk) * wgtiau(it) > salfixmin )  
     560                     pts(:,:,jk,jp_sal,Krhs) = pts(:,:,jk,jp_sal,Krhs) + s_bkginc(:,:,jk) * zincwgt 
    563561                  END WHERE 
    564562               ELSE 
    565                   tsa(:,:,jk,jp_sal) = tsa(:,:,jk,jp_sal) + s_bkginc(:,:,jk) * zincwgt 
     563                  pts(:,:,jk,jp_sal,Krhs) = pts(:,:,jk,jp_sal,Krhs) + s_bkginc(:,:,jk) * zincwgt 
    566564               ENDIF 
    567565            END DO 
     
    584582            IF (ln_temnofreeze) THEN 
    585583               ! Do not apply negative increments if the temperature will fall below freezing 
    586                WHERE( t_bkginc(:,:,:) > 0.0_wp .OR. tsn(:,:,:,jp_tem) + t_bkginc(:,:,:) > fzptnz(:,:,:) )  
    587                   tsn(:,:,:,jp_tem) = t_bkg(:,:,:) + t_bkginc(:,:,:)    
     584               WHERE( t_bkginc(:,:,:) > 0.0_wp .OR. pts(:,:,:,jp_tem,Kmm) + t_bkginc(:,:,:) > fzptnz(:,:,:) )  
     585                  pts(:,:,:,jp_tem,Kmm) = t_bkg(:,:,:) + t_bkginc(:,:,:)    
    588586               END WHERE 
    589587            ELSE 
    590                tsn(:,:,:,jp_tem) = t_bkg(:,:,:) + t_bkginc(:,:,:)    
     588               pts(:,:,:,jp_tem,Kmm) = t_bkg(:,:,:) + t_bkginc(:,:,:)    
    591589            ENDIF 
    592590            IF (ln_salfix) THEN 
    593591               ! Do not apply negative increments if the salinity will fall below a specified 
    594592               ! minimum value salfixmin 
    595                WHERE( s_bkginc(:,:,:) > 0.0_wp .OR. tsn(:,:,:,jp_sal) + s_bkginc(:,:,:) > salfixmin )  
    596                   tsn(:,:,:,jp_sal) = s_bkg(:,:,:) + s_bkginc(:,:,:)    
     593               WHERE( s_bkginc(:,:,:) > 0.0_wp .OR. pts(:,:,:,jp_sal,Kmm) + s_bkginc(:,:,:) > salfixmin )  
     594                  pts(:,:,:,jp_sal,Kmm) = s_bkg(:,:,:) + s_bkginc(:,:,:)    
    597595               END WHERE 
    598596            ELSE 
    599                tsn(:,:,:,jp_sal) = s_bkg(:,:,:) + s_bkginc(:,:,:)    
     597               pts(:,:,:,jp_sal,Kmm) = s_bkg(:,:,:) + s_bkginc(:,:,:)    
    600598            ENDIF 
    601599 
    602             tsb(:,:,:,:) = tsn(:,:,:,:)                 ! Update before fields 
    603  
    604             CALL eos( tsb, rhd, rhop, gdept_0(:,:,:) )  ! Before potential and in situ densities 
     600            pts(:,:,:,:,Kbb) = pts(:,:,:,:,Kmm)                 ! Update before fields 
     601 
     602            CALL eos( pts(:,:,:,:,Kbb), rhd, rhop, gdept_0(:,:,:) )  ! Before potential and in situ densities 
    605603!!gm  fabien 
    606 !            CALL eos( tsb, rhd, rhop )                ! Before potential and in situ densities 
     604!            CALL eos( pts(:,:,:,:,Kbb), rhd, rhop )                ! Before potential and in situ densities 
    607605!!gm 
    608606 
    609             IF( ln_zps .AND. .NOT. lk_c1d .AND. .NOT. ln_isfcav)      & 
    610                &  CALL zps_hde    ( kt, jpts, tsb, gtsu, gtsv,        &  ! Partial steps: before horizontal gradient 
    611                &                              rhd, gru , grv          )  ! of t, s, rd at the last ocean level 
    612             IF( ln_zps .AND. .NOT. lk_c1d .AND.       ln_isfcav)      & 
    613                &  CALL zps_hde_isf( nit000, jpts, tsb, gtsu, gtsv, gtui, gtvi,    &    ! Partial steps for top cell (ISF) 
    614                &                                  rhd, gru , grv , grui, grvi       ! of t, s, rd at the last ocean level 
     607            IF( ln_zps .AND. .NOT. lk_c1d .AND. .NOT. ln_isfcav)           & 
     608               &  CALL zps_hde    ( kt, Kmm, jpts, pts(:,:,:,:,Kbb), gtsu, gtsv,        &  ! Partial steps: before horizontal gradient 
     609               &                              rhd, gru , grv               )  ! of t, s, rd at the last ocean level 
     610            IF( ln_zps .AND. .NOT. lk_c1d .AND.       ln_isfcav)                       & 
     611               &  CALL zps_hde_isf( nit000, Kmm, jpts, pts(:,:,:,:,Kbb), gtsu, gtsv, gtui, gtvi,    &  ! Partial steps for top cell (ISF) 
     612               &                                  rhd, gru , grv , grui, grvi          )  ! of t, s, rd at the last ocean level 
    615613 
    616614            DEALLOCATE( t_bkginc ) 
     
    627625 
    628626 
    629    SUBROUTINE dyn_asm_inc( kt ) 
     627   SUBROUTINE dyn_asm_inc( kt, Kbb, Kmm, puu, pvv, Krhs ) 
    630628      !!---------------------------------------------------------------------- 
    631629      !!                    ***  ROUTINE dyn_asm_inc  *** 
     
    637635      !! ** Action  :  
    638636      !!---------------------------------------------------------------------- 
    639       INTEGER, INTENT(IN) :: kt   ! Current time step 
     637      INTEGER                             , INTENT( in )  ::  kt             ! ocean time-step index 
     638      INTEGER                             , INTENT( in )  ::  Kbb, Kmm, Krhs ! ocean time level indices 
     639      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) ::  puu, pvv       ! ocean velocities and RHS of momentum equation 
    640640      ! 
    641641      INTEGER :: jk 
     
    661661            ! Update the dynamic tendencies 
    662662            DO jk = 1, jpkm1 
    663                ua(:,:,jk) = ua(:,:,jk) + u_bkginc(:,:,jk) * zincwgt 
    664                va(:,:,jk) = va(:,:,jk) + v_bkginc(:,:,jk) * zincwgt 
     663               puu(:,:,jk,Krhs) = puu(:,:,jk,Krhs) + u_bkginc(:,:,jk) * zincwgt 
     664               pvv(:,:,jk,Krhs) = pvv(:,:,jk,Krhs) + v_bkginc(:,:,jk) * zincwgt 
    665665            END DO 
    666666            ! 
     
    680680            ! 
    681681            ! Initialize the now fields with the background + increment 
    682             un(:,:,:) = u_bkg(:,:,:) + u_bkginc(:,:,:) 
    683             vn(:,:,:) = v_bkg(:,:,:) + v_bkginc(:,:,:)   
    684             ! 
    685             ub(:,:,:) = un(:,:,:)         ! Update before fields 
    686             vb(:,:,:) = vn(:,:,:) 
     682            puu(:,:,:,Kmm) = u_bkg(:,:,:) + u_bkginc(:,:,:) 
     683            pvv(:,:,:,Kmm) = v_bkg(:,:,:) + v_bkginc(:,:,:)   
     684            ! 
     685            puu(:,:,:,Kbb) = puu(:,:,:,Kmm)         ! Update before fields 
     686            pvv(:,:,:,Kbb) = pvv(:,:,:,Kmm) 
    687687            ! 
    688688            DEALLOCATE( u_bkg    ) 
     
    697697 
    698698 
    699    SUBROUTINE ssh_asm_inc( kt ) 
     699   SUBROUTINE ssh_asm_inc( kt, Kbb, Kmm ) 
    700700      !!---------------------------------------------------------------------- 
    701701      !!                    ***  ROUTINE ssh_asm_inc  *** 
     
    707707      !! ** Action  :  
    708708      !!---------------------------------------------------------------------- 
    709       INTEGER, INTENT(IN) :: kt   ! Current time step 
     709      INTEGER, INTENT(IN) :: kt         ! Current time step 
     710      INTEGER, INTENT(IN) :: Kbb, Kmm   ! Current time step 
    710711      ! 
    711712      INTEGER :: it 
     
    754755            neuler = 0                                   ! Force Euler forward step 
    755756            ! 
    756             sshn(:,:) = ssh_bkg(:,:) + ssh_bkginc(:,:)   ! Initialize the now fields the background + increment 
    757             ! 
    758             sshb(:,:) = sshn(:,:)                        ! Update before fields 
    759             e3t_b(:,:,:) = e3t_n(:,:,:) 
    760 !!gm why not e3u_b, e3v_b, gdept_b ???? 
     757            ssh(:,:,Kmm) = ssh_bkg(:,:) + ssh_bkginc(:,:)   ! Initialize the now fields the background + increment 
     758            ! 
     759            ssh(:,:,Kbb) = ssh(:,:,Kmm)                        ! Update before fields 
     760            e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm) 
     761!!gm why not e3u(:,:,:,Kbb), e3v(:,:,:,Kbb), gdept(:,:,:,Kbb) ???? 
    761762            ! 
    762763            DEALLOCATE( ssh_bkg    ) 
     
    770771 
    771772 
    772    SUBROUTINE ssh_asm_div( kt, phdivn ) 
     773   SUBROUTINE ssh_asm_div( kt, Kbb, Kmm, phdivn ) 
    773774      !!---------------------------------------------------------------------- 
    774775      !!                  ***  ROUTINE ssh_asm_div  *** 
     
    784785      !!---------------------------------------------------------------------- 
    785786      INTEGER, INTENT(IN) :: kt                               ! ocean time-step index 
     787      INTEGER, INTENT(IN) :: Kbb, Kmm                         ! time level indices 
    786788      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   phdivn   ! horizontal divergence 
    787789      !! 
     
    791793      !  
    792794#if defined key_asminc 
    793       CALL ssh_asm_inc( kt ) !==   (calculate increments) 
     795      CALL ssh_asm_inc( kt, Kbb, Kmm ) !==   (calculate increments) 
    794796      ! 
    795797      IF( ln_linssh ) THEN  
    796          phdivn(:,:,1) = phdivn(:,:,1) - ssh_iau(:,:) / e3t_n(:,:,1) * tmask(:,:,1) 
     798         phdivn(:,:,1) = phdivn(:,:,1) - ssh_iau(:,:) / e3t(:,:,1,Kmm) * tmask(:,:,1) 
    797799      ELSE  
    798800         ALLOCATE( ztim(jpi,jpj) ) 
    799          ztim(:,:) = ssh_iau(:,:) / ( ht_n(:,:) + 1.0 - ssmask(:,:) ) 
     801         ztim(:,:) = ssh_iau(:,:) / ( ht(:,:) + 1.0 - ssmask(:,:) ) 
    800802         DO jk = 1, jpkm1                                  
    801803            phdivn(:,:,jk) = phdivn(:,:,jk) - ztim(:,:) * tmask(:,:,jk)  
  • NEMO/trunk/src/OCE/BDY/bdy_oce.F90

    r11536 r12377  
    141141   LOGICAL, ALLOCATABLE, DIMENSION(:,:,:,:) ::   lsend_bdyext   !: mark needed communication for given boundary, grid and neighbour 
    142142   LOGICAL, ALLOCATABLE, DIMENSION(:,:,:,:) ::   lrecv_bdyext   !:  when searching towards the exterior of the computational domain 
     143   !! * Substitutions 
     144#  include "do_loop_substitute.h90" 
    143145   !!---------------------------------------------------------------------- 
    144146   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
  • NEMO/trunk/src/OCE/BDY/bdydta.F90

    r12049 r12377  
    2323   USE phycst         ! physical constants 
    2424   USE sbcapr         ! atmospheric pressure forcing 
    25    USE sbctide        ! Tidal forcing or not 
     25   USE tide_mod, ONLY: ln_tide ! tidal forcing 
    2626   USE bdy_oce        ! ocean open boundary conditions   
    2727   USE bdytides       ! tidal forcing at boundaries 
     
    6868!$AGRIF_END_DO_NOT_TREAT 
    6969 
     70   !! * Substitutions 
     71#  include "do_loop_substitute.h90" 
    7072   !!---------------------------------------------------------------------- 
    7173   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
     
    7577CONTAINS 
    7678 
    77    SUBROUTINE bdy_dta( kt, kit, kt_offset ) 
     79   SUBROUTINE bdy_dta( kt, Kmm ) 
    7880      !!---------------------------------------------------------------------- 
    7981      !!                   ***  SUBROUTINE bdy_dta  *** 
     
    8587      !!---------------------------------------------------------------------- 
    8688      INTEGER, INTENT(in)           ::   kt           ! ocean time-step index  
    87       INTEGER, INTENT(in), OPTIONAL ::   kit          ! subcycle time-step index (for timesplitting option) 
    88       INTEGER, INTENT(in), OPTIONAL ::   kt_offset    ! time offset in units of timesteps. NB. if kit 
    89       !                                               ! is present then units = subcycle timesteps. 
    90       !                                               ! kt_offset = 0 => get data at "now" time level 
    91       !                                               ! kt_offset = -1 => get data at "before" time level 
    92       !                                               ! kt_offset = +1 => get data at "after" time level 
    93       !                                               ! etc. 
     89      INTEGER, INTENT(in)           ::   Kmm          ! ocean time level index 
    9490      ! 
    9591      INTEGER ::  jbdy, jfld, jstart, jend, ib, jl    ! dummy loop indices 
     
    105101      ! Initialise data arrays once for all from initial conditions where required 
    106102      !--------------------------------------------------------------------------- 
    107       IF( kt == nit000 .AND. .NOT.PRESENT(kit) ) THEN 
     103      IF( kt == nit000 ) THEN 
    108104 
    109105         ! Calculate depth-mean currents 
     
    122118                     ii = idx_bdy(jbdy)%nbi(ib,igrd) 
    123119                     ij = idx_bdy(jbdy)%nbj(ib,igrd) 
    124                      dta_bdy(jbdy)%ssh(ib) = sshn(ii,ij) * tmask(ii,ij,1)          
     120                     dta_bdy(jbdy)%ssh(ib) = ssh(ii,ij,Kmm) * tmask(ii,ij,1)          
    125121                  END DO 
    126122               ENDIF 
     
    130126                     ii = idx_bdy(jbdy)%nbi(ib,igrd) 
    131127                     ij = idx_bdy(jbdy)%nbj(ib,igrd) 
    132                      dta_bdy(jbdy)%u2d(ib) = un_b(ii,ij) * umask(ii,ij,1)          
     128                     dta_bdy(jbdy)%u2d(ib) = uu_b(ii,ij,Kmm) * umask(ii,ij,1)          
    133129                  END DO 
    134130                  igrd = 3 
     
    136132                     ii = idx_bdy(jbdy)%nbi(ib,igrd) 
    137133                     ij = idx_bdy(jbdy)%nbj(ib,igrd) 
    138                      dta_bdy(jbdy)%v2d(ib) = vn_b(ii,ij) * vmask(ii,ij,1)          
     134                     dta_bdy(jbdy)%v2d(ib) = vv_b(ii,ij,Kmm) * vmask(ii,ij,1)          
    139135                  END DO 
    140136               ENDIF 
     
    149145                        ii = idx_bdy(jbdy)%nbi(ib,igrd) 
    150146                        ij = idx_bdy(jbdy)%nbj(ib,igrd) 
    151                         dta_bdy(jbdy)%u3d(ib,ik) =  ( un(ii,ij,ik) - un_b(ii,ij) ) * umask(ii,ij,ik)          
     147                        dta_bdy(jbdy)%u3d(ib,ik) =  ( uu(ii,ij,ik,Kmm) - uu_b(ii,ij,Kmm) ) * umask(ii,ij,ik)          
    152148                     END DO 
    153149                  END DO 
     
    157153                        ii = idx_bdy(jbdy)%nbi(ib,igrd) 
    158154                        ij = idx_bdy(jbdy)%nbj(ib,igrd) 
    159                         dta_bdy(jbdy)%v3d(ib,ik) =  ( vn(ii,ij,ik) - vn_b(ii,ij) ) * vmask(ii,ij,ik)          
     155                        dta_bdy(jbdy)%v3d(ib,ik) =  ( vv(ii,ij,ik,Kmm) - vv_b(ii,ij,Kmm) ) * vmask(ii,ij,ik)          
    160156                     END DO 
    161157                  END DO 
     
    171167                        ii = idx_bdy(jbdy)%nbi(ib,igrd) 
    172168                        ij = idx_bdy(jbdy)%nbj(ib,igrd) 
    173                         dta_bdy(jbdy)%tem(ib,ik) = tsn(ii,ij,ik,jp_tem) * tmask(ii,ij,ik)          
    174                         dta_bdy(jbdy)%sal(ib,ik) = tsn(ii,ij,ik,jp_sal) * tmask(ii,ij,ik)          
     169                        dta_bdy(jbdy)%tem(ib,ik) = ts(ii,ij,ik,jp_tem,Kmm) * tmask(ii,ij,ik)          
     170                        dta_bdy(jbdy)%sal(ib,ik) = ts(ii,ij,ik,jp_sal,Kmm) * tmask(ii,ij,ik)          
    175171                     END DO 
    176172                  END DO 
     
    216212         ! read/update all bdy data 
    217213         ! ------------------------ 
    218          CALL fld_read( kt, 1, bf_alias, kit = kit, kt_offset = kt_offset ) 
    219  
     214         ! BDY: use pt_offset=0.5 as applied at the end of the step and fldread is referenced at the middle of the step 
     215         CALL fld_read( kt, 1, bf_alias, pt_offset = 0.5_wp, Kmm = Kmm ) 
    220216         ! apply some corrections in some specific cases... 
    221217         ! -------------------------------------------------- 
     
    254250               ij   = idx_bdy(jbdy)%nbj(ib,igrd) 
    255251               DO ik = 1, jpkm1 
    256                   dta_alias%u2d(ib) = dta_alias%u2d(ib) + e3u_n(ii,ij,ik) * umask(ii,ij,ik) * dta_alias%u3d(ib,ik) 
     252                  dta_alias%u2d(ib) = dta_alias%u2d(ib) + e3u(ii,ij,ik,Kmm) * umask(ii,ij,ik) * dta_alias%u3d(ib,ik) 
    257253               END DO 
    258                dta_alias%u2d(ib) =  dta_alias%u2d(ib) * r1_hu_n(ii,ij) 
     254               dta_alias%u2d(ib) =  dta_alias%u2d(ib) * r1_hu(ii,ij,Kmm) 
    259255               DO ik = 1, jpkm1            ! compute barocline zonal velocity and put it in u3d 
    260256                  dta_alias%u3d(ib,ik) = dta_alias%u3d(ib,ik) - dta_alias%u2d(ib) 
     
    267263               ij   = idx_bdy(jbdy)%nbj(ib,igrd) 
    268264               DO ik = 1, jpkm1 
    269                   dta_alias%v2d(ib) = dta_alias%v2d(ib) + e3v_n(ii,ij,ik) * vmask(ii,ij,ik) * dta_alias%v3d(ib,ik) 
     265                  dta_alias%v2d(ib) = dta_alias%v2d(ib) + e3v(ii,ij,ik,Kmm) * vmask(ii,ij,ik) * dta_alias%v3d(ib,ik) 
    270266               END DO 
    271                dta_alias%v2d(ib) =  dta_alias%v2d(ib) * r1_hv_n(ii,ij) 
     267               dta_alias%v2d(ib) =  dta_alias%v2d(ib) * r1_hv(ii,ij,Kmm) 
    272268               DO ik = 1, jpkm1            ! compute barocline meridional velocity and put it in v3d 
    273269                  dta_alias%v3d(ib,ik) = dta_alias%v3d(ib,ik) - dta_alias%v2d(ib) 
     
    275271            END DO 
    276272         ENDIF   ! ltotvel 
    277  
    278          ! update tidal harmonic forcing 
    279          IF( PRESENT(kit) .AND. nn_dyn2d_dta(jbdy) .GE. 2 ) THEN 
    280             CALL bdytide_update( kt = kt, idx = idx_bdy(jbdy), dta = dta_alias, td = tides(jbdy),   &  
    281                &                 kit = kit, kt_offset = kt_offset ) 
    282          ENDIF 
    283273 
    284274         !  atm surface pressure : add inverted barometer effect to ssh if it was read 
     
    343333                  nblen => idx_bdy(jbdy)%nblen 
    344334                  nblenrim => idx_bdy(jbdy)%nblenrim 
    345                   IF( cn_dyn2d(jbdy) == 'frs' ) THEN ; ilen1(:)=nblen(:) ; ELSE ; ilen1(:)=nblenrim(:) ; ENDIF  
    346                      IF ( dta_bdy(jbdy)%lneed_ssh   ) dta_bdy_s(jbdy)%ssh(1:ilen1(1)) = dta_bdy(jbdy)%ssh(1:ilen1(1)) 
    347                      IF ( dta_bdy(jbdy)%lneed_dyn2d ) dta_bdy_s(jbdy)%u2d(1:ilen1(2)) = dta_bdy(jbdy)%u2d(1:ilen1(2)) 
    348                      IF ( dta_bdy(jbdy)%lneed_dyn2d ) dta_bdy_s(jbdy)%v2d(1:ilen1(3)) = dta_bdy(jbdy)%v2d(1:ilen1(3)) 
    349                   ENDIF 
    350                END DO 
    351             ELSE ! Add tides if not split-explicit free surface else this is done in ts loop 
    352                ! 
    353                CALL bdy_dta_tides( kt=kt, kt_offset=kt_offset ) 
    354             ENDIF 
     335                  IF( cn_dyn2d(jbdy) == 'frs' ) THEN   ;   ilen1(:)=nblen(:) 
     336                  ELSE                                 ;   ilen1(:)=nblenrim(:) 
     337                  ENDIF 
     338                  IF ( dta_bdy(jbdy)%lneed_ssh   ) dta_bdy_s(jbdy)%ssh(1:ilen1(1)) = dta_bdy(jbdy)%ssh(1:ilen1(1)) 
     339                  IF ( dta_bdy(jbdy)%lneed_dyn2d ) dta_bdy_s(jbdy)%u2d(1:ilen1(2)) = dta_bdy(jbdy)%u2d(1:ilen1(2)) 
     340                  IF ( dta_bdy(jbdy)%lneed_dyn2d ) dta_bdy_s(jbdy)%v2d(1:ilen1(3)) = dta_bdy(jbdy)%v2d(1:ilen1(3)) 
     341               ENDIF 
     342            END DO 
     343         ELSE ! Add tides if not split-explicit free surface else this is done in ts loop 
     344            ! 
     345            ! BDY: use pt_offset=1.0 as applied at the end of the step and bdy_dta_tides is referenced at the middle of the step 
     346            CALL bdy_dta_tides( kt=kt, pt_offset = 1._wp ) 
    355347         ENDIF 
    356          ! 
    357          IF( ln_timing )   CALL timing_stop('bdy_dta') 
    358          ! 
    359       END SUBROUTINE bdy_dta 
     348      ENDIF 
     349      ! 
     350      IF( ln_timing )   CALL timing_stop('bdy_dta') 
     351      ! 
     352   END SUBROUTINE bdy_dta 
    360353 
    361354 
     
    373366      INTEGER ::   ierror, ios     !  
    374367      ! 
     368      INTEGER ::   nbdy_rdstart, nbdy_loc 
     369      CHARACTER(LEN=50)                      ::   cerrmsg       ! error string 
    375370      CHARACTER(len=3)                       ::   cl3           !  
    376371      CHARACTER(len=100)                     ::   cn_dir        ! Root directory for location of data files 
     
    415410      ! Read namelists 
    416411      ! -------------- 
    417       REWIND(numnam_cfg) 
     412      nbdy_rdstart = 1 
    418413      DO jbdy = 1, nb_bdy 
    419414 
     
    421416         WRITE(ctmp2, '(a,i2)') 'block nambdy_dta number ', jbdy 
    422417 
    423          ! There is only one nambdy_dta block in namelist_ref -> use it for each bdy so we do a rewind  
    424          REWIND(numnam_ref) 
     418         ! There is only one nambdy_dta block in namelist_ref -> use it for each bdy so we read from the beginning 
    425419         READ  ( numnam_ref, nambdy_dta, IOSTAT = ios, ERR = 901) 
    426420901      IF( ios /= 0 )   CALL ctl_nam ( ios , 'nambdy_dta in reference namelist' ) 
     
    431425            & .OR. ( dta_bdy(jbdy)%lneed_tra   .AND.       nn_tra_dta(jbdy)    == 1 )   & 
    432426            & .OR. ( dta_bdy(jbdy)%lneed_ice   .AND.       nn_ice_dta(jbdy)    == 1 )   )   THEN 
    433             ! WARNING: we don't do a rewind here, each bdy reads its own nambdy_dta block one after another 
    434             READ  ( numnam_cfg, nambdy_dta, IOSTAT = ios, ERR = 902 ) 
     427            ! 
     428            ! Need to support possibility of reading more than one 
     429            ! nambdy_dta from the namelist_cfg internal file. 
     430            ! Do this by finding the jbdy'th occurence of nambdy_dta in the 
     431            ! character buffer as the starting point. 
     432            ! 
     433            nbdy_loc = INDEX( numnam_cfg( nbdy_rdstart: ), 'nambdy_dta' ) 
     434            IF( nbdy_loc .GT. 0 ) THEN 
     435               nbdy_rdstart = nbdy_rdstart + nbdy_loc 
     436            ELSE 
     437               WRITE(cerrmsg,'(A,I4,A)') 'Error: entry number ',jbdy,' of nambdy_dta not found' 
     438               ios = -1 
     439               CALL ctl_nam ( ios , cerrmsg ) 
     440            ENDIF 
     441            READ  ( numnam_cfg( MAX( 1, nbdy_rdstart - 2 ): ), nambdy_dta, IOSTAT = ios, ERR = 902) 
    435442902         IF( ios >  0 )   CALL ctl_nam ( ios , 'nambdy_dta in configuration namelist' ) 
    436443            IF(lwm) WRITE( numond, nambdy_dta )            
     
    442449            IF( nn_ice_dta(jbdy) == 1 ) THEN   ! if we get ice bdy data from netcdf file 
    443450               CALL fld_fill(  bf(jp_bdya_i,jbdy:jbdy), bn_a_i, cn_dir, 'bdy_dta', 'a_i'//' '//ctmp1, ctmp2 )   ! use namelist info 
    444                CALL fld_clopn( bf(jp_bdya_i,jbdy), nyear, nmonth, nday )   ! not a problem when we call it again after 
     451               CALL fld_def( bf(jp_bdya_i,jbdy) ) 
     452               CALL iom_open( bf(jp_bdya_i,jbdy)%clname, bf(jp_bdya_i,jbdy)%num ) 
    445453               idvar = iom_varid( bf(jp_bdya_i,jbdy)%num, bf(jp_bdya_i,jbdy)%clvar, kndims=indims, kdimsz=i4dimsz, lduld=lluld ) 
    446454               IF( indims == 4 .OR. ( indims == 3 .AND. .NOT. lluld ) ) THEN   ;   ipl = i4dimsz(3)   ! xylt or xyl 
    447455               ELSE                                                            ;   ipl = 1            ! xy or xyt 
    448456               ENDIF 
     457               CALL iom_close( bf(jp_bdya_i,jbdy)%num ) 
    449458               bf(jp_bdya_i,jbdy)%clrootname = 'NOT USED'   ! reset to default value as this subdomain may not need to read this bdy 
    450459            ENDIF 
  • NEMO/trunk/src/OCE/BDY/bdydyn.F90

    r10068 r12377  
    3737CONTAINS 
    3838 
    39    SUBROUTINE bdy_dyn( kt, dyn3d_only ) 
     39   SUBROUTINE bdy_dyn( kt, Kbb, puu, pvv, Kaa, dyn3d_only ) 
    4040      !!---------------------------------------------------------------------- 
    4141      !!                  ***  SUBROUTINE bdy_dyn  *** 
     
    4444      !! 
    4545      !!---------------------------------------------------------------------- 
    46       INTEGER, INTENT(in)           ::   kt           ! Main time step counter 
    47       LOGICAL, INTENT(in), OPTIONAL ::   dyn3d_only   ! T => only update baroclinic velocities 
     46      INTEGER                             , INTENT(in)    ::   kt           ! Main time step counter 
     47      INTEGER                             , INTENT(in)    ::   Kbb, Kaa     ! Ocean time level indices 
     48      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) ::   puu, pvv     ! Ocean velocities (to be updated at open boundaries) 
     49      LOGICAL, OPTIONAL                   , INTENT(in)    ::   dyn3d_only   ! T => only update baroclinic velocities 
    4850      ! 
    4951      INTEGER ::   jk, ii, ij, ib_bdy, ib, igrd     ! Loop counter 
    5052      LOGICAL ::   ll_dyn2d, ll_dyn3d, ll_orlanski 
    51       REAL(wp), DIMENSION(jpi,jpj) :: pua2d, pva2d     ! after barotropic velocities 
     53      REAL(wp), DIMENSION(jpi,jpj) :: zua2d, zva2d     ! after barotropic velocities 
    5254      !!---------------------------------------------------------------------- 
    5355      ! 
     
    7072 
    7173      !                          ! "After" velocities:  
    72       pua2d(:,:) = 0._wp 
    73       pva2d(:,:) = 0._wp      
     74      zua2d(:,:) = 0._wp 
     75      zva2d(:,:) = 0._wp      
    7476      DO jk = 1, jpkm1 
    75          pua2d(:,:) = pua2d(:,:) + e3u_a(:,:,jk) * ua(:,:,jk) * umask(:,:,jk) 
    76          pva2d(:,:) = pva2d(:,:) + e3v_a(:,:,jk) * va(:,:,jk) * vmask(:,:,jk) 
     77         zua2d(:,:) = zua2d(:,:) + e3u(:,:,jk,Kaa) * puu(:,:,jk,Kaa) * umask(:,:,jk) 
     78         zva2d(:,:) = zva2d(:,:) + e3v(:,:,jk,Kaa) * pvv(:,:,jk,Kaa) * vmask(:,:,jk) 
    7779      END DO 
    78       pua2d(:,:) = pua2d(:,:) * r1_hu_a(:,:) 
    79       pva2d(:,:) = pva2d(:,:) * r1_hv_a(:,:) 
     80      zua2d(:,:) = zua2d(:,:) * r1_hu(:,:,Kaa) 
     81      zva2d(:,:) = zva2d(:,:) * r1_hv(:,:,Kaa) 
    8082 
    8183      DO jk = 1 , jpkm1 
    82          ua(:,:,jk) = ( ua(:,:,jk) - pua2d(:,:) ) * umask(:,:,jk) 
    83          va(:,:,jk) = ( va(:,:,jk) - pva2d(:,:) ) * vmask(:,:,jk) 
     84         puu(:,:,jk,Kaa) = ( puu(:,:,jk,Kaa) - zua2d(:,:) ) * umask(:,:,jk) 
     85         pvv(:,:,jk,Kaa) = ( pvv(:,:,jk,Kaa) - zva2d(:,:) ) * vmask(:,:,jk) 
    8486      END DO 
    8587 
     
    8789      IF( ll_orlanski ) THEN     ! "Before" velocities (Orlanski condition only)  
    8890         DO jk = 1 , jpkm1 
    89             ub(:,:,jk) = ( ub(:,:,jk) - ub_b(:,:) ) * umask(:,:,jk) 
    90             vb(:,:,jk) = ( vb(:,:,jk) - vb_b(:,:) ) * vmask(:,:,jk) 
     91            puu(:,:,jk,Kbb) = ( puu(:,:,jk,Kbb) - uu_b(:,:,Kbb) ) * umask(:,:,jk) 
     92            pvv(:,:,jk,Kbb) = ( pvv(:,:,jk,Kbb) - vv_b(:,:,Kbb) ) * vmask(:,:,jk) 
    9193         END DO 
    9294      ENDIF 
     
    9799      !------------------------------------------------------- 
    98100 
    99       IF( ll_dyn2d )   CALL bdy_dyn2d( kt, pua2d, pva2d, ub_b, vb_b, r1_hu_a(:,:), r1_hv_a(:,:), ssha ) 
     101      IF( ll_dyn2d )   CALL bdy_dyn2d( kt, zua2d, zva2d, uu_b(:,:,Kbb), vv_b(:,:,Kbb), r1_hu(:,:,Kaa), r1_hv(:,:,Kaa), ssh(:,:,Kaa) ) 
    100102 
    101       IF( ll_dyn3d )   CALL bdy_dyn3d( kt ) 
     103      IF( ll_dyn3d )   CALL bdy_dyn3d( kt, Kbb, puu, pvv, Kaa ) 
    102104 
    103105      !------------------------------------------------------- 
     
    106108      ! 
    107109      DO jk = 1 , jpkm1 
    108          ua(:,:,jk) = ( ua(:,:,jk) + pua2d(:,:) ) * umask(:,:,jk) 
    109          va(:,:,jk) = ( va(:,:,jk) + pva2d(:,:) ) * vmask(:,:,jk) 
     110         puu(:,:,jk,Kaa) = ( puu(:,:,jk,Kaa) + zua2d(:,:) ) * umask(:,:,jk) 
     111         pvv(:,:,jk,Kaa) = ( pvv(:,:,jk,Kaa) + zva2d(:,:) ) * vmask(:,:,jk) 
    110112      END DO 
    111113      ! 
    112114      IF ( ll_orlanski ) THEN 
    113115         DO jk = 1 , jpkm1 
    114             ub(:,:,jk) = ( ub(:,:,jk) + ub_b(:,:) ) * umask(:,:,jk) 
    115             vb(:,:,jk) = ( vb(:,:,jk) + vb_b(:,:) ) * vmask(:,:,jk) 
     116            puu(:,:,jk,Kbb) = ( puu(:,:,jk,Kbb) + uu_b(:,:,Kbb) ) * umask(:,:,jk) 
     117            pvv(:,:,jk,Kbb) = ( pvv(:,:,jk,Kbb) + vv_b(:,:,Kbb) ) * vmask(:,:,jk) 
    116118         END DO 
    117119      END IF 
  • NEMO/trunk/src/OCE/BDY/bdydyn3d.F90

    r11536 r12377  
    3333CONTAINS 
    3434 
    35    SUBROUTINE bdy_dyn3d( kt ) 
     35   SUBROUTINE bdy_dyn3d( kt, Kbb, puu, pvv, Kaa ) 
    3636      !!---------------------------------------------------------------------- 
    3737      !!                  ***  SUBROUTINE bdy_dyn3d  *** 
     
    4040      !! 
    4141      !!---------------------------------------------------------------------- 
    42       INTEGER, INTENT(in) ::   kt   ! Main time step counter 
     42      INTEGER                             , INTENT( in    ) ::   kt        ! Main time step counter 
     43      INTEGER                             , INTENT( in    ) ::   Kbb, Kaa  ! Time level indices 
     44      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT( inout ) ::   puu, pvv  ! Ocean velocities (to be updated at open boundaries) 
    4345      ! 
    4446      INTEGER  ::   ib_bdy, ir     ! BDY set index, rim index 
     
    5860            CASE('none')        ;   CYCLE 
    5961            CASE('frs' )        ! treat the whole boundary at once 
    60                IF( ir == 0) CALL bdy_dyn3d_frs( idx_bdy(ib_bdy), dta_bdy(ib_bdy), kt, ib_bdy ) 
     62                       IF( ir == 0) CALL bdy_dyn3d_frs( puu, pvv, Kaa, idx_bdy(ib_bdy), dta_bdy(ib_bdy), kt, ib_bdy ) 
    6163            CASE('specified')   ! treat the whole rim      at once 
    62                IF( ir == 0) CALL bdy_dyn3d_spe( idx_bdy(ib_bdy), dta_bdy(ib_bdy), kt, ib_bdy ) 
     64                       IF( ir == 0) CALL bdy_dyn3d_spe( puu, pvv, Kaa, idx_bdy(ib_bdy), dta_bdy(ib_bdy), kt, ib_bdy ) 
    6365            CASE('zero')        ! treat the whole rim      at once 
    64                IF( ir == 0) CALL bdy_dyn3d_zro( idx_bdy(ib_bdy), dta_bdy(ib_bdy), kt, ib_bdy ) 
    65             CASE('orlanski' )   ;   CALL bdy_dyn3d_orlanski( idx_bdy(ib_bdy), dta_bdy(ib_bdy), ib_bdy, llrim0, ll_npo=.false. ) 
    66             CASE('orlanski_npo');   CALL bdy_dyn3d_orlanski( idx_bdy(ib_bdy), dta_bdy(ib_bdy), ib_bdy, llrim0, ll_npo=.true.  ) 
    67             CASE('zerograd')    ;   CALL bdy_dyn3d_zgrad( idx_bdy(ib_bdy), dta_bdy(ib_bdy), kt, ib_bdy, llrim0 ) 
    68             CASE('neumann')     ;   CALL bdy_dyn3d_nmn( idx_bdy(ib_bdy), ib_bdy, llrim0 ) 
     66                       IF( ir == 0) CALL bdy_dyn3d_zro( puu, pvv, Kaa, idx_bdy(ib_bdy), dta_bdy(ib_bdy), kt, ib_bdy ) 
     67            CASE('orlanski' )   ;   CALL bdy_dyn3d_orlanski( Kbb, puu, pvv, Kaa, idx_bdy(ib_bdy), dta_bdy(ib_bdy), ib_bdy, llrim0, ll_npo=.false. ) 
     68            CASE('orlanski_npo');   CALL bdy_dyn3d_orlanski( Kbb, puu, pvv, Kaa, idx_bdy(ib_bdy), dta_bdy(ib_bdy), ib_bdy, llrim0, ll_npo=.true.  ) 
     69            CASE('zerograd')    ;   CALL bdy_dyn3d_zgrad( puu, pvv, Kaa, idx_bdy(ib_bdy), dta_bdy(ib_bdy), kt, ib_bdy, llrim0 ) 
     70            CASE('neumann')     ;   CALL bdy_dyn3d_nmn( puu, pvv, Kaa, idx_bdy(ib_bdy), ib_bdy, llrim0 ) 
    6971            CASE DEFAULT        ;   CALL ctl_stop( 'bdy_dyn3d : unrecognised option for open boundaries for baroclinic velocities' ) 
    7072            END SELECT 
     
    9799         ! 
    98100         IF( ANY(llsend2) .OR. ANY(llrecv2) ) THEN   ! if need to send/recv in at least one direction 
    99             CALL lbc_lnk( 'bdydyn2d', ua, 'U', -1., kfillmode=jpfillnothing ,lsend=llsend2, lrecv=llrecv2 ) 
     101            CALL lbc_lnk( 'bdydyn2d', puu(:,:,:,Kaa), 'U', -1., kfillmode=jpfillnothing ,lsend=llsend2, lrecv=llrecv2 ) 
    100102         END IF 
    101103         IF( ANY(llsend3) .OR. ANY(llrecv3) ) THEN   ! if need to send/recv in at least one direction 
    102             CALL lbc_lnk( 'bdydyn2d', va, 'V', -1., kfillmode=jpfillnothing ,lsend=llsend3, lrecv=llrecv3 ) 
     104            CALL lbc_lnk( 'bdydyn2d', pvv(:,:,:,Kaa), 'V', -1., kfillmode=jpfillnothing ,lsend=llsend3, lrecv=llrecv3 ) 
    103105         END IF 
    104106      END DO   ! ir 
     
    107109 
    108110 
    109    SUBROUTINE bdy_dyn3d_spe( idx, dta, kt , ib_bdy ) 
     111   SUBROUTINE bdy_dyn3d_spe( puu, pvv, Kaa, idx, dta, kt, ib_bdy ) 
    110112      !!---------------------------------------------------------------------- 
    111113      !!                  ***  SUBROUTINE bdy_dyn3d_spe  *** 
     
    115117      !! 
    116118      !!---------------------------------------------------------------------- 
    117       INTEGER        , INTENT(in) ::   kt      ! time step index 
    118       TYPE(OBC_INDEX), INTENT(in) ::   idx     ! OBC indices 
    119       TYPE(OBC_DATA) , INTENT(in) ::   dta     ! OBC external data 
    120       INTEGER        , INTENT(in) ::   ib_bdy  ! BDY set index 
     119      INTEGER                             , INTENT( in    ) ::   Kaa       ! Time level index 
     120      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT( inout ) ::   puu, pvv  ! Ocean velocities (to be updated at open boundaries) 
     121      TYPE(OBC_INDEX)                     , INTENT( in    ) ::   idx       ! OBC indices 
     122      TYPE(OBC_DATA)                      , INTENT( in    ) ::   dta       ! OBC external data 
     123      INTEGER                             , INTENT( in    ) ::   kt        ! Time step 
     124      INTEGER                             , INTENT( in    ) ::   ib_bdy    ! BDY set index 
    121125      ! 
    122126      INTEGER  ::   jb, jk         ! dummy loop indices 
     
    129133            ii   = idx%nbi(jb,igrd) 
    130134            ij   = idx%nbj(jb,igrd) 
    131             ua(ii,ij,jk) = dta%u3d(jb,jk) * umask(ii,ij,jk) 
     135            puu(ii,ij,jk,Kaa) = dta%u3d(jb,jk) * umask(ii,ij,jk) 
    132136         END DO 
    133137      END DO 
     
    138142            ii   = idx%nbi(jb,igrd) 
    139143            ij   = idx%nbj(jb,igrd) 
    140             va(ii,ij,jk) = dta%v3d(jb,jk) * vmask(ii,ij,jk) 
     144            pvv(ii,ij,jk,Kaa) = dta%v3d(jb,jk) * vmask(ii,ij,jk) 
    141145         END DO 
    142146      END DO 
     
    145149 
    146150 
    147    SUBROUTINE bdy_dyn3d_zgrad( idx, dta, kt, ib_bdy, llrim0 ) 
     151   SUBROUTINE bdy_dyn3d_zgrad( puu, pvv, Kaa, idx, dta, kt, ib_bdy, llrim0 ) 
    148152      !!---------------------------------------------------------------------- 
    149153      !!                  ***  SUBROUTINE bdy_dyn3d_zgrad  *** 
     
    152156      !! 
    153157      !!---------------------------------------------------------------------- 
    154       INTEGER                     ::   kt 
    155       TYPE(OBC_INDEX), INTENT(in) ::   idx      ! OBC indices 
    156       TYPE(OBC_DATA),  INTENT(in) ::   dta      ! OBC external data 
    157       INTEGER,         INTENT(in) ::   ib_bdy   ! BDY set index 
    158       LOGICAL,         INTENT(in) ::   llrim0   ! indicate if rim 0 is treated 
     158      INTEGER                             , INTENT( in    ) ::   Kaa       ! Time level index 
     159      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT( inout ) ::   puu, pvv  ! Ocean velocities (to be updated at open boundaries) 
     160      TYPE(OBC_INDEX)                     , INTENT( in    ) ::   idx       ! OBC indices 
     161      TYPE(OBC_DATA)                      , INTENT( in    ) ::   dta       ! OBC external data 
     162      INTEGER                             , INTENT( in    ) ::   kt 
     163      INTEGER                             , INTENT( in    ) ::   ib_bdy    ! BDY set index 
     164      LOGICAL                             , INTENT( in    ) ::   llrim0   ! indicate if rim 0 is treated 
    159165      !! 
    160166      INTEGER  ::   jb, jk         ! dummy loop indices 
     
    178184            ! 
    179185            DO jk = 1, jpkm1 
    180                ua(ii,ij,jk) = ua(ii,ij+flagv,jk) * umask(ii,ij+flagv,jk) 
     186               puu(ii,ij,jk,Kaa) = puu(ii,ij+flagv,jk,Kaa) * umask(ii,ij+flagv,jk) 
    181187            END DO 
    182188            ! 
     
    198204            ! 
    199205            DO jk = 1, jpkm1 
    200                va(ii,ij,jk) = va(ii+flagu,ij,jk) * vmask(ii+flagu,ij,jk) 
     206               pvv(ii,ij,jk,Kaa) = pvv(ii+flagu,ij,jk,Kaa) * vmask(ii+flagu,ij,jk) 
    201207            END DO 
    202208            ! 
     
    207213 
    208214 
    209    SUBROUTINE bdy_dyn3d_zro( idx, dta, kt, ib_bdy ) 
     215   SUBROUTINE bdy_dyn3d_zro( puu, pvv, Kaa, idx, dta, kt, ib_bdy ) 
    210216      !!---------------------------------------------------------------------- 
    211217      !!                  ***  SUBROUTINE bdy_dyn3d_zro  *** 
     
    214220      !! 
    215221      !!---------------------------------------------------------------------- 
    216       INTEGER        , INTENT(in) ::   kt      ! time step index 
    217       TYPE(OBC_INDEX), INTENT(in) ::   idx     ! OBC indices 
    218       TYPE(OBC_DATA) , INTENT(in) ::   dta     ! OBC external data 
    219       INTEGER,         INTENT(in) ::   ib_bdy  ! BDY set index 
     222      INTEGER                             , INTENT( in    ) ::   kt        ! time step index 
     223      INTEGER                             , INTENT( in    ) ::   Kaa       ! Time level index 
     224      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT( inout ) ::   puu, pvv  ! Ocean velocities (to be updated at open boundaries) 
     225      TYPE(OBC_INDEX)                     , INTENT( in    ) ::   idx       ! OBC indices 
     226      TYPE(OBC_DATA)                      , INTENT( in    ) ::   dta       ! OBC external data 
     227      INTEGER                             , INTENT( in    ) ::   ib_bdy    ! BDY set index 
    220228      ! 
    221229      INTEGER  ::   ib, ik         ! dummy loop indices 
     
    228236         ij = idx%nbj(ib,igrd) 
    229237         DO ik = 1, jpkm1 
    230             ua(ii,ij,ik) = 0._wp 
     238            puu(ii,ij,ik,Kaa) = 0._wp 
    231239         END DO 
    232240      END DO 
     
    237245         ij = idx%nbj(ib,igrd) 
    238246         DO ik = 1, jpkm1 
    239             va(ii,ij,ik) = 0._wp 
     247            pvv(ii,ij,ik,Kaa) = 0._wp 
    240248         END DO 
    241249      END DO 
     
    244252 
    245253 
    246    SUBROUTINE bdy_dyn3d_frs( idx, dta, kt, ib_bdy ) 
     254   SUBROUTINE bdy_dyn3d_frs( puu, pvv, Kaa, idx, dta, kt, ib_bdy ) 
    247255      !!---------------------------------------------------------------------- 
    248256      !!                  ***  SUBROUTINE bdy_dyn3d_frs  *** 
     
    255263      !!               topography. Tellus, 365-382. 
    256264      !!---------------------------------------------------------------------- 
    257       INTEGER        , INTENT(in) ::   kt      ! time step index 
    258       TYPE(OBC_INDEX), INTENT(in) ::   idx     ! OBC indices 
    259       TYPE(OBC_DATA) , INTENT(in) ::   dta     ! OBC external data 
    260       INTEGER,         INTENT(in) ::   ib_bdy  ! BDY set index 
     265      INTEGER                             , INTENT( in    ) ::   kt        ! time step index 
     266      INTEGER                             , INTENT( in    ) ::   Kaa       ! Time level index 
     267      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT( inout ) ::   puu, pvv  ! Ocean velocities (to be updated at open boundaries) 
     268      TYPE(OBC_INDEX)                     , INTENT( in    ) ::   idx       ! OBC indices 
     269      TYPE(OBC_DATA)                      , INTENT( in    ) ::   dta       ! OBC external data 
     270      INTEGER                             , INTENT( in    ) ::   ib_bdy    ! BDY set index 
    261271      ! 
    262272      INTEGER  ::   jb, jk         ! dummy loop indices 
     
    271281            ij   = idx%nbj(jb,igrd) 
    272282            zwgt = idx%nbw(jb,igrd) 
    273             ua(ii,ij,jk) = ( ua(ii,ij,jk) + zwgt * ( dta%u3d(jb,jk) - ua(ii,ij,jk) ) ) * umask(ii,ij,jk) 
     283            puu(ii,ij,jk,Kaa) = ( puu(ii,ij,jk,Kaa) + zwgt * ( dta%u3d(jb,jk) - puu(ii,ij,jk,Kaa) ) ) * umask(ii,ij,jk) 
    274284         END DO 
    275285      END DO 
     
    281291            ij   = idx%nbj(jb,igrd) 
    282292            zwgt = idx%nbw(jb,igrd) 
    283             va(ii,ij,jk) = ( va(ii,ij,jk) + zwgt * ( dta%v3d(jb,jk) - va(ii,ij,jk) ) ) * vmask(ii,ij,jk) 
     293            pvv(ii,ij,jk,Kaa) = ( pvv(ii,ij,jk,Kaa) + zwgt * ( dta%v3d(jb,jk) - pvv(ii,ij,jk,Kaa) ) ) * vmask(ii,ij,jk) 
    284294         END DO 
    285295      END DO    
     
    288298 
    289299 
    290    SUBROUTINE bdy_dyn3d_orlanski( idx, dta, ib_bdy, llrim0, ll_npo ) 
     300   SUBROUTINE bdy_dyn3d_orlanski( Kbb, puu, pvv, Kaa, idx, dta, ib_bdy, llrim0, ll_npo ) 
    291301      !!---------------------------------------------------------------------- 
    292302      !!                 ***  SUBROUTINE bdy_dyn3d_orlanski  *** 
     
    298308      !! References:  Marchesiello, McWilliams and Shchepetkin, Ocean Modelling vol. 3 (2001)     
    299309      !!---------------------------------------------------------------------- 
    300       TYPE(OBC_INDEX),              INTENT(in) ::   idx  ! OBC indices 
    301       TYPE(OBC_DATA),               INTENT(in) ::   dta  ! OBC external data 
    302       INTEGER,                      INTENT(in) ::   ib_bdy   ! BDY set index 
    303       LOGICAL,                      INTENT(in) ::   llrim0   ! indicate if rim 0 is treated 
    304       LOGICAL,                      INTENT(in) ::   ll_npo   ! switch for NPO version 
     310      INTEGER                             , INTENT( in    ) ::   Kbb, Kaa  ! Time level indices 
     311      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT( inout ) ::   puu, pvv  ! Ocean velocities (to be updated at open boundaries) 
     312      TYPE(OBC_INDEX)                     , INTENT( in    ) ::   idx       ! OBC indices 
     313      TYPE(OBC_DATA)                      , INTENT( in    ) ::   dta       ! OBC external data 
     314      INTEGER                             , INTENT( in    ) ::   ib_bdy    ! BDY set index 
     315      LOGICAL                             , INTENT( in    ) ::   llrim0    ! indicate if rim 0 is treated 
     316      LOGICAL                             , INTENT( in    ) ::   ll_npo    ! switch for NPO version 
    305317 
    306318      INTEGER  ::   jb, igrd                               ! dummy loop indices 
    307319      !!---------------------------------------------------------------------- 
    308320      ! 
    309       !! Note that at this stage the ub and ua arrays contain the baroclinic velocities.  
     321      !! Note that at this stage the puu(:,:,:,Kbb) and puu(:,:,:,Kaa) arrays contain the baroclinic velocities.  
    310322      ! 
    311323      igrd = 2      ! Orlanski bc on u-velocity;  
    312324      !             
    313       CALL bdy_orlanski_3d( idx, igrd, ub, ua, dta%u3d, ll_npo, llrim0 ) 
     325      CALL bdy_orlanski_3d( idx, igrd, puu(:,:,:,Kbb), puu(:,:,:,Kaa), dta%u3d, ll_npo, llrim0 ) 
    314326 
    315327      igrd = 3      ! Orlanski bc on v-velocity 
    316328      !   
    317       CALL bdy_orlanski_3d( idx, igrd, vb, va, dta%v3d, ll_npo, llrim0 ) 
     329      CALL bdy_orlanski_3d( idx, igrd, pvv(:,:,:,Kbb), pvv(:,:,:,Kaa), dta%v3d, ll_npo, llrim0 ) 
    318330      ! 
    319331   END SUBROUTINE bdy_dyn3d_orlanski 
    320332 
    321333 
    322    SUBROUTINE bdy_dyn3d_dmp( kt ) 
     334   SUBROUTINE bdy_dyn3d_dmp( kt, Kbb, puu, pvv, Krhs ) 
    323335      !!---------------------------------------------------------------------- 
    324336      !!                  ***  SUBROUTINE bdy_dyn3d_dmp  *** 
     
    327339      !! 
    328340      !!---------------------------------------------------------------------- 
    329       INTEGER, INTENT(in) ::   kt   ! time step index 
     341      INTEGER                             , INTENT( in    ) ::   kt        ! time step 
     342      INTEGER                             , INTENT( in    ) ::   Kbb, Krhs ! Time level indices 
     343      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT( inout ) ::   puu, pvv  ! Ocean velocities and trends (to be updated at open boundaries) 
    330344      ! 
    331345      INTEGER  ::   jb, jk         ! dummy loop indices 
     
    345359               zwgt = idx_bdy(ib_bdy)%nbd(jb,igrd) 
    346360               DO jk = 1, jpkm1 
    347                   ua(ii,ij,jk) = ( ua(ii,ij,jk) + zwgt * ( dta_bdy(ib_bdy)%u3d(jb,jk) - & 
    348                                    ub(ii,ij,jk) + ub_b(ii,ij)) ) * umask(ii,ij,jk) 
     361                  puu(ii,ij,jk,Krhs) = ( puu(ii,ij,jk,Krhs) + zwgt * ( dta_bdy(ib_bdy)%u3d(jb,jk) - & 
     362                                   puu(ii,ij,jk,Kbb) + uu_b(ii,ij,Kbb)) ) * umask(ii,ij,jk) 
    349363               END DO 
    350364            END DO 
     
    356370               zwgt = idx_bdy(ib_bdy)%nbd(jb,igrd) 
    357371               DO jk = 1, jpkm1 
    358                   va(ii,ij,jk) = ( va(ii,ij,jk) + zwgt * ( dta_bdy(ib_bdy)%v3d(jb,jk) -  & 
    359                                    vb(ii,ij,jk) + vb_b(ii,ij)) ) * vmask(ii,ij,jk) 
     372                  pvv(ii,ij,jk,Krhs) = ( pvv(ii,ij,jk,Krhs) + zwgt * ( dta_bdy(ib_bdy)%v3d(jb,jk) -  & 
     373                                   pvv(ii,ij,jk,Kbb) + vv_b(ii,ij,Kbb)) ) * vmask(ii,ij,jk) 
    360374               END DO 
    361375            END DO 
     
    368382 
    369383 
    370    SUBROUTINE bdy_dyn3d_nmn( idx, ib_bdy, llrim0 ) 
     384   SUBROUTINE bdy_dyn3d_nmn( puu, pvv, Kaa, idx, ib_bdy, llrim0 ) 
    371385      !!---------------------------------------------------------------------- 
    372386      !!                 ***  SUBROUTINE bdy_dyn3d_nmn  *** 
     
    377391      !! 
    378392      !!---------------------------------------------------------------------- 
    379       TYPE(OBC_INDEX), INTENT(in) ::   idx      ! OBC indices 
    380       INTEGER,         INTENT(in) ::   ib_bdy   ! BDY set index 
    381       LOGICAL,         INTENT(in) ::   llrim0   ! indicate if rim 0 is treated 
     393      INTEGER                             , INTENT( in    ) ::   Kaa       ! Time level index 
     394      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT( inout ) ::   puu, pvv  ! Ocean velocities (to be updated at open boundaries) 
     395      TYPE(OBC_INDEX)                     , INTENT( in    ) ::   idx       ! OBC indices 
     396      INTEGER                             , INTENT( in    ) ::   ib_bdy    ! BDY set index 
     397      LOGICAL                             , INTENT( in    ) ::   llrim0    ! indicate if rim 0 is treated 
    382398      INTEGER  ::   igrd                        ! dummy indice 
    383399      !!---------------------------------------------------------------------- 
    384400      ! 
    385       !! Note that at this stage the ub and ua arrays contain the baroclinic velocities.  
     401      !! Note that at this stage the puu(:,:,:,Kbb) and puu(:,:,:,Kaa) arrays contain the baroclinic velocities.  
    386402      ! 
    387403      igrd = 2      ! Neumann bc on u-velocity;  
    388404      !             
    389       CALL bdy_nmn( idx, igrd, ua, llrim0 )   ! ua is masked 
     405      CALL bdy_nmn( idx, igrd, puu(:,:,:,Kaa), llrim0 ) 
    390406 
    391407      igrd = 3      ! Neumann bc on v-velocity 
    392408      !   
    393       CALL bdy_nmn( idx, igrd, va, llrim0 )   ! va is masked 
     409      CALL bdy_nmn( idx, igrd, pvv(:,:,:,Kaa), llrim0 ) 
    394410      ! 
    395411   END SUBROUTINE bdy_dyn3d_nmn 
  • NEMO/trunk/src/OCE/BDY/bdyini.F90

    r12142 r12377  
    2222   USE bdydta         ! open boundary cond. setting   (bdy_dta_init routine) 
    2323   USE bdytides       ! open boundary cond. setting   (bdytide_init routine) 
    24    USE sbctide        ! Tidal forcing or not 
     24   USE tide_mod, ONLY: ln_tide ! tidal forcing 
    2525   USE phycst   , ONLY: rday 
    2626   ! 
     
    7575      ! Read namelist parameters 
    7676      ! ------------------------ 
    77       REWIND( numnam_ref )              ! Namelist nambdy in reference namelist :Unstructured open boundaries 
    7877      READ  ( numnam_ref, nambdy, IOSTAT = ios, ERR = 901) 
    7978901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nambdy in reference namelist' ) 
     
    9392      cn_ice         (2:jp_bdy) = cn_ice         (1) 
    9493      nn_ice_dta     (2:jp_bdy) = nn_ice_dta     (1) 
    95       REWIND( numnam_cfg )              ! Namelist nambdy in configuration namelist :Unstructured open boundaries 
    9694      READ  ( numnam_cfg, nambdy, IOSTAT = ios, ERR = 902 ) 
    9795902   IF( ios >  0 )   CALL ctl_nam ( ios , 'nambdy in configuration namelist' ) 
     
    364362      ! ------------------------------------------------- 
    365363 
    366       REWIND( numnam_cfg )      
    367364      nblendta(:,:) = 0 
    368365      nbdysege = 0 
     
    10801077      INTEGER          ::   ios                 ! Local integer output status for namelist read 
    10811078      INTEGER          ::   nbdyind, nbdybeg, nbdyend 
     1079      INTEGER          ::   nbdy_count, nbdy_rdstart, nbdy_loc 
    10821080      CHARACTER(LEN=1) ::   ctypebdy   !     -        -  
     1081      CHARACTER(LEN=50)::   cerrmsg    !     -        -  
    10831082      NAMELIST/nambdy_index/ ctypebdy, nbdyind, nbdybeg, nbdyend 
    10841083      !!---------------------------------------------------------------------- 
    1085  
    1086       ! No REWIND here because may need to read more than one nambdy_index namelist. 
    1087       ! Read only namelist_cfg to avoid unseccessfull overwrite  
    1088       ! keep full control of the configuration namelist 
    1089       READ  ( numnam_cfg, nambdy_index, IOSTAT = ios, ERR = 904 ) 
     1084      ! Need to support possibility of reading more than one nambdy_index from 
     1085      ! the namelist_cfg internal file. 
     1086      ! Do this by finding the kb_bdy'th occurence of nambdy_index in the 
     1087      ! character buffer as the starting point. 
     1088      nbdy_rdstart = 1 
     1089      DO nbdy_count = 1, kb_bdy 
     1090       nbdy_loc = INDEX( numnam_cfg( nbdy_rdstart: ), 'nambdy_index' ) 
     1091       IF( nbdy_loc .GT. 0 ) THEN 
     1092          nbdy_rdstart = nbdy_rdstart + nbdy_loc 
     1093       ELSE 
     1094          WRITE(cerrmsg,'(A,I4,A)') 'Error: entry number ',kb_bdy,' of nambdy_index not found' 
     1095          ios = -1 
     1096          CALL ctl_nam ( ios , cerrmsg ) 
     1097       ENDIF 
     1098      END DO 
     1099      nbdy_rdstart = MAX( 1, nbdy_rdstart - 2 ) 
     1100      READ  ( numnam_cfg( nbdy_rdstart: ), nambdy_index, IOSTAT = ios, ERR = 904) 
    10901101904   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nambdy_index in configuration namelist' ) 
    10911102      IF(lwm) WRITE ( numond, nambdy_index ) 
  • NEMO/trunk/src/OCE/BDY/bdylib.F90

    r11536 r12377  
    3535CONTAINS 
    3636 
    37    SUBROUTINE bdy_frs( idx, pta, dta ) 
     37   SUBROUTINE bdy_frs( idx, phia, dta ) 
    3838      !!---------------------------------------------------------------------- 
    3939      !!                 ***  SUBROUTINE bdy_frs  *** 
     
    4545      TYPE(OBC_INDEX),                     INTENT(in) ::   idx  ! OBC indices 
    4646      REAL(wp), DIMENSION(:,:),            INTENT(in) ::   dta  ! OBC external data 
    47       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pta  ! tracer trend 
     47      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   phia  ! tracer trend 
    4848      !! 
    4949      REAL(wp) ::   zwgt           ! boundary weight 
     
    5858            ij = idx%nbj(ib,igrd) 
    5959            zwgt = idx%nbw(ib,igrd) 
    60             pta(ii,ij,ik) = ( pta(ii,ij,ik) + zwgt * (dta(ib,ik) - pta(ii,ij,ik) ) ) * tmask(ii,ij,ik) 
     60            phia(ii,ij,ik) = ( phia(ii,ij,ik) + zwgt * (dta(ib,ik) - phia(ii,ij,ik) ) ) * tmask(ii,ij,ik) 
    6161         END DO 
    6262      END DO 
     
    6565 
    6666 
    67    SUBROUTINE bdy_spe( idx, pta, dta ) 
     67   SUBROUTINE bdy_spe( idx, phia, dta ) 
    6868      !!---------------------------------------------------------------------- 
    6969      !!                 ***  SUBROUTINE bdy_spe  *** 
     
    7474      TYPE(OBC_INDEX),                     INTENT(in) ::   idx  ! OBC indices 
    7575      REAL(wp), DIMENSION(:,:),            INTENT(in) ::   dta  ! OBC external data 
    76       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pta  ! tracer trend 
     76      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   phia  ! tracer trend 
    7777      !! 
    7878      INTEGER  ::   ib, ik, igrd   ! dummy loop indices 
     
    8585         ij = idx%nbj(ib,igrd) 
    8686         DO ik = 1, jpkm1 
    87             pta(ii,ij,ik) = dta(ib,ik) * tmask(ii,ij,ik) 
     87            phia(ii,ij,ik) = dta(ib,ik) * tmask(ii,ij,ik) 
    8888         END DO 
    8989      END DO 
     
    9292 
    9393 
    94    SUBROUTINE bdy_orl( idx, ptb, pta, dta, lrim0, ll_npo ) 
     94   SUBROUTINE bdy_orl( idx, phib, phia, dta, lrim0, ll_npo ) 
    9595      !!---------------------------------------------------------------------- 
    9696      !!                 ***  SUBROUTINE bdy_orl  *** 
     
    102102      TYPE(OBC_INDEX),                     INTENT(in) ::   idx  ! OBC indices 
    103103      REAL(wp), DIMENSION(:,:),            INTENT(in) ::   dta  ! OBC external data 
    104       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   ptb  ! before tracer field 
    105       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pta  ! tracer trend 
     104      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   phib  ! before tracer field 
     105      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   phia  ! tracer trend 
    106106      LOGICAL                 , OPTIONAL,  INTENT(in) ::   lrim0   ! indicate if rim 0 is treated 
    107107      LOGICAL,                             INTENT(in) ::   ll_npo  ! switch for NPO version 
     
    112112      igrd = 1                       ! Everything is at T-points here 
    113113      ! 
    114       CALL bdy_orlanski_3d( idx, igrd, ptb(:,:,:), pta(:,:,:), dta, lrim0, ll_npo ) 
     114      CALL bdy_orlanski_3d( idx, igrd, phib(:,:,:), phia(:,:,:), dta, lrim0, ll_npo ) 
    115115      ! 
    116116   END SUBROUTINE bdy_orl 
  • NEMO/trunk/src/OCE/BDY/bdytides.F90

    r11536 r12377  
    1818   USE phycst         ! physical constants 
    1919   USE bdy_oce        ! ocean open boundary conditions 
    20    USE tideini        !  
     20   USE tide_mod       !  
    2121   USE daymod         ! calendar 
    2222   ! 
     
    3030 
    3131   PUBLIC   bdytide_init     ! routine called in bdy_init 
    32    PUBLIC   bdytide_update   ! routine called in bdy_dta 
    3332   PUBLIC   bdy_dta_tides    ! routine called in dyn_spg_ts 
    3433 
     
    4544   TYPE(OBC_DATA)  , PUBLIC, DIMENSION(jp_bdy) :: dta_bdy_s  !: bdy external data (slow component) 
    4645 
     46   INTEGER ::   kt_tide 
     47 
     48   !! * Substitutions 
     49#  include "do_loop_substitute.h90" 
    4750   !!---------------------------------------------------------------------- 
    4851   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
     
    6467      CHARACTER(len=80)                         ::   filtide             !: Filename root for tidal input files 
    6568      LOGICAL                                   ::   ln_bdytide_2ddta    !: If true, read 2d harmonic data 
    66       LOGICAL                                   ::   ln_bdytide_conj     !: If true, assume complex conjugate tidal data 
    6769      !! 
    6870      INTEGER                                   ::   ib_bdy, itide, ib   !: dummy loop indices 
     
    7173      INTEGER, DIMENSION(3)                     ::   ilen0       !: length of boundary data (from OBC arrays) 
    7274      INTEGER                                   ::   ios                 ! Local integer output status for namelist read 
     75      INTEGER                                   ::   nbdy_rdstart, nbdy_loc 
     76      CHARACTER(LEN=50)                         ::   cerrmsg             !: error string 
    7377      CHARACTER(len=80)                         ::   clfile              !: full file name for tidal input file  
    7478      REAL(wp),ALLOCATABLE, DIMENSION(:,:,:)    ::   dta_read            !: work space to read in tidal harmonics data 
     
    7781      TYPE(TIDES_DATA),  POINTER                ::   td                  !: local short cut    
    7882      !! 
    79       NAMELIST/nambdy_tide/filtide, ln_bdytide_2ddta, ln_bdytide_conj 
     83      NAMELIST/nambdy_tide/filtide, ln_bdytide_2ddta 
    8084      !!---------------------------------------------------------------------- 
    8185      ! 
     
    8488      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~' 
    8589 
    86       REWIND(numnam_cfg) 
    87  
     90 
     91      nbdy_rdstart = 1 
    8892      DO ib_bdy = 1, nb_bdy 
    8993         IF( nn_dyn2d_dta(ib_bdy) >= 2 ) THEN 
     
    9498            filtide(:) = '' 
    9599 
    96             REWIND( numnam_ref ) 
    97100            READ  ( numnam_ref, nambdy_tide, IOSTAT = ios, ERR = 901) 
    98101901         IF( ios /= 0 )   CALL ctl_nam ( ios , 'nambdy_tide in reference namelist' ) 
    99             ! Don't REWIND here - may need to read more than one of these namelists.  
    100             READ  ( numnam_cfg, nambdy_tide, IOSTAT = ios, ERR = 902 ) 
     102            ! 
     103            ! Need to support possibility of reading more than one 
     104            ! nambdy_tide from the namelist_cfg internal file. 
     105            ! Do this by finding the ib_bdy'th occurence of nambdy_tide in the 
     106            ! character buffer as the starting point. 
     107            ! 
     108            nbdy_loc = INDEX( numnam_cfg( nbdy_rdstart: ), 'nambdy_tide' ) 
     109            IF( nbdy_loc .GT. 0 ) THEN 
     110               nbdy_rdstart = nbdy_rdstart + nbdy_loc 
     111            ELSE 
     112               WRITE(cerrmsg,'(A,I4,A)') 'Error: entry number ',ib_bdy,' of nambdy_tide not found' 
     113               ios = -1 
     114               CALL ctl_nam ( ios , cerrmsg ) 
     115            ENDIF 
     116            READ  ( numnam_cfg( MAX( 1, nbdy_rdstart - 2 ): ), nambdy_tide, IOSTAT = ios, ERR = 902) 
    101117902         IF( ios >  0 )   CALL ctl_nam ( ios , 'nambdy_tide in configuration namelist' ) 
    102118            IF(lwm) WRITE ( numond, nambdy_tide ) 
     
    105121            IF(lwp) WRITE(numout,*) '          Namelist nambdy_tide : tidal harmonic forcing at open boundaries' 
    106122            IF(lwp) WRITE(numout,*) '             read tidal data in 2d files: ', ln_bdytide_2ddta 
    107             IF(lwp) WRITE(numout,*) '             assume complex conjugate   : ', ln_bdytide_conj 
    108123            IF(lwp) WRITE(numout,*) '             Number of tidal components to read: ', nb_harmo 
    109124            IF(lwp) THEN  
    110125                    WRITE(numout,*) '             Tidal components: '  
    111126               DO itide = 1, nb_harmo 
    112                   WRITE(numout,*)  '                 ', Wave(ntide(itide))%cname_tide  
     127                  WRITE(numout,*)  '                 ', tide_harmonics(itide)%cname_tide  
    113128               END DO 
    114129            ENDIF  
     
    151166               igrd = 1                       ! Everything is at T-points here 
    152167               DO itide = 1, nb_harmo 
    153                   CALL iom_get( inum, jpdom_autoglo, TRIM(Wave(ntide(itide))%cname_tide)//'_z1', ztr(:,:) ) 
    154                   CALL iom_get( inum, jpdom_autoglo, TRIM(Wave(ntide(itide))%cname_tide)//'_z2', zti(:,:) )  
     168                  CALL iom_get( inum, jpdom_autoglo, TRIM(tide_harmonics(itide)%cname_tide)//'_z1', ztr(:,:) ) 
     169                  CALL iom_get( inum, jpdom_autoglo, TRIM(tide_harmonics(itide)%cname_tide)//'_z2', zti(:,:) )  
    155170                  DO ib = 1, ilen0(igrd) 
    156171                     ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 
     
    168183               igrd = 2                       ! Everything is at U-points here 
    169184               DO itide = 1, nb_harmo 
    170                   CALL iom_get  ( inum, jpdom_autoglo, TRIM(Wave(ntide(itide))%cname_tide)//'_u1', ztr(:,:) ) 
    171                   CALL iom_get  ( inum, jpdom_autoglo, TRIM(Wave(ntide(itide))%cname_tide)//'_u2', zti(:,:) ) 
     185                  CALL iom_get  ( inum, jpdom_autoglo, TRIM(tide_harmonics(itide)%cname_tide)//'_u1', ztr(:,:) ) 
     186                  CALL iom_get  ( inum, jpdom_autoglo, TRIM(tide_harmonics(itide)%cname_tide)//'_u2', zti(:,:) ) 
    172187                  DO ib = 1, ilen0(igrd) 
    173188                     ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 
     
    185200               igrd = 3                       ! Everything is at V-points here 
    186201               DO itide = 1, nb_harmo 
    187                   CALL iom_get  ( inum, jpdom_autoglo, TRIM(Wave(ntide(itide))%cname_tide)//'_v1', ztr(:,:) ) 
    188                   CALL iom_get  ( inum, jpdom_autoglo, TRIM(Wave(ntide(itide))%cname_tide)//'_v2', zti(:,:) ) 
     202                  CALL iom_get  ( inum, jpdom_autoglo, TRIM(tide_harmonics(itide)%cname_tide)//'_v1', ztr(:,:) ) 
     203                  CALL iom_get  ( inum, jpdom_autoglo, TRIM(tide_harmonics(itide)%cname_tide)//'_v2', zti(:,:) ) 
    189204                  DO ib = 1, ilen0(igrd) 
    190205                     ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 
     
    210225               DO itide = 1, nb_harmo 
    211226                  !                                                              ! SSH fields 
    212                   clfile = TRIM(filtide)//TRIM(Wave(ntide(itide))%cname_tide)//'_grid_T.nc' 
     227                  clfile = TRIM(filtide)//TRIM(tide_harmonics(itide)%cname_tide)//'_grid_T.nc' 
    213228                  CALL iom_open( clfile, inum ) 
    214229                  CALL fld_map( inum, 'z1' , dta_read(1:ilen0(1),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,1) ) 
     
    218233                  CALL iom_close( inum ) 
    219234                  !                                                              ! U fields 
    220                   clfile = TRIM(filtide)//TRIM(Wave(ntide(itide))%cname_tide)//'_grid_U.nc' 
     235                  clfile = TRIM(filtide)//TRIM(tide_harmonics(itide)%cname_tide)//'_grid_U.nc' 
    221236                  CALL iom_open( clfile, inum ) 
    222237                  CALL fld_map( inum, 'u1' , dta_read(1:ilen0(2),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,2) ) 
     
    226241                  CALL iom_close( inum ) 
    227242                  !                                                              ! V fields 
    228                   clfile = TRIM(filtide)//TRIM(Wave(ntide(itide))%cname_tide)//'_grid_V.nc' 
     243                  clfile = TRIM(filtide)//TRIM(tide_harmonics(itide)%cname_tide)//'_grid_V.nc' 
    229244                  CALL iom_open( clfile, inum ) 
    230245                  CALL fld_map( inum, 'v1' , dta_read(1:ilen0(3),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,3) ) 
     
    240255            ENDIF ! ln_bdytide_2ddta=.true. 
    241256            ! 
    242             IF( ln_bdytide_conj ) THEN    ! assume complex conjugate in data files 
    243                td%ssh0(:,:,2) = - td%ssh0(:,:,2) 
    244                td%u0  (:,:,2) = - td%u0  (:,:,2) 
    245                td%v0  (:,:,2) = - td%v0  (:,:,2) 
    246             ENDIF 
    247             ! 
    248257            ! Allocate slow varying data in the case of time splitting: 
    249258            ! Do it anyway because at this stage knowledge of free surface scheme is unknown 
     
    262271 
    263272 
    264    SUBROUTINE bdytide_update( kt, idx, dta, td, kit, kt_offset ) 
    265       !!---------------------------------------------------------------------- 
    266       !!                 ***  SUBROUTINE bdytide_update  *** 
    267       !!                 
    268       !! ** Purpose : - Add tidal forcing to ssh, u2d and v2d OBC data arrays.  
    269       !!                 
    270       !!---------------------------------------------------------------------- 
    271       INTEGER          , INTENT(in   ) ::   kt          ! Main timestep counter 
    272       TYPE(OBC_INDEX)  , INTENT(in   ) ::   idx         ! OBC indices 
    273       TYPE(OBC_DATA)   , INTENT(inout) ::   dta         ! OBC external data 
    274       TYPE(TIDES_DATA) , INTENT(inout) ::   td          ! tidal harmonics data 
    275       INTEGER, OPTIONAL, INTENT(in   ) ::   kit         ! Barotropic timestep counter (for timesplitting option) 
    276       INTEGER, OPTIONAL, INTENT(in   ) ::   kt_offset   ! time offset in units of timesteps. NB. if kit 
    277       !                                                 ! is present then units = subcycle timesteps. 
    278       !                                                 ! kt_offset = 0  => get data at "now"    time level 
    279       !                                                 ! kt_offset = -1 => get data at "before" time level 
    280       !                                                 ! kt_offset = +1 => get data at "after"  time level 
    281       !                                                 ! etc. 
    282       ! 
    283       INTEGER  ::   itide, igrd, ib       ! dummy loop indices 
    284       INTEGER  ::   time_add              ! time offset in units of timesteps 
    285       INTEGER, DIMENSION(3) ::   ilen0    ! length of boundary data (from OBC arrays) 
    286       REAL(wp) ::   z_arg, z_sarg, zflag, zramp   ! local scalars     
    287       REAL(wp), DIMENSION(jpmax_harmo) :: z_sist, z_cost 
    288       !!---------------------------------------------------------------------- 
    289       ! 
    290       ilen0(1) =  SIZE(td%ssh(:,1,1)) 
    291       ilen0(2) =  SIZE(td%u(:,1,1)) 
    292       ilen0(3) =  SIZE(td%v(:,1,1)) 
    293  
    294       zflag=1 
    295       IF ( PRESENT(kit) ) THEN 
    296         IF ( kit /= 1 ) zflag=0 
    297       ENDIF 
    298  
    299       IF ( (nsec_day == NINT(0.5_wp * rdt) .OR. kt==nit000) .AND. zflag==1 ) THEN 
    300         ! 
    301         kt_tide = kt - (nsec_day - 0.5_wp * rdt)/rdt 
    302         ! 
    303         IF(lwp) THEN 
    304            WRITE(numout,*) 
    305            WRITE(numout,*) 'bdytide_update : (re)Initialization of the tidal bdy forcing at kt=',kt 
    306            WRITE(numout,*) '~~~~~~~~~~~~~~ ' 
    307         ENDIF 
    308         ! 
    309         CALL tide_init_elevation ( idx, td ) 
    310         CALL tide_init_velocities( idx, td ) 
    311         ! 
    312       ENDIF  
    313  
    314       time_add = 0 
    315       IF( PRESENT(kt_offset) ) THEN 
    316          time_add = kt_offset 
    317       ENDIF 
    318           
    319       IF( PRESENT(kit) ) THEN   
    320          z_arg = ((kt-kt_tide) * rdt + (kit+0.5_wp*(time_add-1)) * rdt / REAL(nn_baro,wp) ) 
    321       ELSE                               
    322          z_arg = ((kt-kt_tide)+time_add) * rdt 
    323       ENDIF 
    324  
    325       ! Linear ramp on tidal component at open boundaries  
    326       zramp = 1._wp 
    327       IF (ln_tide_ramp) zramp = MIN(MAX( (z_arg + (kt_tide-nit000)*rdt)/(rdttideramp*rday),0._wp),1._wp) 
    328  
    329       DO itide = 1, nb_harmo 
    330          z_sarg = z_arg * omega_tide(itide) 
    331          z_cost(itide) = COS( z_sarg ) 
    332          z_sist(itide) = SIN( z_sarg ) 
    333       END DO 
    334  
    335       DO itide = 1, nb_harmo 
    336          igrd=1                              ! SSH on tracer grid 
    337          DO ib = 1, ilen0(igrd) 
    338             dta%ssh(ib) = dta%ssh(ib) + zramp*(td%ssh(ib,itide,1)*z_cost(itide) + td%ssh(ib,itide,2)*z_sist(itide)) 
    339          END DO 
    340          igrd=2                              ! U grid 
    341          DO ib = 1, ilen0(igrd) 
    342             dta%u2d(ib) = dta%u2d(ib) + zramp*(td%u  (ib,itide,1)*z_cost(itide) + td%u  (ib,itide,2)*z_sist(itide)) 
    343          END DO 
    344          igrd=3                              ! V grid 
    345          DO ib = 1, ilen0(igrd)  
    346             dta%v2d(ib) = dta%v2d(ib) + zramp*(td%v  (ib,itide,1)*z_cost(itide) + td%v  (ib,itide,2)*z_sist(itide)) 
    347          END DO 
    348       END DO 
    349       ! 
    350    END SUBROUTINE bdytide_update 
    351  
    352  
    353    SUBROUTINE bdy_dta_tides( kt, kit, kt_offset ) 
     273   SUBROUTINE bdy_dta_tides( kt, kit, pt_offset ) 
    354274      !!---------------------------------------------------------------------- 
    355275      !!                 ***  SUBROUTINE bdy_dta_tides  *** 
     
    360280      INTEGER,           INTENT(in) ::   kt          ! Main timestep counter 
    361281      INTEGER, OPTIONAL, INTENT(in) ::   kit         ! Barotropic timestep counter (for timesplitting option) 
    362       INTEGER, OPTIONAL, INTENT(in) ::   kt_offset   ! time offset in units of timesteps. NB. if kit 
    363       !                                              ! is present then units = subcycle timesteps. 
    364       !                                              ! kt_offset = 0  => get data at "now"    time level 
    365       !                                              ! kt_offset = -1 => get data at "before" time level 
    366       !                                              ! kt_offset = +1 => get data at "after"  time level 
    367       !                                              ! etc. 
     282      REAL(wp),OPTIONAL, INTENT(in) ::   pt_offset   ! time offset in units of timesteps 
    368283      ! 
    369284      LOGICAL  ::   lk_first_btstp            ! =.TRUE. if time splitting and first barotropic step 
    370285      INTEGER  ::   itide, ib_bdy, ib, igrd   ! loop indices 
    371       INTEGER  ::   time_add                  ! time offset in units of timesteps 
    372286      INTEGER, DIMENSION(jpbgrd)   ::   ilen0  
    373287      INTEGER, DIMENSION(1:jpbgrd) ::   nblen, nblenrim  ! short cuts 
    374       REAL(wp) ::   z_arg, z_sarg, zramp, zoff, z_cost, z_sist       
     288      REAL(wp) ::   z_arg, z_sarg, zramp, zoff, z_cost, z_sist, zt_offset    
    375289      !!---------------------------------------------------------------------- 
    376290      ! 
     
    378292      IF ( PRESENT(kit).AND.( kit /= 1 ) ) THEN ; lk_first_btstp=.FALSE. ; ENDIF 
    379293 
    380       time_add = 0 
    381       IF( PRESENT(kt_offset) ) THEN 
    382          time_add = kt_offset 
    383       ENDIF 
     294      zt_offset = 0._wp 
     295      IF( PRESENT(pt_offset) )   zt_offset = pt_offset 
    384296       
    385297      ! Absolute time from model initialization:    
    386298      IF( PRESENT(kit) ) THEN   
    387          z_arg = ( kt + (kit+time_add-1) / REAL(nn_baro,wp) ) * rdt 
     299         z_arg = ( REAL(kt, wp) + ( REAL(kit, wp) + zt_offset - 1. ) / REAL(nn_baro, wp) ) * rdt 
    388300      ELSE                               
    389          z_arg = ( kt + time_add ) * rdt 
     301         z_arg = ( REAL(kt, wp) + zt_offset ) * rdt 
    390302      ENDIF 
    391303 
    392304      ! Linear ramp on tidal component at open boundaries  
    393305      zramp = 1. 
    394       IF (ln_tide_ramp) zramp = MIN(MAX( (z_arg - nit000*rdt)/(rdttideramp*rday),0.),1.) 
     306      IF (ln_tide_ramp) zramp = MIN(MAX( (z_arg - REAL(nit000,wp)*rdt)/(rn_tide_ramp_dt*rday),0.),1.) 
    395307 
    396308      DO ib_bdy = 1,nb_bdy 
     
    409321            IF ( ( nsec_day == NINT(0.5_wp * rdt) .OR. kt==nit000 ) .AND. lk_first_btstp ) THEN 
    410322               ! 
    411                kt_tide = kt - (nsec_day - 0.5_wp * rdt)/rdt 
     323               kt_tide = kt - NINT((REAL(nsec_day,wp) - 0.5_wp * rdt)/rdt) 
    412324               ! 
    413325               IF(lwp) THEN 
     
    421333               ! 
    422334            ENDIF 
    423             zoff = -kt_tide * rdt ! time offset relative to nodal factor computation time 
     335            zoff = REAL(-kt_tide,wp) * rdt ! time offset relative to nodal factor computation time 
    424336            ! 
    425337            ! If time splitting, initialize arrays from slow varying open boundary data: 
     
    433345            DO itide = 1, nb_harmo 
    434346               ! 
    435                z_sarg = (z_arg + zoff) * omega_tide(itide) 
     347               z_sarg = (z_arg + zoff) * tide_harmonics(itide)%omega 
    436348               z_cost = zramp * COS( z_sarg ) 
    437349               z_sist = zramp * SIN( z_sarg ) 
     
    491403         END DO 
    492404         DO ib = 1 , ilen0(igrd) 
    493             mod_tide(ib)=mod_tide(ib)*ftide(itide) 
    494             phi_tide(ib)=phi_tide(ib)+v0tide(itide)+utide(itide) 
     405            mod_tide(ib)=mod_tide(ib)*tide_harmonics(itide)%f 
     406            phi_tide(ib)=phi_tide(ib)+tide_harmonics(itide)%v0+tide_harmonics(itide)%u 
    495407         ENDDO 
    496408         DO ib = 1 , ilen0(igrd) 
     
    530442         END DO 
    531443         DO ib = 1, ilen0(igrd) 
    532             mod_tide(ib)=mod_tide(ib)*ftide(itide) 
    533             phi_tide(ib)=phi_tide(ib)+v0tide(itide)+utide(itide) 
     444            mod_tide(ib)=mod_tide(ib)*tide_harmonics(itide)%f 
     445            phi_tide(ib)=phi_tide(ib)+tide_harmonics(itide)%v0 + tide_harmonics(itide)%u 
    534446         ENDDO 
    535447         DO ib = 1, ilen0(igrd) 
     
    551463         END DO 
    552464         DO ib = 1, ilen0(igrd) 
    553             mod_tide(ib)=mod_tide(ib)*ftide(itide) 
    554             phi_tide(ib)=phi_tide(ib)+v0tide(itide)+utide(itide) 
     465            mod_tide(ib)=mod_tide(ib)*tide_harmonics(itide)%f 
     466            phi_tide(ib)=phi_tide(ib)+tide_harmonics(itide)%v0 + tide_harmonics(itide)%u 
    555467         ENDDO 
    556468         DO ib = 1, ilen0(igrd) 
  • NEMO/trunk/src/OCE/BDY/bdytra.F90

    r11536 r12377  
    4040CONTAINS 
    4141 
    42    SUBROUTINE bdy_tra( kt ) 
     42   SUBROUTINE bdy_tra( kt, Kbb, pts, Kaa ) 
    4343      !!---------------------------------------------------------------------- 
    4444      !!                  ***  SUBROUTINE bdy_tra  *** 
     
    4747      !! 
    4848      !!---------------------------------------------------------------------- 
    49       INTEGER, INTENT(in) ::   kt   ! Main time step counter 
     49      INTEGER                                  , INTENT(in)    :: kt        ! Main time step counter 
     50      INTEGER                                  , INTENT(in)    :: Kbb, Kaa  ! time level indices 
     51      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts,jpt), INTENT(inout) :: pts       ! tracer fields 
    5052      ! 
    5153      INTEGER                        :: ib_bdy, jn, igrd, ir   ! Loop indeces 
     
    7072               CASE('none'        )   ;   CYCLE 
    7173               CASE('frs'         )   ! treat the whole boundary at once 
    72                   IF( ir == 0 ) CALL bdy_frs ( idx_bdy(ib_bdy),                tsa(:,:,:,jn), zdta(jn)%tra ) 
     74                  IF( ir == 0 ) CALL bdy_frs ( idx_bdy(ib_bdy),                pts(:,:,:,jn,Kaa), zdta(jn)%tra ) 
    7375               CASE('specified'   )   ! treat the whole rim      at once 
    74                   IF( ir == 0 ) CALL bdy_spe ( idx_bdy(ib_bdy),                tsa(:,:,:,jn), zdta(jn)%tra ) 
    75                CASE('neumann'     )   ;   CALL bdy_nmn ( idx_bdy(ib_bdy), igrd         , tsa(:,:,:,jn), llrim0 )   ! tsa masked 
    76                CASE('orlanski'    )   ;   CALL bdy_orl ( idx_bdy(ib_bdy), tsb(:,:,:,jn), tsa(:,:,:,jn), & 
     76                  IF( ir == 0 ) CALL bdy_spe ( idx_bdy(ib_bdy),                pts(:,:,:,jn,Kaa), zdta(jn)%tra ) 
     77               CASE('neumann'     )   ;   CALL bdy_nmn ( idx_bdy(ib_bdy), igrd         , pts(:,:,:,jn,Kaa), llrim0 )   ! tsa masked 
     78               CASE('orlanski'    )   ;   CALL bdy_orl ( idx_bdy(ib_bdy), pts(:,:,:,jn,Kbb), pts(:,:,:,jn,Kaa), & 
    7779                    & zdta(jn)%tra, llrim0, ll_npo=.false. ) 
    78                CASE('orlanski_npo')   ;   CALL bdy_orl ( idx_bdy(ib_bdy), tsb(:,:,:,jn), tsa(:,:,:,jn), & 
     80               CASE('orlanski_npo')   ;   CALL bdy_orl ( idx_bdy(ib_bdy), pts(:,:,:,jn,Kbb), pts(:,:,:,jn,Kaa), & 
    7981                    & zdta(jn)%tra, llrim0, ll_npo=.true.  ) 
    80                CASE('runoff'      )   ;   CALL bdy_rnf ( idx_bdy(ib_bdy),                tsa(:,:,:,jn), jn, llrim0 ) 
     82               CASE('runoff'      )   ;   CALL bdy_rnf ( idx_bdy(ib_bdy),                pts(:,:,:,jn,Kaa), jn, llrim0 ) 
    8183               CASE DEFAULT           ;   CALL ctl_stop( 'bdy_tra : unrecognised option for open boundaries for T and S' ) 
    8284               END SELECT 
     
    98100         END DO 
    99101         IF( ANY(llsend1) .OR. ANY(llrecv1) ) THEN   ! if need to send/recv in at least one direction 
    100             CALL lbc_lnk( 'bdytra', tsa, 'T',  1., kfillmode=jpfillnothing ,lsend=llsend1, lrecv=llrecv1 ) 
     102            CALL lbc_lnk( 'bdytra', pts(:,:,:,jn,Kaa), 'T',  1., kfillmode=jpfillnothing ,lsend=llsend1, lrecv=llrecv1 ) 
    101103         END IF 
    102104         ! 
     
    106108 
    107109 
    108    SUBROUTINE bdy_rnf( idx, pta, jpa, llrim0 ) 
     110   SUBROUTINE bdy_rnf( idx, pt, jpa, llrim0 ) 
    109111      !!---------------------------------------------------------------------- 
    110112      !!                 ***  SUBROUTINE bdy_rnf  *** 
     
    116118      !!---------------------------------------------------------------------- 
    117119      TYPE(OBC_INDEX),                     INTENT(in) ::   idx      ! OBC indices 
    118       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pta      ! tracer trend 
     120      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pt       ! tracer trend 
    119121      INTEGER,                             INTENT(in) ::   jpa      ! TRA index 
    120122      LOGICAL,                             INTENT(in) ::   llrim0   ! indicate if rim 0 is treated 
    121123      ! 
    122124      INTEGER  ::   ib, ii, ij, igrd   ! dummy loop indices 
    123       INTEGER  ::   ik, ip, jp ! 2D addresses 
    124125      !!---------------------------------------------------------------------- 
    125126      ! 
    126127      igrd = 1                       ! Everything is at T-points here 
    127128      IF(      jpa == jp_tem ) THEN 
    128          CALL bdy_nmn( idx, igrd, pta, llrim0 ) 
     129         CALL bdy_nmn( idx, igrd, pt, llrim0 ) 
    129130      ELSE IF( jpa == jp_sal ) THEN 
    130131         IF( .NOT. llrim0 )   RETURN 
     
    132133            ii = idx%nbi(ib,igrd) 
    133134            ij = idx%nbj(ib,igrd) 
    134             pta(ii,ij,1:jpkm1) = 0.1 * tmask(ii,ij,1:jpkm1) 
     135            pt(ii,ij,1:jpkm1) = 0.1 * tmask(ii,ij,1:jpkm1) 
    135136         END DO 
    136137      END IF 
     
    139140 
    140141 
    141    SUBROUTINE bdy_tra_dmp( kt ) 
     142   SUBROUTINE bdy_tra_dmp( kt, Kbb, pts, Krhs ) 
    142143      !!---------------------------------------------------------------------- 
    143144      !!                 ***  SUBROUTINE bdy_tra_dmp  *** 
     
    146147      !!  
    147148      !!---------------------------------------------------------------------- 
    148       INTEGER, INTENT(in) ::   kt   ! 
     149      INTEGER                                  , INTENT(in)    :: kt        ! time step 
     150      INTEGER                                  , INTENT(in)    :: Kbb, Krhs ! time level indices 
     151      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts,jpt), INTENT(inout) :: pts       ! active tracers and RHS of tracer equation 
    149152      ! 
    150153      REAL(wp) ::   zwgt           ! boundary weight 
     
    165168               zwgt = idx_bdy(ib_bdy)%nbd(ib,igrd) 
    166169               DO ik = 1, jpkm1 
    167                   zta = zwgt * ( dta_bdy(ib_bdy)%tem(ib,ik) - tsb(ii,ij,ik,jp_tem) ) * tmask(ii,ij,ik) 
    168                   zsa = zwgt * ( dta_bdy(ib_bdy)%sal(ib,ik) - tsb(ii,ij,ik,jp_sal) ) * tmask(ii,ij,ik) 
    169                   tsa(ii,ij,ik,jp_tem) = tsa(ii,ij,ik,jp_tem) + zta 
    170                   tsa(ii,ij,ik,jp_sal) = tsa(ii,ij,ik,jp_sal) + zsa 
     170                  zta = zwgt * ( dta_bdy(ib_bdy)%tem(ib,ik) - pts(ii,ij,ik,jp_tem,Kbb) ) * tmask(ii,ij,ik) 
     171                  zsa = zwgt * ( dta_bdy(ib_bdy)%sal(ib,ik) - pts(ii,ij,ik,jp_sal,Kbb) ) * tmask(ii,ij,ik) 
     172                  pts(ii,ij,ik,jp_tem,Krhs) = pts(ii,ij,ik,jp_tem,Krhs) + zta 
     173                  pts(ii,ij,ik,jp_sal,Krhs) = pts(ii,ij,ik,jp_sal,Krhs) + zsa 
    171174               END DO 
    172175            END DO 
  • NEMO/trunk/src/OCE/BDY/bdyvol.F90

    r12148 r12377  
    1414   USE bdy_oce        ! ocean open boundary conditions 
    1515   USE sbc_oce        ! ocean surface boundary conditions 
     16   USE isf_oce, ONLY : fwfisf_cav, fwfisf_par  ! ice shelf 
    1617   USE dom_oce        ! ocean space and time domain  
    1718   USE phycst         ! physical constants 
    18    USE sbcisf         ! ice shelf 
    1919   ! 
    2020   USE in_out_manager ! I/O manager 
     
    7777      ! Calculate the cumulate surface Flux z_cflxemp (m3/s) over all the domain 
    7878      ! ----------------------------------------------------------------------- 
    79       IF ( kc == 1 ) z_cflxemp = glob_sum( 'bdyvol', ( emp(:,:) - rnf(:,:) + fwfisf(:,:) ) * bdytmask(:,:) * e1e2t(:,:)  ) / rau0 
     79      IF ( kc == 1 ) z_cflxemp = glob_sum( 'bdyvol', ( emp(:,:) - rnf(:,:) + fwfisf_cav(:,:) + fwfisf_par(:,:) ) * bdytmask(:,:) * e1e2t(:,:)  ) / rau0 
    8080 
    8181      ! Compute bdy surface each cycle if non linear free surface 
  • NEMO/trunk/src/OCE/C1D/c1d.F90

    r11536 r12377  
    5050      !!---------------------------------------------------------------------- 
    5151      ! 
    52       REWIND( numnam_ref )              ! Namelist namc1d in reference namelist : Tracer advection scheme 
    5352      READ  ( numnam_ref, namc1d, IOSTAT = ios, ERR = 901) 
    5453901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namc1d in reference namelist' ) 
    5554      ! 
    56       REWIND( numnam_cfg )              ! Namelist namtra_adv in configuration namelist : Tracer advection scheme 
    5755      READ  ( numnam_cfg, namc1d, IOSTAT = ios, ERR = 902 ) 
    5856902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namc1d in configuration namelist' ) 
  • NEMO/trunk/src/OCE/C1D/dtauvd.F90

    r11536 r12377  
    3131   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_uvd   ! structure for input U & V current (file information and data) 
    3232 
     33   !! * Substitutions 
     34#  include "do_loop_substitute.h90" 
    3335   !!---------------------------------------------------------------------- 
    3436   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
     
    6062      ierr0 = 0   ;   ierr1 = 0   ;   ierr2 = 0  ;   ierr3 = 0 
    6163 
    62       REWIND( numnam_ref )              ! Namelist namc1d_uvd in reference namelist :  
    6364      READ  ( numnam_ref, namc1d_uvd, IOSTAT = ios, ERR = 901) 
    6465901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namc1d_uvd in reference namelist' ) 
    6566      ! 
    66       REWIND( numnam_cfg )              ! Namelist namc1d_uvd in configuration namelist : Parameters of the run 
    6767      READ  ( numnam_cfg, namc1d_uvd, IOSTAT = ios, ERR = 902 ) 
    6868902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namc1d_uvd in configuration namelist' ) 
     
    117117 
    118118 
    119    SUBROUTINE dta_uvd( kt, puvd ) 
     119   SUBROUTINE dta_uvd( kt, Kmm, puvd ) 
    120120      !!---------------------------------------------------------------------- 
    121121      !!                   ***  ROUTINE dta_uvd  *** 
     
    133133      !!---------------------------------------------------------------------- 
    134134      INTEGER                           , INTENT(in   ) ::   kt     ! ocean time-step 
     135      INTEGER                           , INTENT(in   ) ::   Kmm    ! time level index 
    135136      REAL(wp), DIMENSION(jpi,jpj,jpk,2), INTENT(  out) ::   puvd   ! U & V current data 
    136137      ! 
     
    157158         ENDIF 
    158159         ! 
    159          DO jj = 1, jpj                   ! vertical interpolation of U & V current: 
    160             DO ji = 1, jpi                ! determines the interpolated U & V current profiles at each (i,j) point 
    161                DO jk = 1, jpk 
    162                   zl = gdept_n(ji,jj,jk) 
    163                   IF    ( zl < gdept_1d(1  ) ) THEN          ! extrapolate above the first level of data 
    164                      zup(jk) =  puvd(ji,jj,1    ,1) 
    165                      zvp(jk) =  puvd(ji,jj,1    ,2) 
    166                   ELSEIF( zl > gdept_1d(jpk) ) THEN          ! extrapolate below the last level of data 
    167                      zup(jk) =  puvd(ji,jj,jpkm1,1) 
    168                      zvp(jk) =  puvd(ji,jj,jpkm1,2) 
    169                   ELSE                                      ! inbetween : vertical interpolation between jkk & jkk+1 
    170                      DO jkk = 1, jpkm1                      ! when  gdept(jkk) < zl < gdept(jkk+1) 
    171                         IF( (zl-gdept_1d(jkk)) * (zl-gdept_1d(jkk+1)) <= 0._wp ) THEN 
    172                            zi = ( zl - gdept_1d(jkk) ) / (gdept_1d(jkk+1)-gdept_1d(jkk)) 
    173                            zup(jk) = puvd(ji,jj,jkk,1) + ( puvd(ji,jj,jkk+1,1 ) - puvd(ji,jj,jkk,1) ) * zi  
    174                            zvp(jk) = puvd(ji,jj,jkk,2) + ( puvd(ji,jj,jkk+1,2 ) - puvd(ji,jj,jkk,2) ) * zi 
    175                         ENDIF 
    176                      END DO 
    177                   ENDIF 
    178                END DO 
    179                DO jk = 1, jpkm1           ! apply mask 
    180                   puvd(ji,jj,jk,1) = zup(jk) * umask(ji,jj,jk) 
    181                   puvd(ji,jj,jk,2) = zvp(jk) * vmask(ji,jj,jk) 
    182                END DO 
    183                puvd(ji,jj,jpk,1) = 0._wp 
    184                puvd(ji,jj,jpk,2) = 0._wp 
     160         DO_2D_11_11 
     161            DO jk = 1, jpk 
     162               zl = gdept(ji,jj,jk,Kmm) 
     163               IF    ( zl < gdept_1d(1  ) ) THEN          ! extrapolate above the first level of data 
     164                  zup(jk) =  puvd(ji,jj,1    ,1) 
     165                  zvp(jk) =  puvd(ji,jj,1    ,2) 
     166               ELSEIF( zl > gdept_1d(jpk) ) THEN          ! extrapolate below the last level of data 
     167                  zup(jk) =  puvd(ji,jj,jpkm1,1) 
     168                  zvp(jk) =  puvd(ji,jj,jpkm1,2) 
     169               ELSE                                      ! inbetween : vertical interpolation between jkk & jkk+1 
     170                  DO jkk = 1, jpkm1                      ! when  gdept(jkk) < zl < gdept(jkk+1) 
     171                     IF( (zl-gdept_1d(jkk)) * (zl-gdept_1d(jkk+1)) <= 0._wp ) THEN 
     172                        zi = ( zl - gdept_1d(jkk) ) / (gdept_1d(jkk+1)-gdept_1d(jkk)) 
     173                        zup(jk) = puvd(ji,jj,jkk,1) + ( puvd(ji,jj,jkk+1,1 ) - puvd(ji,jj,jkk,1) ) * zi  
     174                        zvp(jk) = puvd(ji,jj,jkk,2) + ( puvd(ji,jj,jkk+1,2 ) - puvd(ji,jj,jkk,2) ) * zi 
     175                     ENDIF 
     176                  END DO 
     177               ENDIF 
    185178            END DO 
    186          END DO 
     179            DO jk = 1, jpkm1           ! apply mask 
     180               puvd(ji,jj,jk,1) = zup(jk) * umask(ji,jj,jk) 
     181               puvd(ji,jj,jk,2) = zvp(jk) * vmask(ji,jj,jk) 
     182            END DO 
     183            puvd(ji,jj,jpk,1) = 0._wp 
     184            puvd(ji,jj,jpk,2) = 0._wp 
     185         END_2D 
    187186         !  
    188187         DEALLOCATE( zup, zvp ) 
     
    194193         ! 
    195194         IF( ln_zps ) THEN                ! zps-coordinate (partial steps) interpolation at the last ocean level 
    196             DO jj = 1, jpj 
    197                DO ji = 1, jpi 
    198                   ik = mbkt(ji,jj)  
    199                   IF( ik > 1 ) THEN 
    200                      zl = ( gdept_1d(ik) - gdept_0(ji,jj,ik) ) / ( gdept_1d(ik) - gdept_1d(ik-1) ) 
    201                      puvd(ji,jj,ik,1) = (1.-zl) * puvd(ji,jj,ik,1) + zl * puvd(ji,jj,ik-1,1) 
    202                      puvd(ji,jj,ik,2) = (1.-zl) * puvd(ji,jj,ik,2) + zl * puvd(ji,jj,ik-1,2) 
    203                   ENDIF 
    204                END DO 
    205             END DO 
     195            DO_2D_11_11 
     196               ik = mbkt(ji,jj)  
     197               IF( ik > 1 ) THEN 
     198                  zl = ( gdept_1d(ik) - gdept_0(ji,jj,ik) ) / ( gdept_1d(ik) - gdept_1d(ik-1) ) 
     199                  puvd(ji,jj,ik,1) = (1.-zl) * puvd(ji,jj,ik,1) + zl * puvd(ji,jj,ik-1,1) 
     200                  puvd(ji,jj,ik,2) = (1.-zl) * puvd(ji,jj,ik,2) + zl * puvd(ji,jj,ik-1,2) 
     201               ENDIF 
     202            END_2D 
    206203         ENDIF 
    207204         ! 
  • NEMO/trunk/src/OCE/C1D/dyncor_c1d.F90

    r10068 r12377  
    3131 
    3232   !! * Substitutions 
    33 #  include "vectopt_loop_substitute.h90" 
     33#  include "do_loop_substitute.h90" 
    3434   !!---------------------------------------------------------------------- 
    3535   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
     
    5656 
    5757 
    58    SUBROUTINE dyn_cor_c1d( kt ) 
     58   SUBROUTINE dyn_cor_c1d( kt, Kmm, puu, pvv, Krhs ) 
    5959      !!---------------------------------------------------------------------- 
    6060      !!                   ***  ROUTINE dyn_cor_c1d  *** 
     
    6363      !!               the general trend of the momentum equation in 1D case. 
    6464      !!---------------------------------------------------------------------- 
    65       INTEGER, INTENT( in ) ::   kt   ! ocean time-step index 
     65      INTEGER                             , INTENT(in   ) ::   kt        ! ocean time-step index 
     66      INTEGER                             , INTENT(in   ) ::   Kmm, Krhs ! ocean time level indices 
     67      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) ::   puu, pvv  ! ocean velocities and RHS of momentum equation 
    6668      !! 
    6769      INTEGER ::   ji, jj, jk   ! dummy loop indices 
     
    7577      ! 
    7678      IF( ln_stcor ) THEN 
    77          DO jk = 1, jpkm1 
    78             DO jj = 2, jpjm1 
    79                DO ji = fs_2, fs_jpim1   ! vector opt. 
    80                   ua(ji,jj,jk) = ua(ji,jj,jk) + ff_t(ji,jj) * (vn(ji,jj,jk) + vsd(ji,jj,jk)) 
    81                   va(ji,jj,jk) = va(ji,jj,jk) - ff_t(ji,jj) * (un(ji,jj,jk) + usd(ji,jj,jk)) 
    82                END DO 
    83             END DO 
    84          END DO 
     79         DO_3D_00_00( 1, jpkm1 ) 
     80            puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + ff_t(ji,jj) * (pvv(ji,jj,jk,Kmm) + vsd(ji,jj,jk)) 
     81            pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) - ff_t(ji,jj) * (puu(ji,jj,jk,Kmm) + usd(ji,jj,jk)) 
     82         END_3D 
    8583      ELSE 
    86          DO jk = 1, jpkm1 
    87             DO jj = 2, jpjm1 
    88                DO ji = fs_2, fs_jpim1   ! vector opt. 
    89                   ua(ji,jj,jk) = ua(ji,jj,jk) + ff_t(ji,jj) * vn(ji,jj,jk) 
    90                   va(ji,jj,jk) = va(ji,jj,jk) - ff_t(ji,jj) * un(ji,jj,jk) 
    91                END DO 
    92             END DO 
    93          END DO 
     84         DO_3D_00_00( 1, jpkm1 ) 
     85            puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + ff_t(ji,jj) * pvv(ji,jj,jk,Kmm) 
     86            pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) - ff_t(ji,jj) * puu(ji,jj,jk,Kmm) 
     87         END_3D 
    9488      END IF 
    9589       
    9690      ! 
    97       IF(ln_ctl)   CALL prt_ctl( tab3d_1=ua, clinfo1=' cor  - Ua: ', mask1=umask,  & 
    98          &                       tab3d_2=va, clinfo2=' Va: '       , mask2=vmask ) 
     91      IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab3d_1=puu(:,:,:,Krhs), clinfo1=' cor  - Ua: ', mask1=umask,  & 
     92         &                                  tab3d_2=pvv(:,:,:,Krhs), clinfo2=' Va: '       , mask2=vmask ) 
    9993      ! 
    10094   END SUBROUTINE dyn_cor_c1d 
  • NEMO/trunk/src/OCE/C1D/dyndmp.F90

    r11536 r12377  
    4343 
    4444   !! * Substitutions 
    45 #  include "vectopt_loop_substitute.h90" 
     45#  include "do_loop_substitute.h90" 
    4646   !!---------------------------------------------------------------------- 
    4747   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
     
    7979      !!---------------------------------------------------------------------- 
    8080      ! 
    81       REWIND( numnam_ref )              ! Namelist namc1d_dyndmp in reference namelist :  
    8281      READ  ( numnam_ref, namc1d_dyndmp, IOSTAT = ios, ERR = 901) 
    8382901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namc1d_dyndmp in reference namelist' ) 
    84       REWIND( numnam_cfg )              ! Namelist namc1d_dyndmp in configuration namelist : Parameters of the run 
    8583      READ  ( numnam_cfg, namc1d_dyndmp, IOSTAT = ios, ERR = 902 ) 
    8684902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namc1d_dyndmp in configuration namelist' ) 
     
    130128 
    131129 
    132    SUBROUTINE dyn_dmp( kt ) 
     130   SUBROUTINE dyn_dmp( kt, Kbb, Kmm, puu, pvv, Krhs ) 
    133131      !!---------------------------------------------------------------------- 
    134132      !!                   ***  ROUTINE dyn_dmp  *** 
     
    140138      !! ** Method  :   Compute Newtonian damping towards u_dta and v_dta  
    141139      !!      and add to the general momentum trends: 
    142       !!                     ua = ua + resto_uv * (u_dta - ub) 
    143       !!                     va = va + resto_uv * (v_dta - vb) 
     140      !!                     puu(Krhs) = puu(Krhs) + resto_uv * (u_dta - puu(Kbb)) 
     141      !!                     pvv(Krhs) = pvv(Krhs) + resto_uv * (v_dta - pvv(Kbb)) 
    144142      !!      The trend is computed either throughout the water column 
    145143      !!      (nn_zdmp=0), where the vertical mixing is weak (nn_zdmp=1) or 
    146144      !!      below the well mixed layer (nn_zdmp=2) 
    147145      !! 
    148       !! ** Action  : - (ua,va)   momentum trends updated with the damping trend 
    149       !!---------------------------------------------------------------------- 
    150       INTEGER, INTENT(in) ::   kt   ! ocean time-step index 
     146      !! ** Action  : - (puu(:,:,:,Krhs),pvv(:,:,:,Krhs))   momentum trends updated with the damping trend 
     147      !!---------------------------------------------------------------------- 
     148      INTEGER                             , INTENT(in   ) ::   kt             ! ocean time-step index 
     149      INTEGER                             , INTENT(in   ) ::   Kbb, Kmm, Krhs ! ocean time level indices 
     150      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) ::   puu, pvv       ! ocean velocities and RHS of momentum equation 
    151151      !! 
    152152      INTEGER  ::   ji, jj, jk   ! dummy loop indices 
     
    159159      ! 
    160160      !                           !==   read and interpolate U & V current data at kt   ==! 
    161       CALL dta_uvd( kt, zuv_dta ) !!! NOTE: This subroutine must be altered for use outside 
     161      CALL dta_uvd( kt, Kmm, zuv_dta ) !!! NOTE: This subroutine must be altered for use outside 
    162162                                  !!!       the C1D context (use of U,V grid variables) 
    163163      ! 
     
    165165      ! 
    166166      CASE( 0 )                   ! Newtonian damping throughout the water column 
    167          DO jk = 1, jpkm1 
    168             DO jj = 2, jpjm1 
    169                DO ji = fs_2, fs_jpim1   ! vector opt. 
    170                   zua = resto_uv(ji,jj,jk) * ( zuv_dta(ji,jj,jk,1) - ub(ji,jj,jk) ) 
    171                   zva = resto_uv(ji,jj,jk) * ( zuv_dta(ji,jj,jk,2) - vb(ji,jj,jk) ) 
    172                   ua(ji,jj,jk) = ua(ji,jj,jk) + zua 
    173                   va(ji,jj,jk) = va(ji,jj,jk) + zva 
    174                   utrdmp(ji,jj,jk) = zua           ! save the trends 
    175                   vtrdmp(ji,jj,jk) = zva       
    176                END DO 
    177             END DO 
    178          END DO 
     167         DO_3D_00_00( 1, jpkm1 ) 
     168            zua = resto_uv(ji,jj,jk) * ( zuv_dta(ji,jj,jk,1) - puu(ji,jj,jk,Kbb) ) 
     169            zva = resto_uv(ji,jj,jk) * ( zuv_dta(ji,jj,jk,2) - pvv(ji,jj,jk,Kbb) ) 
     170            puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + zua 
     171            pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + zva 
     172            utrdmp(ji,jj,jk) = zua           ! save the trends 
     173            vtrdmp(ji,jj,jk) = zva       
     174         END_3D 
    179175         ! 
    180176      CASE ( 1 )                  ! no damping above the turbocline (avt > 5 cm2/s) 
    181          DO jk = 1, jpkm1 
    182             DO jj = 2, jpjm1 
    183                DO ji = fs_2, fs_jpim1   ! vector opt. 
    184                   IF( avt(ji,jj,jk) <= avt_c ) THEN 
    185                      zua = resto_uv(ji,jj,jk) * ( zuv_dta(ji,jj,jk,1) - ub(ji,jj,jk) ) 
    186                      zva = resto_uv(ji,jj,jk) * ( zuv_dta(ji,jj,jk,2) - vb(ji,jj,jk) ) 
    187                   ELSE 
    188                      zua = 0._wp 
    189                      zva = 0._wp   
    190                   ENDIF 
    191                   ua(ji,jj,jk) = ua(ji,jj,jk) + zua 
    192                   va(ji,jj,jk) = va(ji,jj,jk) + zva 
    193                   utrdmp(ji,jj,jk) = zua           ! save the trends 
    194                   vtrdmp(ji,jj,jk) = zva 
    195                END DO 
    196             END DO 
    197          END DO 
     177         DO_3D_00_00( 1, jpkm1 ) 
     178            IF( avt(ji,jj,jk) <= avt_c ) THEN 
     179               zua = resto_uv(ji,jj,jk) * ( zuv_dta(ji,jj,jk,1) - puu(ji,jj,jk,Kbb) ) 
     180               zva = resto_uv(ji,jj,jk) * ( zuv_dta(ji,jj,jk,2) - pvv(ji,jj,jk,Kbb) ) 
     181            ELSE 
     182               zua = 0._wp 
     183               zva = 0._wp   
     184            ENDIF 
     185            puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + zua 
     186            pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + zva 
     187            utrdmp(ji,jj,jk) = zua           ! save the trends 
     188            vtrdmp(ji,jj,jk) = zva 
     189         END_3D 
    198190         ! 
    199191      CASE ( 2 )                  ! no damping in the mixed layer 
    200          DO jk = 1, jpkm1 
    201             DO jj = 2, jpjm1 
    202                DO ji = fs_2, fs_jpim1   ! vector opt. 
    203                   IF( gdept_n(ji,jj,jk) >= hmlp (ji,jj) ) THEN 
    204                      zua = resto_uv(ji,jj,jk) * ( zuv_dta(ji,jj,jk,1) - ub(ji,jj,jk) ) 
    205                      zva = resto_uv(ji,jj,jk) * ( zuv_dta(ji,jj,jk,2) - vb(ji,jj,jk) ) 
    206                   ELSE 
    207                      zua = 0._wp 
    208                      zva = 0._wp   
    209                   ENDIF 
    210                   ua(ji,jj,jk) = ua(ji,jj,jk) + zua 
    211                   va(ji,jj,jk) = va(ji,jj,jk) + zva 
    212                   utrdmp(ji,jj,jk) = zua           ! save the trends 
    213                   vtrdmp(ji,jj,jk) = zva 
    214                END DO 
    215             END DO 
    216          END DO 
     192         DO_3D_00_00( 1, jpkm1 ) 
     193            IF( gdept(ji,jj,jk,Kmm) >= hmlp (ji,jj) ) THEN 
     194               zua = resto_uv(ji,jj,jk) * ( zuv_dta(ji,jj,jk,1) - puu(ji,jj,jk,Kbb) ) 
     195               zva = resto_uv(ji,jj,jk) * ( zuv_dta(ji,jj,jk,2) - pvv(ji,jj,jk,Kbb) ) 
     196            ELSE 
     197               zua = 0._wp 
     198               zva = 0._wp   
     199            ENDIF 
     200            puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + zua 
     201            pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + zva 
     202            utrdmp(ji,jj,jk) = zua           ! save the trends 
     203            vtrdmp(ji,jj,jk) = zva 
     204         END_3D 
    217205         ! 
    218206      END SELECT 
    219207      ! 
    220208      !                           ! Control print 
    221       IF( ln_ctl   )   CALL prt_ctl( tab3d_1=ua(:,:,:), clinfo1=' dmp  - Ua: ', mask1=umask,   & 
    222          &                           tab3d_2=va(:,:,:), clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' ) 
     209      IF( sn_cfctl%l_prtctl   )   CALL prt_ctl( tab3d_1=puu(:,:,:,Krhs), clinfo1=' dmp  - Ua: ', mask1=umask,   & 
     210         &                                      tab3d_2=pvv(:,:,:,Krhs), clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' ) 
    223211      ! 
    224212      ! 
  • NEMO/trunk/src/OCE/C1D/step_c1d.F90

    r10068 r12377  
    66   !! History :   2.0  !  2004-04  (C. Ethe)  adapted from step.F90 for C1D 
    77   !!             3.0  !  2008-04  (G. Madec)  redo the adaptation to include SBC 
     8   !!             4.1  !  2019-08  (A. Coward, D. Storkey) rewrite in preparation for new timestepping scheme 
    89   !!---------------------------------------------------------------------- 
    910#if defined key_c1d 
     
    1415   !!---------------------------------------------------------------------- 
    1516   USE step_oce        ! time stepping definition modules  
     17   USE step, ONLY : Nbb, Nnn, Naa, Nrhs ! time level indices 
    1618#if defined key_top 
    1719   USE trcstp          ! passive tracer time-stepping      (trc_stp routine) 
    1820#endif 
    1921   USE dyncor_c1d      ! Coriolis term (c1d case)         (dyn_cor_1d     ) 
    20    USE dynnxt          ! time-stepping                    (dyn_nxt routine) 
     22   USE dynatf          ! time filtering                   (dyn_atf routine) 
    2123   USE dyndmp          ! U & V momentum damping           (dyn_dmp routine) 
    2224   USE restart         ! restart  
     
    6567      ! Update data, open boundaries, surface boundary condition (including sea-ice) 
    6668      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    67                          CALL sbc    ( kstp )         ! Sea Boundary Condition (including sea-ice) 
     69                         CALL sbc    ( kstp, Nbb, Nnn )  ! Sea Boundary Condition (including sea-ice) 
    6870 
    6971      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    70       ! Ocean physics update                (ua, va, ta, sa used as workspace) 
     72      ! Ocean physics update         
    7173      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    72                          CALL eos_rab( tsb, rab_b )   ! before local thermal/haline expension ratio at T-points 
    73                          CALL eos_rab( tsn, rab_n )   ! now    local thermal/haline expension ratio at T-points 
    74                          CALL bn2( tsb, rab_b, rn2b ) ! before Brunt-Vaisala frequency 
    75                          CALL bn2( tsn, rab_n, rn2 ) ! now    Brunt-Vaisala frequency 
     74                         CALL eos_rab( ts(:,:,:,:,Nbb), rab_b, Nnn )  ! before local thermal/haline expension ratio at T-points 
     75                         CALL eos_rab( ts(:,:,:,:,Nnn), rab_n, Nnn )  ! now    local thermal/haline expension ratio at T-points 
     76                         CALL bn2( ts(:,:,:,:,Nbb), rab_b, rn2b, Nnn ) ! before Brunt-Vaisala frequency 
     77                         CALL bn2( ts(:,:,:,:,Nnn), rab_n, rn2 , Nnn ) ! now    Brunt-Vaisala frequency 
    7678       
    7779      !  VERTICAL PHYSICS 
    78                          CALL zdf_phy( kstp )         ! vertical physics update (bfr, avt, avs, avm + MLD) 
     80                         CALL zdf_phy( kstp, Nbb, Nnn, Nrhs  )    ! vertical physics update (bfr, avt, avs, avm + MLD) 
    7981 
    80       IF(.NOT.ln_linssh )   CALL ssh_nxt       ( kstp )  ! after ssh (includes call to div_hor) 
    81       IF(.NOT.ln_linssh )   CALL dom_vvl_sf_nxt( kstp )  ! after vertical scale factors  
     82      IF(.NOT.ln_linssh )   CALL ssh_nxt       ( kstp, Nbb, Nnn, ssh, Naa )  ! after ssh (includes call to div_hor) 
     83      IF(.NOT.ln_linssh )   CALL dom_vvl_sf_nxt( kstp, Nbb, Nnn,      Naa )  ! after vertical scale factors  
    8284 
    83       IF(.NOT.ln_linssh )   CALL wzv           ( kstp )  ! now cross-level velocity  
     85      IF(.NOT.ln_linssh )   CALL wzv           ( kstp, Nbb, Nnn, ww,  Naa )  ! now cross-level velocity  
    8486      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    85       ! diagnostics and outputs             (ua, va, ta, sa used as workspace) 
     87      ! diagnostics and outputs        
    8688      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    87                          CALL dia_wri( kstp )       ! ocean model: outputs 
    88       IF( lk_diahth  )   CALL dia_hth( kstp )       ! Thermocline depth (20°C) 
     89                         CALL dia_wri( kstp, Nnn )  ! ocean model: outputs 
     90      IF( lk_diahth  )   CALL dia_hth( kstp, Nnn )  ! Thermocline depth (20°C) 
    8991 
    9092 
     
    9395      ! Passive Tracer Model 
    9496      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    95                         CALL trc_stp( kstp )       ! time-stepping 
     97                        CALL trc_stp( kstp, Nbb, Nnn, Nrhs, Naa  )   ! time-stepping 
    9698#endif 
    9799 
    98100      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    99       ! Active tracers                              (ua, va used as workspace) 
     101      ! Active tracers                              (uu(:,:,:,Nrhs), vv(:,:,:,Nrhs) used as workspace) 
    100102      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    101                         tsa(:,:,:,:) = 0._wp       ! set tracer trends to zero 
     103                        ts(:,:,:,:,Nrhs) = 0._wp       ! set tracer trends to zero 
    102104 
    103                         CALL tra_sbc( kstp )       ! surface boundary condition 
    104       IF( ln_traqsr )   CALL tra_qsr( kstp )       ! penetrative solar radiation qsr 
    105       IF( ln_tradmp )   CALL tra_dmp( kstp )       ! internal damping trends- tracers 
    106       IF(.NOT.ln_linssh)CALL tra_adv( kstp )       ! horizontal & vertical advection 
    107       IF( ln_zdfosm  )  CALL tra_osm( kstp )       ! OSMOSIS non-local tracer fluxes 
    108                         CALL tra_zdf( kstp )       ! vertical mixing 
    109                         CALL eos( tsn, rhd, rhop, gdept_0(:,:,:) )   ! now potential density for zdfmxl 
    110       IF( ln_zdfnpc )   CALL tra_npc( kstp )       ! applied non penetrative convective adjustment on (t,s) 
    111                         CALL tra_nxt( kstp )       ! tracer fields at next time step 
     105                        CALL tra_sbc( kstp,      Nnn, ts, Nrhs  )  ! surface boundary condition 
     106      IF( ln_traqsr )   CALL tra_qsr( kstp,      Nnn, ts, Nrhs  )  ! penetrative solar radiation qsr 
     107      IF( ln_tradmp )   CALL tra_dmp( kstp, Nbb, Nnn, ts, Nrhs  )  ! internal damping trends- tracers 
     108      IF(.NOT.ln_linssh)CALL tra_adv( kstp, Nbb, Nnn, ts, Nrhs  )  ! horizontal & vertical advection 
     109      IF( ln_zdfosm  )  CALL tra_osm( kstp, Nnn     , ts, Nrhs  )  ! OSMOSIS non-local tracer fluxes 
     110                        CALL tra_zdf( kstp, Nbb, Nnn, Nrhs, ts, Naa   )         ! vertical mixing 
     111                        CALL eos( ts(:,:,:,:,Nnn), rhd, rhop, gdept_0(:,:,:) )  ! now potential density for zdfmxl 
     112      IF( ln_zdfnpc )   CALL tra_npc( kstp,      Nnn, Nrhs, ts, Naa   )         ! applied non penetrative convective adjustment on (t,s) 
     113                        CALL tra_atf( kstp, Nbb, Nnn, Nrhs,     Naa, ts   )     ! time filtering of "now" tracer fields 
    112114 
    113115 
     
    116118      ! Dynamics                                    (ta, sa used as workspace) 
    117119      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    118                         ua(:,:,:) = 0._wp          ! set dynamics trends to zero 
    119                         va(:,:,:) = 0._wp 
     120                        uu(:,:,:,Nrhs) = 0._wp          ! set dynamics trends to zero 
     121                        vv(:,:,:,Nrhs) = 0._wp 
    120122 
    121       IF( ln_dyndmp )   CALL dyn_dmp    ( kstp )   ! internal damping trends- momentum 
    122                         CALL dyn_cor_c1d( kstp )   ! vorticity term including Coriolis 
    123       IF( ln_zdfosm  )  CALL dyn_osm    ( kstp )   ! OSMOSIS non-local velocity fluxes 
    124                         CALL dyn_zdf    ( kstp )   ! vertical diffusion 
    125                         CALL dyn_nxt    ( kstp )   ! lateral velocity at next time step 
    126       IF(.NOT.ln_linssh)CALL ssh_swp    ( kstp )   ! swap of sea surface height 
    127  
    128       IF(.NOT.ln_linssh)CALL dom_vvl_sf_swp( kstp )! swap of vertical scale factors 
     123      IF( ln_dyndmp )   CALL dyn_dmp    ( kstp, Nbb, Nnn      , uu, vv, Nrhs )  ! internal damping trends- momentum 
     124                        CALL dyn_cor_c1d( kstp,      Nnn      , uu, vv, Nrhs )  ! vorticity term including Coriolis 
     125      IF( ln_zdfosm  )  CALL dyn_osm    ( kstp,      Nnn      , uu, vv, Nrhs )  ! OSMOSIS non-local velocity fluxes 
     126                        CALL dyn_zdf    ( kstp, Nbb, Nnn, Nrhs, uu, vv, Naa  )  ! vertical diffusion 
     127                        CALL dyn_atf    ( kstp, Nbb, Nnn, Naa , uu, vv, e3t, e3u, e3v )  ! time filtering of "now" fields 
     128      IF(.NOT.ln_linssh)CALL ssh_atf    ( kstp, Nbb, Nnn, Naa , ssh )                    ! time filtering of "now" sea surface height 
     129      ! 
     130      ! Swap time levels 
     131      Nrhs = Nbb 
     132      Nbb = Nnn 
     133      Nnn = Naa 
     134      Naa = Nrhs 
     135      ! 
     136      IF(.NOT.ln_linssh)CALL dom_vvl_sf_update( kstp, Nbb, Nnn, Naa )                    ! update of vertical scale factors 
    129137 
    130138      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    131139      ! Control and restarts 
    132140      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    133                              CALL stp_ctl( kstp, indic ) 
    134       IF( kstp == nit000 )   CALL iom_close( numror )      ! close input  ocean restart file 
    135       IF( lrst_oce       )   CALL rst_write( kstp )        ! write output ocean restart file 
     141                             CALL stp_ctl( kstp, Nnn, indic ) 
     142      IF( kstp == nit000 )   CALL iom_close( numror )          ! close input  ocean restart file 
     143      IF( lrst_oce       )   CALL rst_write( kstp, Nbb, Nnn )  ! write output ocean restart file 
    136144      ! 
    137145#if defined key_iomput 
  • NEMO/trunk/src/OCE/CRS/crsdomwri.F90

    r10425 r12377  
    155155      DO jk = 1,jpk    
    156156         DO jj = 1, jpj_crsm1    
    157             DO ji = 1, jpi_crsm1  ! jes what to do for fs_jpim1??vector opt. 
     157            DO ji = 1, jpi_crsm1  ! jes what to do for jpim1??vector opt. 
    158158               zdepu(ji,jj,jk) = MIN( gdept_crs(ji,jj,jk) , gdept_crs(ji+1,jj  ,jk) ) * umask_crs(ji,jj,jk) 
    159159               zdepv(ji,jj,jk) = MIN( gdept_crs(ji,jj,jk) , gdept_crs(ji  ,jj+1,jk) ) * vmask_crs(ji,jj,jk) 
  • NEMO/trunk/src/OCE/CRS/crsfld.F90

    r10425 r12377  
    3232 
    3333   !! * Substitutions 
    34 #  include "vectopt_loop_substitute.h90" 
     34#  include "do_loop_substitute.h90" 
    3535   !!---------------------------------------------------------------------- 
    3636   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
     
    4040CONTAINS 
    4141 
    42    SUBROUTINE crs_fld( kt ) 
     42   SUBROUTINE crs_fld( kt, Kmm ) 
    4343      !!--------------------------------------------------------------------- 
    4444      !!                  ***  ROUTINE crs_fld  *** 
     
    5454      !!---------------------------------------------------------------------- 
    5555      INTEGER, INTENT(in) ::   kt   ! ocean time-step index 
     56      INTEGER, INTENT(in) ::   Kmm  ! time level index 
    5657      ! 
    5758      INTEGER  ::   ji, jj, jk        ! dummy loop indices 
     
    6768 
    6869      ! Depth work arrrays 
    69       ze3t(:,:,:) = e3t_n(:,:,:) 
    70       ze3u(:,:,:) = e3u_n(:,:,:) 
    71       ze3v(:,:,:) = e3v_n(:,:,:) 
    72       ze3w(:,:,:) = e3w_n(:,:,:) 
     70      ze3t(:,:,:) = e3t(:,:,:,Kmm) 
     71      ze3u(:,:,:) = e3u(:,:,:,Kmm) 
     72      ze3v(:,:,:) = e3v(:,:,:,Kmm) 
     73      ze3w(:,:,:) = e3w(:,:,:,Kmm) 
    7374 
    7475      IF( kt == nit000  ) THEN 
     
    9697 
    9798      !  Temperature 
    98       zt(:,:,:) = tsn(:,:,:,jp_tem)  ;      zt_crs(:,:,:) = 0._wp 
     99      zt(:,:,:) = ts(:,:,:,jp_tem,Kmm)  ;      zt_crs(:,:,:) = 0._wp 
    99100      CALL crs_dom_ope( zt, 'VOL', 'T', tmask, zt_crs, p_e12=e1e2t, p_e3=ze3t, psgn=1.0 ) 
    100101      tsn_crs(:,:,:,jp_tem) = zt_crs(:,:,:) 
     
    105106       
    106107      !  Salinity 
    107       zs(:,:,:) = tsn(:,:,:,jp_sal)  ;      zs_crs(:,:,:) = 0._wp 
     108      zs(:,:,:) = ts(:,:,:,jp_sal,Kmm)  ;      zs_crs(:,:,:) = 0._wp 
    108109      CALL crs_dom_ope( zs, 'VOL', 'T', tmask, zs_crs, p_e12=e1e2t, p_e3=ze3t, psgn=1.0 ) 
    109110      tsn_crs(:,:,:,jp_sal) = zt_crs(:,:,:) 
     
    113114 
    114115      !  U-velocity 
    115       CALL crs_dom_ope( un, 'SUM', 'U', umask, un_crs, p_e12=e2u, p_e3=ze3u, p_surf_crs=e2e3u_msk, psgn=-1.0 ) 
     116      CALL crs_dom_ope( uu(:,:,:,Kmm), 'SUM', 'U', umask, un_crs, p_e12=e2u, p_e3=ze3u, p_surf_crs=e2e3u_msk, psgn=-1.0 ) 
    116117      ! 
    117118      zt(:,:,:) = 0._wp     ;    zs(:,:,:) = 0._wp  ;   zt_crs(:,:,:) = 0._wp   ;    zs_crs(:,:,:) = 0._wp 
    118       DO jk = 1, jpkm1 
    119          DO jj = 2, jpjm1 
    120             DO ji = 2, jpim1    
    121                zt(ji,jj,jk)  = un(ji,jj,jk) * 0.5 * ( tsn(ji,jj,jk,jp_tem) + tsn(ji+1,jj,jk,jp_tem) )  
    122                zs(ji,jj,jk)  = un(ji,jj,jk) * 0.5 * ( tsn(ji,jj,jk,jp_sal) + tsn(ji+1,jj,jk,jp_sal) )  
    123             END DO 
    124          END DO 
    125       END DO 
     119      DO_3D_00_00( 1, jpkm1 ) 
     120         zt(ji,jj,jk)  = uu(ji,jj,jk,Kmm) * 0.5 * ( ts(ji,jj,jk,jp_tem,Kmm) + ts(ji+1,jj,jk,jp_tem,Kmm) )  
     121         zs(ji,jj,jk)  = uu(ji,jj,jk,Kmm) * 0.5 * ( ts(ji,jj,jk,jp_sal,Kmm) + ts(ji+1,jj,jk,jp_sal,Kmm) )  
     122      END_3D 
    126123      CALL crs_dom_ope( zt, 'SUM', 'U', umask, zt_crs, p_e12=e2u, p_e3=ze3u, p_surf_crs=e2e3u_msk, psgn=-1.0 ) 
    127124      CALL crs_dom_ope( zs, 'SUM', 'U', umask, zs_crs, p_e12=e2u, p_e3=ze3u, p_surf_crs=e2e3u_msk, psgn=-1.0 ) 
     
    132129 
    133130      !  V-velocity 
    134       CALL crs_dom_ope( vn, 'SUM', 'V', vmask, vn_crs, p_e12=e1v, p_e3=ze3v, p_surf_crs=e1e3v_msk, psgn=-1.0 ) 
     131      CALL crs_dom_ope( vv(:,:,:,Kmm), 'SUM', 'V', vmask, vn_crs, p_e12=e1v, p_e3=ze3v, p_surf_crs=e1e3v_msk, psgn=-1.0 ) 
    135132      !                                                                                  
    136133      zt(:,:,:) = 0._wp     ;    zs(:,:,:) = 0._wp  ;   zt_crs(:,:,:) = 0._wp   ;    zs_crs(:,:,:) = 0._wp 
    137       DO jk = 1, jpkm1 
    138          DO jj = 2, jpjm1 
    139             DO ji = 2, jpim1    
    140                zt(ji,jj,jk)  = vn(ji,jj,jk) * 0.5 * ( tsn(ji,jj,jk,jp_tem) + tsn(ji,jj+1,jk,jp_tem) )  
    141                zs(ji,jj,jk)  = vn(ji,jj,jk) * 0.5 * ( tsn(ji,jj,jk,jp_sal) + tsn(ji,jj+1,jk,jp_sal) )  
    142             END DO 
    143          END DO 
    144       END DO 
     134      DO_3D_00_00( 1, jpkm1 ) 
     135         zt(ji,jj,jk)  = vv(ji,jj,jk,Kmm) * 0.5 * ( ts(ji,jj,jk,jp_tem,Kmm) + ts(ji,jj+1,jk,jp_tem,Kmm) )  
     136         zs(ji,jj,jk)  = vv(ji,jj,jk,Kmm) * 0.5 * ( ts(ji,jj,jk,jp_sal,Kmm) + ts(ji,jj+1,jk,jp_sal,Kmm) )  
     137      END_3D 
    145138      CALL crs_dom_ope( zt, 'SUM', 'V', vmask, zt_crs, p_e12=e1v, p_e3=ze3v, p_surf_crs=e1e3v_msk, psgn=-1.0 ) 
    146139      CALL crs_dom_ope( zs, 'SUM', 'V', vmask, zs_crs, p_e12=e1v, p_e3=ze3v, p_surf_crs=e1e3v_msk, psgn=-1.0 ) 
     
    152145      IF( iom_use( "eken") ) THEN     !      kinetic energy 
    153146         z3d(:,:,jk) = 0._wp  
    154          DO jk = 1, jpkm1 
    155             DO jj = 2, jpjm1 
    156                DO ji = fs_2, fs_jpim1   ! vector opt. 
    157                   zztmp  = r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) 
    158                   z3d(ji,jj,jk) = 0.25_wp * zztmp * (                                    & 
    159                      &            un(ji-1,jj,jk)**2 * e2u(ji-1,jj) * e3u_n(ji-1,jj,jk)   & 
    160                      &          + un(ji  ,jj,jk)**2 * e2u(ji  ,jj) * e3u_n(ji  ,jj,jk)   & 
    161                      &          + vn(ji,jj-1,jk)**2 * e1v(ji,jj-1) * e3v_n(ji,jj-1,jk)   & 
    162                      &          + vn(ji,jj  ,jk)**2 * e1v(ji,jj  ) * e3v_n(ji,jj  ,jk)   ) 
    163                END DO 
    164             END DO 
    165          END DO 
     147         DO_3D_00_00( 1, jpkm1 ) 
     148            zztmp  = r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm) 
     149            z3d(ji,jj,jk) = 0.25_wp * zztmp * (                                    & 
     150               &            uu(ji-1,jj,jk,Kmm)**2 * e2u(ji-1,jj) * e3u(ji-1,jj,jk,Kmm)   & 
     151               &          + uu(ji  ,jj,jk,Kmm)**2 * e2u(ji  ,jj) * e3u(ji  ,jj,jk,Kmm)   & 
     152               &          + vv(ji,jj-1,jk,Kmm)**2 * e1v(ji,jj-1) * e3v(ji,jj-1,jk,Kmm)   & 
     153               &          + vv(ji,jj  ,jk,Kmm)**2 * e1v(ji,jj  ) * e3v(ji,jj  ,jk,Kmm)   ) 
     154         END_3D 
    166155         CALL lbc_lnk( 'crsfld', z3d, 'T', 1. ) 
    167156         ! 
     
    191180      !  W-velocity 
    192181      IF( ln_crs_wn ) THEN 
    193          CALL crs_dom_ope( wn, 'SUM', 'W', tmask, wn_crs, p_e12=e1e2t, p_surf_crs=e1e2w_msk, psgn=1.0 ) 
    194        !  CALL crs_dom_ope( wn, 'VOL', 'W', tmask, wn_crs, p_e12=e1e2t, p_e3=ze3w ) 
     182         CALL crs_dom_ope( ww, 'SUM', 'W', tmask, wn_crs, p_e12=e1e2t, p_surf_crs=e1e2w_msk, psgn=1.0 ) 
     183       !  CALL crs_dom_ope( ww, 'VOL', 'W', tmask, wn_crs, p_e12=e1e2t, p_e3=ze3w ) 
    195184      ELSE 
    196185        wn_crs(:,:,jpk) = 0._wp 
     
    219208       
    220209      !  sbc fields   
    221       CALL crs_dom_ope( sshn , 'VOL', 'T', tmask, sshn_crs , p_e12=e1e2t, p_e3=ze3t           , psgn=1.0 )   
     210      CALL crs_dom_ope( ssh(:,:,Kmm) , 'VOL', 'T', tmask, sshn_crs , p_e12=e1e2t, p_e3=ze3t           , psgn=1.0 )   
    222211      CALL crs_dom_ope( utau , 'SUM', 'U', umask, utau_crs , p_e12=e2u  , p_surf_crs=e2u_crs  , psgn=1.0 ) 
    223212      CALL crs_dom_ope( vtau , 'SUM', 'V', vmask, vtau_crs , p_e12=e1v  , p_surf_crs=e1v_crs  , psgn=1.0 ) 
  • NEMO/trunk/src/OCE/CRS/crsini.F90

    r11536 r12377  
    3535CONTAINS 
    3636    
    37    SUBROUTINE crs_init  
     37   SUBROUTINE crs_init( Kmm ) 
    3838      !!------------------------------------------------------------------- 
    3939      !!                     *** SUBROUTINE crs_init 
     
    6868      !!               - Read in pertinent data ? 
    6969      !!------------------------------------------------------------------- 
     70      INTEGER, INTENT(in) :: Kmm   ! time level index 
     71      ! 
    7072      INTEGER  :: ji,jj,jk      ! dummy indices 
    7173      INTEGER  :: ierr                                ! allocation error status 
     
    8082     !--------------------------------------------------------- 
    8183     ! 
    82       REWIND( numnam_ref )              ! Namelist namrun in reference namelist : Parameters of the run 
    8384      READ  ( numnam_ref, namcrs, IOSTAT = ios, ERR = 901) 
    8485901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namcrs in reference namelist' ) 
    85       REWIND( numnam_cfg )              ! Namelist namrun in configuration namelist : Parameters of the run 
    8686      READ  ( numnam_cfg, namcrs, IOSTAT = ios, ERR = 902 ) 
    8787902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namcrs in configuration namelist' ) 
     
    9898        WRITE(numout,*) '      create a mesh file (=T)               ln_msh_crs = ', ln_msh_crs 
    9999        WRITE(numout,*) '      type of Kz coarsening (0,1,2)         nn_crs_kz  = ', nn_crs_kz 
    100         WRITE(numout,*) '      wn coarsened or computed using hdivn  ln_crs_wn  = ', ln_crs_wn 
     100        WRITE(numout,*) '      ww coarsened or computed using hdiv  ln_crs_wn  = ', ln_crs_wn 
    101101     ENDIF 
    102102               
     
    174174      
    175175     ! 
    176      ze3t(:,:,:) = e3t_n(:,:,:) 
    177      ze3u(:,:,:) = e3u_n(:,:,:) 
    178      ze3v(:,:,:) = e3v_n(:,:,:) 
    179      ze3w(:,:,:) = e3w_n(:,:,:) 
     176     ze3t(:,:,:) = e3t(:,:,:,Kmm) 
     177     ze3u(:,:,:) = e3u(:,:,:,Kmm) 
     178     ze3v(:,:,:) = e3v(:,:,:,Kmm) 
     179     ze3w(:,:,:) = e3w(:,:,:,Kmm) 
    180180 
    181181     !    3.d.2   Surfaces  
  • NEMO/trunk/src/OCE/DIA/dia25h.F90

    r11536 r12377  
    3939CONTAINS 
    4040 
    41    SUBROUTINE dia_25h_init  
     41   SUBROUTINE dia_25h_init( Kbb ) 
    4242      !!--------------------------------------------------------------------------- 
    4343      !!                  ***  ROUTINE dia_25h_init  *** 
     
    4747      !! ** Method : Read namelist 
    4848      !!--------------------------------------------------------------------------- 
     49      INTEGER, INTENT(in) :: Kbb       ! Time level index 
     50      ! 
    4951      INTEGER ::   ios                 ! Local integer output status for namelist read 
    5052      INTEGER ::   ierror              ! Local integer for memory allocation 
     
    5355      !!---------------------------------------------------------------------- 
    5456      ! 
    55       REWIND ( numnam_ref )              ! Read Namelist nam_dia25h in reference namelist : 25hour mean diagnostics 
    5657      READ   ( numnam_ref, nam_dia25h, IOSTAT=ios, ERR= 901 ) 
    5758901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nam_dia25h in reference namelist' ) 
    58       REWIND( numnam_cfg )              ! Namelist nam_dia25h in configuration namelist  25hour diagnostics 
    5959      READ  ( numnam_cfg, nam_dia25h, IOSTAT = ios, ERR = 902 ) 
    6060902   IF( ios >  0 )   CALL ctl_nam ( ios , 'nam_dia25h in configuration namelist' ) 
     
    9595      ! ------------------------- ! 
    9696      cnt_25h = 1  ! sets the first value of sum at timestep 1 (note - should strictly be at timestep zero so before values used where possible)  
    97       tn_25h  (:,:,:) = tsb (:,:,:,jp_tem) 
    98       sn_25h  (:,:,:) = tsb (:,:,:,jp_sal) 
    99       sshn_25h(:,:)   = sshb(:,:) 
    100       un_25h  (:,:,:) = ub  (:,:,:) 
    101       vn_25h  (:,:,:) = vb  (:,:,:) 
     97      tn_25h  (:,:,:) = ts (:,:,:,jp_tem,Kbb) 
     98      sn_25h  (:,:,:) = ts (:,:,:,jp_sal,Kbb) 
     99      sshn_25h(:,:)   = ssh(:,:,Kbb) 
     100      un_25h  (:,:,:) = uu  (:,:,:,Kbb) 
     101      vn_25h  (:,:,:) = vv  (:,:,:,Kbb) 
    102102      avt_25h (:,:,:) = avt (:,:,:) 
    103103      avm_25h (:,:,:) = avm (:,:,:) 
     
    116116 
    117117 
    118    SUBROUTINE dia_25h( kt 
     118   SUBROUTINE dia_25h( kt, Kmm 
    119119      !!---------------------------------------------------------------------- 
    120120      !!                 ***  ROUTINE dia_25h  *** 
     
    125125      !!---------------------------------------------------------------------- 
    126126      INTEGER, INTENT(in) ::   kt   ! ocean time-step index 
     127      INTEGER, INTENT(in) ::   Kmm  ! ocean time level index 
    127128      !! 
    128129      INTEGER ::   ji, jj, jk 
     
    150151      ! wn_25h could not be initialised in dia_25h_init, so we do it here instead 
    151152      IF( kt == nn_it000 ) THEN 
    152          wn_25h(:,:,:) = wn(:,:,:) 
     153         wn_25h(:,:,:) = ww(:,:,:) 
    153154      ENDIF 
    154155 
     
    161162         ENDIF 
    162163 
    163          tn_25h  (:,:,:)     = tn_25h  (:,:,:) + tsn (:,:,:,jp_tem) 
    164          sn_25h  (:,:,:)     = sn_25h  (:,:,:) + tsn (:,:,:,jp_sal) 
    165          sshn_25h(:,:)       = sshn_25h(:,:)   + sshn(:,:) 
    166          un_25h  (:,:,:)     = un_25h  (:,:,:) + un  (:,:,:) 
    167          vn_25h  (:,:,:)     = vn_25h  (:,:,:) + vn  (:,:,:) 
    168          wn_25h  (:,:,:)     = wn_25h  (:,:,:) + wn  (:,:,:) 
     164         tn_25h  (:,:,:)     = tn_25h  (:,:,:) + ts (:,:,:,jp_tem,Kmm) 
     165         sn_25h  (:,:,:)     = sn_25h  (:,:,:) + ts (:,:,:,jp_sal,Kmm) 
     166         sshn_25h(:,:)       = sshn_25h(:,:)   + ssh(:,:,Kmm) 
     167         un_25h  (:,:,:)     = un_25h  (:,:,:) + uu  (:,:,:,Kmm) 
     168         vn_25h  (:,:,:)     = vn_25h  (:,:,:) + vv  (:,:,:,Kmm) 
     169         wn_25h  (:,:,:)     = wn_25h  (:,:,:) + ww  (:,:,:) 
    169170         avt_25h (:,:,:)     = avt_25h (:,:,:) + avt (:,:,:) 
    170171         avm_25h (:,:,:)     = avm_25h (:,:,:) + avm (:,:,:) 
     
    245246         ! 
    246247         ! After the write reset the values to cnt=1 and sum values equal current value  
    247          tn_25h  (:,:,:) = tsn (:,:,:,jp_tem) 
    248          sn_25h  (:,:,:) = tsn (:,:,:,jp_sal) 
    249          sshn_25h(:,:)   = sshn(:,:) 
    250          un_25h  (:,:,:) = un  (:,:,:) 
    251          vn_25h  (:,:,:) = vn  (:,:,:) 
    252          wn_25h  (:,:,:) = wn  (:,:,:) 
     248         tn_25h  (:,:,:) = ts (:,:,:,jp_tem,Kmm) 
     249         sn_25h  (:,:,:) = ts (:,:,:,jp_sal,Kmm) 
     250         sshn_25h(:,:)   = ssh(:,:,Kmm) 
     251         un_25h  (:,:,:) = uu  (:,:,:,Kmm) 
     252         vn_25h  (:,:,:) = vv  (:,:,:,Kmm) 
     253         wn_25h  (:,:,:) = ww  (:,:,:) 
    253254         avt_25h (:,:,:) = avt (:,:,:) 
    254255         avm_25h (:,:,:) = avm (:,:,:) 
  • NEMO/trunk/src/OCE/DIA/diaar5.F90

    r12276 r12377  
    3939       
    4040   !! * Substitutions 
    41 #  include "vectopt_loop_substitute.h90" 
     41#  include "do_loop_substitute.h90" 
    4242   !!---------------------------------------------------------------------- 
    4343   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
     
    6262 
    6363 
    64    SUBROUTINE dia_ar5( kt ) 
     64   SUBROUTINE dia_ar5( kt, Kmm ) 
    6565      !!---------------------------------------------------------------------- 
    6666      !!                    ***  ROUTINE dia_ar5  *** 
     
    7070      ! 
    7171      INTEGER, INTENT( in ) ::   kt   ! ocean time-step index 
     72      INTEGER, INTENT( in ) ::   Kmm  ! ocean time level index 
    7273      ! 
    7374      INTEGER  ::   ji, jj, jk, iks, ikb                      ! dummy loop arguments 
     
    8990         ALLOCATE( zrhd(jpi,jpj,jpk) , zrhop(jpi,jpj,jpk) ) 
    9091         ALLOCATE( ztsn(jpi,jpj,jpk,jpts) ) 
    91          zarea_ssh(:,:) = area(:,:) * sshn(:,:) 
     92         zarea_ssh(:,:) = area(:,:) * ssh(:,:,Kmm) 
    9293      ENDIF 
    9394      ! 
     
    99100         zrhd(:,:,jpk) = 0._wp        ! ocean volume ; rhd is used as workspace 
    100101         DO jk = 1, jpkm1 
    101             zrhd(:,:,jk) = area(:,:) * e3t_n(:,:,jk) * tmask(:,:,jk) 
     102            zrhd(:,:,jk) = area(:,:) * e3t(:,:,jk,Kmm) * tmask(:,:,jk) 
    102103         END DO 
    103104         CALL iom_put( 'volcello'  , zrhd(:,:,:)  )  ! WARNING not consistent with CMIP DR where volcello is at ca. 2000 
    104          CALL iom_put( 'masscello' , rau0 * e3t_n(:,:,:) * tmask(:,:,:) )  ! ocean mass 
     105         CALL iom_put( 'masscello' , rau0 * e3t(:,:,:,Kmm) * tmask(:,:,:) )  ! ocean mass 
    105106      ENDIF  
    106107      ! 
    107108      IF( iom_use( 'e3tb' ) )  THEN    ! bottom layer thickness 
    108          DO jj = 1, jpj 
    109             DO ji = 1, jpi 
    110                ikb = mbkt(ji,jj) 
    111                z2d(ji,jj) = e3t_n(ji,jj,ikb) 
    112             END DO 
    113          END DO 
     109         DO_2D_11_11 
     110            ikb = mbkt(ji,jj) 
     111            z2d(ji,jj) = e3t(ji,jj,ikb,Kmm) 
     112         END_2D 
    114113         CALL iom_put( 'e3tb', z2d ) 
    115114      ENDIF  
     
    122121         CALL iom_put( 'voltot', zvol               ) 
    123122         CALL iom_put( 'sshtot', zvolssh / area_tot ) 
    124          CALL iom_put( 'sshdyn', sshn(:,:) - (zvolssh / area_tot) ) 
     123         CALL iom_put( 'sshdyn', ssh(:,:,Kmm) - (zvolssh / area_tot) ) 
    125124         ! 
    126125      ENDIF 
     
    128127      IF( iom_use( 'botpres' ) .OR. iom_use( 'sshthster' )  .OR. iom_use( 'sshsteric' )  ) THEN     
    129128         !                      
    130          ztsn(:,:,:,jp_tem) = tsn(:,:,:,jp_tem)                    ! thermosteric ssh 
     129         ztsn(:,:,:,jp_tem) = ts(:,:,:,jp_tem,Kmm)                    ! thermosteric ssh 
    131130         ztsn(:,:,:,jp_sal) = sn0(:,:,:) 
    132          CALL eos( ztsn, zrhd, gdept_n(:,:,:) )                       ! now in situ density using initial salinity 
     131         CALL eos( ztsn, zrhd, gdept(:,:,:,Kmm) )                       ! now in situ density using initial salinity 
    133132         ! 
    134133         zbotpres(:,:) = 0._wp                        ! no atmospheric surface pressure, levitating sea-ice 
    135134         DO jk = 1, jpkm1 
    136             zbotpres(:,:) = zbotpres(:,:) + e3t_n(:,:,jk) * zrhd(:,:,jk) 
     135            zbotpres(:,:) = zbotpres(:,:) + e3t(:,:,jk,Kmm) * zrhd(:,:,jk) 
    137136         END DO 
    138137         IF( ln_linssh ) THEN 
     
    141140                  DO jj = 1, jpj 
    142141                     iks = mikt(ji,jj) 
    143                      zbotpres(ji,jj) = zbotpres(ji,jj) + sshn(ji,jj) * zrhd(ji,jj,iks) + riceload(ji,jj) 
     142                     zbotpres(ji,jj) = zbotpres(ji,jj) + ssh(ji,jj,Kmm) * zrhd(ji,jj,iks) + riceload(ji,jj) 
    144143                  END DO 
    145144               END DO 
    146145            ELSE 
    147                zbotpres(:,:) = zbotpres(:,:) + sshn(:,:) * zrhd(:,:,1) 
     146               zbotpres(:,:) = zbotpres(:,:) + ssh(:,:,Kmm) * zrhd(:,:,1) 
    148147            END IF 
    149148!!gm 
     
    157156       
    158157         !                                         ! steric sea surface height 
    159          CALL eos( tsn, zrhd, zrhop, gdept_n(:,:,:) )                 ! now in situ and potential density 
     158         CALL eos( ts(:,:,:,:,Kmm), zrhd, zrhop, gdept(:,:,:,Kmm) )                 ! now in situ and potential density 
    160159         zrhop(:,:,jpk) = 0._wp 
    161160         CALL iom_put( 'rhop', zrhop ) 
     
    163162         zbotpres(:,:) = 0._wp                        ! no atmospheric surface pressure, levitating sea-ice 
    164163         DO jk = 1, jpkm1 
    165             zbotpres(:,:) = zbotpres(:,:) + e3t_n(:,:,jk) * zrhd(:,:,jk) 
     164            zbotpres(:,:) = zbotpres(:,:) + e3t(:,:,jk,Kmm) * zrhd(:,:,jk) 
    166165         END DO 
    167166         IF( ln_linssh ) THEN 
     
    170169                  DO jj = 1,jpj 
    171170                     iks = mikt(ji,jj) 
    172                      zbotpres(ji,jj) = zbotpres(ji,jj) + sshn(ji,jj) * zrhd(ji,jj,iks) + riceload(ji,jj) 
     171                     zbotpres(ji,jj) = zbotpres(ji,jj) + ssh(ji,jj,Kmm) * zrhd(ji,jj,iks) + riceload(ji,jj) 
    173172                  END DO 
    174173               END DO 
    175174            ELSE 
    176                zbotpres(:,:) = zbotpres(:,:) + sshn(:,:) * zrhd(:,:,1) 
     175               zbotpres(:,:) = zbotpres(:,:) + ssh(:,:,Kmm) * zrhd(:,:,1) 
    177176            END IF 
    178177         END IF 
     
    183182         !                                         ! ocean bottom pressure 
    184183         zztmp = rau0 * grav * 1.e-4_wp               ! recover pressure from pressure anomaly and cover to dbar = 1.e4 Pa 
    185          zbotpres(:,:) = zztmp * ( zbotpres(:,:) + sshn(:,:) + thick0(:,:) ) 
     184         zbotpres(:,:) = zztmp * ( zbotpres(:,:) + ssh(:,:,Kmm) + thick0(:,:) ) 
    186185         CALL iom_put( 'botpres', zbotpres ) 
    187186         ! 
     
    191190          !                                         ! Mean density anomalie, temperature and salinity 
    192191          ztsn(:,:,:,:) = 0._wp                    ! ztsn(:,:,1,jp_tem/sal) is used here as 2D Workspace for temperature & salinity 
    193           DO jk = 1, jpkm1 
    194              DO jj = 1, jpj 
    195                 DO ji = 1, jpi 
    196                    zztmp = area(ji,jj) * e3t_n(ji,jj,jk) 
    197                    ztsn(ji,jj,1,jp_tem) = ztsn(ji,jj,1,jp_tem) + zztmp * tsn(ji,jj,jk,jp_tem) 
    198                    ztsn(ji,jj,1,jp_sal) = ztsn(ji,jj,1,jp_sal) + zztmp * tsn(ji,jj,jk,jp_sal) 
    199                 ENDDO 
    200              ENDDO 
    201           ENDDO 
     192          DO_3D_11_11( 1, jpkm1 ) 
     193             zztmp = area(ji,jj) * e3t(ji,jj,jk,Kmm) 
     194             ztsn(ji,jj,1,jp_tem) = ztsn(ji,jj,1,jp_tem) + zztmp * ts(ji,jj,jk,jp_tem,Kmm) 
     195             ztsn(ji,jj,1,jp_sal) = ztsn(ji,jj,1,jp_sal) + zztmp * ts(ji,jj,jk,jp_sal,Kmm) 
     196          END_3D 
    202197 
    203198          IF( ln_linssh ) THEN 
     
    206201                  DO jj = 1, jpj 
    207202                     iks = mikt(ji,jj) 
    208                      ztsn(ji,jj,1,jp_tem) = ztsn(ji,jj,1,jp_tem) + zarea_ssh(ji,jj) * tsn(ji,jj,iks,jp_tem)  
    209                      ztsn(ji,jj,1,jp_sal) = ztsn(ji,jj,1,jp_sal) + zarea_ssh(ji,jj) * tsn(ji,jj,iks,jp_sal)  
     203                     ztsn(ji,jj,1,jp_tem) = ztsn(ji,jj,1,jp_tem) + zarea_ssh(ji,jj) * ts(ji,jj,iks,jp_tem,Kmm)  
     204                     ztsn(ji,jj,1,jp_sal) = ztsn(ji,jj,1,jp_sal) + zarea_ssh(ji,jj) * ts(ji,jj,iks,jp_sal,Kmm)  
    210205                  END DO 
    211206               END DO 
    212207            ELSE 
    213                ztsn(:,:,1,jp_tem) = ztsn(:,:,1,jp_tem) + zarea_ssh(:,:) * tsn(:,:,1,jp_tem)  
    214                ztsn(:,:,1,jp_sal) = ztsn(:,:,1,jp_sal) + zarea_ssh(:,:) * tsn(:,:,1,jp_sal)  
     208               ztsn(:,:,1,jp_tem) = ztsn(:,:,1,jp_tem) + zarea_ssh(:,:) * ts(:,:,1,jp_tem,Kmm)  
     209               ztsn(:,:,1,jp_sal) = ztsn(:,:,1,jp_sal) + zarea_ssh(:,:) * ts(:,:,1,jp_sal,Kmm)  
    215210            END IF 
    216211         ENDIF 
     
    233228            ztpot(:,:,jpk) = 0._wp 
    234229            DO jk = 1, jpkm1 
    235                ztpot(:,:,jk) = eos_pt_from_ct( tsn(:,:,jk,jp_tem), tsn(:,:,jk,jp_sal) ) 
     230               ztpot(:,:,jk) = eos_pt_from_ct( ts(:,:,jk,jp_tem,Kmm), ts(:,:,jk,jp_sal,Kmm) ) 
    236231            END DO 
    237232            ! 
     
    242237               z2d(:,:) = 0._wp 
    243238               DO jk = 1, jpkm1 
    244                  z2d(:,:) = z2d(:,:) + area(:,:) * e3t_n(:,:,jk) * ztpot(:,:,jk) 
     239                 z2d(:,:) = z2d(:,:) + area(:,:) * e3t(:,:,jk,Kmm) * ztpot(:,:,jk) 
    245240               END DO 
    246241               ztemp = glob_sum( 'diaar5', z2d(:,:)  )  
     
    255250             IF( iom_use( 'tosmint_pot') ) THEN 
    256251               z2d(:,:) = 0._wp 
    257                DO jk = 1, jpkm1 
    258                   DO jj = 1, jpj 
    259                      DO ji = 1, jpi   ! vector opt. 
    260                         z2d(ji,jj) = z2d(ji,jj) + rau0 * e3t_n(ji,jj,jk) *  ztpot(ji,jj,jk) 
    261                      END DO 
    262                   END DO 
    263                END DO 
     252               DO_3D_11_11( 1, jpkm1 ) 
     253                  z2d(ji,jj) = z2d(ji,jj) + rau0 * e3t(ji,jj,jk,Kmm) *  ztpot(ji,jj,jk) 
     254               END_3D 
    264255               CALL iom_put( 'tosmint_pot', z2d )  
    265256            ENDIF 
     
    268259      ELSE        
    269260         IF( iom_use('ssttot') ) THEN   ! Output sst in case we use EOS-80 
    270             zsst  = glob_sum( 'diaar5', area(:,:) * tsn(:,:,1,jp_tem) ) 
     261            zsst  = glob_sum( 'diaar5', area(:,:) * ts(:,:,1,jp_tem,Kmm) ) 
    271262            CALL iom_put('ssttot', zsst / area_tot ) 
    272263         ENDIF 
     
    280271         zpe(:,:) = 0._wp 
    281272         IF( ln_zdfddm ) THEN 
    282             DO jk = 2, jpk 
    283                DO jj = 1, jpj 
    284                   DO ji = 1, jpi 
    285                      IF( rn2(ji,jj,jk) > 0._wp ) THEN 
    286                         zrw = ( gdept_n(ji,jj,jk) - gdepw_n(ji,jj,jk) ) / e3w_n(ji,jj,jk) 
    287                         ! 
    288                         zaw = rab_n(ji,jj,jk,jp_tem) * (1. - zrw) + rab_n(ji,jj,jk-1,jp_tem)* zrw 
    289                         zbw = rab_n(ji,jj,jk,jp_sal) * (1. - zrw) + rab_n(ji,jj,jk-1,jp_sal)* zrw 
    290                         ! 
    291                         zpe(ji, jj) = zpe(ji,jj)   & 
    292                            &        -  grav * (  avt(ji,jj,jk) * zaw * (tsn(ji,jj,jk-1,jp_tem) - tsn(ji,jj,jk,jp_tem) )  & 
    293                            &                   - avs(ji,jj,jk) * zbw * (tsn(ji,jj,jk-1,jp_sal) - tsn(ji,jj,jk,jp_sal) ) ) 
    294                      ENDIF 
    295                   END DO 
    296                END DO 
    297              END DO 
     273            DO_3D_11_11( 2, jpk ) 
     274               IF( rn2(ji,jj,jk) > 0._wp ) THEN 
     275                  zrw = ( gdept(ji,jj,jk,Kmm) - gdepw(ji,jj,jk,Kmm) ) / e3w(ji,jj,jk,Kmm) 
     276                  ! 
     277                  zaw = rab_n(ji,jj,jk,jp_tem) * (1. - zrw) + rab_n(ji,jj,jk-1,jp_tem)* zrw 
     278                  zbw = rab_n(ji,jj,jk,jp_sal) * (1. - zrw) + rab_n(ji,jj,jk-1,jp_sal)* zrw 
     279                  ! 
     280                  zpe(ji, jj) = zpe(ji,jj)   & 
     281                     &        -  grav * (  avt(ji,jj,jk) * zaw * (ts(ji,jj,jk-1,jp_tem,Kmm) - ts(ji,jj,jk,jp_tem,Kmm) )  & 
     282                     &                   - avs(ji,jj,jk) * zbw * (ts(ji,jj,jk-1,jp_sal,Kmm) - ts(ji,jj,jk,jp_sal,Kmm) ) ) 
     283               ENDIF 
     284            END_3D 
    298285          ELSE 
    299             DO jk = 1, jpk 
    300                DO ji = 1, jpi 
    301                   DO jj = 1, jpj 
    302                      zpe(ji,jj) = zpe(ji,jj) + avt(ji,jj,jk) * MIN(0._wp,rn2(ji,jj,jk)) * rau0 * e3w_n(ji,jj,jk) 
    303                   END DO 
    304                END DO 
    305             END DO 
     286            DO_3D_11_11( 1, jpk ) 
     287               zpe(ji,jj) = zpe(ji,jj) + avt(ji,jj,jk) * MIN(0._wp,rn2(ji,jj,jk)) * rau0 * e3w(ji,jj,jk,Kmm) 
     288            END_3D 
    306289         ENDIF 
    307290          CALL iom_put( 'tnpeo', zpe ) 
     
    320303 
    321304 
    322    SUBROUTINE dia_ar5_hst( ktra, cptr, pua, pva )  
     305   SUBROUTINE dia_ar5_hst( ktra, cptr, puflx, pvflx )  
    323306      !!---------------------------------------------------------------------- 
    324307      !!                    ***  ROUTINE dia_ar5_htr *** 
     
    329312      INTEGER                         , INTENT(in )  :: ktra  ! tracer index 
    330313      CHARACTER(len=3)                , INTENT(in)   :: cptr  ! transport type  'adv'/'ldf' 
    331       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in)   :: pua   ! 3D input array of advection/diffusion 
    332       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in)   :: pva   ! 3D input array of advection/diffusion 
     314      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in)   :: puflx  ! u-flux of advection/diffusion 
     315      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in)   :: pvflx  ! v-flux of advection/diffusion 
    333316      ! 
    334317      INTEGER    ::  ji, jj, jk 
     
    336319 
    337320     
    338       z2d(:,:) = pua(:,:,1)  
    339       DO jk = 1, jpkm1 
    340          DO jj = 2, jpjm1 
    341             DO ji = fs_2, fs_jpim1   ! vector opt. 
    342                z2d(ji,jj) = z2d(ji,jj) + pua(ji,jj,jk)  
    343             END DO 
    344          END DO 
    345        END DO 
     321      z2d(:,:) = puflx(:,:,1)  
     322      DO_3D_00_00( 1, jpkm1 ) 
     323         z2d(ji,jj) = z2d(ji,jj) + puflx(ji,jj,jk)  
     324      END_3D 
    346325       CALL lbc_lnk( 'diaar5', z2d, 'U', -1. ) 
    347326       IF( cptr == 'adv' ) THEN 
     
    354333       ENDIF 
    355334       ! 
    356        z2d(:,:) = pva(:,:,1)  
    357        DO jk = 1, jpkm1 
    358           DO jj = 2, jpjm1 
    359              DO ji = fs_2, fs_jpim1   ! vector opt. 
    360                 z2d(ji,jj) = z2d(ji,jj) + pva(ji,jj,jk)  
    361              END DO 
    362           END DO 
    363        END DO 
     335       z2d(:,:) = pvflx(:,:,1)  
     336       DO_3D_00_00( 1, jpkm1 ) 
     337          z2d(ji,jj) = z2d(ji,jj) + pvflx(ji,jj,jk)  
     338       END_3D 
    364339       CALL lbc_lnk( 'diaar5', z2d, 'V', -1. ) 
    365340       IF( cptr == 'adv' ) THEN 
     
    406381         zvol0 (:,:) = 0._wp 
    407382         thick0(:,:) = 0._wp 
    408          DO jk = 1, jpkm1 
    409             DO jj = 1, jpj               ! interpolation of salinity at the last ocean level (i.e. the partial step) 
    410                DO ji = 1, jpi 
    411                   idep = tmask(ji,jj,jk) * e3t_0(ji,jj,jk) 
    412                   zvol0 (ji,jj) = zvol0 (ji,jj) +  idep * area(ji,jj) 
    413                   thick0(ji,jj) = thick0(ji,jj) +  idep     
    414                END DO 
    415             END DO 
    416          END DO 
     383         DO_3D_11_11( 1, jpkm1 ) 
     384            idep = tmask(ji,jj,jk) * e3t_0(ji,jj,jk) 
     385            zvol0 (ji,jj) = zvol0 (ji,jj) +  idep * area(ji,jj) 
     386            thick0(ji,jj) = thick0(ji,jj) +  idep     
     387         END_3D 
    417388         vol0 = glob_sum( 'diaar5', zvol0 ) 
    418389         DEALLOCATE( zvol0 ) 
     
    428399            sn0(:,:,:) = sn0(:,:,:) * tmask(:,:,:) 
    429400            IF( ln_zps ) THEN               ! z-coord. partial steps 
    430                DO jj = 1, jpj               ! interpolation of salinity at the last ocean level (i.e. the partial step) 
    431                   DO ji = 1, jpi 
    432                      ik = mbkt(ji,jj) 
    433                      IF( ik > 1 ) THEN 
    434                         zztmp = ( gdept_1d(ik) - gdept_0(ji,jj,ik) ) / ( gdept_1d(ik) - gdept_1d(ik-1) ) 
    435                         sn0(ji,jj,ik) = ( 1._wp - zztmp ) * sn0(ji,jj,ik) + zztmp * sn0(ji,jj,ik-1) 
    436                      ENDIF 
    437                   END DO 
    438                END DO 
     401               DO_2D_11_11 
     402                  ik = mbkt(ji,jj) 
     403                  IF( ik > 1 ) THEN 
     404                     zztmp = ( gdept_1d(ik) - gdept_0(ji,jj,ik) ) / ( gdept_1d(ik) - gdept_1d(ik-1) ) 
     405                     sn0(ji,jj,ik) = ( 1._wp - zztmp ) * sn0(ji,jj,ik) + zztmp * sn0(ji,jj,ik-1) 
     406                  ENDIF 
     407               END_2D 
    439408            ENDIF 
    440409            ! 
  • NEMO/trunk/src/OCE/DIA/diacfl.F90

    r11532 r12377  
    3333 
    3434   !! * Substitutions 
    35 #  include "vectopt_loop_substitute.h90" 
     35#  include "do_loop_substitute.h90" 
    3636   !!---------------------------------------------------------------------- 
    3737   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
     
    4141CONTAINS 
    4242 
    43    SUBROUTINE dia_cfl ( kt ) 
     43   SUBROUTINE dia_cfl ( kt, Kmm ) 
    4444      !!---------------------------------------------------------------------- 
    4545      !!                  ***  ROUTINE dia_cfl  *** 
     
    4949      !!---------------------------------------------------------------------- 
    5050      INTEGER, INTENT(in) ::   kt   ! ocean time-step index 
     51      INTEGER, INTENT(in) ::   Kmm  ! ocean time level index 
    5152      ! 
    5253      INTEGER                          ::   ji, jj, jk                       ! dummy loop indices 
     
    6465      ! 
    6566      !                 
    66       DO jk = 1, jpk       ! calculate Courant numbers 
    67          DO jj = 1, jpj 
    68             DO ji = 1, jpi 
    69                zCu_cfl(ji,jj,jk) = ABS( un(ji,jj,jk) ) * z2dt / e1u  (ji,jj)      ! for i-direction 
    70                zCv_cfl(ji,jj,jk) = ABS( vn(ji,jj,jk) ) * z2dt / e2v  (ji,jj)      ! for j-direction 
    71                zCw_cfl(ji,jj,jk) = ABS( wn(ji,jj,jk) ) * z2dt / e3w_n(ji,jj,jk)   ! for k-direction 
    72             END DO 
    73          END DO          
    74       END DO 
     67      DO_3D_11_11( 1, jpk ) 
     68         zCu_cfl(ji,jj,jk) = ABS( uu(ji,jj,jk,Kmm) ) * z2dt / e1u  (ji,jj)      ! for i-direction 
     69         zCv_cfl(ji,jj,jk) = ABS( vv(ji,jj,jk,Kmm) ) * z2dt / e2v  (ji,jj)      ! for j-direction 
     70         zCw_cfl(ji,jj,jk) = ABS( ww(ji,jj,jk) ) * z2dt / e3w(ji,jj,jk,Kmm)   ! for k-direction 
     71      END_3D 
    7572      ! 
    7673      ! write outputs 
  • NEMO/trunk/src/OCE/DIA/diadct.F90

    r11536 r12377  
    123123      !!--------------------------------------------------------------------- 
    124124 
    125      REWIND( numnam_ref )              ! Namelist nam_diadct in reference namelist : Diagnostic: transport through sections 
    126125     READ  ( numnam_ref, nam_diadct, IOSTAT = ios, ERR = 901) 
    127126901  IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_diadct in reference namelist' ) 
    128127 
    129      REWIND( numnam_cfg )              ! Namelist nam_diadct in configuration namelist : Diagnostic: transport through sections 
    130128     READ  ( numnam_cfg, nam_diadct, IOSTAT = ios, ERR = 902 ) 
    131129902  IF( ios >  0 ) CALL ctl_nam ( ios , 'nam_diadct in configuration namelist' ) 
     
    175173  
    176174  
    177   SUBROUTINE dia_dct( kt ) 
     175  SUBROUTINE dia_dct( kt, Kmm ) 
    178176     !!--------------------------------------------------------------------- 
    179177     !!               ***  ROUTINE diadct  ***   
     
    192190     !!               Reinitialise all relevant arrays to zero  
    193191     !!--------------------------------------------------------------------- 
    194      INTEGER, INTENT(in) ::   kt 
     192     INTEGER, INTENT(in) ::   kt    ! ocean time step 
     193     INTEGER, INTENT(in) ::   Kmm   ! time level index 
    195194     ! 
    196195     INTEGER ::   jsec              ! loop on sections 
     
    232231 
    233232           !Compute transport through section   
    234            CALL transport(secs(jsec),lldebug,jsec)  
     233           CALL transport(Kmm,secs(jsec),lldebug,jsec)  
    235234 
    236235        ENDDO 
     
    246245           ! Sum over each class  
    247246           DO jsec=1,nb_sec  
    248               CALL dia_dct_sum(secs(jsec),jsec)  
     247              CALL dia_dct_sum(Kmm,secs(jsec),jsec)  
    249248           ENDDO  
    250249 
     
    558557 
    559558 
    560    SUBROUTINE transport(sec,ld_debug,jsec) 
     559   SUBROUTINE transport(Kmm,sec,ld_debug,jsec) 
    561560     !!------------------------------------------------------------------------------------------- 
    562561     !!                     ***  ROUTINE transport  *** 
     
    578577     !! 
    579578     !!------------------------------------------------------------------------------------------- 
     579     INTEGER      ,INTENT(IN)    :: Kmm         ! time level index 
    580580     TYPE(SECTION),INTENT(INOUT) :: sec 
    581581     LOGICAL      ,INTENT(IN)    :: ld_debug 
     
    673673            SELECT CASE( sec%direction(jseg) ) 
    674674               CASE(0,1)  
    675                   ztn   = interp(k%I,k%J,jk,'V',tsn(:,:,:,jp_tem) )  
    676                   zsn   = interp(k%I,k%J,jk,'V',tsn(:,:,:,jp_sal) )  
    677                   zrhop = interp(k%I,k%J,jk,'V',rhop)  
    678                   zrhoi = interp(k%I,k%J,jk,'V',rhd*rau0+rau0)  
    679                   zsshn =  0.5*( sshn(k%I,k%J) + sshn(k%I,k%J+1)    ) * vmask(k%I,k%J,1)  
     675                  ztn   = interp(Kmm,k%I,k%J,jk,'V',ts(:,:,:,jp_tem,Kmm) )  
     676                  zsn   = interp(Kmm,k%I,k%J,jk,'V',ts(:,:,:,jp_sal,Kmm) )  
     677                  zrhop = interp(Kmm,k%I,k%J,jk,'V',rhop)  
     678                  zrhoi = interp(Kmm,k%I,k%J,jk,'V',rhd*rau0+rau0)  
     679                  zsshn =  0.5*( ssh(k%I,k%J,Kmm) + ssh(k%I,k%J+1,Kmm)    ) * vmask(k%I,k%J,1)  
    680680               CASE(2,3)  
    681                   ztn   = interp(k%I,k%J,jk,'U',tsn(:,:,:,jp_tem) )  
    682                   zsn   = interp(k%I,k%J,jk,'U',tsn(:,:,:,jp_sal) )  
    683                   zrhop = interp(k%I,k%J,jk,'U',rhop)  
    684                   zrhoi = interp(k%I,k%J,jk,'U',rhd*rau0+rau0)  
    685                   zsshn =  0.5*( sshn(k%I,k%J) + sshn(k%I+1,k%J)    ) * umask(k%I,k%J,1)   
     681                  ztn   = interp(Kmm,k%I,k%J,jk,'U',ts(:,:,:,jp_tem,Kmm) )  
     682                  zsn   = interp(Kmm,k%I,k%J,jk,'U',ts(:,:,:,jp_sal,Kmm) )  
     683                  zrhop = interp(Kmm,k%I,k%J,jk,'U',rhop)  
     684                  zrhoi = interp(Kmm,k%I,k%J,jk,'U',rhd*rau0+rau0)  
     685                  zsshn =  0.5*( ssh(k%I,k%J,Kmm) + ssh(k%I+1,k%J,Kmm)    ) * umask(k%I,k%J,1)   
    686686               END SELECT  
    687687               ! 
    688                zdep= gdept_n(k%I,k%J,jk)  
     688               zdep= gdept(k%I,k%J,jk,Kmm)  
    689689   
    690690               SELECT CASE( sec%direction(jseg) )                !compute velocity with the correct direction  
    691691               CASE(0,1)    
    692692                  zumid=0._wp 
    693                   zvmid=isgnv*vn(k%I,k%J,jk)*vmask(k%I,k%J,jk)  
     693                  zvmid=isgnv*vv(k%I,k%J,jk,Kmm)*vmask(k%I,k%J,jk)  
    694694               CASE(2,3)  
    695                   zumid=isgnu*un(k%I,k%J,jk)*umask(k%I,k%J,jk)  
     695                  zumid=isgnu*uu(k%I,k%J,jk,Kmm)*umask(k%I,k%J,jk)  
    696696                  zvmid=0._wp 
    697697               END SELECT  
     
    699699               !zTnorm=transport through one cell;  
    700700               !velocity* cell's length * cell's thickness  
    701                zTnorm = zumid*e2u(k%I,k%J) * e3u_n(k%I,k%J,jk)     &  
    702                   &   + zvmid*e1v(k%I,k%J) * e3v_n(k%I,k%J,jk)  
     701               zTnorm = zumid*e2u(k%I,k%J) * e3u(k%I,k%J,jk,Kmm)     &  
     702                  &   + zvmid*e1v(k%I,k%J) * e3v(k%I,k%J,jk,Kmm)  
    703703 
    704704!!gm  THIS is WRONG  no transport due to ssh in linear free surface case !!!!! 
     
    765765 
    766766 
    767   SUBROUTINE dia_dct_sum(sec,jsec)  
     767  SUBROUTINE dia_dct_sum(Kmm,sec,jsec)  
    768768     !!-------------------------------------------------------------  
    769769     !! Purpose: Average the transport over nn_dctwri time steps   
     
    784784     !!  
    785785     !!-------------------------------------------------------------  
     786     INTEGER      ,INTENT(IN)    :: Kmm         ! time level index 
    786787     TYPE(SECTION),INTENT(INOUT) :: sec  
    787788     INTEGER      ,INTENT(IN)    :: jsec        ! numeric identifier of section  
     
    845846              SELECT CASE( sec%direction(jseg) )  
    846847              CASE(0,1)  
    847                  ztn   = interp(k%I,k%J,jk,'V',tsn(:,:,:,jp_tem) )  
    848                  zsn   = interp(k%I,k%J,jk,'V',tsn(:,:,:,jp_sal) )  
    849                  zrhop = interp(k%I,k%J,jk,'V',rhop)  
    850                  zrhoi = interp(k%I,k%J,jk,'V',rhd*rau0+rau0)  
     848                 ztn   = interp(Kmm,k%I,k%J,jk,'V',ts(:,:,:,jp_tem,Kmm) )  
     849                 zsn   = interp(Kmm,k%I,k%J,jk,'V',ts(:,:,:,jp_sal,Kmm) )  
     850                 zrhop = interp(Kmm,k%I,k%J,jk,'V',rhop)  
     851                 zrhoi = interp(Kmm,k%I,k%J,jk,'V',rhd*rau0+rau0)  
    851852 
    852853              CASE(2,3)  
    853                  ztn   = interp(k%I,k%J,jk,'U',tsn(:,:,:,jp_tem) )  
    854                  zsn   = interp(k%I,k%J,jk,'U',tsn(:,:,:,jp_sal) )  
    855                  zrhop = interp(k%I,k%J,jk,'U',rhop)  
    856                  zrhoi = interp(k%I,k%J,jk,'U',rhd*rau0+rau0)  
    857                  zsshn =  0.5*( sshn(k%I,k%J)    + sshn(k%I+1,k%J)    ) * umask(k%I,k%J,1)   
     854                 ztn   = interp(Kmm,k%I,k%J,jk,'U',ts(:,:,:,jp_tem,Kmm) )  
     855                 zsn   = interp(Kmm,k%I,k%J,jk,'U',ts(:,:,:,jp_sal,Kmm) )  
     856                 zrhop = interp(Kmm,k%I,k%J,jk,'U',rhop)  
     857                 zrhoi = interp(Kmm,k%I,k%J,jk,'U',rhd*rau0+rau0)  
     858                 zsshn =  0.5*( ssh(k%I,k%J,Kmm)    + ssh(k%I+1,k%J,Kmm)    ) * umask(k%I,k%J,1)   
    858859              END SELECT  
    859860  
    860               zdep= gdept_n(k%I,k%J,jk)  
     861              zdep= gdept(k%I,k%J,jk,Kmm)  
    861862   
    862863              !-------------------------------  
     
    11011102 
    11021103 
    1103    FUNCTION interp(ki, kj, kk, cd_point, ptab) 
     1104   FUNCTION interp(Kmm, ki, kj, kk, cd_point, ptab) 
    11041105  !!---------------------------------------------------------------------- 
    11051106  !! 
     
    11621163  !!---------------------------------------------------------------------- 
    11631164  !*arguments 
     1165  INTEGER, INTENT(IN)                          :: Kmm          ! time level index 
    11641166  INTEGER, INTENT(IN)                          :: ki, kj, kk   ! coordinate of point 
    11651167  CHARACTER(len=1), INTENT(IN)                 :: cd_point     ! type of point (U, V) 
     
    11961198  IF( ln_sco )THEN   ! s-coordinate case 
    11971199 
    1198      zdepu = ( gdept_n(ii1,ij1,kk) +  gdept_n(ii2,ij2,kk) ) * 0.5_wp  
    1199      zdep1 = gdept_n(ii1,ij1,kk) - zdepu 
    1200      zdep2 = gdept_n(ii2,ij2,kk) - zdepu 
     1200     zdepu = ( gdept(ii1,ij1,kk,Kmm) +  gdept(ii2,ij2,kk,Kmm) ) * 0.5_wp  
     1201     zdep1 = gdept(ii1,ij1,kk,Kmm) - zdepu 
     1202     zdep2 = gdept(ii2,ij2,kk,Kmm) - zdepu 
    12011203 
    12021204     ! weights 
     
    12101212  ELSE       ! full step or partial step case  
    12111213 
    1212      ze3t  = e3t_n(ii2,ij2,kk) - e3t_n(ii1,ij1,kk)  
    1213      zwgt1 = ( e3w_n(ii2,ij2,kk) - e3w_n(ii1,ij1,kk) ) / e3w_n(ii2,ij2,kk) 
    1214      zwgt2 = ( e3w_n(ii1,ij1,kk) - e3w_n(ii2,ij2,kk) ) / e3w_n(ii1,ij1,kk) 
     1214     ze3t  = e3t(ii2,ij2,kk,Kmm) - e3t(ii1,ij1,kk,Kmm)  
     1215     zwgt1 = ( e3w(ii2,ij2,kk,Kmm) - e3w(ii1,ij1,kk,Kmm) ) / e3w(ii2,ij2,kk,Kmm) 
     1216     zwgt2 = ( e3w(ii1,ij1,kk,Kmm) - e3w(ii2,ij2,kk,Kmm) ) / e3w(ii1,ij1,kk,Kmm) 
    12151217 
    12161218     IF(kk .NE. 1)THEN 
     
    12451247      IMPLICIT NONE 
    12461248   END SUBROUTINE dia_dct_init 
    1247    SUBROUTINE dia_dct( kt ) 
     1249 
     1250   SUBROUTINE dia_dct( kt, Kmm )         ! Dummy routine 
    12481251      IMPLICIT NONE 
    1249       INTEGER, INTENT(in) ::   kt 
     1252      INTEGER, INTENT( in ) :: kt   ! ocean time-step index 
     1253      INTEGER, INTENT( in ) :: Kmm  ! ocean time level index 
     1254      WRITE(*,*) 'dia_dct: You should not have seen this print! error?', kt 
    12501255   END SUBROUTINE dia_dct 
    12511256   ! 
  • NEMO/trunk/src/OCE/DIA/diahsb.F90

    r11536 r12377  
    1717   USE phycst         ! physical constants 
    1818   USE sbc_oce        ! surface thermohaline fluxes 
     19   USE isf_oce        ! ice shelf fluxes 
    1920   USE sbcrnf         ! river runoff 
    20    USE sbcisf         ! ice shelves 
    2121   USE domvvl         ! vertical scale factors 
    2222   USE traqsr         ! penetrative solar radiation 
     
    4848   REAL(wp), DIMENSION(:,:)  , ALLOCATABLE ::   ssh_hc_loc_ini, ssh_sc_loc_ini   ! 
    4949   REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   hc_loc_ini, sc_loc_ini, e3t_ini  ! 
    50  
    51    !! * Substitutions 
    52 #  include "vectopt_loop_substitute.h90" 
     50   REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   tmask_ini 
     51 
    5352   !!---------------------------------------------------------------------- 
    5453   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
     
    5857CONTAINS 
    5958 
    60    SUBROUTINE dia_hsb( kt ) 
     59   SUBROUTINE dia_hsb( kt, Kbb, Kmm ) 
    6160      !!--------------------------------------------------------------------------- 
    6261      !!                  ***  ROUTINE dia_hsb  *** 
     
    6968      !! 
    7069      !!--------------------------------------------------------------------------- 
    71       INTEGER, INTENT(in) ::   kt   ! ocean time-step index 
     70      INTEGER, INTENT(in) ::   kt         ! ocean time-step index 
     71      INTEGER, INTENT(in) ::   Kbb, Kmm   ! ocean time level indices 
    7272      ! 
    7373      INTEGER    ::   ji, jj, jk                  ! dummy loop indice 
     
    8686      IF( ln_timing )   CALL timing_start('dia_hsb')       
    8787      ! 
    88       tsn(:,:,:,1) = tsn(:,:,:,1) * tmask(:,:,:) ; tsb(:,:,:,1) = tsb(:,:,:,1) * tmask(:,:,:) ; 
    89       tsn(:,:,:,2) = tsn(:,:,:,2) * tmask(:,:,:) ; tsb(:,:,:,2) = tsb(:,:,:,2) * tmask(:,:,:) ; 
     88      ts(:,:,:,1,Kmm) = ts(:,:,:,1,Kmm) * tmask(:,:,:) ; ts(:,:,:,1,Kbb) = ts(:,:,:,1,Kbb) * tmask(:,:,:) ; 
     89      ts(:,:,:,2,Kmm) = ts(:,:,:,2,Kmm) * tmask(:,:,:) ; ts(:,:,:,2,Kbb) = ts(:,:,:,2,Kbb) * tmask(:,:,:) ; 
    9090      ! ------------------------- ! 
    9191      ! 1 - Trends due to forcing ! 
    9292      ! ------------------------- ! 
    93       z_frc_trd_v = r1_rau0 * glob_sum( 'diahsb', - ( emp(:,:) - rnf(:,:) + fwfisf(:,:) ) * surf(:,:) )   ! volume fluxes 
     93      z_frc_trd_v = r1_rau0 * glob_sum( 'diahsb', - ( emp(:,:) - rnf(:,:) + fwfisf_cav(:,:) + fwfisf_par(:,:) ) * surf(:,:) )   ! volume fluxes 
    9494      z_frc_trd_t =           glob_sum( 'diahsb', sbc_tsc(:,:,jp_tem) * surf(:,:) )                       ! heat fluxes 
    9595      z_frc_trd_s =           glob_sum( 'diahsb', sbc_tsc(:,:,jp_sal) * surf(:,:) )                       ! salt fluxes 
     
    9898      IF( ln_rnf_sal)   z_frc_trd_s = z_frc_trd_s + glob_sum( 'diahsb', rnf_tsc(:,:,jp_sal) * surf(:,:) ) 
    9999      !                    ! Add ice shelf heat & salt input 
    100       IF( ln_isf    )   z_frc_trd_t = z_frc_trd_t + glob_sum( 'diahsb', risf_tsc(:,:,jp_tem) * surf(:,:) ) 
     100      IF( ln_isf    )   z_frc_trd_t = z_frc_trd_t & 
     101         &                          + glob_sum( 'diahsb', ( risf_cav_tsc(:,:,jp_tem) + risf_par_tsc(:,:,jp_tem) ) * surf(:,:) ) 
    101102      !                    ! Add penetrative solar radiation 
    102103      IF( ln_traqsr )   z_frc_trd_t = z_frc_trd_t + r1_rau0_rcp * glob_sum( 'diahsb', qsr     (:,:) * surf(:,:) ) 
     
    108109            DO ji=1,jpi 
    109110               DO jj=1,jpj 
    110                   z2d0(ji,jj) = surf(ji,jj) * wn(ji,jj,mikt(ji,jj)) * tsb(ji,jj,mikt(ji,jj),jp_tem) 
    111                   z2d1(ji,jj) = surf(ji,jj) * wn(ji,jj,mikt(ji,jj)) * tsb(ji,jj,mikt(ji,jj),jp_sal) 
     111                  z2d0(ji,jj) = surf(ji,jj) * ww(ji,jj,mikt(ji,jj)) * ts(ji,jj,mikt(ji,jj),jp_tem,Kbb) 
     112                  z2d1(ji,jj) = surf(ji,jj) * ww(ji,jj,mikt(ji,jj)) * ts(ji,jj,mikt(ji,jj),jp_sal,Kbb) 
    112113               END DO 
    113114            END DO 
    114115         ELSE 
    115             z2d0(:,:) = surf(:,:) * wn(:,:,1) * tsb(:,:,1,jp_tem) 
    116             z2d1(:,:) = surf(:,:) * wn(:,:,1) * tsb(:,:,1,jp_sal) 
     116            z2d0(:,:) = surf(:,:) * ww(:,:,1) * ts(:,:,1,jp_tem,Kbb) 
     117            z2d1(:,:) = surf(:,:) * ww(:,:,1) * ts(:,:,1,jp_sal,Kbb) 
    117118         END IF 
    118119         z_wn_trd_t = - glob_sum( 'diahsb', z2d0 )  
     
    135136 
    136137      !                    ! volume variation (calculated with ssh) 
    137       zdiff_v1 = glob_sum_full( 'diahsb', surf(:,:)*sshn(:,:) - surf_ini(:,:)*ssh_ini(:,:) ) 
     138      zdiff_v1 = glob_sum_full( 'diahsb', surf(:,:)*ssh(:,:,Kmm) - surf_ini(:,:)*ssh_ini(:,:) ) 
    138139 
    139140      !                    ! heat & salt content variation (associated with ssh) 
     
    142143            DO ji = 1, jpi 
    143144               DO jj = 1, jpj 
    144                   z2d0(ji,jj) = surf(ji,jj) * ( tsn(ji,jj,mikt(ji,jj),jp_tem) * sshn(ji,jj) - ssh_hc_loc_ini(ji,jj) )  
    145                   z2d1(ji,jj) = surf(ji,jj) * ( tsn(ji,jj,mikt(ji,jj),jp_sal) * sshn(ji,jj) - ssh_sc_loc_ini(ji,jj) )  
     145                  z2d0(ji,jj) = surf(ji,jj) * ( ts(ji,jj,mikt(ji,jj),jp_tem,Kmm) * ssh(ji,jj,Kmm) - ssh_hc_loc_ini(ji,jj) )  
     146                  z2d1(ji,jj) = surf(ji,jj) * ( ts(ji,jj,mikt(ji,jj),jp_sal,Kmm) * ssh(ji,jj,Kmm) - ssh_sc_loc_ini(ji,jj) )  
    146147               END DO 
    147148            END DO 
    148149         ELSE                          ! no under ice-shelf seas 
    149             z2d0(:,:) = surf(:,:) * ( tsn(:,:,1,jp_tem) * sshn(:,:) - ssh_hc_loc_ini(:,:) )  
    150             z2d1(:,:) = surf(:,:) * ( tsn(:,:,1,jp_sal) * sshn(:,:) - ssh_sc_loc_ini(:,:) )  
     150            z2d0(:,:) = surf(:,:) * ( ts(:,:,1,jp_tem,Kmm) * ssh(:,:,Kmm) - ssh_hc_loc_ini(:,:) )  
     151            z2d1(:,:) = surf(:,:) * ( ts(:,:,1,jp_sal,Kmm) * ssh(:,:,Kmm) - ssh_sc_loc_ini(:,:) )  
    151152         END IF 
    152153         z_ssh_hc = glob_sum_full( 'diahsb', z2d0 )  
     
    155156      ! 
    156157      DO jk = 1, jpkm1           ! volume variation (calculated with scale factors) 
    157          zwrk(:,:,jk) = ( surf(:,:)*e3t_n(:,:,jk) - surf_ini(:,:)*e3t_ini(:,:,jk) ) * tmask(:,:,jk) 
     158         zwrk(:,:,jk) = surf(:,:)*e3t(:,:,jk,Kmm)*tmask(:,:,jk) - surf_ini(:,:)*e3t_ini(:,:,jk)*tmask_ini(:,:,jk) 
    158159      END DO 
    159       zdiff_v2 = glob_sum_full( 'diahsb', zwrk(:,:,:) ) 
     160      zdiff_v2 = glob_sum_full( 'diahsb', zwrk(:,:,:) )     ! glob_sum_full needed as tmask and tmask_ini could be different 
    160161      DO jk = 1, jpkm1           ! heat content variation 
    161          zwrk(:,:,jk) = ( surf(:,:)*e3t_n(:,:,jk)*tsn(:,:,jk,jp_tem) - surf_ini(:,:)*hc_loc_ini(:,:,jk) ) * tmask(:,:,jk) 
     162         zwrk(:,:,jk) = ( surf(:,:)*e3t(:,:,jk,Kmm)*ts(:,:,jk,jp_tem,Kmm) - surf_ini(:,:)*hc_loc_ini(:,:,jk) ) 
    162163      END DO 
    163164      zdiff_hc = glob_sum_full( 'diahsb', zwrk(:,:,:) ) 
    164165      DO jk = 1, jpkm1           ! salt content variation 
    165          zwrk(:,:,jk) = ( surf(:,:)*e3t_n(:,:,jk)*tsn(:,:,jk,jp_sal) - surf_ini(:,:)*sc_loc_ini(:,:,jk) ) * tmask(:,:,jk) 
     166         zwrk(:,:,jk) = ( surf(:,:)*e3t(:,:,jk,Kmm)*ts(:,:,jk,jp_sal,Kmm) - surf_ini(:,:)*sc_loc_ini(:,:,jk) ) 
    166167      END DO 
    167168      zdiff_sc = glob_sum_full( 'diahsb', zwrk(:,:,:) ) 
     
    185186      ! ----------------------- ! 
    186187      DO jk = 1, jpkm1           ! total ocean volume (calculated with scale factors) 
    187          zwrk(:,:,jk) = surf(:,:) * e3t_n(:,:,jk) * tmask(:,:,jk) 
     188         zwrk(:,:,jk) = surf(:,:) * e3t(:,:,jk,Kmm) * tmask(:,:,jk) 
    188189      END DO 
    189       zvol_tot = glob_sum_full( 'diahsb', zwrk(:,:,:) ) 
     190      zvol_tot = glob_sum( 'diahsb', zwrk(:,:,:) ) 
    190191 
    191192!!gm to be added ? 
    192193!      IF( ln_linssh ) THEN            ! fixed volume, add the ssh contribution 
    193 !        zvol_tot = zvol_tot + glob_sum( 'diahsb', surf(:,:) * sshn(:,:) ) 
     194!        zvol_tot = zvol_tot + glob_sum( 'diahsb', surf(:,:) * ssh(:,:,Kmm) ) 
    194195!      ENDIF 
    195196!!gm end 
     
    233234      ENDIF 
    234235      ! 
    235       IF( lrst_oce )   CALL dia_hsb_rst( kt, 'WRITE' ) 
     236      IF( lrst_oce )   CALL dia_hsb_rst( kt, Kmm, 'WRITE' ) 
    236237      ! 
    237238      IF( ln_timing )   CALL timing_stop('dia_hsb') 
     
    240241 
    241242 
    242    SUBROUTINE dia_hsb_rst( kt, cdrw ) 
     243   SUBROUTINE dia_hsb_rst( kt, Kmm, cdrw ) 
    243244      !!--------------------------------------------------------------------- 
    244245      !!                   ***  ROUTINE dia_hsb_rst  *** 
     
    249250      !!---------------------------------------------------------------------- 
    250251      INTEGER         , INTENT(in) ::   kt     ! ocean time-step 
     252      INTEGER         , INTENT(in) ::   Kmm    ! ocean time level index 
    251253      CHARACTER(len=*), INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag 
    252254      ! 
     
    270272            CALL iom_get( numror, jpdom_autoglo, 'ssh_ini'   , ssh_ini   , ldxios = lrxios ) 
    271273            CALL iom_get( numror, jpdom_autoglo, 'e3t_ini'   , e3t_ini   , ldxios = lrxios ) 
     274            CALL iom_get( numror, jpdom_autoglo, 'tmask_ini' , tmask_ini , ldxios = lrxios ) 
    272275            CALL iom_get( numror, jpdom_autoglo, 'hc_loc_ini', hc_loc_ini, ldxios = lrxios ) 
    273276            CALL iom_get( numror, jpdom_autoglo, 'sc_loc_ini', sc_loc_ini, ldxios = lrxios ) 
     
    281284            IF(lwp) WRITE(numout,*) 
    282285            surf_ini(:,:) = e1e2t(:,:) * tmask_i(:,:)         ! initial ocean surface 
    283             ssh_ini(:,:) = sshn(:,:)                          ! initial ssh 
     286            ssh_ini(:,:) = ssh(:,:,Kmm)                          ! initial ssh 
    284287            DO jk = 1, jpk 
    285288              ! if ice sheet/oceqn coupling, need to mask ini variables here (mask could change at the next NEMO instance). 
    286                e3t_ini   (:,:,jk) = e3t_n(:,:,jk)                      * tmask(:,:,jk)  ! initial vertical scale factors 
    287                hc_loc_ini(:,:,jk) = tsn(:,:,jk,jp_tem) * e3t_n(:,:,jk) * tmask(:,:,jk)  ! initial heat content 
    288                sc_loc_ini(:,:,jk) = tsn(:,:,jk,jp_sal) * e3t_n(:,:,jk) * tmask(:,:,jk)  ! initial salt content 
     289               e3t_ini   (:,:,jk) = e3t(:,:,jk,Kmm)                      * tmask(:,:,jk)  ! initial vertical scale factors 
     290               tmask_ini (:,:,jk) = tmask(:,:,jk)                                       ! initial mask 
     291               hc_loc_ini(:,:,jk) = ts(:,:,jk,jp_tem,Kmm) * e3t(:,:,jk,Kmm) * tmask(:,:,jk)  ! initial heat content 
     292               sc_loc_ini(:,:,jk) = ts(:,:,jk,jp_sal,Kmm) * e3t(:,:,jk,Kmm) * tmask(:,:,jk)  ! initial salt content 
    289293            END DO 
    290294            frc_v = 0._wp                                           ! volume       trend due to forcing 
     
    295299                  DO ji = 1, jpi 
    296300                     DO jj = 1, jpj 
    297                         ssh_hc_loc_ini(ji,jj) = tsn(ji,jj,mikt(ji,jj),jp_tem) * sshn(ji,jj)   ! initial heat content in ssh 
    298                         ssh_sc_loc_ini(ji,jj) = tsn(ji,jj,mikt(ji,jj),jp_sal) * sshn(ji,jj)   ! initial salt content in ssh 
     301                        ssh_hc_loc_ini(ji,jj) = ts(ji,jj,mikt(ji,jj),jp_tem,Kmm) * ssh(ji,jj,Kmm)   ! initial heat content in ssh 
     302                        ssh_sc_loc_ini(ji,jj) = ts(ji,jj,mikt(ji,jj),jp_sal,Kmm) * ssh(ji,jj,Kmm)   ! initial salt content in ssh 
    299303                     END DO 
    300304                   END DO 
    301305                ELSE 
    302                   ssh_hc_loc_ini(:,:) = tsn(:,:,1,jp_tem) * sshn(:,:)   ! initial heat content in ssh 
    303                   ssh_sc_loc_ini(:,:) = tsn(:,:,1,jp_sal) * sshn(:,:)   ! initial salt content in ssh 
     306                  ssh_hc_loc_ini(:,:) = ts(:,:,1,jp_tem,Kmm) * ssh(:,:,Kmm)   ! initial heat content in ssh 
     307                  ssh_sc_loc_ini(:,:) = ts(:,:,1,jp_sal,Kmm) * ssh(:,:,Kmm)   ! initial salt content in ssh 
    304308               END IF 
    305309               frc_wn_t = 0._wp                                       ! initial heat content misfit due to free surface 
     
    325329         CALL iom_rstput( kt, nitrst, numrow, 'ssh_ini'   , ssh_ini   , ldxios = lwxios ) 
    326330         CALL iom_rstput( kt, nitrst, numrow, 'e3t_ini'   , e3t_ini   , ldxios = lwxios ) 
     331         CALL iom_rstput( kt, nitrst, numrow, 'tmask_ini' , tmask_ini , ldxios = lwxios ) 
    327332         CALL iom_rstput( kt, nitrst, numrow, 'hc_loc_ini', hc_loc_ini, ldxios = lwxios ) 
    328333         CALL iom_rstput( kt, nitrst, numrow, 'sc_loc_ini', sc_loc_ini, ldxios = lwxios ) 
     
    338343 
    339344 
    340    SUBROUTINE dia_hsb_init 
     345   SUBROUTINE dia_hsb_init( Kmm ) 
    341346      !!--------------------------------------------------------------------------- 
    342347      !!                  ***  ROUTINE dia_hsb  *** 
     
    350355      !!             - Compute coefficients for conversion 
    351356      !!--------------------------------------------------------------------------- 
     357      INTEGER, INTENT(in) :: Kmm ! time level index 
     358      ! 
    352359      INTEGER ::   ierror, ios   ! local integer 
    353360      !! 
     
    360367         WRITE(numout,*) '~~~~~~~~~~~~ ' 
    361368      ENDIF 
    362       REWIND( numnam_ref )              ! Namelist namhsb in reference namelist 
    363369      READ  ( numnam_ref, namhsb, IOSTAT = ios, ERR = 901) 
    364370901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namhsb in reference namelist' ) 
    365       REWIND( numnam_cfg )              ! Namelist namhsb in configuration namelist 
    366371      READ  ( numnam_cfg, namhsb, IOSTAT = ios, ERR = 902 ) 
    367372902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namhsb in configuration namelist' ) 
     
    396401      ! ------------------- ! 
    397402      ALLOCATE( hc_loc_ini(jpi,jpj,jpk), sc_loc_ini(jpi,jpj,jpk), surf_ini(jpi,jpj), & 
    398          &      e3t_ini(jpi,jpj,jpk), surf(jpi,jpj),  ssh_ini(jpi,jpj), STAT=ierror  ) 
     403         &      e3t_ini(jpi,jpj,jpk), surf(jpi,jpj),  ssh_ini(jpi,jpj), tmask_ini(jpi,jpj,jpk),STAT=ierror  ) 
    399404      IF( ierror > 0 ) THEN 
    400405         CALL ctl_stop( 'dia_hsb_init: unable to allocate hc_loc_ini' )   ;   RETURN 
     
    417422      ! 4 - initial conservation variables ! 
    418423      ! ---------------------------------- ! 
    419       CALL dia_hsb_rst( nit000, 'READ' )  !* read or initialize all required files 
     424      CALL dia_hsb_rst( nit000, Kmm, 'READ' )  !* read or initialize all required files 
    420425      ! 
    421426   END SUBROUTINE dia_hsb_init 
  • NEMO/trunk/src/OCE/DIA/diahth.F90

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

    r12276 r12377  
    4646   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   hstr_ove, hstr_btr, hstr_vtr   !: heat Salt TRansports(overturn, baro, merional) 
    4747 
    48    LOGICAL, PUBLIC ::   ln_diaptr   !  Poleward transport flag (T) or not (F) 
    49    LOGICAL, PUBLIC ::   ln_subbas   !  Atlantic/Pacific/Indian basins calculation 
     48   LOGICAL , PUBLIC ::   l_diaptr        !: tracers  trend flag (set from namelist in trdini) 
    5049   INTEGER, PARAMETER, PUBLIC ::   nptr = 5  ! (glo, atl, pac, ind, ipc) 
    5150 
     
    6059   REAL(wp), TARGET, ALLOCATABLE, SAVE, DIMENSION(:,:) :: p_fval2d 
    6160 
     61   LOGICAL ::   ll_init = .TRUE.        !: tracers  trend flag (set from namelist in trdini) 
    6262   !! * Substitutions 
    63 #  include "vectopt_loop_substitute.h90" 
     63#  include "do_loop_substitute.h90" 
    6464   !!---------------------------------------------------------------------- 
    6565   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
     
    6969CONTAINS 
    7070 
    71    SUBROUTINE dia_ptr( pvtr ) 
     71   SUBROUTINE dia_ptr( kt, Kmm, pvtr ) 
    7272      !!---------------------------------------------------------------------- 
    7373      !!                  ***  ROUTINE dia_ptr  *** 
    7474      !!---------------------------------------------------------------------- 
     75      INTEGER                         , INTENT(in)           ::   kt     ! ocean time-step index      
     76      INTEGER                         , INTENT(in)           ::   Kmm    ! time level index 
    7577      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in), OPTIONAL ::   pvtr   ! j-effective transport 
    7678      ! 
     
    9294      ! 
    9395      IF( ln_timing )   CALL timing_start('dia_ptr') 
    94       ! 
     96 
     97      IF( kt == nit000 .AND. ll_init )   CALL dia_ptr_init 
     98      ! 
     99      IF( .NOT. l_diaptr )   RETURN 
     100 
    95101      IF( PRESENT( pvtr ) ) THEN 
    96102         IF( iom_use( 'zomsf' ) ) THEN    ! effective MSF 
     
    111117            zmask(:,:,:) = 0._wp 
    112118            zts(:,:,:,:) = 0._wp 
    113             DO jk = 1, jpkm1 
    114                DO jj = 1, jpjm1 
    115                   DO ji = 1, jpi 
    116                      zvfc = e1v(ji,jj) * e3v_n(ji,jj,jk) 
    117                      zmask(ji,jj,jk)      = vmask(ji,jj,jk)      * zvfc 
    118                      zts(ji,jj,jk,jp_tem) = (tsn(ji,jj,jk,jp_tem)+tsn(ji,jj+1,jk,jp_tem)) * 0.5 * zvfc  !Tracers averaged onto V grid 
    119                      zts(ji,jj,jk,jp_sal) = (tsn(ji,jj,jk,jp_sal)+tsn(ji,jj+1,jk,jp_sal)) * 0.5 * zvfc 
    120                   ENDDO 
    121                ENDDO 
    122              ENDDO 
     119            DO_3D_10_11( 1, jpkm1 ) 
     120               zvfc = e1v(ji,jj) * e3v(ji,jj,jk,Kmm) 
     121               zmask(ji,jj,jk)      = vmask(ji,jj,jk)      * zvfc 
     122               zts(ji,jj,jk,jp_tem) = (ts(ji,jj,jk,jp_tem,Kmm)+ts(ji,jj+1,jk,jp_tem,Kmm)) * 0.5 * zvfc  !Tracers averaged onto V grid 
     123               zts(ji,jj,jk,jp_sal) = (ts(ji,jj,jk,jp_sal,Kmm)+ts(ji,jj+1,jk,jp_sal,Kmm)) * 0.5 * zvfc 
     124            END_3D 
    123125         ENDIF 
    124126         IF( iom_use( 'sopstove' ) .OR. iom_use( 'sophtove' ) ) THEN 
     
    186188         zts(:,:,:,:) = 0._wp 
    187189         IF( iom_use( 'zotem' ) .OR. iom_use( 'zosal' ) .OR. iom_use( 'zosrf' )  ) THEN    ! i-mean i-k-surface  
    188             DO jk = 1, jpkm1 
    189                DO jj = 1, jpj 
    190                   DO ji = 1, jpi 
    191                      zsfc = e1t(ji,jj) * e3t_n(ji,jj,jk) 
    192                      zmask(ji,jj,jk)      = tmask(ji,jj,jk)      * zsfc 
    193                      zts(ji,jj,jk,jp_tem) = tsn(ji,jj,jk,jp_tem) * zsfc 
    194                      zts(ji,jj,jk,jp_sal) = tsn(ji,jj,jk,jp_sal) * zsfc 
    195                   END DO 
    196                END DO 
    197             END DO 
     190            DO_3D_11_11( 1, jpkm1 ) 
     191               zsfc = e1t(ji,jj) * e3t(ji,jj,jk,Kmm) 
     192               zmask(ji,jj,jk)      = tmask(ji,jj,jk)      * zsfc 
     193               zts(ji,jj,jk,jp_tem) = ts(ji,jj,jk,jp_tem,Kmm) * zsfc 
     194               zts(ji,jj,jk,jp_sal) = ts(ji,jj,jk,jp_sal,Kmm) * zsfc 
     195            END_3D 
    198196            ! 
    199197            DO jn = 1, nptr 
     
    280278         IF( iom_use( 'sopstvtr' ) .OR. iom_use( 'sophtvtr' ) ) THEN 
    281279            zts(:,:,:,:) = 0._wp 
    282             DO jk = 1, jpkm1 
    283                DO jj = 1, jpjm1 
    284                   DO ji = 1, jpi 
    285                      zvfc = e1v(ji,jj) * e3v_n(ji,jj,jk) 
    286                      zts(ji,jj,jk,jp_tem) = (tsn(ji,jj,jk,jp_tem)+tsn(ji,jj+1,jk,jp_tem)) * 0.5 * zvfc  !Tracers averaged onto V grid 
    287                      zts(ji,jj,jk,jp_sal) = (tsn(ji,jj,jk,jp_sal)+tsn(ji,jj+1,jk,jp_sal)) * 0.5 * zvfc 
    288                   ENDDO 
    289                ENDDO 
    290              ENDDO 
     280            DO_3D_10_11( 1, jpkm1 ) 
     281               zvfc = e1v(ji,jj) * e3v(ji,jj,jk,Kmm) 
     282               zts(ji,jj,jk,jp_tem) = (ts(ji,jj,jk,jp_tem,Kmm)+ts(ji,jj+1,jk,jp_tem,Kmm)) * 0.5 * zvfc  !Tracers averaged onto V grid 
     283               zts(ji,jj,jk,jp_sal) = (ts(ji,jj,jk,jp_sal,Kmm)+ts(ji,jj+1,jk,jp_sal,Kmm)) * 0.5 * zvfc 
     284            END_3D 
    291285             CALL dia_ptr_hst( jp_tem, 'vtr', zts(:,:,:,jp_tem) ) 
    292286             CALL dia_ptr_hst( jp_sal, 'vtr', zts(:,:,:,jp_sal) ) 
     
    326320      !! ** Purpose :   Initialization, namelist read 
    327321      !!---------------------------------------------------------------------- 
    328       INTEGER ::  inum, jn, ios, ierr           ! local integers 
    329       !! 
    330       NAMELIST/namptr/ ln_diaptr, ln_subbas 
     322      INTEGER ::  inum, jn           ! local integers 
     323      !! 
    331324      REAL(wp), DIMENSION(jpi,jpj) :: zmsk 
    332325      !!---------------------------------------------------------------------- 
    333326 
    334  
    335       REWIND( numnam_ref )              ! Namelist namptr in reference namelist : Poleward transport 
    336       READ  ( numnam_ref, namptr, IOSTAT = ios, ERR = 901) 
    337 901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namptr in reference namelist' ) 
    338  
    339       REWIND( numnam_cfg )              ! Namelist namptr in configuration namelist : Poleward transport 
    340       READ  ( numnam_cfg, namptr, IOSTAT = ios, ERR = 902 ) 
    341 902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namptr in configuration namelist' ) 
    342       IF(lwm) WRITE ( numond, namptr ) 
    343  
     327      l_diaptr = .FALSE. 
     328      IF(   iom_use( 'zomsf'    ) .OR. iom_use( 'zotem'    ) .OR. iom_use( 'zosal'    ) .OR.  & 
     329         &  iom_use( 'zosrf'    ) .OR. iom_use( 'sopstove' ) .OR. iom_use( 'sophtove' ) .OR.  & 
     330         &  iom_use( 'sopstbtr' ) .OR. iom_use( 'sophtbtr' ) .OR. iom_use( 'sophtadv' ) .OR.  & 
     331         &  iom_use( 'sopstadv' ) .OR. iom_use( 'sophtldf' ) .OR. iom_use( 'sopstldf' ) .OR.  &  
     332         &  iom_use( 'sophteiv' ) .OR. iom_use( 'sopsteiv' ) .OR. iom_use( 'sopstvtr' ) .OR.  & 
     333         &  iom_use( 'sophtvtr' ) .OR. iom_use( 'uocetr_vsum_cumul' ) )  l_diaptr  = .TRUE. 
     334 
     335  
    344336      IF(lwp) THEN                     ! Control print 
    345337         WRITE(numout,*) 
     
    347339         WRITE(numout,*) '~~~~~~~~~~~~' 
    348340         WRITE(numout,*) '   Namelist namptr : set ptr parameters' 
    349          WRITE(numout,*) '      Poleward heat & salt transport (T) or not (F)      ln_diaptr  = ', ln_diaptr 
     341         WRITE(numout,*) '      Poleward heat & salt transport (T) or not (F)      l_diaptr  = ', l_diaptr 
    350342      ENDIF 
    351343 
    352       IF( ln_diaptr ) THEN   
     344      IF( l_diaptr ) THEN   
    353345         ! 
    354346         IF( dia_ptr_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'dia_ptr_init : unable to allocate arrays' ) 
     
    389381         hstr_vtr(:,:,:) = 0._wp           ! 
    390382         ! 
     383         ll_init = .FALSE. 
     384         ! 
    391385      ENDIF  
    392386      !  
     
    394388 
    395389 
    396    SUBROUTINE dia_ptr_hst( ktra, cptr, pva )  
     390   SUBROUTINE dia_ptr_hst( ktra, cptr, pvflx )  
    397391      !!---------------------------------------------------------------------- 
    398392      !!                    ***  ROUTINE dia_ptr_hst *** 
     
    403397      INTEGER                         , INTENT(in )  :: ktra  ! tracer index 
    404398      CHARACTER(len=3)                , INTENT(in)   :: cptr  ! transport type  'adv'/'ldf'/'eiv' 
    405       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in)   :: pva   ! 3D input array of advection/diffusion 
     399      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in)   :: pvflx   ! 3D input array of advection/diffusion 
    406400      INTEGER                                        :: jn    ! 
    407401 
     
    410404         IF( ktra == jp_tem )  THEN 
    411405             DO jn = 1, nptr 
    412                 hstr_adv(:,jp_tem,jn) = ptr_sj( pva(:,:,:), btmsk(:,:,jn) ) 
     406                hstr_adv(:,jp_tem,jn) = ptr_sj( pvflx(:,:,:), btmsk(:,:,jn) ) 
    413407             ENDDO 
    414408         ENDIF 
    415409         IF( ktra == jp_sal )  THEN 
    416410             DO jn = 1, nptr 
    417                 hstr_adv(:,jp_sal,jn) = ptr_sj( pva(:,:,:), btmsk(:,:,jn) ) 
     411                hstr_adv(:,jp_sal,jn) = ptr_sj( pvflx(:,:,:), btmsk(:,:,jn) ) 
    418412             ENDDO 
    419413         ENDIF 
     
    423417         IF( ktra == jp_tem )  THEN 
    424418             DO jn = 1, nptr 
    425                 hstr_ldf(:,jp_tem,jn) = ptr_sj( pva(:,:,:), btmsk(:,:,jn) ) 
     419                hstr_ldf(:,jp_tem,jn) = ptr_sj( pvflx(:,:,:), btmsk(:,:,jn) ) 
    426420             ENDDO 
    427421         ENDIF 
    428422         IF( ktra == jp_sal )  THEN 
    429423             DO jn = 1, nptr 
    430                 hstr_ldf(:,jp_sal,jn) = ptr_sj( pva(:,:,:), btmsk(:,:,jn) ) 
     424                hstr_ldf(:,jp_sal,jn) = ptr_sj( pvflx(:,:,:), btmsk(:,:,jn) ) 
    431425             ENDDO 
    432426         ENDIF 
     
    436430         IF( ktra == jp_tem )  THEN 
    437431             DO jn = 1, nptr 
    438                 hstr_eiv(:,jp_tem,jn) = ptr_sj( pva(:,:,:), btmsk(:,:,jn) ) 
     432                hstr_eiv(:,jp_tem,jn) = ptr_sj( pvflx(:,:,:), btmsk(:,:,jn) ) 
    439433             ENDDO 
    440434         ENDIF 
    441435         IF( ktra == jp_sal )  THEN 
    442436             DO jn = 1, nptr 
    443                 hstr_eiv(:,jp_sal,jn) = ptr_sj( pva(:,:,:), btmsk(:,:,jn) ) 
     437                hstr_eiv(:,jp_sal,jn) = ptr_sj( pvflx(:,:,:), btmsk(:,:,jn) ) 
    444438             ENDDO 
    445439         ENDIF 
     
    449443         IF( ktra == jp_tem )  THEN 
    450444             DO jn = 1, nptr 
    451                 hstr_vtr(:,jp_tem,jn) = ptr_sj( pva(:,:,:), btmsk(:,:,jn) ) 
     445                hstr_vtr(:,jp_tem,jn) = ptr_sj( pvflx(:,:,:), btmsk(:,:,jn) ) 
    452446             ENDDO 
    453447         ENDIF 
    454448         IF( ktra == jp_sal )  THEN 
    455449             DO jn = 1, nptr 
    456                 hstr_vtr(:,jp_sal,jn) = ptr_sj( pva(:,:,:), btmsk(:,:,jn) ) 
     450                hstr_vtr(:,jp_sal,jn) = ptr_sj( pvflx(:,:,:), btmsk(:,:,jn) ) 
    457451             ENDDO 
    458452         ENDIF 
     
    486480 
    487481 
    488    FUNCTION ptr_sj_3d( pva, pmsk )   RESULT ( p_fval ) 
     482   FUNCTION ptr_sj_3d( pvflx, pmsk )   RESULT ( p_fval ) 
    489483      !!---------------------------------------------------------------------- 
    490484      !!                    ***  ROUTINE ptr_sj_3d  *** 
     
    492486      !! ** Purpose :   i-k sum computation of a j-flux array 
    493487      !! 
    494       !! ** Method  : - i-k sum of pva using the interior 2D vmask (vmask_i). 
    495       !!              pva is supposed to be a masked flux (i.e. * vmask*e1v*e3v) 
    496       !! 
    497       !! ** Action  : - p_fval: i-k-mean poleward flux of pva 
    498       !!---------------------------------------------------------------------- 
    499       REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk)  ::   pva   ! mask flux array at V-point 
     488      !! ** Method  : - i-k sum of pvflx using the interior 2D vmask (vmask_i). 
     489      !!              pvflx is supposed to be a masked flux (i.e. * vmask*e1v*e3v) 
     490      !! 
     491      !! ** Action  : - p_fval: i-k-mean poleward flux of pvflx 
     492      !!---------------------------------------------------------------------- 
     493      REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk)  ::   pvflx  ! mask flux array at V-point 
    500494      REAL(wp), INTENT(in), DIMENSION(jpi,jpj)      ::   pmsk   ! Optional 2D basin mask 
    501495      ! 
     
    509503      ijpj = jpj 
    510504      p_fval(:) = 0._wp 
    511       DO jk = 1, jpkm1 
    512          DO jj = 2, jpjm1 
    513             DO ji = fs_2, fs_jpim1   ! Vector opt. 
    514                p_fval(jj) = p_fval(jj) + pva(ji,jj,jk) * pmsk(ji,jj) * tmask_i(ji,jj) 
    515             END DO 
    516          END DO 
    517       END DO 
     505      DO_3D_00_00( 1, jpkm1 ) 
     506         p_fval(jj) = p_fval(jj) + pvflx(ji,jj,jk) * pmsk(ji,jj) * tmask_i(ji,jj) 
     507      END_3D 
    518508#if defined key_mpp_mpi 
    519509      CALL mpp_sum( 'diaptr', p_fval, ijpj, ncomm_znl) 
     
    523513 
    524514 
    525    FUNCTION ptr_sj_2d( pva, pmsk )   RESULT ( p_fval ) 
     515   FUNCTION ptr_sj_2d( pvflx, pmsk )   RESULT ( p_fval ) 
    526516      !!---------------------------------------------------------------------- 
    527517      !!                    ***  ROUTINE ptr_sj_2d  *** 
    528518      !! 
    529       !! ** Purpose :   "zonal" and vertical sum computation of a i-flux array 
    530       !! 
    531       !! ** Method  : - i-k sum of pva using the interior 2D vmask (vmask_i). 
    532       !!      pva is supposed to be a masked flux (i.e. * vmask*e1v*e3v) 
    533       !! 
    534       !! ** Action  : - p_fval: i-k-mean poleward flux of pva 
    535       !!---------------------------------------------------------------------- 
    536       REAL(wp) , INTENT(in), DIMENSION(jpi,jpj) ::   pva   ! mask flux array at V-point 
     519      !! ** Purpose :   "zonal" and vertical sum computation of a j-flux array 
     520      !! 
     521      !! ** Method  : - i-k sum of pvflx using the interior 2D vmask (vmask_i). 
     522      !!      pvflx is supposed to be a masked flux (i.e. * vmask*e1v*e3v) 
     523      !! 
     524      !! ** Action  : - p_fval: i-k-mean poleward flux of pvflx 
     525      !!---------------------------------------------------------------------- 
     526      REAL(wp) , INTENT(in), DIMENSION(jpi,jpj) ::   pvflx  ! mask flux array at V-point 
    537527      REAL(wp) , INTENT(in), DIMENSION(jpi,jpj) ::   pmsk   ! Optional 2D basin mask 
    538528      ! 
     
    546536      ijpj = jpj 
    547537      p_fval(:) = 0._wp 
    548       DO jj = 2, jpjm1 
    549          DO ji = fs_2, fs_jpim1   ! Vector opt. 
    550             p_fval(jj) = p_fval(jj) + pva(ji,jj) * pmsk(ji,jj) * tmask_i(ji,jj) 
    551          END DO 
    552       END DO 
     538      DO_2D_00_00 
     539         p_fval(jj) = p_fval(jj) + pvflx(ji,jj) * pmsk(ji,jj) * tmask_i(ji,jj) 
     540      END_2D 
    553541#if defined key_mpp_mpi 
    554542      CALL mpp_sum( 'diaptr', p_fval, ijpj, ncomm_znl ) 
     
    577565      p_fval(:,:) = 0._wp 
    578566      DO jc = 1, jpnj ! looping over all processors in j axis 
    579          DO jj = 2, jpjm1 
    580             DO ji = fs_2, fs_jpim1   ! Vector opt. 
    581                p_fval(ji,jj) = p_fval(ji,jj-1) + pva(ji,jj) * tmask_i(ji,jj) 
    582             END DO 
    583          END DO 
     567         DO_2D_00_00 
     568            p_fval(ji,jj) = p_fval(ji,jj-1) + pva(ji,jj) * tmask_i(ji,jj) 
     569         END_2D 
    584570         CALL lbc_lnk( 'diaptr', p_fval, 'U', -1. ) 
    585571      END DO 
     
    595581      !! ** Purpose :   i-sum computation of an array 
    596582      !! 
    597       !! ** Method  : - i-sum of pva using the interior 2D vmask (vmask_i). 
    598       !! 
    599       !! ** Action  : - p_fval: i-mean poleward flux of pva 
     583      !! ** Method  : - i-sum of field using the interior 2D vmask (pmsk). 
     584      !! 
     585      !! ** Action  : - p_fval: i-sum of masked field 
    600586      !!---------------------------------------------------------------------- 
    601587      !! 
     
    618604      p_fval(:,:) = 0._wp 
    619605      ! 
    620       DO jk = 1, jpkm1 
    621          DO jj = 2, jpjm1 
    622             DO ji = fs_2, fs_jpim1   ! Vector opt. 
    623                p_fval(jj,jk) = p_fval(jj,jk) + pta(ji,jj,jk) * pmsk(ji,jj) * tmask_i(ji,jj) 
    624             END DO 
    625          END DO 
    626       END DO 
     606      DO_3D_00_00( 1, jpkm1 ) 
     607         p_fval(jj,jk) = p_fval(jj,jk) + pta(ji,jj,jk) * pmsk(ji,jj) * tmask_i(ji,jj) 
     608      END_3D 
    627609      ! 
    628610#if defined key_mpp_mpi 
  • NEMO/trunk/src/OCE/DIA/diawri.F90

    r12206 r12377  
    2626   !!---------------------------------------------------------------------- 
    2727   USE oce            ! ocean dynamics and tracers  
     28   USE isf_oce 
     29   USE isfcpl 
     30   USE abl            ! abl variables in case ln_abl = .true. 
    2831   USE dom_oce        ! ocean space and time domain 
    2932   USE phycst         ! physical constants 
     
    5659   USE lib_mpp         ! MPP library 
    5760   USE timing          ! preformance summary 
    58    USE diurnal_bulk    ! diurnal warm layer 
    59    USE cool_skin       ! Cool skin 
     61   USE diu_bulk        ! diurnal warm layer 
     62   USE diu_coolskin    ! Cool skin 
    6063 
    6164   IMPLICIT NONE 
     
    6568   PUBLIC   dia_wri_state 
    6669   PUBLIC   dia_wri_alloc           ! Called by nemogcm module 
    67  
     70#if ! defined key_iomput    
     71   PUBLIC   dia_wri_alloc_abl       ! Called by sbcabl  module (if ln_abl = .true.) 
     72#endif 
    6873   INTEGER ::   nid_T, nz_T, nh_T, ndim_T, ndim_hT   ! grid_T file 
    6974   INTEGER ::          nb_T              , ndim_bT   ! grid_T file 
     
    7176   INTEGER ::   nid_V, nz_V, nh_V, ndim_V, ndim_hV   ! grid_V file 
    7277   INTEGER ::   nid_W, nz_W, nh_W                    ! grid_W file 
     78   INTEGER ::   nid_A, nz_A, nh_A, ndim_A, ndim_hA   ! grid_ABL file    
    7379   INTEGER ::   ndex(1)                              ! ??? 
    7480   INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_hT, ndex_hU, ndex_hV 
     81   INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_hA, ndex_A ! ABL 
    7582   INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_T, ndex_U, ndex_V 
    7683   INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_bT 
    7784 
    7885   !! * Substitutions 
    79 #  include "vectopt_loop_substitute.h90" 
     86#  include "do_loop_substitute.h90" 
    8087   !!---------------------------------------------------------------------- 
    8188   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
     
    96103 
    97104    
    98    SUBROUTINE dia_wri( kt ) 
     105   SUBROUTINE dia_wri( kt, Kmm ) 
    99106      !!--------------------------------------------------------------------- 
    100107      !!                  ***  ROUTINE dia_wri  *** 
     
    106113      !!---------------------------------------------------------------------- 
    107114      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index 
     115      INTEGER, INTENT( in ) ::   Kmm     ! ocean time level index 
    108116      !! 
    109117      INTEGER ::   ji, jj, jk       ! dummy loop indices 
     
    119127      ! Output the initial state and forcings 
    120128      IF( ninist == 1 ) THEN                        
    121          CALL dia_wri_state( 'output.init' ) 
     129         CALL dia_wri_state( Kmm, 'output.init' ) 
    122130         ninist = 0 
    123131      ENDIF 
     
    128136      CALL iom_put("e3v_0", e3v_0(:,:,:) ) 
    129137      ! 
    130       CALL iom_put( "e3t" , e3t_n(:,:,:) ) 
    131       CALL iom_put( "e3u" , e3u_n(:,:,:) ) 
    132       CALL iom_put( "e3v" , e3v_n(:,:,:) ) 
    133       CALL iom_put( "e3w" , e3w_n(:,:,:) ) 
     138      CALL iom_put( "e3t" , e3t(:,:,:,Kmm) ) 
     139      CALL iom_put( "e3u" , e3u(:,:,:,Kmm) ) 
     140      CALL iom_put( "e3v" , e3v(:,:,:,Kmm) ) 
     141      CALL iom_put( "e3w" , e3w(:,:,:,Kmm) ) 
    134142      IF( iom_use("e3tdef") )   & 
    135          CALL iom_put( "e3tdef"  , ( ( e3t_n(:,:,:) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100 * tmask(:,:,:) ) ** 2 ) 
     143         CALL iom_put( "e3tdef"  , ( ( e3t(:,:,:,Kmm) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100 * tmask(:,:,:) ) ** 2 ) 
    136144 
    137145      IF( ll_wd ) THEN 
    138          CALL iom_put( "ssh" , (sshn+ssh_ref)*tmask(:,:,1) )   ! sea surface height (brought back to the reference used for wetting and drying) 
     146         CALL iom_put( "ssh" , (ssh(:,:,Kmm)+ssh_ref)*tmask(:,:,1) )   ! sea surface height (brought back to the reference used for wetting and drying) 
    139147      ELSE 
    140          CALL iom_put( "ssh" , sshn )              ! sea surface height 
     148         CALL iom_put( "ssh" , ssh(:,:,Kmm) )              ! sea surface height 
    141149      ENDIF 
    142150 
    143151      IF( iom_use("wetdep") )   &                  ! wet depth 
    144          CALL iom_put( "wetdep" , ht_0(:,:) + sshn(:,:) ) 
     152         CALL iom_put( "wetdep" , ht_0(:,:) + ssh(:,:,Kmm) ) 
    145153       
    146       CALL iom_put( "toce", tsn(:,:,:,jp_tem) )    ! 3D temperature 
    147       CALL iom_put(  "sst", tsn(:,:,1,jp_tem) )    ! surface temperature 
     154      CALL iom_put( "toce", ts(:,:,:,jp_tem,Kmm) )    ! 3D temperature 
     155      CALL iom_put(  "sst", ts(:,:,1,jp_tem,Kmm) )    ! surface temperature 
    148156      IF ( iom_use("sbt") ) THEN 
    149          DO jj = 1, jpj 
    150             DO ji = 1, jpi 
    151                ikbot = mbkt(ji,jj) 
    152                z2d(ji,jj) = tsn(ji,jj,ikbot,jp_tem) 
    153             END DO 
    154          END DO 
     157         DO_2D_11_11 
     158            ikbot = mbkt(ji,jj) 
     159            z2d(ji,jj) = ts(ji,jj,ikbot,jp_tem,Kmm) 
     160         END_2D 
    155161         CALL iom_put( "sbt", z2d )                ! bottom temperature 
    156162      ENDIF 
    157163       
    158       CALL iom_put( "soce", tsn(:,:,:,jp_sal) )    ! 3D salinity 
    159       CALL iom_put(  "sss", tsn(:,:,1,jp_sal) )    ! surface salinity 
     164      CALL iom_put( "soce", ts(:,:,:,jp_sal,Kmm) )    ! 3D salinity 
     165      CALL iom_put(  "sss", ts(:,:,1,jp_sal,Kmm) )    ! surface salinity 
    160166      IF ( iom_use("sbs") ) THEN 
    161          DO jj = 1, jpj 
    162             DO ji = 1, jpi 
    163                ikbot = mbkt(ji,jj) 
    164                z2d(ji,jj) = tsn(ji,jj,ikbot,jp_sal) 
    165             END DO 
    166          END DO 
     167         DO_2D_11_11 
     168            ikbot = mbkt(ji,jj) 
     169            z2d(ji,jj) = ts(ji,jj,ikbot,jp_sal,Kmm) 
     170         END_2D 
    167171         CALL iom_put( "sbs", z2d )                ! bottom salinity 
    168172      ENDIF 
     
    171175         zztmp = rau0 * 0.25 
    172176         z2d(:,:) = 0._wp 
    173          DO jj = 2, jpjm1 
    174             DO ji = fs_2, fs_jpim1   ! vector opt. 
    175                zztmp2 = (  ( rCdU_bot(ji+1,jj)+rCdU_bot(ji  ,jj) ) * un(ji  ,jj,mbku(ji  ,jj))  )**2   & 
    176                   &   + (  ( rCdU_bot(ji  ,jj)+rCdU_bot(ji-1,jj) ) * un(ji-1,jj,mbku(ji-1,jj))  )**2   & 
    177                   &   + (  ( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj  ) ) * vn(ji,jj  ,mbkv(ji,jj  ))  )**2   & 
    178                   &   + (  ( rCdU_bot(ji,jj  )+rCdU_bot(ji,jj-1) ) * vn(ji,jj-1,mbkv(ji,jj-1))  )**2 
    179                z2d(ji,jj) = zztmp * SQRT( zztmp2 ) * tmask(ji,jj,1)  
    180                ! 
    181             END DO 
    182          END DO 
     177         DO_2D_00_00 
     178            zztmp2 = (  ( rCdU_bot(ji+1,jj)+rCdU_bot(ji  ,jj) ) * uu(ji  ,jj,mbku(ji  ,jj),Kmm)  )**2   & 
     179               &   + (  ( rCdU_bot(ji  ,jj)+rCdU_bot(ji-1,jj) ) * uu(ji-1,jj,mbku(ji-1,jj),Kmm)  )**2   & 
     180               &   + (  ( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj  ) ) * vv(ji,jj  ,mbkv(ji,jj  ),Kmm)  )**2   & 
     181               &   + (  ( rCdU_bot(ji,jj  )+rCdU_bot(ji,jj-1) ) * vv(ji,jj-1,mbkv(ji,jj-1),Kmm)  )**2 
     182            z2d(ji,jj) = zztmp * SQRT( zztmp2 ) * tmask(ji,jj,1)  
     183            ! 
     184         END_2D 
    183185         CALL lbc_lnk( 'diawri', z2d, 'T', 1. ) 
    184186         CALL iom_put( "taubot", z2d )            
    185187      ENDIF 
    186188          
    187       CALL iom_put( "uoce", un(:,:,:) )            ! 3D i-current 
    188       CALL iom_put(  "ssu", un(:,:,1) )            ! surface i-current 
     189      CALL iom_put( "uoce", uu(:,:,:,Kmm) )            ! 3D i-current 
     190      CALL iom_put(  "ssu", uu(:,:,1,Kmm) )            ! surface i-current 
    189191      IF ( iom_use("sbu") ) THEN 
    190          DO jj = 1, jpj 
    191             DO ji = 1, jpi 
    192                ikbot = mbku(ji,jj) 
    193                z2d(ji,jj) = un(ji,jj,ikbot) 
    194             END DO 
    195          END DO 
     192         DO_2D_11_11 
     193            ikbot = mbku(ji,jj) 
     194            z2d(ji,jj) = uu(ji,jj,ikbot,Kmm) 
     195