Changeset 13427


Ignore:
Timestamp:
2020-08-21T18:26:57+02:00 (5 months ago)
Author:
techene
Message:

debug in order to pass non linear SETTE test when agrif in not activated see #2385

Location:
NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/src
Files:
16 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/src/OCE/DOM/dom_oce.F90

    r13286 r13427  
    287287         &       e3w_0(jpi,jpj,jpk) , e3uw_0(jpi,jpj,jpk) , e3vw_0(jpi,jpj,jpk)                      ,  STAT=ierr(ii) ) 
    288288         ! 
    289 #if ! defined key_qco 
     289#if defined key_qco 
     290      ii = ii+1 
     291      ALLOCATE( r3t  (jpi,jpj,jpt)   , r3u  (jpi,jpj,jpt)    , r3v  (jpi,jpj,jpt)    , r3f  (jpi,jpj) ,      & 
     292         &      r3t_f(jpi,jpj)       , r3u_f(jpi,jpj)        , r3v_f(jpi,jpj)                         ,  STAT=ierr(ii) )              
     293#else 
    290294      ii = ii+1 
    291295      ALLOCATE( e3t(jpi,jpj,jpk,jpt) , e3u (jpi,jpj,jpk,jpt) , e3v (jpi,jpj,jpk,jpt) , e3f(jpi,jpj,jpk) ,      & 
     
    294298         ! 
    295299      ii = ii+1 
    296       ALLOCATE( r3t  (jpi,jpj,jpt)   , r3u  (jpi,jpj,jpt)    , r3v  (jpi,jpj,jpt)    , r3f  (jpi,jpj) ,  & 
    297          &      r3t_f(jpi,jpj)       , r3u_f(jpi,jpj)        , r3v_f(jpi,jpj)                         ,  STAT=ierr(ii) )        
    298          ! 
    299       ii = ii+1 
    300300      ALLOCATE( ht_0(jpi,jpj) ,    hu_0(jpi,jpj)    ,    hv_0(jpi,jpj)     , hf_0(jpi,jpj) ,       & 
    301301         &   r1_ht_0(jpi,jpj) , r1_hu_0(jpi,jpj) ,    r1_hv_0(jpi,jpj),   r1_hf_0(jpi,jpj) ,   STAT=ierr(ii)  ) 
     
    304304      ii = ii+1 
    305305      ALLOCATE( ht  (jpi,jpj) ,    hu  (jpi,jpj,jpt),    hv  (jpi,jpj,jpt)                 ,       & 
    306          &                      r1_hu  (jpi,jpj,jpt), r1_hv  (jpi,jpj,jpt)                 ,   STAT=ierr(ii)  ) 
    307 #else 
    308       ii = ii+1 
    309       ALLOCATE(                    hu  (jpi,jpj,jpt),    hv  (jpi,jpj,jpt)                 ,       & 
    310306         &                      r1_hu  (jpi,jpj,jpt), r1_hv  (jpi,jpj,jpt)                 ,   STAT=ierr(ii)  ) 
    311307#endif 
  • NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/src/OCE/DOM/domain.F90

    r13286 r13427  
    177177      ! 
    178178      IF( ln_linssh ) THEN       != Fix in time : set to the reference one for all 
    179       ! 
     179         ! 
    180180         DO jt = 1, jpt                         ! depth of t- and w-grid-points 
    181181            gdept(:,:,:,jt) = gdept_0(:,:,:) 
     
    185185         ! 
    186186         DO jt = 1, jpt                         ! vertical scale factors 
    187             e3t(:,:,:,jt) =  e3t_0(:,:,:) 
    188             e3u(:,:,:,jt) =  e3u_0(:,:,:) 
    189             e3v(:,:,:,jt) =  e3v_0(:,:,:) 
    190             e3w(:,:,:,jt) =  e3w_0(:,:,:) 
     187            e3t (:,:,:,jt) =  e3t_0(:,:,:) 
     188            e3u (:,:,:,jt) =  e3u_0(:,:,:) 
     189            e3v (:,:,:,jt) =  e3v_0(:,:,:) 
     190            e3w (:,:,:,jt) =  e3w_0(:,:,:) 
    191191            e3uw(:,:,:,jt) = e3uw_0(:,:,:) 
    192192            e3vw(:,:,:,jt) = e3vw_0(:,:,:) 
    193193         END DO 
    194             e3f(:,:,:)    =  e3f_0(:,:,:) 
     194            e3f (:,:,:)    =  e3f_0(:,:,:) 
    195195         ! 
    196196         DO jt = 1, jpt                         ! water column thickness and its inverse 
    197             hu(:,:,jt)    =    hu_0(:,:) 
    198             hv(:,:,jt)    =    hv_0(:,:) 
     197               hu(:,:,jt) =    hu_0(:,:) 
     198               hv(:,:,jt) =    hv_0(:,:) 
    199199            r1_hu(:,:,jt) = r1_hu_0(:,:) 
    200200            r1_hv(:,:,jt) = r1_hv_0(:,:) 
    201201         END DO 
    202             ht(:,:) =    ht_0(:,:) 
     202               ht   (:,:) =    ht_0(:,:) 
    203203         ! 
    204204      ELSE                       != time varying : initialize before/now/after variables 
  • NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/src/OCE/DOM/domqco.F90

    r13295 r13427  
    88   !!            3.3  !  2011-10  (M. Leclair) totally rewrote domvvl: vvl option includes z_star and z_tilde coordinates 
    99   !!            3.6  !  2014-11  (P. Mathiot) add ice shelf capability 
    10    !!            4.1  !  2019-08  (A. Coward, D. Storkey) rename dom_vvl_sf_swp -> dom_vvl_sf_update for new timestepping 
    11    !!            4.x  !  2020-02  (G. Madec, S. Techene) pure z* (quasi-eulerian) coordinate 
     10   !!            4.1  !  2019-08  (A. Coward, D. Storkey) add time level indices for prognostic variables 
     11   !!            4.x  !  2020-02  (S. Techene, G. Madec) quasi-eulerian coordinate (z* or s*) from domvvl 
    1212   !!---------------------------------------------------------------------- 
    1313 
    1414   !!---------------------------------------------------------------------- 
    15    !!   dom_qe_init   : define initial vertical scale factors, depths and column thickness 
    16    !!   dom_qe_r3c    : Compute ssh/h_0 ratioat t-, u-, v-, and optionally f-points 
    17    !!       qe_rst_read : read/write restart file 
    18    !!   dom_qe_ctl    : Check the vvl options 
     15   !!   dom_qco_init   : define initial vertical scale factors, depths and column thickness 
     16   !!   dom_qco_zgr    : Set ssh/h_0 ratio at t 
     17   !!   dom_qco_r3c    : Compute ssh/h_0 ratio at t-, u-, v-, and optionally f-points 
     18   !!       qco_rst_read : read/write restart file 
     19   !!       qco_ctl    : Check the vvl options 
    1920   !!---------------------------------------------------------------------- 
    2021   USE oce            ! ocean dynamics and tracers 
     
    7980      !! 
    8081      !!---------------------------------------------------------------------- 
    81       INTEGER, INTENT(in) :: Kbb, Kmm, Kaa 
     82      INTEGER, INTENT(in) ::   Kbb, Kmm, Kaa   ! time level indices 
     83      !!---------------------------------------------------------------------- 
    8284      ! 
    8385      IF(lwp) WRITE(numout,*) 
     
    8587      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~' 
    8688      ! 
    87       CALL dom_qco_ctl     ! choose vertical coordinate (z_star, z_tilde or layer) 
    88       ! 
    89       !                    ! Read or initialize e3t_(b/n), tilde_e3t_(b/n) and hdiv_lf 
    90       CALL qe_rst_read( nit000, Kbb, Kmm ) 
    91       ! 
    92       CALL dom_qco_zgr(Kbb, Kmm, Kaa) ! interpolation scale factor, depth and water column 
    93       ! 
    94       ! IF(lwxios) THEN   ! define variables in restart file when writing with XIOS 
    95       !    CALL iom_set_rstw_var_active('e3t_b') 
    96       !    CALL iom_set_rstw_var_active('e3t_n') 
    97       ! ENDIF 
     89      CALL qco_ctl                            ! choose vertical coordinate (z_star, z_tilde or layer) 
     90      !              ! CAUTION COM A METTRE !!! 
     91!!st      CALL qco_rst_read2( nit000, Kbb, Kmm )  ! Read or initialize ssh_(Kbb/Kmm) and r3 
     92!!st CAUTION if read2 removed change restart.F90 !  
     93      ! 
     94      CALL qco_rst_read( nit000, Kbb, Kmm )   ! Read or initialize ssh_(Kbb/Kmm) 
     95      ! 
     96      CALL dom_qco_zgr( Kbb, Kmm, Kaa )       ! interpolation scale factor, depth and water column 
     97      ! 
     98      IF(lwxios) THEN   ! define variables in restart file when writing with XIOS 
     99         CALL iom_set_rstw_var_active('sshb') 
     100         CALL iom_set_rstw_var_active('sshn') 
     101      ENDIF 
    98102      ! 
    99103   END SUBROUTINE dom_qco_init 
    100104 
    101105 
    102    SUBROUTINE dom_qco_zgr(Kbb, Kmm, Kaa) 
     106   SUBROUTINE dom_qco_zgr( Kbb, Kmm, Kaa ) 
    103107      !!---------------------------------------------------------------------- 
    104108      !!                ***  ROUTINE dom_qco_init  *** 
    105109      !! 
    106       !! ** Purpose :  Initialization of all ssh. to h._0 ratio 
    107       !! 
    108       !! ** Method  :  - interpolate scale factors 
     110      !! ** Purpose :  Initialization of all ssh./h._0 ratio 
     111      !! 
     112      !! ** Method  :  - call domqco using Kbb and Kmm 
    109113      !! 
    110114      !! ** Action  : - r3(t/u/v)_b 
    111115      !!              - r3(t/u/v/f)_n 
    112       !! 
    113       !! Reference  : Leclair, M., and G. Madec, 2011, Ocean Modelling. 
    114       !!---------------------------------------------------------------------- 
    115       INTEGER, INTENT(in) :: Kbb, Kmm, Kaa 
     116      !!---------------------------------------------------------------------- 
     117      INTEGER, INTENT(in) ::   Kbb, Kmm, Kaa   ! time level indices 
    116118      !!---------------------------------------------------------------------- 
    117119      ! 
     
    121123      CALL dom_qco_r3c( ssh(:,:,Kmm), r3t(:,:,Kmm), r3u(:,:,Kmm), r3v(:,:,Kmm), r3f(:,:) ) 
    122124      ! 
     125#if defined key_agrif 
     126      ! We need to define r3[tuv](Kaa) for AGRIF initialisation (should not be a problem for the restartability...) 
     127      r3t(:,:,Kaa) = r3t(:,:,Kmm) 
     128      r3u(:,:,Kaa) = r3u(:,:,Kmm) 
     129      r3v(:,:,Kaa) = r3v(:,:,Kmm) 
     130#endif 
     131      ! 
    123132   END SUBROUTINE dom_qco_zgr 
    124133 
     
    148157      !                                      !==  ratio at u-,v-point  ==! 
    149158      ! 
    150       IF( ln_dynadv_vec ) THEN                     !- Vector Form   (thickness weighted averaging) 
     159!!st      IF( ln_dynadv_vec ) THEN                     !- Vector Form   (thickness weighted averaging) 
     160#if ! defined key_qcoTest_FluxForm 
     161      !                                ! no 'key_qcoTest_FluxForm' : surface weighted ssh average 
    151162         DO_2D( 0, 0, 0, 0 ) 
    152163            pr3u(ji,jj) = 0.5_wp * (  e1e2t(ji  ,jj) * pssh(ji  ,jj)  & 
     
    155166               &                    + e1e2t(ji,jj+1) * pssh(ji,jj+1)  ) * r1_hv_0(ji,jj) * r1_e1e2v(ji,jj) 
    156167         END_2D 
    157       ELSE                                         !- Flux Form   (simple averaging) 
     168!!st      ELSE                                         !- Flux Form   (simple averaging) 
     169#else 
    158170         DO_2D( 0, 0, 0, 0 ) 
    159             pr3u(ji,jj) = 0.5_wp * (  pssh(ji  ,jj) + pssh(ji+1,jj)  ) * r1_hu_0(ji,jj) 
    160             pr3v(ji,jj) = 0.5_wp * (  pssh(ji,jj  ) + pssh(ji,jj+1)  ) * r1_hv_0(ji,jj) 
     171            pr3u(ji,jj) = 0.5_wp * (  pssh(ji,jj) + pssh(ji+1,jj  )  ) * r1_hu_0(ji,jj) 
     172            pr3v(ji,jj) = 0.5_wp * (  pssh(ji,jj) + pssh(ji  ,jj+1)  ) * r1_hv_0(ji,jj) 
    161173         END_2D 
    162       ENDIF 
     174!!st      ENDIF 
     175#endif          
    163176      ! 
    164177      IF( .NOT.PRESENT( pr3f ) ) THEN              !- lbc on ratio at u-, v-points only 
     
    168181      ELSE                                   !==  ratio at f-point  ==! 
    169182         ! 
    170          IF( ln_dynadv_vec )   THEN                !- Vector Form   (thickness weighted averaging) 
     183!!st         IF( ln_dynadv_vec )   THEN                !- Vector Form   (thickness weighted averaging) 
     184#if ! defined key_qcoTest_FluxForm 
     185         !                                ! no 'key_qcoTest_FluxForm' : surface weighted ssh average 
     186 
    171187            DO_2D( 1, 0, 1, 0 )                               ! start from 1 since lbc_lnk('F') doesn't update the 1st row/line 
    172188               pr3f(ji,jj) = 0.25_wp * (  e1e2t(ji  ,jj  ) * pssh(ji  ,jj  )  & 
     
    175191                  &                     + e1e2t(ji+1,jj+1) * pssh(ji+1,jj+1)  ) * r1_hf_0(ji,jj) * r1_e1e2f(ji,jj) 
    176192            END_2D 
    177          ELSE                                      !- Flux Form   (simple averaging) 
     193!!st         ELSE                                      !- Flux Form   (simple averaging) 
     194#else 
    178195            DO_2D( 1, 0, 1, 0 )                               ! start from 1 since lbc_lnk('F') doesn't update the 1st row/line 
    179                pr3f(ji,jj) = 0.25_wp * (  pssh(ji  ,jj  ) + pssh(ji+1,jj  )  & 
    180                   &                     + pssh(ji  ,jj+1) + pssh(ji+1,jj+1)  ) * r1_hf_0(ji,jj) 
     196               pr3f(ji,jj) = 0.25_wp * (  pssh(ji,jj  ) + pssh(ji+1,jj  )  & 
     197                  &                     + pssh(ji,jj+1) + pssh(ji+1,jj+1)  ) * r1_hf_0(ji,jj) 
    181198            END_2D 
    182          ENDIF 
     199!!st         ENDIF 
     200#endif 
    183201         !                                                 ! lbc on ratio at u-,v-,f-points 
    184202         CALL lbc_lnk_multi( 'dom_qco_r3c', pr3u, 'U', 1._wp, pr3v, 'V', 1._wp, pr3f, 'F', 1._wp ) 
     
    189207 
    190208 
    191    SUBROUTINE qe_rst_read( kt, Kbb, Kmm ) 
     209   SUBROUTINE qco_rst_read( kt, Kbb, Kmm ) 
    192210      !!--------------------------------------------------------------------- 
    193       !!                   ***  ROUTINE qe_rst_read  *** 
     211      !!                   ***  ROUTINE qco_rst_read  *** 
    194212      !! 
    195213      !! ** Purpose :   Read ssh in restart file 
     
    199217      !!                it is set to the _0 values. 
    200218      !!---------------------------------------------------------------------- 
    201       INTEGER         , INTENT(in) ::   kt        ! ocean time-step 
    202       INTEGER         , INTENT(in) ::   Kbb, Kmm  ! ocean time level indices 
     219      INTEGER, INTENT(in) ::   kt         ! ocean time-step 
     220      INTEGER, INTENT(in) ::   Kbb, Kmm   ! ocean time level indices 
    203221      ! 
    204222      INTEGER ::   ji, jj, jk 
     
    206224      !!---------------------------------------------------------------------- 
    207225      ! 
    208          IF( ln_rstart ) THEN                   !* Read the restart file 
    209             CALL rst_read_open                  !  open the restart file if necessary 
    210             ! 
    211             id1 = iom_varid( numror, 'sshb', ldstop = .FALSE. ) 
    212             id2 = iom_varid( numror, 'sshn', ldstop = .FALSE. ) 
    213             ! 
    214             !                             ! --------- ! 
    215             !                             ! all cases ! 
    216             !                             ! --------- ! 
    217             ! 
    218             IF( MIN( id1, id2 ) > 0 ) THEN       ! all required arrays exist 
    219                CALL iom_get( numror, jpdom_auto, 'sshb'   , ssh(:,:,Kbb), ldxios = lrxios    ) 
    220                CALL iom_get( numror, jpdom_auto, 'sshn'   , ssh(:,:,Kmm), ldxios = lrxios    ) 
    221                ! needed to restart if land processor not computed 
    222                IF(lwp) write(numout,*) 'qe_rst_read : ssh(:,:,Kbb) and ssh(:,:,Kmm) found in restart files' 
    223                WHERE ( ssmask(:,:) == 0.0_wp )   !!gm/st ==> sm should not be necessary on ssh when it was required on e3 
    224                   ssh(:,:,Kmm) = 0._wp 
    225                   ssh(:,:,Kbb) = 0._wp 
    226                END WHERE 
    227                IF( l_1st_euler ) THEN 
    228                   ssh(:,:,Kbb) = ssh(:,:,Kmm) 
    229                ENDIF 
    230             ELSE IF( id1 > 0 ) THEN 
    231                IF(lwp) write(numout,*) 'qe_rst_read WARNING : ssh(:,:,Kmm) not found in restart files' 
    232                IF(lwp) write(numout,*) 'sshn set equal to sshb.' 
    233                IF(lwp) write(numout,*) 'neuler is forced to 0' 
    234                CALL iom_get( numror, jpdom_auto, 'sshb', ssh(:,:,Kbb), ldxios = lrxios ) 
    235                ssh(:,:,Kmm) = ssh(:,:,Kbb) 
    236                l_1st_euler = .TRUE. 
    237             ELSE IF( id2 > 0 ) THEN 
    238                IF(lwp) write(numout,*) 'qe_rst_read WARNING : ssh(:,:,Kbb) not found in restart files' 
    239                IF(lwp) write(numout,*) 'sshb set equal to sshn.' 
    240                IF(lwp) write(numout,*) 'neuler is forced to 0' 
    241                CALL iom_get( numror, jpdom_auto, 'sshn', ssh(:,:,Kmm), ldxios = lrxios ) 
     226      IF( ln_rstart ) THEN                   !* Read the restart file 
     227         CALL rst_read_open                  !  open the restart file if necessary 
     228         ! 
     229         id1 = iom_varid( numror, 'sshb', ldstop = .FALSE. ) 
     230         id2 = iom_varid( numror, 'sshn', ldstop = .FALSE. ) 
     231         ! 
     232         !                             ! --------- ! 
     233         !                             ! all cases ! 
     234         !                             ! --------- ! 
     235         ! 
     236         IF( MIN( id1, id2 ) > 0 ) THEN       ! all required arrays exist 
     237            CALL iom_get( numror, jpdom_auto, 'sshb'   , ssh(:,:,Kbb), ldxios = lrxios    ) 
     238            CALL iom_get( numror, jpdom_auto, 'sshn'   , ssh(:,:,Kmm), ldxios = lrxios    ) 
     239            ! needed to restart if land processor not computed 
     240            IF(lwp) write(numout,*) 'qco_rst_read : ssh(:,:,Kbb) and ssh(:,:,Kmm) found in restart files' 
     241            !!WHERE ( ssmask(:,:) == 0.0_wp )   !!gm/st ==> sm should not be necessary on ssh while it was required on e3 
     242            !!   ssh(:,:,Kmm) = 0._wp 
     243            !!   ssh(:,:,Kbb) = 0._wp 
     244            !!END WHERE 
     245            IF( l_1st_euler ) THEN 
    242246               ssh(:,:,Kbb) = ssh(:,:,Kmm) 
    243                l_1st_euler = .TRUE. 
    244             ELSE 
    245                IF(lwp) write(numout,*) 'qe_rst_read WARNING : ssh(:,:,Kmm) not found in restart file' 
    246                IF(lwp) write(numout,*) 'ssh_b and ssh_n set to zero' 
    247                IF(lwp) write(numout,*) 'neuler is forced to 0' 
    248                ssh(:,:,:) = 0._wp 
    249                l_1st_euler = .TRUE. 
    250247            ENDIF 
    251             ! 
    252          ELSE                                   !* Initialize at "rest" 
    253             ! 
    254             IF( ll_wd ) THEN   ! MJB ll_wd edits start here - these are essential 
     248         ELSE IF( id1 > 0 ) THEN 
     249            IF(lwp) write(numout,*) 'qco_rst_read WARNING : ssh(:,:,Kmm) not found in restart files' 
     250            IF(lwp) write(numout,*) 'sshn set equal to sshb.' 
     251            IF(lwp) write(numout,*) 'neuler is forced to 0' 
     252            CALL iom_get( numror, jpdom_auto, 'sshb', ssh(:,:,Kbb), ldxios = lrxios ) 
     253            ssh(:,:,Kmm) = ssh(:,:,Kbb) 
     254            l_1st_euler = .TRUE. 
     255         ELSE IF( id2 > 0 ) THEN 
     256            IF(lwp) write(numout,*) 'qco_rst_read WARNING : ssh(:,:,Kbb) not found in restart files' 
     257            IF(lwp) write(numout,*) 'sshb set equal to sshn.' 
     258            IF(lwp) write(numout,*) 'neuler is forced to 0' 
     259            CALL iom_get( numror, jpdom_auto, 'sshn', ssh(:,:,Kmm), ldxios = lrxios ) 
     260            ssh(:,:,Kbb) = ssh(:,:,Kmm) 
     261            l_1st_euler = .TRUE. 
     262         ELSE 
     263            IF(lwp) write(numout,*) 'qco_rst_read WARNING : ssh(:,:,Kmm) not found in restart file' 
     264            IF(lwp) write(numout,*) 'ssh_b and ssh_n set to zero' 
     265            IF(lwp) write(numout,*) 'neuler is forced to 0' 
     266            ssh(:,:,:) = 0._wp 
     267            l_1st_euler = .TRUE. 
     268         ENDIF 
     269         ! 
     270      ELSE                                   !* Initialize at "rest" 
     271         ! 
     272         IF( ll_wd ) THEN   ! MJB ll_wd edits start here - these are essential 
     273            ! 
     274            IF( cn_cfg == 'wad' ) THEN            ! Wetting and drying test case 
     275               CALL usr_def_istate( gdept(:,:,:,Kbb), tmask, ts(:,:,:,:,Kbb), uu(:,:,:,Kbb), vv(:,:,:,Kbb), ssh(:,:,Kbb)  ) 
     276               ts (:,:,:,:,Kmm) = ts (:,:,:,:,Kbb)       ! set now values from to before ones 
     277               ssh(:,:    ,Kmm) = ssh(:,:    ,Kbb) 
     278               uu (:,:,:  ,Kmm) = uu (:,:,:  ,Kbb) 
     279               vv (:,:,:  ,Kmm) = vv (:,:,:  ,Kbb) 
     280            ELSE                                  ! if not test case 
     281               ssh(:,:,Kmm) = -ssh_ref 
     282               ssh(:,:,Kbb) = -ssh_ref 
    255283               ! 
    256                IF( cn_cfg == 'wad' ) THEN            ! Wetting and drying test case 
    257                   CALL usr_def_istate( gdept(:,:,:,Kbb), tmask, ts(:,:,:,:,Kbb), uu(:,:,:,Kbb), vv(:,:,:,Kbb), ssh(:,:,Kbb)  ) 
    258                   ts (:,:,:,:,Kmm) = ts (:,:,:,:,Kbb)       ! set now values from to before ones 
    259                   ssh(:,:    ,Kmm) = ssh(:,:    ,Kbb) 
    260                   uu (:,:,:  ,Kmm) = uu (:,:,:  ,Kbb) 
    261                   vv (:,:,:  ,Kmm) = vv (:,:,:  ,Kbb) 
    262                ELSE                                  ! if not test case 
    263                   ssh(:,:,Kmm) = -ssh_ref 
    264                   ssh(:,:,Kbb) = -ssh_ref 
    265                   ! 
    266                   DO_2D( 1, 1, 1, 1 ) 
    267                      IF( ht_0(ji,jj)-ssh_ref <  rn_wdmin1 ) THEN ! if total depth is less than min depth 
    268                         ssh(ji,jj,Kbb) = rn_wdmin1 - (ht_0(ji,jj) ) 
    269                         ssh(ji,jj,Kmm) = rn_wdmin1 - (ht_0(ji,jj) ) 
    270                      ENDIF 
    271                   END_2D 
    272                ENDIF 
    273  
    274                DO ji = 1, jpi 
    275                   DO jj = 1, jpj 
    276                      IF ( ht_0(ji,jj) .LE. 0.0 .AND. NINT( ssmask(ji,jj) ) .EQ. 1) THEN 
    277                        CALL ctl_stop( 'qe_rst_read: ht_0 must be positive at potentially wet points' ) 
    278                      ENDIF 
    279                   END DO 
     284               DO_2D( 1, 1, 1, 1 ) 
     285                  IF( ht_0(ji,jj)-ssh_ref <  rn_wdmin1 ) THEN ! if total depth is less than min depth 
     286                     ssh(ji,jj,Kbb) = rn_wdmin1 - (ht_0(ji,jj) ) 
     287                     ssh(ji,jj,Kmm) = rn_wdmin1 - (ht_0(ji,jj) ) 
     288                  ENDIF 
     289               END_2D 
     290            ENDIF 
     291            ! 
     292            DO ji = 1, jpi 
     293               DO jj = 1, jpj 
     294                  IF ( ht_0(ji,jj) .LE. 0.0 .AND. NINT( ssmask(ji,jj) ) .EQ. 1) THEN 
     295                    CALL ctl_stop( 'qco_rst_read: ht_0 must be positive at potentially wet points' ) 
     296                  ENDIF 
    280297               END DO 
     298            END DO 
     299            ! 
     300         ELSE 
     301            ! 
     302            ! Just to read set ssh in fact, called latter once vertical grid 
     303            ! is set up: 
     304!           CALL usr_def_istate( gdept_0, tmask, ts(:,:,:,:,Kbb), uu(:,:,:,Kbb), vv(:,:,:,Kbb), ssh(:,:,Kbb)  ) 
     305!           ! 
     306            ssh(:,:,:) = 0._wp 
     307            ! 
     308         ENDIF           ! end of ll_wd edits 
     309         ! 
     310      ENDIF 
     311      ! 
     312   END SUBROUTINE qco_rst_read 
     313 
     314    
     315   SUBROUTINE qco_rst_read2( kt, Kbb, Kmm ) 
     316      !!--------------------------------------------------------------------- 
     317      !!                   ***  ROUTINE qco_rst_read  *** 
     318      !! 
     319      !! ** Purpose :   Read ssh in restart file 
     320      !! 
     321      !! ** Method  :   use of IOM library 
     322      !!                if the restart does not contain ssh, 
     323      !!                it is set to the _0 values. 
     324      !!---------------------------------------------------------------------- 
     325      INTEGER, INTENT(in) ::   kt         ! ocean time-step 
     326      INTEGER, INTENT(in) ::   Kbb, Kmm   ! ocean time level indices 
     327      ! 
     328      INTEGER ::   ji, jj, jk 
     329      INTEGER ::   id1, id2     ! local integers 
     330      !!---------------------------------------------------------------------- 
     331      ! 
     332      IF( ln_rstart ) THEN                   !* Read the restart file 
     333         CALL rst_read_open                  !  open the restart file if necessary 
     334         ! 
     335         id1 = iom_varid( numror, 'sshb', ldstop = .FALSE. ) 
     336         id2 = iom_varid( numror, 'sshn', ldstop = .FALSE. ) 
     337         ! 
     338         !                             ! --------- ! 
     339         !                             ! all cases ! 
     340         !                             ! --------- ! 
     341         ! 
     342         IF( MIN( id1, id2 ) > 0 ) THEN       ! all required arrays exist 
     343            CALL iom_get( numror, jpdom_auto, 'sshb'   , ssh(:,:,Kbb), ldxios = lrxios    ) 
     344            CALL iom_get( numror, jpdom_auto, 'sshn'   , ssh(:,:,Kmm), ldxios = lrxios    ) 
     345            CALL iom_get( numror, jpdom_auto, 'r3tb'   , r3t(:,:,Kbb), ldxios = lrxios    ) 
     346            CALL iom_get( numror, jpdom_auto, 'r3tn'   , r3t(:,:,Kmm), ldxios = lrxios    ) 
     347            CALL iom_get( numror, jpdom_auto, 'r3ub'   , r3u(:,:,Kbb), ldxios = lrxios, cd_type = 'U'    ) 
     348            CALL iom_get( numror, jpdom_auto, 'r3un'   , r3u(:,:,Kmm), ldxios = lrxios, cd_type = 'U'    ) 
     349            CALL iom_get( numror, jpdom_auto, 'r3vb'   , r3v(:,:,Kbb), ldxios = lrxios, cd_type = 'V'    ) 
     350            CALL iom_get( numror, jpdom_auto, 'r3vn'   , r3v(:,:,Kmm), ldxios = lrxios, cd_type = 'V'    ) 
     351            CALL iom_get( numror, jpdom_auto, 'r3f'    , r3f(:,:)    , ldxios = lrxios, cd_type = 'F'    ) 
     352             
     353            ! needed to restart if land processor not computed 
     354            IF(lwp) write(numout,*) 'qco_rst_read : ssh(:,:,Kbb) and ssh(:,:,Kmm) found in restart files' 
     355            !!WHERE ( ssmask(:,:) == 0.0_wp )   !!gm/st ==> sm should not be necessary on ssh while it was required on e3 
     356            !!   ssh(:,:,Kmm) = 0._wp 
     357            !!   ssh(:,:,Kbb) = 0._wp 
     358            !!END WHERE 
     359            IF( l_1st_euler ) THEN 
     360               ssh(:,:,Kbb) = ssh(:,:,Kmm) 
     361            ENDIF 
     362         ELSE IF( id1 > 0 ) THEN 
     363            IF(lwp) write(numout,*) 'qco_rst_read WARNING : ssh(:,:,Kmm) not found in restart files' 
     364            IF(lwp) write(numout,*) 'sshn set equal to sshb.' 
     365            IF(lwp) write(numout,*) 'neuler is forced to 0' 
     366            CALL iom_get( numror, jpdom_auto, 'sshb', ssh(:,:,Kbb), ldxios = lrxios ) 
     367            ssh(:,:,Kmm) = ssh(:,:,Kbb) 
     368            l_1st_euler = .TRUE. 
     369         ELSE IF( id2 > 0 ) THEN 
     370            IF(lwp) write(numout,*) 'qco_rst_read WARNING : ssh(:,:,Kbb) not found in restart files' 
     371            IF(lwp) write(numout,*) 'sshb set equal to sshn.' 
     372            IF(lwp) write(numout,*) 'neuler is forced to 0' 
     373            CALL iom_get( numror, jpdom_auto, 'sshn', ssh(:,:,Kmm), ldxios = lrxios ) 
     374            ssh(:,:,Kbb) = ssh(:,:,Kmm) 
     375            l_1st_euler = .TRUE. 
     376         ELSE 
     377            IF(lwp) write(numout,*) 'qco_rst_read WARNING : ssh(:,:,Kmm) not found in restart file' 
     378            IF(lwp) write(numout,*) 'ssh_b and ssh_n set to zero' 
     379            IF(lwp) write(numout,*) 'neuler is forced to 0' 
     380            ssh(:,:,:) = 0._wp 
     381            l_1st_euler = .TRUE. 
     382         ENDIF 
     383         ! 
     384      ELSE                                   !* Initialize at "rest" 
     385         ! 
     386         IF( ll_wd ) THEN   ! MJB ll_wd edits start here - these are essential 
     387            ! 
     388            IF( cn_cfg == 'wad' ) THEN            ! Wetting and drying test case 
     389               CALL usr_def_istate( gdept(:,:,:,Kbb), tmask, ts(:,:,:,:,Kbb), uu(:,:,:,Kbb), vv(:,:,:,Kbb), ssh(:,:,Kbb)  ) 
     390               ts (:,:,:,:,Kmm) = ts (:,:,:,:,Kbb)       ! set now values from to before ones 
     391               ssh(:,:    ,Kmm) = ssh(:,:    ,Kbb) 
     392               uu (:,:,:  ,Kmm) = uu (:,:,:  ,Kbb) 
     393               vv (:,:,:  ,Kmm) = vv (:,:,:  ,Kbb) 
     394            ELSE                                  ! if not test case 
     395               ssh(:,:,Kmm) = -ssh_ref 
     396               ssh(:,:,Kbb) = -ssh_ref 
    281397               ! 
    282             ELSE 
    283                ! 
    284                ! Just to read set ssh in fact, called latter once vertical grid 
    285                ! is set up: 
    286 !               CALL usr_def_istate( gdept_0, tmask, ts(:,:,:,:,Kbb), uu(:,:,:,Kbb), vv(:,:,:,Kbb), ssh(:,:,Kbb)  ) 
    287 !               ! 
    288                 ssh(:,:,:) = 0._wp 
    289                ! 
    290             ENDIF           ! end of ll_wd edits 
    291             ! 
    292          ENDIF 
    293       ! 
    294    END SUBROUTINE qe_rst_read 
    295  
    296  
    297    SUBROUTINE dom_qco_ctl 
     398               DO_2D( 1, 1, 1, 1 ) 
     399                  IF( ht_0(ji,jj)-ssh_ref <  rn_wdmin1 ) THEN ! if total depth is less than min depth 
     400                     ssh(ji,jj,Kbb) = rn_wdmin1 - (ht_0(ji,jj) ) 
     401                     ssh(ji,jj,Kmm) = rn_wdmin1 - (ht_0(ji,jj) ) 
     402                  ENDIF 
     403               END_2D 
     404            ENDIF 
     405            ! 
     406            DO ji = 1, jpi 
     407               DO jj = 1, jpj 
     408                  IF ( ht_0(ji,jj) .LE. 0.0 .AND. NINT( ssmask(ji,jj) ) .EQ. 1) THEN 
     409                    CALL ctl_stop( 'qco_rst_read: ht_0 must be positive at potentially wet points' ) 
     410                  ENDIF 
     411               END DO 
     412            END DO 
     413            ! 
     414         ELSE 
     415            ! 
     416            ! Just to read set ssh in fact, called latter once vertical grid 
     417            ! is set up: 
     418!           CALL usr_def_istate( gdept_0, tmask, ts(:,:,:,:,Kbb), uu(:,:,:,Kbb), vv(:,:,:,Kbb), ssh(:,:,Kbb)  ) 
     419!           ! 
     420            ssh(:,:,:) = 0._wp 
     421            r3t(:,:,:) = 0._wp 
     422            r3u(:,:,:) = 0._wp 
     423            r3v(:,:,:) = 0._wp 
     424            r3f(:,:  ) = 0._wp 
     425            ! 
     426         ENDIF           ! end of ll_wd edits 
     427         ! 
     428      ENDIF 
     429      ! 
     430   END SUBROUTINE qco_rst_read2 
     431 
     432 
     433   SUBROUTINE qco_ctl 
    298434      !!--------------------------------------------------------------------- 
    299       !!                  ***  ROUTINE dom_qco_ctl  *** 
     435      !!                  ***  ROUTINE qco_ctl  *** 
    300436      !! 
    301437      !! ** Purpose :   Control the consistency between namelist options 
     
    317453      IF(lwp) THEN                    ! Namelist print 
    318454         WRITE(numout,*) 
    319          WRITE(numout,*) 'dom_qco_ctl : choice/control of the variable vertical coordinate' 
    320          WRITE(numout,*) '~~~~~~~~~~~' 
     455         WRITE(numout,*) 'qco_ctl : choice/control of the variable vertical coordinate' 
     456         WRITE(numout,*) '~~~~~~~~' 
    321457         WRITE(numout,*) '   Namelist nam_vvl : chose a vertical coordinate' 
    322458         WRITE(numout,*) '      zstar                      ln_vvl_zstar   = ', ln_vvl_zstar 
     
    362498#endif 
    363499      ! 
    364    END SUBROUTINE dom_qco_ctl 
     500   END SUBROUTINE qco_ctl 
    365501 
    366502   !!====================================================================== 
  • NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/src/OCE/DOM/domzgr_substitute.h90

    r13237 r13427  
    1919#   define  e3uw(i,j,k,t)  (e3uw_0(i,j,k)*(1._wp+r3u(i,j,t))) 
    2020#   define  e3vw(i,j,k,t)  (e3vw_0(i,j,k)*(1._wp+r3v(i,j,t))) 
    21 #   define  ht(i,j)        (ht_0(i,j)+ssh(i,j,Kmm)) 
     21#   define  ht(i,j)        (ht_0(i,j)*(1._wp+r3t(i,j,Kmm))) 
    2222#   define  hu(i,j,t)      (hu_0(i,j)*(1._wp+r3u(i,j,t))) 
    2323#   define  hv(i,j,t)      (hv_0(i,j)*(1._wp+r3v(i,j,t))) 
     
    2929#endif 
    3030!!---------------------------------------------------------------------- 
     31!!#   define  e3t_f(i,j,k)   (e3t_0(i,j,k)*(1._wp+r3t_f(i,j)*tmask(i,j,k))) 
     32!!#   define  e3u_f(i,j,k)   (e3u_0(i,j,k)*(1._wp+r3u_f(i,j)*umask(i,j,k))) 
     33!!#   define  e3v_f(i,j,k)   (e3v_0(i,j,k)*(1._wp+r3v_f(i,j)*vmask(i,j,k))) 
  • NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/src/OCE/DOM/istate.F90

    r13295 r13427  
    4242   PRIVATE 
    4343 
    44    PUBLIC   istate_init   ! routine called by step.F90 
     44   PUBLIC   istate_init   ! routine called by nemogcm.F90 
    4545 
    4646   !! * Substitutions 
     
    8080!!gm 
    8181 
    82       rhd  (:,:,:  ) = 0._wp   ;   rhop (:,:,:  ) = 0._wp      ! set one for all to 0 at level jpk 
    83       rn2b (:,:,:  ) = 0._wp   ;   rn2  (:,:,:  ) = 0._wp      ! set one for all to 0 at levels 1 and jpk 
    84       ts  (:,:,:,:,Kaa) = 0._wp                                   ! set one for all to 0 at level jpk 
    85       rab_b(:,:,:,:) = 0._wp   ;   rab_n(:,:,:,:) = 0._wp      ! set one for all to 0 at level jpk 
     82      rhd  (:,:,:      ) = 0._wp   ;   rhop (:,:,:  ) = 0._wp      ! set one for all to 0 at level jpk 
     83      rn2b (:,:,:      ) = 0._wp   ;   rn2  (:,:,:  ) = 0._wp      ! set one for all to 0 at levels 1 and jpk 
     84      ts   (:,:,:,:,Kaa) = 0._wp                                   ! set one for all to 0 at level jpk 
     85      rab_b(:,:,:,:    ) = 0._wp   ;   rab_n(:,:,:,:) = 0._wp      ! set one for all to 0 at level jpk 
    8686#if defined key_agrif 
    8787      uu   (:,:,:  ,Kaa) = 0._wp   ! used in agrif_oce_sponge at initialization 
     
    9696         CALL agrif_istate( Kbb, Kmm, Kaa )   ! Interp from parent 
    9797         ! 
    98          ts  (:,:,:,:,Kmm) = ts (:,:,:,:,Kbb)  
    99          ssh (:,:,Kmm)     = ssh(:,:,Kbb) 
    100          uu   (:,:,:,Kmm)   = uu  (:,:,:,Kbb) 
    101          vv   (:,:,:,Kmm)   = vv  (:,:,:,Kbb) 
     98         ts (:,:,:,:,Kmm) = ts (:,:,:,:,Kbb)  
     99         ssh(:,:,    Kmm) = ssh(:,:    ,Kbb) 
     100         uu (:,:,:  ,Kmm) = uu (:,:,:  ,Kbb) 
     101         vv (:,:,:  ,Kmm) = vv (:,:,:  ,Kbb) 
    102102      ELSE 
    103103#endif 
     
    117117            CALL dta_tsd( nit000, ts(:,:,:,:,Kbb) )       ! read 3D T and S data at nit000 
    118118            ! 
    119             ssh(:,:,Kbb)  = 0._wp               ! set the ocean at rest 
    120             uu  (:,:,:,Kbb) = 0._wp 
    121             vv  (:,:,:,Kbb) = 0._wp   
     119            ssh(:,:  ,Kbb) = 0._wp               ! set the ocean at rest 
     120            uu (:,:,:,Kbb) = 0._wp 
     121            vv (:,:,:,Kbb) = 0._wp   
    122122            ! 
    123123            IF( ll_wd ) THEN 
     
    139139            CALL usr_def_istate( zgdept, tmask, ts(:,:,:,:,Kbb), uu(:,:,:,Kbb), vv(:,:,:,Kbb), ssh(:,:,Kbb)  )          
    140140         ENDIF 
    141          ts  (:,:,:,:,Kmm) = ts (:,:,:,:,Kbb)       ! set now values from to before ones 
    142          ssh (:,:,Kmm)     = ssh(:,:,Kbb)    
    143          uu   (:,:,:,Kmm)   = uu  (:,:,:,Kbb) 
    144          vv   (:,:,:,Kmm)   = vv  (:,:,:,Kbb) 
     141         ts (:,:,:,:,Kmm) = ts (:,:,:,:,Kbb)       ! set now values from to before ones 
     142         ssh(:,:    ,Kmm) = ssh(:,:    ,Kbb)    
     143         uu (:,:,:  ,Kmm) = uu (:,:,:  ,Kbb) 
     144         vv (:,:,:  ,Kmm) = vv (:,:,:  ,Kbb) 
    145145 
    146146!!gm POTENTIAL BUG : 
  • NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/src/OCE/DYN/dynadv.F90

    r12377 r13427  
    127127      IF( ioptio /= 1 )   CALL ctl_stop( 'choose ONE and only ONE advection scheme' ) 
    128128      IF( nn_dynkeg /= nkeg_C2 .AND. nn_dynkeg /= nkeg_HW )   CALL ctl_stop( 'KEG scheme wrong value of nn_dynkeg' ) 
    129  
     129#if defined key_qcoTest_FluxForm 
     130      IF( ln_dynadv_vec  ) THEN CALL ctl_stop( 'STOP', 'key_qcoTest_FluxForm requires flux form advection' ) 
     131#endif 
    130132 
    131133      IF(lwp) THEN                    ! Print the choice 
  • NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/src/OCE/DYN/dynatf_qco.F90

    r13295 r13427  
    199199      ! JC: Would be more clever to swap variables than to make a full vertical 
    200200      ! integration 
    201       ! 
     201      ! CAUTION : calculation need to be done in the same way than see GM   
    202202      uu_b(:,:,Kaa) = e3u(:,:,1,Kaa) * puu(:,:,1,Kaa) * umask(:,:,1) 
    203       uu_b(:,:,Kmm) = e3u(:,:,1,Kmm) * puu(:,:,1,Kmm) * umask(:,:,1) 
     203      uu_b(:,:,Kmm) = (e3u_0(:,:,1) * ( 1._wp + r3u_f(:,:) * umask(:,:,1) )) * puu(:,:,1,Kmm) * umask(:,:,1) 
    204204      vv_b(:,:,Kaa) = e3v(:,:,1,Kaa) * pvv(:,:,1,Kaa) * vmask(:,:,1) 
    205       vv_b(:,:,Kmm) = e3v(:,:,1,Kmm) * pvv(:,:,1,Kmm) * vmask(:,:,1) 
     205      vv_b(:,:,Kmm) = (e3v_0(:,:,1) * ( 1._wp + r3v_f(:,:) * vmask(:,:,1))) * pvv(:,:,1,Kmm) * vmask(:,:,1) 
    206206      DO jk = 2, jpkm1 
    207207         uu_b(:,:,Kaa) = uu_b(:,:,Kaa) + e3u(:,:,jk,Kaa) * puu(:,:,jk,Kaa) * umask(:,:,jk) 
    208          uu_b(:,:,Kmm) = uu_b(:,:,Kmm) + e3u(:,:,jk,Kmm) * puu(:,:,jk,Kmm) * umask(:,:,jk) 
     208         uu_b(:,:,Kmm) = uu_b(:,:,Kmm) + (e3u_0(:,:,jk) * ( 1._wp + r3u_f(:,:) * umask(:,:,jk) )) * puu(:,:,jk,Kmm) * umask(:,:,jk) 
    209209         vv_b(:,:,Kaa) = vv_b(:,:,Kaa) + e3v(:,:,jk,Kaa) * pvv(:,:,jk,Kaa) * vmask(:,:,jk) 
    210          vv_b(:,:,Kmm) = vv_b(:,:,Kmm) + e3v(:,:,jk,Kmm) * pvv(:,:,jk,Kmm) * vmask(:,:,jk) 
     210         vv_b(:,:,Kmm) = vv_b(:,:,Kmm) + (e3v_0(:,:,jk) * ( 1._wp + r3v_f(:,:) * vmask(:,:,jk) )) * pvv(:,:,jk,Kmm) * vmask(:,:,jk) 
    211211      END DO 
    212212      uu_b(:,:,Kaa) = uu_b(:,:,Kaa) * r1_hu(:,:,Kaa) 
    213213      vv_b(:,:,Kaa) = vv_b(:,:,Kaa) * r1_hv(:,:,Kaa) 
    214       uu_b(:,:,Kmm) = uu_b(:,:,Kmm) * r1_hu(:,:,Kmm) 
    215       vv_b(:,:,Kmm) = vv_b(:,:,Kmm) * r1_hv(:,:,Kmm) 
     214      uu_b(:,:,Kmm) = uu_b(:,:,Kmm) * (r1_hu_0(:,:)/( 1._wp + r3u_f(:,:) )) 
     215      vv_b(:,:,Kmm) = vv_b(:,:,Kmm) * (r1_hv_0(:,:)/( 1._wp + r3v_f(:,:) )) 
    216216      ! 
    217217      IF( .NOT.ln_dynspg_ts ) THEN        ! output the barotropic currents 
  • NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/src/OCE/DYN/dynspg_ts.F90

    r13295 r13427  
    306306      ENDIF 
    307307      ! 
    308       !                                   !=  Add atmospheric pressure forcing  =! 
    309       !                                   !  ----------------------------------  ! 
    310       IF( ln_bt_fw ) THEN                        ! Add wind forcing 
     308      !                                   !=  Add wind forcing  =! 
     309      !                                   !  ------------------  ! 
     310      IF( ln_bt_fw ) THEN 
    311311         DO_2D( 0, 0, 0, 0 ) 
    312312            zu_frc(ji,jj) =  zu_frc(ji,jj) + r1_rho0 * utau(ji,jj) * r1_hu(ji,jj,Kmm) 
     
    386386      ! 
    387387      IF( ln_linssh ) THEN    ! mid-step ocean depth is fixed (hup2_e=hu_n=hu_0) 
    388          zhup2_e(:,:) = hu(:,:,Kmm) 
    389          zhvp2_e(:,:) = hv(:,:,Kmm) 
    390          zhtp2_e(:,:) = ht(:,:) 
    391       ENDIF 
    392       ! 
    393       IF (ln_bt_fw) THEN                  ! FORWARD integration: start from NOW fields                     
    394          sshn_e(:,:) =    pssh(:,:,Kmm)             
     388         zhup2_e(:,:) = hu_0(:,:) 
     389         zhvp2_e(:,:) = hv_0(:,:) 
     390         zhtp2_e(:,:) = ht_0(:,:) 
     391      ENDIF 
     392      ! 
     393      IF( ln_bt_fw ) THEN                 ! FORWARD integration: start from NOW fields                     
     394         sshn_e(:,:) =    pssh (:,:,Kmm)             
    395395         un_e  (:,:) =    puu_b(:,:,Kmm)             
    396396         vn_e  (:,:) =    pvv_b(:,:,Kmm) 
     
    401401         hvr_e (:,:) = r1_hv(:,:,Kmm) 
    402402      ELSE                                ! CENTRED integration: start from BEFORE fields 
    403          sshn_e(:,:) =    pssh(:,:,Kbb) 
     403         sshn_e(:,:) =    pssh (:,:,Kbb) 
    404404         un_e  (:,:) =    puu_b(:,:,Kbb)          
    405405         vn_e  (:,:) =    pvv_b(:,:,Kbb) 
     
    412412      ! 
    413413      ! Initialize sums: 
    414       puu_b  (:,:,Kaa) = 0._wp       ! After barotropic velocities (or transport if flux form)           
    415       pvv_b  (:,:,Kaa) = 0._wp 
     414      puu_b (:,:,Kaa) = 0._wp       ! After barotropic velocities (or transport if flux form)           
     415      pvv_b (:,:,Kaa) = 0._wp 
    416416      pssh  (:,:,Kaa) = 0._wp       ! Sum for after averaged sea level 
    417       un_adv(:,:) = 0._wp       ! Sum for now transport issued from ts loop 
    418       vn_adv(:,:) = 0._wp 
     417      un_adv(:,:)     = 0._wp       ! Sum for now transport issued from ts loop 
     418      vn_adv(:,:)     = 0._wp 
    419419      ! 
    420420      IF( ln_wd_dl ) THEN 
     
    475475            ! 
    476476            !                          ! ocean u- and v-depth at mid-step   (separate DO-loops remove the need of a lbc_lnk) 
     477#if defined key_qcoTest_FluxForm 
     478            !                                ! 'key_qcoTest_FluxForm' : simple ssh average 
     479            DO_2D( 1, 1, 1, 0 ) 
     480               zhup2_e(ji,jj) = hu_0(ji,jj) + r1_2 * (  zsshp2_e(ji,jj) + zsshp2_e(ji+1,jj  )  ) * ssumask(ji,jj) 
     481            END_2D 
     482            DO_2D( 1, 0, 1, 1 ) 
     483               zhvp2_e(ji,jj) = hv_0(ji,jj) + r1_2 * (  zsshp2_e(ji,jj) + zsshp2_e(ji  ,jj+1)  ) * ssvmask(ji,jj) 
     484            END_2D 
     485#else 
     486            !                                ! no 'key_qcoTest_FluxForm' : surface weighted ssh average 
    477487            DO_2D( 1, 1, 1, 0 ) 
    478488               zhup2_e(ji,jj) = hu_0(ji,jj) + r1_2 * r1_e1e2u(ji,jj)                        & 
     
    485495                    &                                 + e1e2t(ji,jj+1) * zsshp2_e(ji,jj+1)  ) * ssvmask(ji,jj) 
    486496            END_2D 
     497#endif                
    487498            ! 
    488499         ENDIF 
     
    540551         !   
    541552         ! Sea Surface Height at u-,v-points (vvl case only) 
    542          IF( .NOT.ln_linssh ) THEN                                 
     553         IF( .NOT.ln_linssh ) THEN 
     554#if defined key_qcoTest_FluxForm 
     555            !                                ! 'key_qcoTest_FluxForm' : simple ssh average 
     556            DO_2D( 1, 1, 1, 0 ) 
     557               zsshu_a(ji,jj) = r1_2 * (  ssha_e(ji,jj) + ssha_e(ji+1,jj  )  ) * ssumask(ji,jj) 
     558            END_2D 
     559            DO_2D( 1, 0, 1, 1 ) 
     560               zsshv_a(ji,jj) = r1_2 * (  ssha_e(ji,jj) + ssha_e(ji  ,jj+1)  ) * ssvmask(ji,jj) 
     561            END_2D 
     562#else 
    543563            DO_2D( 0, 0, 0, 0 ) 
    544                zsshu_a(ji,jj) = r1_2 * ssumask(ji,jj) * r1_e1e2u(ji,jj)    & 
    545                   &              * ( e1e2t(ji  ,jj  )  * ssha_e(ji  ,jj  ) & 
    546                   &              +   e1e2t(ji+1,jj  )  * ssha_e(ji+1,jj  ) ) 
    547                zsshv_a(ji,jj) = r1_2 * ssvmask(ji,jj) * r1_e1e2v(ji,jj)    & 
    548                   &              * ( e1e2t(ji  ,jj  )  * ssha_e(ji  ,jj  ) & 
    549                   &              +   e1e2t(ji  ,jj+1)  * ssha_e(ji  ,jj+1) ) 
    550             END_2D 
    551          ENDIF    
     564               zsshu_a(ji,jj) = r1_2 * r1_e1e2u(ji,jj) * ( e1e2t(ji  ,jj  ) * ssha_e(ji  ,jj  )   & 
     565                  &                                      + e1e2t(ji+1,jj  ) * ssha_e(ji+1,jj  ) ) * ssumask(ji,jj) 
     566               zsshv_a(ji,jj) = r1_2 * r1_e1e2v(ji,jj) * ( e1e2t(ji  ,jj  ) * ssha_e(ji  ,jj  )   & 
     567                  &                                      + e1e2t(ji  ,jj+1) * ssha_e(ji  ,jj+1) ) * ssvmask(ji,jj) 
     568            END_2D 
     569#endif 
     570         ENDIF 
    552571         !          
    553572         ! Half-step back interpolation of SSH for surface pressure computation at step jit+1/2 
     
    624643               !                    ! hu_e, hv_e hold depth at jn,  zhup2_e, zhvp2_e hold extrapolated depth at jn+1/2 
    625644               !                    ! backward interpolated depth used in spg terms at jn+1/2 
     645#if defined key_qcoTest_FluxForm 
     646            !                                ! 'key_qcoTest_FluxForm' : simple ssh average 
     647               zhu_bck = hu_0(ji,jj) + r1_2 * (  zsshp2_e(ji,jj) + zsshp2_e(ji+1,jj  )  ) * ssumask(ji,jj) 
     648               zhv_bck = hv_0(ji,jj) + r1_2 * (  zsshp2_e(ji,jj) + zsshp2_e(ji  ,jj+1)  ) * ssvmask(ji,jj) 
     649#else 
    626650               zhu_bck = hu_0(ji,jj) + r1_2*r1_e1e2u(ji,jj) * (  e1e2t(ji  ,jj) * zsshp2_e(ji  ,jj)    & 
    627651                    &                                          + e1e2t(ji+1,jj) * zsshp2_e(ji+1,jj)  ) * ssumask(ji,jj) 
    628652               zhv_bck = hv_0(ji,jj) + r1_2*r1_e1e2v(ji,jj) * (  e1e2t(ji,jj  ) * zsshp2_e(ji,jj  )    & 
    629653                    &                                          + e1e2t(ji,jj+1) * zsshp2_e(ji,jj+1)  ) * ssvmask(ji,jj) 
     654#endif 
    630655               !                    ! inverse depth at jn+1 
    631656               z1_hu = ssumask(ji,jj) / ( hu_0(ji,jj) + zsshu_a(ji,jj) + 1._wp - ssumask(ji,jj) ) 
     
    646671         IF ( ll_wd ) THEN ! revert to explicit for bit comparison tests in non wad runs 
    647672            DO_2D( 0, 0, 0, 0 ) 
    648                   ua_e(ji,jj) =  ua_e(ji,jj) /(1.0 -   rDt_e * zCdU_u(ji,jj) * hur_e(ji,jj)) 
    649                   va_e(ji,jj) =  va_e(ji,jj) /(1.0 -   rDt_e * zCdU_v(ji,jj) * hvr_e(ji,jj)) 
     673               ua_e(ji,jj) =  ua_e(ji,jj) / ( 1._wp - rDt_e * zCdU_u(ji,jj) * hur_e(ji,jj) ) 
     674               va_e(ji,jj) =  va_e(ji,jj) / ( 1._wp - rDt_e * zCdU_v(ji,jj) * hvr_e(ji,jj) ) 
    650675            END_2D 
    651676         ENDIF 
    652677        
    653678         IF( .NOT.ln_linssh ) THEN !* Update ocean depth (variable volume case only) 
    654             hu_e (2:jpim1,2:jpjm1) = hu_0(2:jpim1,2:jpjm1) + zsshu_a(2:jpim1,2:jpjm1) 
    655             hur_e(2:jpim1,2:jpjm1) = ssumask(2:jpim1,2:jpjm1) / ( hu_e(2:jpim1,2:jpjm1) + 1._wp - ssumask(2:jpim1,2:jpjm1) ) 
    656             hv_e (2:jpim1,2:jpjm1) = hv_0(2:jpim1,2:jpjm1) + zsshv_a(2:jpim1,2:jpjm1) 
    657             hvr_e(2:jpim1,2:jpjm1) = ssvmask(2:jpim1,2:jpjm1) / ( hv_e(2:jpim1,2:jpjm1) + 1._wp - ssvmask(2:jpim1,2:jpjm1) ) 
     679            hu_e (2:jpim1,2:jpjm1) =    hu_0(2:jpim1,2:jpjm1) + zsshu_a(2:jpim1,2:jpjm1) 
     680            hur_e(2:jpim1,2:jpjm1) = ssumask(2:jpim1,2:jpjm1) / (  hu_e(2:jpim1,2:jpjm1) + 1._wp - ssumask(2:jpim1,2:jpjm1) ) 
     681            hv_e (2:jpim1,2:jpjm1) =    hv_0(2:jpim1,2:jpjm1) + zsshv_a(2:jpim1,2:jpjm1) 
     682            hvr_e(2:jpim1,2:jpjm1) = ssvmask(2:jpim1,2:jpjm1) / (  hv_e(2:jpim1,2:jpjm1) + 1._wp - ssvmask(2:jpim1,2:jpjm1) ) 
    658683         ENDIF 
    659684         ! 
     
    743768      ELSE 
    744769         ! At this stage, pssh(:,:,:,Krhs) has been corrected: compute new depths at velocity points 
     770#if defined key_qcoTest_FluxForm 
     771         !                                ! 'key_qcoTest_FluxForm' : simple ssh average 
    745772         DO_2D( 1, 0, 1, 0 ) 
    746             zsshu_a(ji,jj) = r1_2 * ssumask(ji,jj)  * r1_e1e2u(ji,jj) & 
    747                &              * ( e1e2t(ji  ,jj) * pssh(ji  ,jj,Kaa)      & 
    748                &              +   e1e2t(ji+1,jj) * pssh(ji+1,jj,Kaa) ) 
    749             zsshv_a(ji,jj) = r1_2 * ssvmask(ji,jj)  * r1_e1e2v(ji,jj) & 
    750                &              * ( e1e2t(ji,jj  ) * pssh(ji,jj  ,Kaa)      & 
    751                &              +   e1e2t(ji,jj+1) * pssh(ji,jj+1,Kaa) ) 
    752          END_2D 
     773            zsshu_a(ji,jj) = r1_2 * ( pssh(ji,jj,Kaa) + pssh(ji+1,jj  ,Kaa) ) * ssumask(ji,jj) 
     774            zsshv_a(ji,jj) = r1_2 * ( pssh(ji,jj,Kaa) + pssh(ji  ,jj+1,Kaa) ) * ssvmask(ji,jj) 
     775         END_2D 
     776#else 
     777         DO_2D( 1, 0, 1, 0 ) 
     778            zsshu_a(ji,jj) = r1_2 * r1_e1e2u(ji,jj) * ( e1e2t(ji  ,jj) * pssh(ji  ,jj,Kaa)   & 
     779               &                                      + e1e2t(ji+1,jj) * pssh(ji+1,jj,Kaa) ) * ssumask(ji,jj) 
     780            zsshv_a(ji,jj) = r1_2 * r1_e1e2v(ji,jj) * ( e1e2t(ji,jj  ) * pssh(ji,jj  ,Kaa)   & 
     781               &                                      + e1e2t(ji,jj+1) * pssh(ji,jj+1,Kaa) ) * ssvmask(ji,jj) 
     782         END_2D 
     783#endif    
    753784         CALL lbc_lnk_multi( 'dynspg_ts', zsshu_a, 'U', 1._wp, zsshv_a, 'V', 1._wp ) ! Boundary conditions 
    754785         ! 
    755786         DO jk=1,jpkm1 
    756             puu(:,:,jk,Krhs) = puu(:,:,jk,Krhs) + r1_hu(:,:,Kmm) * ( puu_b(:,:,Kaa) - puu_b(:,:,Kbb) * hu(:,:,Kbb) ) * r1_Dt_b 
    757             pvv(:,:,jk,Krhs) = pvv(:,:,jk,Krhs) + r1_hv(:,:,Kmm) * ( pvv_b(:,:,Kaa) - pvv_b(:,:,Kbb) * hv(:,:,Kbb) ) * r1_Dt_b 
     787            puu(:,:,jk,Krhs) = puu(:,:,jk,Krhs) + r1_hu(:,:,Kmm)   & 
     788               &             * ( puu_b(:,:,Kaa) - puu_b(:,:,Kbb) * hu(:,:,Kbb) ) * r1_Dt_b 
     789            pvv(:,:,jk,Krhs) = pvv(:,:,jk,Krhs) + r1_hv(:,:,Kmm)   & 
     790               &             * ( pvv_b(:,:,Kaa) - pvv_b(:,:,Kbb) * hv(:,:,Kbb) ) * r1_Dt_b 
    758791         END DO 
    759792         ! Save barotropic velocities not transport: 
     
    11011134         CASE ( 0 )                                   ! original formulation  (masked averaging of e3t divided by 4) 
    11021135            DO_2D( 1, 0, 1, 0 ) 
    1103                zwz(ji,jj) =   ( ht(ji  ,jj+1) + ht(ji+1,jj+1) +                    & 
    1104                     &           ht(ji  ,jj  ) + ht(ji+1,jj  )  ) * 0.25_wp   
     1136               zwz(ji,jj) = ( ht(ji,jj+1) + ht(ji+1,jj+1)   & 
     1137                    &       + ht(ji,jj  ) + ht(ji+1,jj  ) ) * 0.25_wp   
    11051138               IF( zwz(ji,jj) /= 0._wp )   zwz(ji,jj) = ff_f(ji,jj) / zwz(ji,jj) 
    11061139            END_2D 
    11071140         CASE ( 1 )                                   ! new formulation  (masked averaging of e3t divided by the sum of mask) 
    11081141            DO_2D( 1, 0, 1, 0 ) 
    1109                zwz(ji,jj) =             (  ht  (ji  ,jj+1) + ht  (ji+1,jj+1)      & 
    1110                     &                    + ht  (ji  ,jj  ) + ht  (ji+1,jj  )  )   & 
    1111                     &       / ( MAX( 1._wp,  ssmask(ji  ,jj+1) + ssmask(ji+1,jj+1)      & 
    1112                     &                      + ssmask(ji  ,jj  ) + ssmask(ji+1,jj  )  )   ) 
     1142               zwz(ji,jj) =     (    ht(ji,jj+1) +     ht(ji+1,jj+1)      & 
     1143                    &            +   ht(ji,jj  ) +     ht(ji+1,jj  )  )   & 
     1144                    &    / ( MAX(ssmask(ji,jj+1) + ssmask(ji+1,jj+1)      & 
     1145                    &          + ssmask(ji,jj  ) + ssmask(ji+1,jj  ) , 1._wp  )   ) 
    11131146               IF( zwz(ji,jj) /= 0._wp )   zwz(ji,jj) = ff_f(ji,jj) / zwz(ji,jj) 
    11141147            END_2D 
     
    11161149         CALL lbc_lnk( 'dynspg_ts', zwz, 'F', 1._wp ) 
    11171150         ! 
    1118          ftne(1,:) = 0._wp ; ftnw(1,:) = 0._wp ; ftse(1,:) = 0._wp ; ftsw(1,:) = 0._wp 
     1151         ftne(1,:) = 0._wp   ;   ftnw(1,:) = 0._wp   ;   ftse(1,:) = 0._wp   ;  ftsw(1,:) = 0._wp 
    11191152         DO_2D( 0, 1, 0, 1 ) 
    11201153            ftne(ji,jj) = zwz(ji-1,jj  ) + zwz(ji  ,jj  ) + zwz(ji  ,jj-1) 
     
    11251158         ! 
    11261159      CASE( np_EET )                  != EEN scheme using e3t energy conserving scheme 
    1127          ftne(1,:) = 0._wp ; ftnw(1,:) = 0._wp ; ftse(1,:) = 0._wp ; ftsw(1,:) = 0._wp 
     1160         ftne(1,:) = 0._wp   ;   ftnw(1,:) = 0._wp   ;   ftse(1,:) = 0._wp   ;  ftsw(1,:) = 0._wp 
    11281161         DO_2D( 0, 1, 0, 1 ) 
    11291162            z1_ht = ssmask(ji,jj) / ( ht(ji,jj) + 1._wp - ssmask(ji,jj) ) 
     
    11671200         ENDIF 
    11681201         ! 
    1169          DO jj = 1, jpjm1 
     1202         DO jj = 1, jpjm1   ! keep only the value at the coast if sco 
    11701203            zhf(:,jj) = zhf(:,jj) * (1._wp- umask(:,jj,1) * umask(:,jj+1,1)) 
    11711204         END DO 
    11721205         ! 
    1173          DO jk = 1, jpkm1 
     1206         DO jk = 1, jpkm1   ! ocean point : sum of masked e3f 
    11741207            DO jj = 1, jpjm1 
    11751208               zhf(:,jj) = zhf(:,jj) + e3f(:,jj,jk) * umask(:,jj,jk) * umask(:,jj+1,jk) 
  • NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/src/OCE/DYN/sshwzv.F90

    r13295 r13427  
    299299         !                                                  ! filtered "now" field 
    300300         pssh(:,:,Kmm) = pssh(:,:,Kmm) + rn_atfp * ( pssh(:,:,Kbb) - 2 * pssh(:,:,Kmm) + pssh(:,:,Kaa) ) 
     301         ! 
    301302         IF( .NOT.ln_linssh ) THEN                          ! "now" <-- with forcing removed 
    302303            zcoef = rn_atfp * rn_Dt * r1_rho0 
     
    307308 
    308309            ! ice sheet coupling 
    309             IF ( ln_isf .AND. ln_isfcpl .AND. kt == nit000+1) pssh(:,:,Kbb) = pssh(:,:,Kbb) - rn_atfp * rn_Dt * ( risfcpl_ssh(:,:) - 0.0 ) * ssmask(:,:) 
     310            IF( ln_isf .AND. ln_isfcpl .AND. kt == nit000+1 )   & 
     311               &   pssh(:,:,Kbb) = pssh(:,:,Kbb) - rn_atfp * rn_Dt * ( risfcpl_ssh(:,:) - 0.0 ) * ssmask(:,:) 
    310312 
    311313         ENDIF 
    312314      ENDIF 
    313315      ! 
    314       IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab2d_1=pssh(:,:,Kmm), clinfo1=' pssh(:,:,Kmm)  - : ', mask1=tmask ) 
     316      IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab2d_1=pssh(:,:,Kmm), clinfo1=' atf  - pssh(:,:,Kmm): ', mask1=tmask ) 
    315317      ! 
    316318      IF( ln_timing )   CALL timing_stop('ssh_atf') 
  • NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/src/OCE/IOM/restart.F90

    r13286 r13427  
    160160                     CALL iom_rstput( kt, nitrst, numrow, 'sshn'   ,ssh(:,:         ,Kmm), ldxios = lwxios      ) 
    161161                     CALL iom_rstput( kt, nitrst, numrow, 'rhop'   , rhop, ldxios = lwxios      ) 
     162!!st uncomment in case of dom_qco_read2 usage 
     163!!#if defined key_qco 
     164!!                     CALL iom_rstput( kt, nitrst, numrow, 'r3tb'   , r3t(:,:,Kbb) , ldxios = lwxios      ) 
     165!!                     CALL iom_rstput( kt, nitrst, numrow, 'r3ub'   , r3u(:,:,Kbb) , ldxios = lwxios      ) 
     166!!                     CALL iom_rstput( kt, nitrst, numrow, 'r3vb'   , r3v(:,:,Kbb) , ldxios = lwxios      )                      
     167!!                     CALL iom_rstput( kt, nitrst, numrow, 'r3tn'   , r3t(:,:,Kmm) , ldxios = lwxios      ) 
     168!!                     CALL iom_rstput( kt, nitrst, numrow, 'r3un'   , r3u(:,:,Kmm) , ldxios = lwxios      ) 
     169!!                     CALL iom_rstput( kt, nitrst, numrow, 'r3vn'   , r3v(:,:,Kmm) , ldxios = lwxios      ) 
     170!!                     CALL iom_rstput( kt, nitrst, numrow, 'r3f'   , r3f(:,:) , ldxios = lwxios      ) 
     171!!#endif  
    162172      ENDIF 
    163173       
     
    275285         CALL iom_get( numror, jpdom_auto, 'tb'     , ts(:,:,:,jp_tem,Kbb), ldxios = lrxios ) 
    276286         CALL iom_get( numror, jpdom_auto, 'sb'     , ts(:,:,:,jp_sal,Kbb), ldxios = lrxios ) 
     287#if ! defined key_qco 
    277288         CALL iom_get( numror, jpdom_auto, 'sshb'   ,ssh(:,:         ,Kbb), ldxios = lrxios ) 
     289#endif 
    278290      ELSE 
    279291         l_1st_euler =  .TRUE.      ! before field not found, forced euler 1st time-step 
     
    285297      CALL iom_get( numror, jpdom_auto, 'tn'     , ts(:,:,:,jp_tem,Kmm), ldxios = lrxios ) 
    286298      CALL iom_get( numror, jpdom_auto, 'sn'     , ts(:,:,:,jp_sal,Kmm), ldxios = lrxios ) 
     299#if ! defined key_qco 
    287300      CALL iom_get( numror, jpdom_auto, 'sshn'   ,ssh(:,:         ,Kmm), ldxios = lrxios ) 
     301#endif 
    288302      IF( iom_varid( numror, 'rhop', ldstop = .FALSE. ) > 0 ) THEN 
    289303         CALL iom_get( numror, jpdom_auto, 'rhop'   , rhop, ldxios = lrxios )   ! now    potential density 
  • NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/src/OCE/SBC/sbcice_cice.F90

    r13295 r13427  
    1212   USE oce             ! ocean dynamics and tracers 
    1313   USE dom_oce         ! ocean space and time domain 
    14 # if ! defined key_qco 
    15    USE domvvl 
     14# if defined key_qco 
     15   USE domqco         ! Variable volume 
    1616# else 
    17    USE domqco 
     17   USE domvvl         ! Variable volume 
    1818# endif 
    1919   USE phycst, only : rcp, rho0, r1_rho0, rhos, rhoi 
  • NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/src/OCE/ZDF/zdfddm.F90

    r13295 r13427  
    3131   !! * Substitutions 
    3232#  include "do_loop_substitute.h90" 
     33#  include "domzgr_substitute.h90" 
    3334   !!---------------------------------------------------------------------- 
    3435   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
  • NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/src/OCE/nemogcm.F90

    r13286 r13427  
    183183      ! 
    184184      DO WHILE( istp <= nitend .AND. nstop == 0 ) 
    185 #if defined key_qco 
     185#  if defined key_qco 
    186186         CALL stp_MLF 
    187 #else 
     187#  else 
    188188         CALL stp 
    189 #endif 
     189#  endif 
    190190         istp = istp + 1 
    191191      END DO 
     
    204204            ENDIF 
    205205             
    206 #if defined key_qco 
     206#  if defined key_qco 
    207207            CALL stp_MLF      ( istp ) 
    208 #else 
     208#  else 
    209209            CALL stp        ( istp )  
    210 #endif 
     210#  endif 
    211211            istp = istp + 1 
    212212 
     
    416416 
    417417      ! Initialise time level indices 
    418       Nbb = 1; Nnn = 2; Naa = 3; Nrhs = Naa 
     418      Nbb = 1   ;   Nnn = 2   ;   Naa = 3   ;  Nrhs = Naa 
    419419#if defined key_agrif 
    420       Kbb_a = Nbb; Kmm_a = Nnn; Krhs_a = Nrhs   ! agrif_oce module copies of time level indices 
     420      Kbb_a = Nbb   ;   Kmm_a = Nnn   ;  Krhs_a = Nrhs   ! agrif_oce module copies of time level indices 
    421421#endif  
    422422      !                             !-------------------------------! 
  • NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/src/OCE/stpMLF.F90

    r13237 r13427  
    3232   !!            4.0  !  2017-05  (G. Madec)  introduction of the vertical physics manager (zdfphy) 
    3333   !!            4.1  !  2019-08  (A. Coward, D. Storkey) rewrite in preparation for new timestepping scheme 
    34    !!---------------------------------------------------------------------- 
    35  
    36    !!---------------------------------------------------------------------- 
    37    !!   stp_MLF       : OPA system time-stepping 
     34   !!            4.x  !  2020-08  (S. Techene, G. Madec)  quasi eulerian coordinate time stepping  
     35   !!---------------------------------------------------------------------- 
     36 
     37#if defined key_qco 
     38   !!---------------------------------------------------------------------- 
     39   !!   'key_qco'       Quasi-Eulerian vertical coordonate 
     40   !!---------------------------------------------------------------------- 
     41    
     42   !!---------------------------------------------------------------------- 
     43   !!   stp_MLF       : NEMO modified Leap Frog time-stepping with qco 
    3844   !!---------------------------------------------------------------------- 
    3945   USE step_oce       ! time stepping definition modules 
     
    4349   USE traatfqco      ! time filtering                   (tra_atf_qco routine) 
    4450   USE dynatfqco      ! time filtering                   (dyn_atf_qco routine) 
    45    USE dynspg_ts      ! surface pressure gradient: split-explicit scheme (define un_adv) 
     51    
    4652   USE bdydyn         ! ocean open boundary conditions (define bdy_dyn) 
    4753 
     
    4955   PRIVATE 
    5056    
    51 #if ! defined key_qco 
    52    !!---------------------------------------------------------------------- 
    53    !!   'key_qco'      EMPTY MODULE      Quasi-Eulerian vertical coordonate 
    54    !!---------------------------------------------------------------------- 
    55 #else 
    5657   PUBLIC   stp_MLF   ! called by nemogcm.F90 
    5758 
     
    195196         zgdept(:,:,jk) = gdept(:,:,jk,Nnn) 
    196197      END DO 
    197                             CALL ssh_nxt       ( kstp, Nbb, Nnn, ssh,  Naa )   ! after ssh (includes call to div_hor) 
    198       IF( .NOT.ln_linssh )  CALL dom_qco_r3c    ( ssh(:,:,Naa), r3t(:,:,Naa), r3u(:,:,Naa), r3v(:,:,Naa) )   ! "after" ssh./h._0 ratio 
    199                             CALL wzv           ( kstp, Nbb, Nnn, Naa, ww  )    ! Nnn cross-level velocity 
    200       IF( ln_zad_Aimp )     CALL wAimp         ( kstp,      Nnn           )    ! Adaptive-implicit vertical advection partitioning 
    201                             CALL eos    ( ts(:,:,:,:,Nnn), rhd, rhop, zgdept ) ! now in situ density for hpg computation 
     198                            CALL ssh_nxt    ( kstp, Nbb, Nnn, ssh,  Naa )   ! after ssh (includes call to div_hor) 
     199      IF( .NOT.ln_linssh )  THEN 
     200                             CALL dom_qco_r3c( ssh(:,:,Naa), r3t(:,:,Naa), r3u(:,:,Naa), r3v(:,:,Naa)           )   ! "after" ssh/h_0 ratio at t,u,v pts 
     201         IF( ln_dynspg_exp ) CALL dom_qco_r3c( ssh(:,:,Nnn), r3t(:,:,Nnn), r3u(:,:,Nnn), r3v(:,:,Nnn), r3f(:,:) )   ! spg_exp : needed only for "now" ssh/h_0 ratio at f point 
     202      ENDIF 
     203                            CALL wzv        ( kstp, Nbb, Nnn, Naa, ww  )    ! Nnn cross-level velocity 
     204      IF( ln_zad_Aimp )     CALL wAimp      ( kstp,      Nnn           )    ! Adaptive-implicit vertical advection partitioning 
     205                            CALL eos        ( ts(:,:,:,:,Nnn), rhd, rhop, zgdept ) ! now in situ density for hpg computation 
    202206 
    203207 
     
    218222                         CALL dyn_hpg( kstp,      Nnn      , uu, vv, Nrhs )  ! horizontal gradient of Hydrostatic pressure 
    219223                         CALL dyn_spg( kstp, Nbb, Nnn, Nrhs, uu, vv, ssh, uu_b, vv_b, Naa )  ! surface pressure gradient 
    220  
    221224                                                      ! With split-explicit free surface, since now transports have been updated and ssh(:,:,Nrhs) as well 
     225 
    222226      IF( ln_dynspg_ts ) THEN                         ! vertical scale factors and vertical velocity need to be updated 
    223227                            CALL div_hor    ( kstp, Nbb, Nnn )                ! Horizontal divergence  (2nd call in time-split case) 
    224          IF(.NOT.ln_linssh) CALL dom_qco_r3c ( ssh(:,:,Naa), r3t(:,:,Naa), r3u(:,:,Naa), r3v(:,:,Naa), r3f(:,:) )   ! "after" ssh./h._0 ratio 
     228         IF(.NOT.ln_linssh) CALL dom_qco_r3c ( ssh(:,:,Naa), r3t(:,:,Naa), r3u(:,:,Naa), r3v(:,:,Naa), r3f(:,:) )   ! update ssh/h_0 ratio at t,u,v,f pts  
    225229      ENDIF 
    226230                            CALL dyn_zdf    ( kstp, Nbb, Nnn, Nrhs, uu, vv, Naa  )  ! vertical diffusion 
     
    305309!!    place. 
    306310!! 
    307                          CALL mlf_baro_corr ( Nnn, Naa, uu, vv )                    ! barotrope ajustment 
    308                          CALL finalize_sbc  ( kstp, Nbb, Naa, uu, vv, ts )          ! boundary condifions 
    309                          CALL tra_atf_qco   ( kstp, Nbb, Nnn, Naa, ts )             ! time filtering of "now" tracer arrays 
    310                          CALL dyn_atf_qco   ( kstp, Nbb, Nnn, Naa, uu, vv  )        ! time filtering of "now" velocities and scale factors 
    311                          r3t(:,:,Nnn) = r3t_f(:,:) 
     311      IF( ln_dynspg_ts ) CALL mlf_baro_corr (            Nnn, Naa, uu, vv     )   ! barotrope adjustment 
     312                         CALL finalize_lbc  ( kstp, Nbb     , Naa, uu, vv, ts )   ! boundary conditions 
     313                         CALL tra_atf_qco   ( kstp, Nbb, Nnn, Naa        , ts )   ! time filtering of "now" tracer arrays 
     314                         CALL dyn_atf_qco   ( kstp, Nbb, Nnn, Naa, uu, vv     )   ! time filtering of "now" velocities  
     315                         r3t(:,:,Nnn) = r3t_f(:,:)                                ! update now ssh/h_0 with time filtered values 
    312316                         r3u(:,:,Nnn) = r3u_f(:,:) 
    313317                         r3v(:,:,Nnn) = r3v_f(:,:) 
     
    347351      ! Control 
    348352      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    349                          CALL stp_ctl      ( kstp, Nbb, Nnn, indic ) 
     353                         CALL stp_ctl      ( kstp, Nnn ) 
    350354 
    351355      IF( kstp == nit000 ) THEN                          ! 1st time step only 
     
    364368      IF( kstp == nitend .OR. indic < 0 ) THEN 
    365369                      CALL iom_context_finalize(      cxios_context          ) ! needed for XIOS+AGRIF 
    366          IF(lrxios) CALL iom_context_finalize(      crxios_context          ) 
     370         IF( lrxios ) CALL iom_context_finalize(      crxios_context          ) 
    367371         IF( ln_crs ) CALL iom_context_finalize( trim(cxios_context)//"_crs" ) ! 
    368372      ENDIF 
     
    380384 
    381385 
    382    SUBROUTINE mlf_baro_corr (Kmm, Kaa, puu, pvv) 
     386   SUBROUTINE mlf_baro_corr( Kmm, Kaa, puu, pvv ) 
    383387      !!---------------------------------------------------------------------- 
    384388      !!                  ***  ROUTINE mlf_baro_corr  *** 
     
    389393      !!             estimate (ln_dynspg_ts=T) 
    390394      !! 
    391       !! ** Action :   puu(Kmm),pvv(Kmm),puu(Kaa),pvv(Kaa)   now and after horizontal velocity 
    392       !!---------------------------------------------------------------------- 
    393       INTEGER                             , INTENT(in   ) :: Kmm, Kaa    ! before and after time level indices 
    394       REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: puu, pvv         ! velocities 
    395       ! 
    396       REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   zue, zve 
     395      !! ** Action :   puu(Kmm),pvv(Kmm)   updated now horizontal velocity (ln_bt_fw=F) 
     396      !!               puu(Kaa),pvv(Kaa)   after horizontal velocity 
     397      !!---------------------------------------------------------------------- 
     398      USE dynspg_ts, ONLY : un_adv, vn_adv   ! updated Kmm barotropic transport  
     399      !! 
     400      INTEGER                             , INTENT(in   ) ::   Kmm, Kaa   ! before and after time level indices 
     401      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) ::   puu, pvv   ! velocities 
    397402      ! 
    398403      INTEGER  ::   jk   ! dummy loop indices 
    399       !!---------------------------------------------------------------------- 
    400  
    401       IF ( ln_dynspg_ts ) THEN 
    402          ALLOCATE( zue(jpi,jpj)     , zve(jpi,jpj)     ) 
    403          ! Ensure below that barotropic velocities match time splitting estimate 
    404          ! Compute actual transport and replace it with ts estimate at "after" time step 
    405          zue(:,:) = e3u(:,:,1,Kaa) * puu(:,:,1,Kaa) * umask(:,:,1) 
    406          zve(:,:) = e3v(:,:,1,Kaa) * pvv(:,:,1,Kaa) * vmask(:,:,1) 
    407          DO jk = 2, jpkm1 
    408             zue(:,:) = zue(:,:) + e3u(:,:,jk,Kaa) * puu(:,:,jk,Kaa) * umask(:,:,jk) 
    409             zve(:,:) = zve(:,:) + e3v(:,:,jk,Kaa) * pvv(:,:,jk,Kaa) * vmask(:,:,jk) 
     404      REAL(wp), DIMENSION(jpi,jpj) ::   zue, zve 
     405      !!---------------------------------------------------------------------- 
     406 
     407      ! Ensure below that barotropic velocities match time splitting estimate 
     408      ! Compute actual transport and replace it with ts estimate at "after" time step 
     409      zue(:,:) = e3u(:,:,1,Kaa) * puu(:,:,1,Kaa) * umask(:,:,1) 
     410      zve(:,:) = e3v(:,:,1,Kaa) * pvv(:,:,1,Kaa) * vmask(:,:,1) 
     411      DO jk = 2, jpkm1 
     412         zue(:,:) = zue(:,:) + e3u(:,:,jk,Kaa) * puu(:,:,jk,Kaa) * umask(:,:,jk) 
     413         zve(:,:) = zve(:,:) + e3v(:,:,jk,Kaa) * pvv(:,:,jk,Kaa) * vmask(:,:,jk) 
     414      END DO 
     415      DO jk = 1, jpkm1 
     416         puu(:,:,jk,Kaa) = ( puu(:,:,jk,Kaa) - zue(:,:) * r1_hu(:,:,Kaa) + uu_b(:,:,Kaa) ) * umask(:,:,jk) 
     417         pvv(:,:,jk,Kaa) = ( pvv(:,:,jk,Kaa) - zve(:,:) * r1_hv(:,:,Kaa) + vv_b(:,:,Kaa) ) * vmask(:,:,jk) 
     418      END DO 
     419      ! 
     420      IF( .NOT.ln_bt_fw ) THEN 
     421         ! Remove advective velocity from "now velocities" 
     422         ! prior to asselin filtering 
     423         ! In the forward case, this is done below after asselin filtering 
     424         ! so that asselin contribution is removed at the same time 
     425         DO jk = 1, jpkm1 
     426            puu(:,:,jk,Kmm) = ( puu(:,:,jk,Kmm) - un_adv(:,:)*r1_hu(:,:,Kmm) + uu_b(:,:,Kmm) )*umask(:,:,jk) 
     427            pvv(:,:,jk,Kmm) = ( pvv(:,:,jk,Kmm) - vn_adv(:,:)*r1_hv(:,:,Kmm) + vv_b(:,:,Kmm) )*vmask(:,:,jk) 
    410428         END DO 
    411          DO jk = 1, jpkm1 
    412             puu(:,:,jk,Kaa) = ( puu(:,:,jk,Kaa) - zue(:,:) * r1_hu(:,:,Kaa) + uu_b(:,:,Kaa) ) * umask(:,:,jk) 
    413             pvv(:,:,jk,Kaa) = ( pvv(:,:,jk,Kaa) - zve(:,:) * r1_hv(:,:,Kaa) + vv_b(:,:,Kaa) ) * vmask(:,:,jk) 
    414          END DO 
    415          ! 
    416          IF( .NOT.ln_bt_fw ) THEN 
    417             ! Remove advective velocity from "now velocities" 
    418             ! prior to asselin filtering 
    419             ! In the forward case, this is done below after asselin filtering 
    420             ! so that asselin contribution is removed at the same time 
    421             DO jk = 1, jpkm1 
    422                puu(:,:,jk,Kmm) = ( puu(:,:,jk,Kmm) - un_adv(:,:)*r1_hu(:,:,Kmm) + uu_b(:,:,Kmm) )*umask(:,:,jk) 
    423                pvv(:,:,jk,Kmm) = ( pvv(:,:,jk,Kmm) - vn_adv(:,:)*r1_hv(:,:,Kmm) + vv_b(:,:,Kmm) )*vmask(:,:,jk) 
    424             END DO 
    425          ENDIF 
    426          ! 
    427          DEALLOCATE( zue, zve ) 
    428          ! 
    429429      ENDIF 
    430430      ! 
     
    432432 
    433433 
    434    SUBROUTINE finalize_sbc (kt, Kbb, Kaa, puu, pvv, pts) 
    435       !!---------------------------------------------------------------------- 
    436       !!                  ***  ROUTINE finalize_sbc  *** 
     434   SUBROUTINE finalize_lbc( kt, Kbb, Kaa, puu, pvv, pts ) 
     435      !!---------------------------------------------------------------------- 
     436      !!                  ***  ROUTINE finalize_lbc  *** 
    437437      !! 
    438438      !! ** Purpose :   Apply the boundary condition on the after velocity 
     
    445445      !! ** Action :   puu(Kaa),pvv(Kaa)   after horizontal velocity and tracers 
    446446      !!---------------------------------------------------------------------- 
    447       INTEGER                             , INTENT(in   ) :: kt               ! ocean time-step index 
    448       INTEGER                             , INTENT(in   ) :: Kbb, Kaa    ! before and after time level indices 
    449       REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: puu, pvv         ! velocities to be time filtered 
    450       REAL(wp), DIMENSION(jpi,jpj,jpk,jpts,jpt), INTENT(inout) :: pts            ! active tracers 
     447      INTEGER                                  , INTENT(in   ) ::   kt         ! ocean time-step index 
     448      INTEGER                                  , INTENT(in   ) ::   Kbb, Kaa   ! before and after time level indices 
     449      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt)     , INTENT(inout) ::   puu, pvv   ! velocities to be time filtered 
     450      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts,jpt), INTENT(inout) ::   pts        ! active tracers 
     451      !!---------------------------------------------------------------------- 
    451452      ! 
    452453      ! Update after tracer and velocity on domain lateral boundaries 
    453454      ! 
    454 #if defined key_agrif 
    455             CALL Agrif_tra                     ! AGRIF zoom boundaries 
    456             CALL Agrif_dyn( kt )               !* AGRIF zoom boundaries 
    457 #endif 
     455# if defined key_agrif 
     456            CALL Agrif_tra                     !* AGRIF zoom boundaries 
     457            CALL Agrif_dyn( kt ) 
     458# endif 
    458459      !                                        ! local domain boundaries  (T-point, unchanged sign) 
    459       CALL lbc_lnk_multi( 'finalize_sbc', puu(:,:,:,       Kaa), 'U', -1., pvv(:,:,:       ,Kaa), 'V', -1.   & 
    460                        &                , pts(:,:,:,jp_tem,Kaa), 'T',  1., pts(:,:,:,jp_sal,Kaa), 'T',  1. )     !* local domain boundaries 
     460      CALL lbc_lnk_multi( 'finalize_lbc', puu(:,:,:,       Kaa), 'U', -1., pvv(:,:,:       ,Kaa), 'V', -1.   & 
     461                       &                , pts(:,:,:,jp_tem,Kaa), 'T',  1., pts(:,:,:,jp_sal,Kaa), 'T',  1. ) 
    461462      ! 
    462463      !                                        !* BDY open boundaries 
     
    467468      ENDIF 
    468469      ! 
    469    END SUBROUTINE finalize_sbc 
    470 #endif 
    471    ! 
     470   END SUBROUTINE finalize_lbc 
     471 
     472#else 
     473   !!---------------------------------------------------------------------- 
     474   !!   default option             EMPTY MODULE           qco not activated 
     475   !!---------------------------------------------------------------------- 
     476#endif 
     477    
    472478   !!====================================================================== 
    473479END MODULE stepMLF 
  • NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/src/OFF/dtadyn.F90

    r13295 r13427  
    2323   USE c1d             ! 1D configuration: lk_c1d 
    2424   USE dom_oce         ! ocean domain: variables 
    25 #if ! defined key_qco  
    26    USE domvvl          ! variable volume 
     25#if defined key_qco  
     26   USE domqco          ! variable volume 
    2727#else 
    28    USE domqco 
     28   USE domvvl 
    2929#endif 
    3030   USE zdf_oce         ! ocean vertical physics: variables 
     
    9797   !! * Substitutions 
    9898#  include "do_loop_substitute.h90" 
     99#  include "domzgr_substitute.h90" 
     100    
    99101   !!---------------------------------------------------------------------- 
    100102   !! NEMO/OFF 4.0 , NEMO Consortium (2018) 
     
    385387        gdepw(:,:,:,Kbb) = gdepw(:,:,:,Kmm) 
    386388        ! 
    387       ENDIF 
    388389#endif 
     390      ENDIF 
    389391      ! 
    390392      IF( ln_dynrnf .AND. ln_dynrnf_depth ) THEN       ! read depht over which runoffs are distributed 
  • NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/src/OFF/nemogcm.F90

    r13286 r13427  
    6464   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 
    6565   USE lbcnfd  , ONLY : isendto, nsndto   ! Setup of north fold exchanges 
    66    USE step, ONLY : Nbb, Nnn, Naa, Nrhs   ! time level indices 
     66#if defined key_qco 
     67   USE stepMLF , ONLY : Nbb, Nnn, Naa, Nrhs   ! time level indices 
     68#else 
     69   USE step    , ONLY : Nbb, Nnn, Naa, Nrhs   ! time level indices 
     70#endif 
    6771   USE halo_mng 
    6872 
     
    126130                                CALL dta_dyn_atf( istp, Nbb, Nnn, Naa )       ! time filter of sea  surface height and vertical scale factors 
    127131# if defined key_qco 
    128                                 CALL dom_qco_r3c( ssh(:,:,Kmm), r3t_f, r3u_f, r3v_f ) 
     132                                CALL dom_qco_r3c( ssh(:,:,Nnn), r3t_f, r3u_f, r3v_f ) 
    129133# endif 
    130134         ENDIF 
    131135                                CALL trc_stp    ( istp, Nbb, Nnn, Nrhs, Naa ) ! time-stepping 
    132136# if defined key_qco 
    133                                 !r3t(:,:,Kmm) = r3t_f(:,:)                     ! update ssh to h0 ratio 
    134                                 !r3u(:,:,Kmm) = r3u_f(:,:) 
    135                                 !r3v(:,:,Kmm) = r3v_f(:,:) 
     137                                !r3t(:,:,Nnn) = r3t_f(:,:)                     ! update ssh to h0 ratio 
     138                                !r3u(:,:,Nnn) = r3u_f(:,:) 
     139                                !r3v(:,:,Nnn) = r3v_f(:,:) 
    136140# endif 
    137141#endif 
     
    143147         ! 
    144148#if ! defined key_qco 
    145 #if ! defined key_sed_off 
     149# if ! defined key_sed_off 
    146150         IF( .NOT.ln_linssh )   CALL dta_dyn_sf_interp( istp, Nnn )  ! calculate now grid parameters 
    147 #endif 
     151# endif 
    148152#endif          
    149153                                CALL stp_ctl    ( istp )             ! Time loop: control and print 
Note: See TracChangeset for help on using the changeset viewer.