Changeset 13604


Ignore:
Timestamp:
2020-10-14T18:07:31+02:00 (3 months ago)
Author:
techene
Message:

#2385 cleaning, relevant features moved to OCE, bug corrected in dynvor

Location:
NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/src/SWE
Files:
1 added
14 deleted
3 edited
1 copied

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/src/SWE/dynvor.F90

    r13295 r13604  
    447447            IF( ln_dynvor_msk ) THEN                     ! mask the relative vorticity 
    448448               DO_2D( 1, 0, 1, 0 ) 
    449                   zwz(ji,jj) = ff_f(ji,jj) * fmask(ji,jj,jk) 
     449                  zwz(ji,jj) = zwz(ji,jj) * fmask(ji,jj,jk) 
    450450               END_2D 
    451451            ENDIF 
  • NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/src/SWE/stpRK3.F90

    r13295 r13604  
    11MODULE stpRK3 
    22   !!====================================================================== 
    3    !!                       ***  MODULE step  *** 
     3   !!                       ***  MODULE stpRK3  *** 
    44   !! Time-stepping   : manager of the shallow water equation time stepping 
     5   !!                   3rd order Runge-Kutta scheme 
    56   !!====================================================================== 
    67   !! History :  NEMO !  2020-03  (A. Nasser, G. Madec)  Original code from  4.0.2 
    7    !!---------------------------------------------------------------------- 
    8  
    9    !!---------------------------------------------------------------------- 
    10    !!   stpRK3             : Shallow Water time-stepping 
    11    !!---------------------------------------------------------------------- 
    12    USE step_oce         ! time stepping definition modules 
    13    USE phycst           ! physical constants 
    14    USE usrdef_nam 
     8   !!             -   !  2020-10  (S. Techene, G. Madec)  cleanning 
     9   !!---------------------------------------------------------------------- 
     10 
     11   !!---------------------------------------------------------------------- 
     12   !!   stp_RK3       : RK3 Shallow Water Eq. time-stepping 
     13   !!---------------------------------------------------------------------- 
     14   USE stp_oce        ! modules used in nemo_init and stp_RK3 
    1515   ! 
    16    USE iom              ! xIOs server  
    17    USE domqco 
     16   USE domqco         ! quasi-eulerian coordinate 
     17   USE phycst         ! physical constants 
     18   USE usrdef_nam     ! user defined namelist parameters 
    1819 
    1920   IMPLICIT NONE 
     
    2122 
    2223   PUBLIC   stp_RK3   ! called by nemogcm.F90 
    23     
    24    !!---------------------------------------------------------------------- 
    25    !! time level indices 
    26    !!---------------------------------------------------------------------- 
    27    INTEGER, PUBLIC ::   Nbb, Nnn, Naa, Nrhs   !! used by nemo_init 
     24 
     25   !                                          !**  time level indices  **! 
     26   INTEGER, PUBLIC ::   Nbb, Nnn, Naa, Nrhs   !: used by nemo_init 
    2827       
    2928   !! * Substitutions 
     
    3736CONTAINS 
    3837 
    39 #if defined key_agrif 
    40    RECURSIVE SUBROUTINE stp_RK3( ) 
    41       INTEGER             ::   kstp   ! ocean time-step index 
    42 #else 
    4338   SUBROUTINE stp_RK3( kstp ) 
    44       INTEGER, INTENT(in) ::   kstp   ! ocean time-step index 
    45 #endif 
    4639      !!---------------------------------------------------------------------- 
    4740      !!                     ***  ROUTINE stp_RK3  *** 
    4841      !! 
    49       !! ** Purpose : - Time stepping of shallow water (SHW) (momentum and ssh eqs.) 
     42      !! ** Purpose : - RK3 Time stepping scheme for shallow water Eq. 
    5043      !! 
    51       !! ** Method  : -1- Update forcings 
    52       !!              -2- Update the ssh at Naa 
    53       !!              -3- Compute the momentum trends (Nrhs) 
    54       !!              -4- Update the horizontal velocity 
    55       !!              -5- Apply Asselin time filter to uu,vv,ssh 
    56       !!              -6- Outputs and diagnostics 
     44      !! ** Method  : 3rd order time stepping scheme which has 3 stages 
     45      !!       * Update calendar and forcings 
     46      !!       stage 1   : n ==> n+1/3 using variables at n 
     47      !!                 - Compute the rhs of momentum 
     48      !!                 - Time step ssh at Naa (n+1/3) 
     49      !!                 - Time step u,v at Naa (n+1/3) 
     50      !!                 - Swap time indices 
     51      !!       stage 2   : n ==> n+1/2 using variables at n and n+1/3 
     52      !!                 - Compute the rhs of momentum 
     53      !!                 - Time step ssh at Naa (n+1/2) 
     54      !!                 - Time step u,v at Naa (n+1/2) 
     55      !!                 - Swap time indices 
     56      !!       stage 3   : n ==> n+1 using variables at n and n+1/2 
     57      !!                 - Compute the rhs of momentum 
     58      !!                 - Time step ssh at Naa (n+1) 
     59      !!                 - Time step u,v at Naa (n+1) 
     60      !!                 - Swap time indices 
     61      !!       * Outputs and diagnostics 
     62      !! 
     63      !!       NB: in stages 1 and 2 lateral mixing and forcing are not taken 
     64      !!           into account in the momentum RHS execpt if key_RK3all is used 
    5765      !!---------------------------------------------------------------------- 
     66      INTEGER, INTENT(in   ) ::   kstp   ! ocean time-step index 
     67      ! 
    5868      INTEGER ::   ji, jj, jk   ! dummy loop indice 
    5969      INTEGER ::   indic        ! error indicator if < 0 
    60 !!gm kcall can be removed, I guess 
    61       INTEGER ::   kcall        ! optional integer argument (dom_vvl_sf_nxt) 
    6270      REAL(wp)::   z1_2rho0,  z5_6,  z3_4  ! local scalars 
    63        
    64       REAL(wp) ::   zue3a, zue3n, zue3b    ! local scalars 
    65       REAL(wp) ::   zve3a, zve3n, zve3b    !   -      - 
    66       REAL(wp) ::   ze3t_tf, ze3u_tf, ze3v_tf, zua, zva 
     71      REAL(wp)::   zue3a, zue3b, zua    ! local scalars 
     72      REAL(wp)::   zve3a, zve3b, zva    !   -      - 
    6773      !! --------------------------------------------------------------------- 
    68 #if defined key_agrif 
    69       kstp = nit000 + Agrif_Nb_Step() 
    70       Kbb_a = Nbb; Kmm_a = Nnn; Krhs_a = Nrhs   ! agrif_oce module copies of time level indices 
    71       IF( lk_agrif_debug ) THEN 
    72          IF( Agrif_Root() .and. lwp)   WRITE(*,*) '---' 
    73          IF(lwp)   WRITE(*,*) 'Grid Number', Agrif_Fixed(),' time step ', kstp, 'int tstep', Agrif_NbStepint() 
    74       ENDIF 
    75       IF( kstp == nit000 + 1 )   lk_agrif_fstep = .FALSE. 
    76 # if defined key_iomput 
    77       IF( Agrif_Nbstepint() == 0 )   CALL iom_swap( cxios_context ) 
    78 # endif 
    79 #endif 
    8074      ! 
    8175      IF( ln_timing )   CALL timing_start('stp_RK3') 
     
    9387                             indic = 0                ! reset to no error condition 
    9488                              
    95       IF( kstp == nit000 ) THEN                       ! initialize IOM context (must be done after nemo_init for AGRIF+XIOS+OASIS) 
    96                              CALL iom_init( cxios_context, ld_closedef=.FALSE. )   ! for model grid (including passible AGRIF zoom) 
    97          IF( lk_diamlr   )   CALL dia_mlr_iom_init    ! with additional setup for multiple-linear-regression analysis 
     89      IF( kstp == nit000 ) THEN                       ! initialize IOM context 
     90                             CALL iom_init( cxios_context, ld_closedef=.FALSE. )   ! for model grid (including possible AGRIF zoom) 
    9891                             CALL iom_init_closedef 
    99          IF( ln_crs      )   CALL iom_init( TRIM(cxios_context)//"_crs" )  ! for coarse grid 
    10092      ENDIF 
    10193      IF( kstp /= nit000 )   CALL day( kstp )         ! Calendar (day was already called at nit000 in day_init) 
    10294                             CALL iom_setkt( kstp - nit000 + 1,      cxios_context          )   ! tell IOM we are at time step kstp 
    103       IF( ln_crs         )   CALL iom_setkt( kstp - nit000 + 1, TRIM(cxios_context)//"_crs" )   ! tell IOM we are at time step kstp 
    104  
    105       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    106       ! Update external forcing (tides, open boundaries, ice shelf interaction and surface boundary condition (including sea-ice) 
    107       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    108       IF( ln_tide    )   CALL tide_update( kstp )                     ! update tide potential 
    109       IF( ln_apr_dyn )   CALL sbc_apr ( kstp )                        ! atmospheric pressure (NB: call before bdy_dta which needs ssh_ib)  
    110       IF( ln_bdy     )   CALL bdy_dta ( kstp, Nnn )                   ! update dynamic & tracer data at open boundaries 
    111                          CALL sbc     ( kstp, Nbb, Nnn )              ! Sea Boundary Condition 
    112  
    113       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    114       ! Ocean physics update 
    115       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    116       !  LATERAL  PHYSICS 
    117       !                                                                        ! eddy diffusivity coeff. 
     95 
     96      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
     97      ! Update external forcing   (SWE: surface boundary condition only) 
     98      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
     99 
     100                             CALL sbc     ( kstp, Nbb, Nnn )                   ! Sea Boundary Condition 
     101 
     102      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
     103      ! Ocean physics update   (SWE: eddy viscosity only) 
     104      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
     105 
    118106      IF( l_ldfdyn_time )   CALL ldf_dyn( kstp, Nbb )                          ! eddy viscosity coeff.  
    119  
    120107 
    121108      !====================================================================== 
     
    127114       
    128115      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    129       !  RK3 1st stage Ocean dynamics : hdiv, ssh, e3, u, v, w 
     116      !  RK3 1st stage Ocean dynamics : u, v, ssh 
    130117      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    131118      rDt   = rn_Dt / 3._wp   
    132119      r1_Dt = 1._wp / rDt 
    133        
    134                             CALL ssh_nxt       ( kstp, Nbb, Nbb, ssh, Naa )    ! after ssh (includes call to div_hor) 
    135  
    136                          uu(:,:,:,Nrhs) = 0._wp            ! set dynamics trends to zero 
    137                          vv(:,:,:,Nrhs) = 0._wp 
    138  
    139                             CALL dyn_adv( kstp, Nbb, Nbb      , uu, vv, Nrhs )  ! advection (VF or FF)   ==> RHS 
    140   
    141                             CALL dyn_vor( kstp,      Nbb      , uu, vv, Nrhs )  ! vorticity              ==> RHS 
     120      ! 
     121      !                                 !==  RHS of the momentum Eq.  ==! 
     122      ! 
     123      uu(:,:,:,Nrhs) = 0._wp                        ! set dynamics trends to zero 
     124      vv(:,:,:,Nrhs) = 0._wp 
     125 
     126      CALL dyn_adv( kstp, Nbb, Nbb, uu, vv, Nrhs )  ! advection (VF or FF) ==> RHS 
     127      CALL dyn_vor( kstp,      Nbb, uu, vv, Nrhs )  ! vorticity            ==> RHS 
    142128#if defined key_RK3all  
    143                             CALL dyn_ldf( kstp, Nbb, Nbb      , uu, vv, Nrhs )  ! lateral mixing 
    144 #endif 
    145       ! 
    146 !!an - calcul du gradient de pression horizontal (explicit) 
    147       DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 
     129      CALL dyn_ldf( kstp, Nbb, Nbb, uu, vv, Nrhs )  ! lateral mixing 
     130#endif 
     131      ! 
     132      DO_3D( 0,0, 0,0, 1,jpkm1 ) 
     133         !                                          ! horizontal pressure gradient 
    148134         uu(ji,jj,jk,Nrhs) = uu(ji,jj,jk,Nrhs) - grav * ( ssh(ji+1,jj,Nbb) - ssh(ji,jj,Nbb) ) * r1_e1u(ji,jj) 
    149135         vv(ji,jj,jk,Nrhs) = vv(ji,jj,jk,Nrhs) - grav * ( ssh(ji,jj+1,Nbb) - ssh(ji,jj,Nbb) ) * r1_e2v(ji,jj) 
     
    151137      ! 
    152138#if defined key_RK3all  
    153       ! add wind stress forcing and layer linear friction to the RHS  
     139      !                                             ! wind stress and layer friction 
    154140      z5_6 = 5._wp/6._wp 
    155141      DO_3D( 0, 0, 0, 0,1,jpkm1) 
     
    160146      END_3D 
    161147#endif 
    162 !!an 
    163                             CALL dom_qco_r3c   ( ssh(:,:,Naa), r3t(:,:,Naa), r3u(:,:,Naa), r3v(:,:,Naa), r3f(:,:) )   ! "after" ssh./h._0 ratio explicit 
    164       IF( ln_dynadv_vec ) THEN      ! vector invariant form : applied on velocity 
     148!!!!st why not ? 
     149!!      z5_6 = 5._wp/6._wp 
     150!!      DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 
     151!!         !                                                              ! horizontal pressure gradient 
     152!!         zrhs_u =        - grav * ( ssh(ji+1,jj,Nnn) - ssh(ji,jj,Nnn) ) * r1_e1u(ji,jj) 
     153!!         zrhs_v =        - grav * ( ssh(ji,jj+1,Nnn) - ssh(ji,jj,Nnn) ) * r1_e2v(ji,jj) 
     154!!         !                                                              ! wind stress and layer friction 
     155!!#if defined key_RK3all 
     156!!         zrhs_u = zrhs_u + r1_rho0 * ( z5_6*utau_b(ji,jj) + (1._wp - z5_6)*utau(ji,jj) ) / e3u(ji,jj,jk,Nbb)   & 
     157!!            &                                  - rn_rfr * uu(ji,jj,jk,Nbb) 
     158!!         zrhs_v = zrhs_v + r1_rho0 * ( z5_6*vtau_b(ji,jj) + (1._wp - z5_6)*vtau(ji,jj) ) / e3v(ji,jj,jk,Nbb)   & 
     159!!            &                                  - rn_rfr * vv(ji,jj,jk,Nbb) 
     160!!         !                                                              ! ==> RHS 
     161!!#endif 
     162 
     163!!         uu(ji,jj,jk,Nrhs) = uu(ji,jj,jk,Nrhs) + zrhs_u 
     164!!         vv(ji,jj,jk,Nrhs) = vv(ji,jj,jk,Nrhs) + zrhs_v 
     165!!      END_3D 
     166!!!!st end 
     167      ! 
     168      !                                 !==  Time stepping of ssh Eq.  ==!   (and update r3_Naa) 
     169      ! 
     170      CALL ssh_nxt( kstp, Nbb, Nbb, ssh, Naa )      ! after ssh 
     171      !                                                                 ! after ssh/h_0 ratio explicit 
     172      CALL dom_qco_r3c( ssh(:,:,Naa), r3t(:,:,Naa), r3u(:,:,Naa), r3v(:,:,Naa), r3f(:,:) ) 
     173      ! 
     174      !                                 !==  Time stepping of momentum Eq.  ==! 
     175      ! 
     176      IF( ln_dynadv_vec ) THEN                      ! vector invariant form : applied on velocity 
    165177         DO_3D( 0, 0, 0, 0,1,jpkm1) 
    166178            uu(ji,jj,jk,Naa) = uu(ji,jj,jk,Nbb) + rDt * uu(ji,jj,jk,Nrhs) * umask(ji,jj,jk) 
     
    168180         END_3D           
    169181      ELSE 
    170          DO_3D( 0, 0, 0, 0,1,jpkm1)       ! flux form : applied on thickness weighted velocity 
    171             uu(ji,jj,jk,Naa) = (         uu(ji,jj,jk,Nbb )*e3u(ji,jj,jk,Nbb)                              & 
    172                &                 + rDt * uu(ji,jj,jk,Nrhs)*e3t(ji,jj,jk,Nbb) * umask(ji,jj,jk)        )   & 
    173                &               /                           e3t(ji,jj,jk,Naa) 
    174             vv(ji,jj,jk,Naa) = (         vv(ji,jj,jk,Nbb )*e3v(ji,jj,jk,Nbb)                              & 
    175                &                 + rDt * vv(ji,jj,jk,Nrhs)*e3t(ji,jj,jk,Nbb) * vmask(ji,jj,jk)        )   & 
    176                &               /                           e3t(ji,jj,jk,Naa) 
     182         DO_3D( 0, 0, 0, 0,1,jpkm1)                 ! flux form : applied on thickness weighted velocity 
     183            zue3b = e3u(ji,jj,jk,Nbb) * uu(ji,jj,jk,Nbb) 
     184            zve3b = e3v(ji,jj,jk,Nbb) * vv(ji,jj,jk,Nbb) 
     185            zue3a = zue3b + rDt * e3u(ji,jj,jk,Nbb) * uu(ji,jj,jk,Nrhs) * umask(ji,jj,jk) 
     186            zve3a = zve3b + rDt * e3v(ji,jj,jk,Nbb) * vv(ji,jj,jk,Nrhs) * vmask(ji,jj,jk) 
     187            ! 
     188            uu(ji,jj,jk,Naa) = zue3a / e3u(ji,jj,jk,Naa)     
     189            vv(ji,jj,jk,Naa) = zve3a / e3v(ji,jj,jk,Naa) 
    177190         END_3D 
    178191      ENDIF 
    179       ! Swap time levels 
     192      ! 
     193      !                                 !==  Swap time levels  ==! 
     194      ! 
    180195      Nrhs= Nnn 
    181196      Nnn = Naa 
    182197      Naa = Nrhs 
    183198 
    184        
    185199      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    186200      !  RK3 2nd stage Ocean dynamics : hdiv, ssh, e3, u, v, w 
     
    188202      rDt   = rn_Dt / 2._wp   
    189203      r1_Dt = 1._wp / rDt 
    190        
    191                             CALL ssh_nxt       ( kstp, Nbb, Nnn, ssh, Naa )    ! after ssh (includes call to div_hor) 
    192  
    193                          uu(:,:,:,Nrhs) = 0._wp            ! set dynamics trends to zero 
    194                          vv(:,:,:,Nrhs) = 0._wp 
    195 !!st TBC for dyn_adv 
    196                             CALL dyn_adv( kstp, Nbb, Nnn      , uu, vv, Nrhs )  ! advection (VF or FF)   ==> RHS  
    197   
    198                             CALL dyn_vor( kstp,      Nnn      , uu, vv, Nrhs )  ! vorticity              ==> RHS 
     204      ! 
     205      !                                 !==  RHS of the momentum Eq.  ==! 
     206      ! 
     207      uu(:,:,:,Nrhs) = 0._wp                        ! set dynamics trends to zero 
     208      vv(:,:,:,Nrhs) = 0._wp 
     209      CALL dyn_adv( kstp, Nbb, Nnn, uu, vv, Nrhs )  ! advection (VF or FF) ==> RHS  
     210      CALL dyn_vor( kstp,      Nnn, uu, vv, Nrhs )  ! vorticity            ==> RHS 
    199211#if defined key_RK3all   
    200                             CALL dyn_ldf( kstp, Nbb, Nbb      , uu, vv, Nrhs )  ! lateral mixing 
    201 #endif 
    202                              
    203       ! 
    204 !!an - calcul du gradient de pression horizontal (explicit) 
     212      CALL dyn_ldf( kstp, Nbb, Nbb, uu, vv, Nrhs )  ! lateral mixing 
     213#endif 
     214      ! 
    205215      DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 
     216         !                                          ! horizontal pressure gradient 
    206217         uu(ji,jj,jk,Nrhs) = uu(ji,jj,jk,Nrhs) - grav * ( ssh(ji+1,jj,Nnn) - ssh(ji,jj,Nnn) ) * r1_e1u(ji,jj) 
    207218         vv(ji,jj,jk,Nrhs) = vv(ji,jj,jk,Nrhs) - grav * ( ssh(ji,jj+1,Nnn) - ssh(ji,jj,Nnn) ) * r1_e2v(ji,jj) 
    208219      END_3D 
    209220      ! 
    210       ! add wind stress forcing and layer linear friction to the RHS  
    211221#if defined key_RK3all 
     222      !                                             ! wind stress and layer friction 
    212223      z3_4 = 3._wp/4._wp 
    213224      DO_3D( 0, 0, 0, 0,1,jpkm1) 
     
    218229      END_3D 
    219230#endif 
    220 !!an 
    221                            CALL dom_qco_r3c   ( ssh(:,:,Naa), r3t(:,:,Naa), r3u(:,:,Naa), r3v(:,:,Naa), r3f(:,:) )   ! "after" ssh./h._0 ratio explicit 
    222       IF( ln_dynadv_vec ) THEN      ! vector invariant form : applied on velocity 
     231      ! 
     232      !                                 !==  Time stepping of ssh Eq.  ==!   (and update r3_Naa) 
     233      ! 
     234      CALL ssh_nxt( kstp, Nbb, Nnn, ssh, Naa )      ! after ssh 
     235      !                                                                 ! after ssh/h_0 ratio explicit 
     236      CALL dom_qco_r3c( ssh(:,:,Naa), r3t(:,:,Naa), r3u(:,:,Naa), r3v(:,:,Naa), r3f(:,:) ) 
     237      ! 
     238      !                                 !==  Time stepping of momentum Eq.  ==! 
     239      ! 
     240      IF( ln_dynadv_vec ) THEN                      ! vector invariant form : applied on velocity 
    223241         DO_3D( 0, 0, 0, 0,1,jpkm1) 
    224242            uu(ji,jj,jk,Naa) = uu(ji,jj,jk,Nbb) + rDt * uu(ji,jj,jk,Nrhs) * umask(ji,jj,jk) 
     
    226244         END_3D           
    227245      ELSE 
    228          DO_3D( 0, 0, 0, 0,1,jpkm1)       ! flux form : applied on thickness weighted velocity 
    229             uu(ji,jj,jk,Naa) = (         uu(ji,jj,jk,Nbb )*e3u(ji,jj,jk,Nbb)                              & 
    230                &                 + rDt * uu(ji,jj,jk,Nrhs)*e3t(ji,jj,jk,Nnn) * umask(ji,jj,jk)        )   & 
    231                &               /                           e3t(ji,jj,jk,Naa) 
    232             vv(ji,jj,jk,Naa) = (         vv(ji,jj,jk,Nbb )*e3v(ji,jj,jk,Nbb)                              & 
    233                &                 + rDt * vv(ji,jj,jk,Nrhs)*e3t(ji,jj,jk,Nnn) * vmask(ji,jj,jk)        )   & 
    234                &               /                           e3t(ji,jj,jk,Naa) 
     246         DO_3D( 0, 0, 0, 0,1,jpkm1)                 ! flux form : applied on thickness weighted velocity 
     247            zue3b = e3u(ji,jj,jk,Nbb) * uu(ji,jj,jk,Nbb) 
     248            zve3b = e3v(ji,jj,jk,Nbb) * vv(ji,jj,jk,Nbb) 
     249            zue3a = zue3b + rDt * e3u(ji,jj,jk,Nnn) * uu(ji,jj,jk,Nrhs) * umask(ji,jj,jk) 
     250            zve3a = zve3b + rDt * e3v(ji,jj,jk,Nnn) * vv(ji,jj,jk,Nrhs) * vmask(ji,jj,jk) 
     251            ! 
     252            uu(ji,jj,jk,Naa) = zue3a / e3u(ji,jj,jk,Naa)     
     253            vv(ji,jj,jk,Naa) = zve3a / e3v(ji,jj,jk,Naa) 
    235254         END_3D 
    236255      ENDIF 
    237       ! Swap time levels 
     256      ! 
     257      !                                 !==  Swap time levels  ==! 
     258      ! 
    238259      Nrhs= Nnn 
    239260      Nnn = Naa 
     
    245266      rDt   = rn_Dt 
    246267      r1_Dt = 1._wp / rDt 
    247  
    248                             CALL ssh_nxt       ( kstp, Nbb, Nnn, ssh, Naa )    ! after ssh (includes call to div_hor) 
    249        
    250                          uu(:,:,:,Nrhs) = 0._wp            ! set dynamics trends to zero 
    251                          vv(:,:,:,Nrhs) = 0._wp 
    252  
    253       IF( ln_bdy     )      CALL bdy_dyn3d_dmp ( kstp, Nbb,      uu, vv, Nrhs )  ! bdy damping trends 
    254  
    255 #if defined key_agrif 
    256       IF(.NOT. Agrif_Root())  &  
    257                &            CALL Agrif_Sponge_dyn        ! momentum sponge 
    258 #endif 
    259                             CALL dyn_adv( kstp, Nbb, Nnn      , uu, vv, Nrhs )  ! advection (VF or FF)   ==> RHS 
    260   
    261                             CALL dyn_vor( kstp,      Nnn      , uu, vv, Nrhs )  ! vorticity              ==> RHS 
    262   
    263                             CALL dyn_ldf( kstp, Nbb, Nnn      , uu, vv, Nrhs )  ! lateral mixing 
    264  
    265 !!an - calcul du gradient de pression horizontal (explicit) 
     268      ! 
     269      !                                 !==  RHS of the momentum Eq.  ==! 
     270      ! 
     271      uu(:,:,:,Nrhs) = 0._wp                        ! set dynamics trends to zero 
     272      vv(:,:,:,Nrhs) = 0._wp 
     273      ! 
     274      CALL dyn_adv( kstp, Nbb, Nnn, uu, vv, Nrhs )  ! advection (VF or FF) ==> RHS 
     275      CALL dyn_vor( kstp,      Nnn, uu, vv, Nrhs )  ! vorticity            ==> RHS 
     276      CALL dyn_ldf( kstp, Nbb, Nnn, uu, vv, Nrhs )  ! lateral mixing 
     277 
    266278      DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 
     279         !                                          ! horizontal pressure gradient 
    267280         uu(ji,jj,jk,Nrhs) = uu(ji,jj,jk,Nrhs) - grav * ( ssh(ji+1,jj,Nnn) - ssh(ji,jj,Nnn) ) * r1_e1u(ji,jj) 
    268281         vv(ji,jj,jk,Nrhs) = vv(ji,jj,jk,Nrhs) - grav * ( ssh(ji,jj+1,Nnn) - ssh(ji,jj,Nnn) ) * r1_e2v(ji,jj) 
    269282      END_3D 
    270       ! 
    271       ! add wind stress forcing and layer linear friction to the RHS  
     283      !                                             ! wind stress and layer friction 
    272284      z1_2rho0 = 0.5_wp * r1_rho0 
    273285      DO_3D( 0, 0, 0, 0,1,jpkm1) 
     
    276288         vv(ji,jj,jk,Nrhs) = vv(ji,jj,jk,Nrhs) + z1_2rho0 * ( vtau_b(ji,jj) + vtau(ji,jj) ) / e3v(ji,jj,jk,Nnn)   & 
    277289            &                                  - rn_rfr * vv(ji,jj,jk,Nbb) 
    278       END_3D    
    279 !!an          
    280                             CALL dom_qco_r3c   ( ssh(:,:,Naa), r3t(:,:,Naa), r3u(:,:,Naa), r3v(:,:,Naa), r3f(:,:) )   ! "after" ssh./h._0 ratio explicit       
    281       IF( ln_dynadv_vec ) THEN      ! vector invariant form : applied on velocity 
     290      END_3D 
     291      ! 
     292      !                                 !==  Time stepping of ssh Eq.  ==!   (and update r3_Naa) 
     293      ! 
     294      CALL ssh_nxt( kstp, Nbb, Nnn, ssh, Naa )      ! after ssh 
     295      !                                                                 ! after ssh/h_0 ratio explicit 
     296      CALL dom_qco_r3c( ssh(:,:,Naa), r3t(:,:,Naa), r3u(:,:,Naa), r3v(:,:,Naa), r3f(:,:) ) 
     297      ! 
     298      !                                 !==  Time stepping of momentum Eq.  ==! 
     299      ! 
     300      IF( ln_dynadv_vec ) THEN                      ! vector invariant form : applied on velocity 
    282301         DO_3D( 1, 1, 1, 1,1,jpkm1) 
    283302            zua = uu(ji,jj,jk,Nbb) + rDt * uu(ji,jj,jk,Nrhs) * umask(ji,jj,jk) 
     
    291310         END_3D 
    292311         ! 
    293       ELSE                          ! flux form : applied on thickness weighted velocity 
     312      ELSE                                          ! flux form : applied on thickness weighted velocity 
    294313         DO_3D( 1, 1, 1, 1,1,jpkm1) 
    295             zue3n = e3u(ji,jj,jk,Nnn) * uu(ji,jj,jk,Nnn) 
    296             zve3n = e3v(ji,jj,jk,Nnn) * vv(ji,jj,jk,Nnn) 
    297314            zue3b = e3u(ji,jj,jk,Nbb) * uu(ji,jj,jk,Nbb) 
    298315            zve3b = e3v(ji,jj,jk,Nbb) * vv(ji,jj,jk,Nbb) 
    299             !                                                ! LF time stepping 
    300             zue3a = zue3b + rDt * e3t(ji,jj,jk,Nbb) * uu(ji,jj,jk,Nrhs) * umask(ji,jj,jk) 
    301             zve3a = zve3b + rDt * e3t(ji,jj,jk,Nbb) * vv(ji,jj,jk,Nrhs) * vmask(ji,jj,jk) 
     316            zue3a = zue3b + rDt * e3u(ji,jj,jk,Nbb) * uu(ji,jj,jk,Nrhs) * umask(ji,jj,jk) 
     317            zve3a = zve3b + rDt * e3v(ji,jj,jk,Nbb) * vv(ji,jj,jk,Nrhs) * vmask(ji,jj,jk) 
    302318            ! 
    303             uu(ji,jj,jk,Naa) = zue3a / e3t(ji,jj,jk,Naa)     
    304             vv(ji,jj,jk,Naa) = zve3a / e3t(ji,jj,jk,Naa) 
     319            uu(ji,jj,jk,Naa) = zue3a / e3u(ji,jj,jk,Naa)     
     320            vv(ji,jj,jk,Naa) = zve3a / e3v(ji,jj,jk,Naa) 
    305321         END_3D 
    306 !!st je ne comprends pas l'histoire des e3t et du Nbb et pas du Nnn pour le rhs ?         
    307       ENDIF 
    308  
     322      ENDIF 
    309323 
    310324      CALL lbc_lnk_multi( 'stp_RK3', uu(:,:,:,Nnn), 'U', -1., vv(:,:,:,Nnn), 'V', -1.,   &   !* local domain boundaries 
    311          &                       uu(:,:,:,Naa), 'U', -1., vv(:,:,:,Naa), 'V', -1.    )      
    312  
    313 !!an          
    314  
    315  
    316       ! Swap time levels 
     325         &                           uu(:,:,:,Naa), 'U', -1., vv(:,:,:,Naa), 'V', -1.    )      
     326      ! 
     327      !                                 !==  Swap time levels  ==! 
     328      ! 
    317329      Nrhs = Nbb 
    318330      Nbb = Naa 
    319331      Naa = Nrhs 
    320       ! 
    321 !                         CALL dom_vvl_sf_update_st( kstp, Nbb, Nnn, Naa )  ! recompute vertical scale factors 
     332 
    322333      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    323334      ! diagnostics and outputs 
    324335      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    325       IF( ln_floats  )   CALL flo_stp   ( kstp, Nbb, Nnn )      ! drifting Floats 
     336       
    326337      IF( ln_diacfl  )   CALL dia_cfl   ( kstp,      Nnn )      ! Courant number diagnostics 
    327      
    328338                         CALL dia_wri   ( kstp,      Nnn )      ! ocean model: outputs 
    329  
    330339      ! 
    331340      IF( lrst_oce   )   CALL rst_write    ( kstp, Nbb, Nnn )   ! write output ocean restart file 
    332            
    333  
    334 #if defined key_agrif 
    335       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    336       ! AGRIF 
    337       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<       
    338                          Kbb_a = Nbb; Kmm_a = Nnn; Krhs_a = Nrhs      ! agrif_oce module copies of time level indices 
    339                          CALL Agrif_Integrate_ChildGrids( stp_RK3 )       ! allows to finish all the Child Grids before updating 
    340  
    341                          IF( Agrif_NbStepint() == 0 ) THEN 
    342                             CALL Agrif_update_all( )                  ! Update all components 
    343                          ENDIF 
    344 #endif 
    345       IF( ln_diaobs  )   CALL dia_obs      ( kstp, Nnn )      ! obs-minus-model (assimilation) diagnostics (call after dynamics update) 
    346341 
    347342      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    348343      ! Control 
    349344      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    350                          CALL stp_ctl      ( kstp, Nbb, Nnn, indic ) 
    351                 
    352                           
     345                         CALL stp_ctl      ( kstp, Nnn ) 
     346 
    353347      IF( kstp == nit000 ) THEN                          ! 1st time step only 
    354348                                        CALL iom_close( numror )   ! close input  ocean restart file 
     
    360354#if defined key_iomput 
    361355      IF( kstp == nitend .OR. indic < 0 ) THEN  
    362                       CALL iom_context_finalize(      cxios_context          ) ! needed for XIOS+AGRIF 
    363                       IF(lrxios) CALL iom_context_finalize(      crxios_context          ) 
    364       ENDIF 
    365 #endif 
    366       ! 
    367       IF( l_1st_euler ) THEN         ! recover Leap-frog timestep 
    368          rDt = 2._wp * rn_Dt    
    369          r1_Dt = 1._wp / rDt 
    370          l_1st_euler = .FALSE.       
    371       ENDIF 
     356         CALL iom_context_finalize( cxios_context ) 
     357      ENDIF 
     358#endif 
    372359      ! 
    373360      IF( ln_timing )   CALL timing_stop('stp_RK3') 
    374361      ! 
    375362   END SUBROUTINE stp_RK3 
    376    ! 
     363 
    377364   !!====================================================================== 
    378365END MODULE stpRK3 
  • NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/src/SWE/stp_oce.F90

    r13426 r13604  
    1 MODULE step_oce 
     1MODULE stp_oce 
    22   !!====================================================================== 
    3    !!                       ***  MODULE step_oce  *** 
     3   !!                       ***  MODULE stp_oce  *** 
    44   !! Ocean time-stepping : module used in both initialisation phase and time stepping 
     5   !!                                     (i.e. nemo_init and stp or stp_MLF routines) 
    56   !!====================================================================== 
    67   !! History :   3.3  !  2010-08  (C. Ethe)  Original code - reorganisation of the initial phase 
     
    910   USE oce             ! ocean dynamics and tracers variables 
    1011   USE dom_oce         ! ocean space and time domain variables 
    11    USE zdf_oce         ! ocean vertical physics variables 
    12    USE zdfdrg  ,  ONLY : ln_drgimp   ! implicit top/bottom friction 
    1312 
    1413   USE daymod          ! calendar                         (day     routine) 
     
    1918   USE sbccpl          ! surface boundary condition: coupled formulation (call send at end of step) 
    2019   USE sbcapr          ! surface boundary condition: atmospheric pressure 
    21    USE tide_mod, ONLY : ln_tide, tide_update 
    2220   USE sbcwave         ! Wave intialisation 
     21   USE tide_mod        ! tides 
     22 
     23   USE bdy_oce  , ONLY : ln_bdy 
     24   USE bdydta          ! open boundary condition data     (bdy_dta routine) 
     25   USE bdytra          ! bdy cond. for tracers            (bdy_tra routine) 
     26   USE bdydyn3d        ! bdy cond. for baroclinic vel.  (bdy_dyn3d routine) 
    2327 
    2428   USE isf_oce         ! ice shelf boundary condition 
    2529   USE isfstp          ! ice shelf boundary condition     (isf_stp routine) 
     30 
     31   USE sshwzv          ! vertical velocity and ssh        (ssh_nxt routine) 
     32   !                                                      (ssh_swp routine) 
     33   !                                                      (wzv     routine) 
     34   USE domvvl          ! variable vertical scale factors  (dom_vvl_sf_nxt routine) 
     35   !                                                      (dom_vvl_sf_swp routine) 
     36    
     37   USE divhor          ! horizontal divergence            (div_hor routine) 
     38   USE dynadv          ! advection                        (dyn_adv routine) 
     39   USE dynvor          ! vorticity term                   (dyn_vor routine) 
     40   USE dynhpg          ! hydrostatic pressure grad.       (dyn_hpg routine) 
     41   USE dynldf          ! lateral momentum diffusion       (dyn_ldf routine) 
     42   USE dynzdf          ! vertical diffusion               (dyn_zdf routine) 
     43   USE dynspg          ! surface pressure gradient        (dyn_spg routine) 
     44   USE dynatf          ! time-filtering                   (dyn_atf routine) 
    2645 
    2746   USE traqsr          ! solar radiation penetration      (tra_qsr routine) 
     
    3958   USE eosbn2          ! equation of state                (eos_bn2 routine) 
    4059 
    41    USE divhor          ! horizontal divergence            (div_hor routine) 
    42    USE dynadv          ! advection                        (dyn_adv routine) 
    43    USE dynvor          ! vorticity term                   (dyn_vor routine) 
    44    USE dynhpg          ! hydrostatic pressure grad.       (dyn_hpg routine) 
    45    USE dynldf          ! lateral momentum diffusion       (dyn_ldf routine) 
    46    USE dynzdf          ! vertical diffusion               (dyn_zdf routine) 
    47    USE dynspg          ! surface pressure gradient        (dyn_spg routine) 
    48  
    49    USE dynatf          ! time-filtering                   (dyn_atf routine) 
    50  
    5160   USE stopar          ! Stochastic parametrization       (sto_par routine) 
    5261   USE stopts  
    53  
    54    USE bdy_oce  , ONLY : ln_bdy 
    55    USE bdydta          ! open boundary condition data     (bdy_dta routine) 
    56    USE bdytra          ! bdy cond. for tracers            (bdy_tra routine) 
    57    USE bdydyn3d        ! bdy cond. for baroclinic vel.  (bdy_dyn3d routine) 
    58  
    59    USE sshwzv          ! vertical velocity and ssh        (ssh_nxt routine) 
    60    !                                                       (ssh_swp routine) 
    61    !                                                       (wzv     routine) 
    62    USE domvvl          ! variable vertical scale factors  (dom_vvl_sf_nxt routine) 
    63    !                                                       (dom_vvl_sf_swp routine) 
    6462 
    6563   USE ldfslp          ! iso-neutral slopes               (ldf_slp routine) 
     
    6765   USE ldftra          ! lateral eddy diffusive coef.     (ldf_tra routine) 
    6866 
     67   USE zdf_oce         ! ocean vertical physics variables 
    6968   USE zdfphy          ! vertical physics manager      (zdf_phy_init routine) 
    70    USE zdfosm  , ONLY : osm_rst, dyn_osm, tra_osm      ! OSMOSIS routines used in step.F90 
     69   USE zdfdrg   , ONLY : ln_drgimp   ! implicit top/bottom friction 
     70   USE zdfosm   , ONLY : osm_rst, dyn_osm, tra_osm      ! OSMOSIS routines used in step.F90 
    7171 
    7272   USE diu_layers      ! diurnal SST bulk and coolskin routines 
     
    8181   USE diahth          ! thermocline depth                (dia_hth routine) 
    8282   USE diahsb          ! heat, salt and volume budgets    (dia_hsb routine) 
    83    USE diacfl 
    84    USE diaobs          ! Observation operator 
     83   USE diacfl          ! CFL diagnostics                  (dia_cfl routine) 
     84   USE diaobs          ! Observation operator             (dia_obs routine) 
    8585   USE diadetide       ! Weights computation for daily detiding of model diagnostics 
    8686   USE diamlr          ! IOM context management for multiple-linear-regression analysis 
     
    9292   USE asminc          ! assimilation increments      (tra_asm_inc routine) 
    9393   !                                                   (dyn_asm_inc routine) 
    94    USE asmbkg 
     94   USE asmbkg          ! writing out state trajectory 
    9595   USE stpctl          ! time stepping control            (stp_ctl routine) 
    9696   USE restart         ! ocean restart                    (rst_wri routine) 
     
    117117   !! Software governed by the CeCILL license (see ./LICENSE) 
    118118   !!====================================================================== 
    119 END MODULE step_oce 
     119END MODULE stp_oce 
  • NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/src/SWE/stpctl.F90

    r12983 r13604  
    33   !!                       ***  MODULE  stpctl  *** 
    44   !! Ocean run control :  gross check of the ocean time stepping 
     5   !!              *** Shallow Water Equation (SWE) case *** 
     6   !!               ( No test on temperature and salinity ) 
    57   !!====================================================================== 
    6    !! History :  OPA  ! 1991-03  (G. Madec) Original code 
    7    !!            6.0  ! 1992-06  (M. Imbard) 
    8    !!            8.0  ! 1997-06  (A.M. Treguier) 
    9    !!   NEMO     1.0  ! 2002-06  (G. Madec)  F90: Free form and module 
    10    !!            2.0  ! 2009-07  (G. Madec)  Add statistic for time-spliting 
    11    !!            3.7  ! 2016-09  (G. Madec)  Remove solver 
    12    !!            4.0  ! 2017-04  (G. Madec)  regroup global communications 
     8   !! History :  SWE  ! 2020-09  (A. Nasser, S. Techene ) OCE/stpctl adaptated to SWE 
    139   !!---------------------------------------------------------------------- 
    1410 
     
    1915   USE dom_oce         ! ocean space and time domain variables  
    2016   USE c1d             ! 1D vertical configuration 
     17   USE zdf_oce ,  ONLY : ln_zad_Aimp       ! ocean vertical physics variables 
     18   USE wet_dry,   ONLY : ll_wd, ssh_ref    ! reference depth for negative bathy 
     19   ! 
    2120   USE diawri          ! Standard run outputs       (dia_wri_state routine) 
    22    ! 
    2321   USE in_out_manager  ! I/O manager 
    2422   USE lbclnk          ! ocean lateral boundary conditions (or mpp link) 
    2523   USE lib_mpp         ! distributed memory computing 
    26    USE zdf_oce ,  ONLY : ln_zad_Aimp       ! ocean vertical physics variables 
    27    USE wet_dry,   ONLY : ll_wd, ssh_ref    ! reference depth for negative bathy 
    28  
    2924   USE netcdf          ! NetCDF library 
     25 
    3026   IMPLICIT NONE 
    3127   PRIVATE 
     
    3329   PUBLIC stp_ctl           ! routine called by step.F90 
    3430 
    35    INTEGER  ::   idrun, idtime, idssh, idu, ids1, ids2, idt1, idt2, idc1, idw1, istatus 
    36    LOGICAL  ::   lsomeoce 
    37 !!stoops 
     31   INTEGER                ::   nrunid   ! netcdf file id 
     32   INTEGER, DIMENSION(2)  ::   nvarid   ! netcdf variable id 
     33!!SWE   INTEGER, DIMENSION(8)  ::   nvarid   ! netcdf variable id 
     34 
    3835#  include "domzgr_substitute.h90" 
    3936   !!---------------------------------------------------------------------- 
     
    4441CONTAINS 
    4542 
    46    SUBROUTINE stp_ctl( kt, Kbb, Kmm, kindic ) 
     43   SUBROUTINE stp_ctl( kt, Kmm ) 
    4744      !!---------------------------------------------------------------------- 
    4845      !!                    ***  ROUTINE stp_ctl  *** 
     
    5249      !! ** Method  : - Save the time step in numstp 
    5350      !!              - Print it each 50 time steps 
    54       !!              - Stop the run IF problem encountered by setting indic=-3 
    55       !!                Problems checked: |ssh| maximum larger than 10 m 
     51      !!              - Stop the run IF problem encountered by setting nstop > 0 
     52      !!                Problems checked: e3t0+ssh minimum smaller that 0 
    5653      !!                                  |U|   maximum larger than 10 m/s  
    57       !!                                  negative sea surface salinity 
     54      !!                                  ( not for SWE : negative sea surface salinity ) 
    5855      !! 
    5956      !! ** Actions :   "time.step" file = last ocean time-step 
    6057      !!                "run.stat"  file = run statistics 
    61       !!                nstop indicator sheared among all local domain (lk_mpp=T) 
     58      !!                 nstop indicator sheared among all local domain 
    6259      !!---------------------------------------------------------------------- 
    6360      INTEGER, INTENT(in   ) ::   kt       ! ocean time-step index 
    64       INTEGER, INTENT(in   ) ::   Kbb, Kmm      ! ocean time level index 
    65       INTEGER, INTENT(inout) ::   kindic   ! error indicator 
    66       !! 
    67       INTEGER                ::   ji, jj, jk          ! dummy loop indices 
    68       INTEGER, DIMENSION(2)  ::   ih                  ! min/max loc indices 
    69       INTEGER, DIMENSION(3)  ::   iu, is1, is2        ! min/max loc indices 
    70       REAL(wp)               ::   zzz                 ! local real  
    71       REAL(wp), DIMENSION(3) ::   zmax 
    72       LOGICAL                ::   ll_wrtstp, ll_colruns, ll_wrtruns 
    73       CHARACTER(len=20) :: clname 
    74       !!---------------------------------------------------------------------- 
    75       ! 
    76       ll_wrtstp  = ( MOD( kt, sn_cfctl%ptimincr ) == 0 ) .OR. ( kt == nitend ) 
    77       ll_colruns = ll_wrtstp .AND. ( sn_cfctl%l_runstat ) 
    78       ll_wrtruns = ll_colruns .AND. lwm 
    79       IF( kt == nit000 .AND. lwp ) THEN 
    80          WRITE(numout,*) 
    81          WRITE(numout,*) 'stp_ctl : time-stepping control' 
    82          WRITE(numout,*) '~~~~~~~' 
    83          !                                ! open time.step file 
    84          IF( lwm ) CALL ctl_opn( numstp, 'time.step', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp, narea ) 
    85          !                                ! open run.stat file(s) at start whatever 
    86          !                                ! the value of sn_cfctl%ptimincr 
    87          IF( lwm .AND. ( sn_cfctl%l_runstat ) ) THEN 
     61      INTEGER, INTENT(in   ) ::   Kmm      ! ocean time level index 
     62      !! 
     63      INTEGER                         ::   ji                                    ! dummy loop indices 
     64      INTEGER                         ::   idtime, istatus 
     65!!SWE      INTEGER , DIMENSION(9)          ::   iareasum, iareamin, iareamax 
     66      INTEGER , DIMENSION(3)          ::   iareasum, iareamin, iareamax 
     67      INTEGER , DIMENSION(3,4)        ::   iloc                                  ! min/max loc indices 
     68      REAL(wp)                        ::   zzz                                   ! local real  
     69!!SWE      REAL(wp), DIMENSION(9)          ::   zmax, zmaxlocal 
     70      REAL(wp), DIMENSION(3)          ::   zmax, zmaxlocal 
     71      LOGICAL                         ::   ll_wrtstp, ll_colruns, ll_wrtruns 
     72      LOGICAL, DIMENSION(jpi,jpj,jpk) ::   llmsk 
     73      CHARACTER(len=20)               ::   clname 
     74      !!---------------------------------------------------------------------- 
     75      ! 
     76      IF( nstop > 0 .AND. ngrdstop > -1 )   RETURN   !   stpctl was already called by a child grid 
     77      ! 
     78      ll_wrtstp  = ( MOD( kt-nit000, sn_cfctl%ptimincr ) == 0 ) .OR. ( kt == nitend ) 
     79      ll_colruns = ll_wrtstp .AND. sn_cfctl%l_runstat .AND. jpnij > 1  
     80      ll_wrtruns = ( ll_colruns .OR. jpnij == 1 ) .AND. lwm 
     81      ! 
     82      IF( kt == nit000 ) THEN 
     83         ! 
     84         IF( lwp ) THEN 
     85            WRITE(numout,*) 
     86            WRITE(numout,*) 'stp_ctl : time-stepping control' 
     87            WRITE(numout,*) '~~~~~~~' 
     88         ENDIF 
     89         !                                ! open time.step    ascii file, done only by 1st subdomain 
     90         IF( lwm )   CALL ctl_opn( numstp, 'time.step', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp, narea ) 
     91         ! 
     92         IF( ll_wrtruns ) THEN 
     93            !                             ! open run.stat     ascii file, done only by 1st subdomain 
    8894            CALL ctl_opn( numrun, 'run.stat', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp, narea ) 
     95            !                             ! open run.stat.nc netcdf file, done only by 1st subdomain 
    8996            clname = 'run.stat.nc' 
    9097            IF( .NOT. Agrif_Root() )   clname = TRIM(Agrif_CFixed())//"_"//TRIM(clname) 
    91             istatus = NF90_CREATE( TRIM(clname), NF90_CLOBBER, idrun ) 
    92             istatus = NF90_DEF_DIM( idrun, 'time', NF90_UNLIMITED, idtime ) 
    93             istatus = NF90_DEF_VAR( idrun, 'abs_ssh_max', NF90_DOUBLE, (/ idtime /), idssh ) 
    94             istatus = NF90_DEF_VAR( idrun,   'abs_u_max', NF90_DOUBLE, (/ idtime /), idu   ) 
    95             istatus = NF90_ENDDEF(idrun) 
    96          ENDIF 
    97       ENDIF 
    98       IF( kt == nit000 )   lsomeoce = COUNT( ssmask(:,:) == 1._wp ) > 0 
    99       ! 
    100       IF(lwm .AND. ll_wrtstp) THEN        !==  current time step  ==!   ("time.step" file) 
     98            istatus = NF90_CREATE( TRIM(clname), NF90_CLOBBER, nrunid ) 
     99            istatus = NF90_DEF_DIM( nrunid, 'time', NF90_UNLIMITED, idtime ) 
     100            istatus = NF90_DEF_VAR( nrunid, 'abs_ssh_max', NF90_DOUBLE, (/ idtime /), nvarid(1) ) 
     101            istatus = NF90_DEF_VAR( nrunid,   'abs_u_max', NF90_DOUBLE, (/ idtime /), nvarid(2) ) 
     102!!SWE            istatus = NF90_DEF_VAR( nrunid,       's_min', NF90_DOUBLE, (/ idtime /), nvarid(3) ) 
     103!!SWE            istatus = NF90_DEF_VAR( nrunid,       's_max', NF90_DOUBLE, (/ idtime /), nvarid(4) ) 
     104!!SWE            istatus = NF90_DEF_VAR( nrunid,       't_min', NF90_DOUBLE, (/ idtime /), nvarid(5) ) 
     105!!SWE            istatus = NF90_DEF_VAR( nrunid,       't_max', NF90_DOUBLE, (/ idtime /), nvarid(6) ) 
     106!!SWE            IF( ln_zad_Aimp ) THEN 
     107!!SWE               istatus = NF90_DEF_VAR( nrunid,   'Cf_max', NF90_DOUBLE, (/ idtime /), nvarid(7) ) 
     108!!SWE               istatus = NF90_DEF_VAR( nrunid,'abs_wi_max',NF90_DOUBLE, (/ idtime /), nvarid(8) ) 
     109!!SWE            ENDIF 
     110            istatus = NF90_ENDDEF(nrunid) 
     111         ENDIF 
     112         !     
     113      ENDIF 
     114      ! 
     115      !                                   !==              write current time step              ==! 
     116      !                                   !==  done only by 1st subdomain at writting timestep  ==! 
     117      IF( lwm .AND. ll_wrtstp ) THEN 
    101118         WRITE ( numstp, '(1x, i8)' )   kt 
    102119         REWIND( numstp ) 
    103120      ENDIF 
    104       ! 
    105       !                                   !==  test of extrema  ==! 
    106       IF( ll_wd ) THEN 
    107          zmax(1) = MAXVAL(  ABS( ssh(:,:,Kmm) + ssh_ref*tmask(:,:,1) )  )        ! ssh max  
     121      !                                   !==            test of local extrema           ==! 
     122      !                                   !==  done by all processes at every time step  ==! 
     123!!SWE      llmsk(:,:,1) = ssmask(:,:) == 1._wp 
     124!!SWE      IF( ll_wd ) THEN 
     125!!SWE         zmax(1) = MAXVAL( ABS( ssh(:,:,Kmm) + ssh_ref ), mask = llmsk(:,:,1) )   ! ssh max 
     126!!SWE      ELSE 
     127!!SWE         zmax(1) = MAXVAL( ABS( ssh(:,:,Kmm)           ), mask = llmsk(:,:,1) )   ! ssh max 
     128!!SWE      ENDIF 
     129      zmax(1) = MINVAL( e3t_0(:,:,1)+ssh(:,:,Kmm)  )                              ! e3t_Kmm min 
     130!!SWE 
     131      llmsk(:,:,:) = umask(:,:,:) == 1._wp 
     132      zmax(2) = MAXVAL(  ABS( uu(:,:,:,Kmm) ), mask = llmsk )                     ! velocity max (zonal only) 
     133!!SWE      llmsk(:,:,:) = tmask(:,:,:) == 1._wp 
     134!!SWE      zmax(3) = MAXVAL( -ts(:,:,:,jp_sal,Kmm), mask = llmsk )                     ! minus salinity max 
     135!!SWE      zmax(4) = MAXVAL(  ts(:,:,:,jp_sal,Kmm), mask = llmsk )                     !       salinity max 
     136!!SWE      IF( ll_colruns .OR. jpnij == 1 ) THEN     ! following variables are used only in the netcdf file 
     137!!SWE         zmax(5) = MAXVAL( -ts(:,:,:,jp_tem,Kmm), mask = llmsk )                  ! minus temperature max 
     138!!SWE         zmax(6) = MAXVAL(  ts(:,:,:,jp_tem,Kmm), mask = llmsk )                  !       temperature max 
     139!!SWE         IF( ln_zad_Aimp ) THEN 
     140!!SWE            zmax(7) = MAXVAL(   Cu_adv(:,:,:)   , mask = llmsk )                  ! partitioning coeff. max 
     141!!SWE            llmsk(:,:,:) = wmask(:,:,:) == 1._wp 
     142!!SWE            zmax(8) = MAXVAL(  ABS( wi(:,:,:) ) , mask = llmsk )                  ! implicit vertical vel. max 
     143!!SWE         ELSE 
     144!!SWE            zmax(7:8) = 0._wp 
     145!!SWE         ENDIF 
     146!!SWE      ELSE 
     147!!SWE         zmax(5:8) = 0._wp 
     148!!SWE      ENDIF 
     149!!SWE      zmax(9) = REAL( nstop, wp )                                              ! stop indicator 
     150!!SWE 
     151      zmax(3) = REAL( nstop , wp )                                            ! stop indicator 
     152!!SWE 
     153      !                                   !==               get global extrema             ==! 
     154      !                                   !==  done by all processes if writting run.stat  ==! 
     155      IF( ll_colruns ) THEN 
     156         zmaxlocal(:) = zmax(:) 
     157         CALL mpp_max( "stpctl", zmax )          ! max over the global domain 
     158!!SWE         nstop = NINT( zmax(9) )                 ! update nstop indicator (now sheared among all local domains) 
     159         nstop = NINT( zmax(3) )                 ! update nstop indicator (now sheared among all local domains) 
     160      ENDIF 
     161      !                                   !==              write "run.stat" files              ==! 
     162      !                                   !==  done only by 1st subdomain at writting timestep  ==! 
     163      IF( ll_wrtruns ) THEN 
     164!!SWE         WRITE(numrun,9500) kt, zmax(1), zmax(2), -zmax(3), zmax(4) 
     165         WRITE(numrun,9500) kt, zmax(1), zmax(2) 
     166         istatus = NF90_PUT_VAR( nrunid, nvarid(1), (/ zmax(1)/), (/kt/), (/1/) ) 
     167         istatus = NF90_PUT_VAR( nrunid, nvarid(2), (/ zmax(2)/), (/kt/), (/1/) ) 
     168!!SWE         istatus = NF90_PUT_VAR( nrunid, nvarid(3), (/-zmax(3)/), (/kt/), (/1/) ) 
     169!!SWE         istatus = NF90_PUT_VAR( nrunid, nvarid(4), (/ zmax(4)/), (/kt/), (/1/) ) 
     170!!SWE         istatus = NF90_PUT_VAR( nrunid, nvarid(5), (/-zmax(5)/), (/kt/), (/1/) ) 
     171!!SWE         istatus = NF90_PUT_VAR( nrunid, nvarid(6), (/ zmax(6)/), (/kt/), (/1/) ) 
     172!!SWE         IF( ln_zad_Aimp ) THEN 
     173!!SWE            istatus = NF90_PUT_VAR( nrunid, nvarid(7), (/ zmax(7)/), (/kt/), (/1/) ) 
     174!!SWE            istatus = NF90_PUT_VAR( nrunid, nvarid(8), (/ zmax(8)/), (/kt/), (/1/) ) 
     175!!SWE         ENDIF 
     176         IF( kt == nitend )   istatus = NF90_CLOSE(nrunid) 
     177      ENDIF 
     178      !                                   !==               error handling               ==! 
     179      !                                   !==  done by all processes at every time step  ==! 
     180      ! 
     181!!SWE specific : start 
     182      IF(   zmax(1) <=   0._wp .OR.           &               ! negative e3t_Kmm 
     183         &  zmax(2) >   10._wp .OR.           &               ! too large velocity ( > 10 m/s) 
     184         &  ISNAN( zmax(1) + zmax(2) ) .OR.   &               ! NaN encounter in the tests 
     185         &  ABS(   zmax(1) + zmax(2) ) > HUGE(1._wp) ) THEN   ! Infinity encounter in the tests 
     186         ! 
     187         iloc(:,:) = 0 
     188         IF( ll_colruns ) THEN   ! zmax is global, so it is the same on all subdomains -> no dead lock with mpp_maxloc 
     189            ! first: close the netcdf file, so we can read it 
     190            IF( lwm .AND. kt /= nitend )   istatus = NF90_CLOSE(nrunid) 
     191            ! get global loc on the min/max 
     192            CALL mpp_minloc( 'stpctl', e3t_0(:,:,1) + ssh(:,:,Kmm), ssmask(:,:  ), zzz, iloc(1:2,1) )   ! mpp_maxloc ok if mask = F  
     193            CALL mpp_maxloc( 'stpctl', ABS( uu(:,:,:,Kmm))        ,  umask(:,:,:), zzz, iloc(1:3,2) ) 
     194            ! find which subdomain has the max. 
     195            iareamin(:) = jpnij+1   ;   iareamax(:) = 0   ;   iareasum(:) = 0 
     196            DO ji = 1, 3 
     197               IF( zmaxlocal(ji) == zmax(ji) ) THEN 
     198                  iareamin(ji) = narea   ;   iareamax(ji) = narea   ;   iareasum(ji) = 1 
     199               ENDIF 
     200            END DO 
     201            CALL mpp_min( "stpctl", iareamin )         ! min over the global domain 
     202            CALL mpp_max( "stpctl", iareamax )         ! max over the global domain 
     203            CALL mpp_sum( "stpctl", iareasum )         ! sum over the global domain 
     204         ELSE                    ! find local min and max locations: 
     205            ! if we are here, this means that the subdomain contains some oce points -> no need to test the mask used in maxloc 
     206            iloc(1:2,1) = MINLOC( e3t_0(:,:,1) + ssh(:,:,Kmm), mask = ssmask(:,:  ) == 1._wp ) + (/ nimpp - 1, njmpp - 1    /) 
     207            iloc(1:3,2) = MAXLOC( ABS(  uu(:,:,:,       Kmm)), mask =  umask(:,:,:) == 1._wp ) + (/ nimpp - 1, njmpp - 1, 0 /) 
     208            iareamin(:) = narea   ;   iareamax(:) = narea   ;   iareasum(:) = 1         ! this is local information 
     209         ENDIF 
     210         ! 
     211         WRITE(ctmp1,*) ' stp_ctl:  e3t0+ssh < 0 m  or  |U| > 10 m/s  or  NaN encounter in the tests' 
     212         CALL wrt_line( ctmp2, kt, 'e3t0+ssh min',  zmax(1), iloc(:,1), iareasum(1), iareamin(1), iareamax(1) ) 
     213         CALL wrt_line( ctmp3, kt, '|U|   max'   ,  zmax(2), iloc(:,2), iareasum(2), iareamin(2), iareamax(2) ) 
     214         IF( Agrif_Root() ) THEN 
     215            WRITE(ctmp6,*) '      ===> output of last computed fields in output.abort* files' 
     216         ELSE 
     217            WRITE(ctmp6,*) '      ===> output of last computed fields in '//TRIM(Agrif_CFixed())//'_output.abort* files' 
     218         ENDIF 
     219         ! 
     220         CALL dia_wri_state( Kmm, 'output.abort' )     ! create an output.abort file 
     221         ! 
     222         IF( ll_colruns .or. jpnij == 1 ) THEN   ! all processes synchronized -> use lwp to print in opened ocean.output files 
     223            IF(lwp) THEN   ;   CALL ctl_stop( ctmp1, ' ', ctmp2, ctmp3, ' ', ctmp6 ) 
     224            ELSE           ;   nstop = MAX(1, nstop)   ! make sure nstop > 0 (automatically done when calling ctl_stop) 
     225            ENDIF 
     226         ELSE                                    ! only mpi subdomains with errors are here -> STOP now 
     227            CALL ctl_stop( 'STOP', ctmp1, ' ', ctmp2, ctmp3, ' ', ctmp6 ) 
     228         ENDIF 
     229         ! 
     230      ENDIF 
     231!!SWE specific : end 
     232      ! 
     233      IF( nstop > 0 ) THEN                                                  ! an error was detected and we did not abort yet... 
     234         ngrdstop = Agrif_Fixed()                                           ! store which grid got this error 
     235         IF( .NOT. ll_colruns .AND. jpnij > 1 )   CALL ctl_stop( 'STOP' )   ! we must abort here to avoid MPI deadlock 
     236      ENDIF 
     237      ! 
     238!!SWE 9500  FORMAT(' it :', i8, '    |ssh|_max: ', D23.16, ' |U|_max: ', D23.16,' S_min: ', D23.16,' S_max: ', D23.16) 
     2399500  FORMAT(' it :', i8, '      e3t_min: ', D23.16, ' |U|_max: ', D23.16) 
     240      ! 
     241   END SUBROUTINE stp_ctl 
     242 
     243 
     244   SUBROUTINE wrt_line( cdline, kt, cdprefix, pval, kloc, ksum, kmin, kmax ) 
     245      !!---------------------------------------------------------------------- 
     246      !!                     ***  ROUTINE wrt_line  *** 
     247      !! 
     248      !! ** Purpose :   write information line 
     249      !! 
     250      !!---------------------------------------------------------------------- 
     251      CHARACTER(len=*),      INTENT(  out) ::   cdline 
     252      CHARACTER(len=*),      INTENT(in   ) ::   cdprefix 
     253      REAL(wp),              INTENT(in   ) ::   pval 
     254      INTEGER, DIMENSION(3), INTENT(in   ) ::   kloc 
     255      INTEGER,               INTENT(in   ) ::   kt, ksum, kmin, kmax 
     256      ! 
     257      CHARACTER(len=80) ::   clsuff 
     258      CHARACTER(len=9 ) ::   clkt, clsum, clmin, clmax 
     259      CHARACTER(len=9 ) ::   cli, clj, clk 
     260      CHARACTER(len=1 ) ::   clfmt 
     261      CHARACTER(len=4 ) ::   cl4   ! needed to be able to compile with Agrif, I don't know why 
     262      INTEGER           ::   ifmtk 
     263      !!---------------------------------------------------------------------- 
     264      WRITE(clkt , '(i9)') kt 
     265       
     266      WRITE(clfmt, '(i1)') INT(LOG10(REAL(jpnij  ,wp))) + 1     ! how many digits to we need to write ? (we decide max = 9) 
     267      !!! WRITE(clsum, '(i'//clfmt//')') ksum                   ! this is creating a compilation error with AGRIF 
     268      cl4 = '(i'//clfmt//')'   ;   WRITE(clsum, cl4) ksum 
     269      WRITE(clfmt, '(i1)') INT(LOG10(REAL(MAX(1,jpnij-1),wp))) + 1    ! how many digits to we need to write ? (we decide max = 9) 
     270      cl4 = '(i'//clfmt//')'   ;   WRITE(clmin, cl4) kmin-1 
     271                                   WRITE(clmax, cl4) kmax-1 
     272      ! 
     273      WRITE(clfmt, '(i1)') INT(LOG10(REAL(jpiglo,wp))) + 1      ! how many digits to we need to write jpiglo? (we decide max = 9) 
     274      cl4 = '(i'//clfmt//')'   ;   WRITE(cli, cl4) kloc(1)      ! this is ok with AGRIF 
     275      WRITE(clfmt, '(i1)') INT(LOG10(REAL(jpjglo,wp))) + 1      ! how many digits to we need to write jpjglo? (we decide max = 9) 
     276      cl4 = '(i'//clfmt//')'   ;   WRITE(clj, cl4) kloc(2)      ! this is ok with AGRIF 
     277      ! 
     278      IF( ksum == 1 ) THEN   ;   WRITE(clsuff,9100) TRIM(clmin) 
     279      ELSE                   ;   WRITE(clsuff,9200) TRIM(clsum), TRIM(clmin), TRIM(clmax) 
     280      ENDIF 
     281      IF(kloc(3) == 0) THEN 
     282         ifmtk = INT(LOG10(REAL(jpk,wp))) + 1                   ! how many digits to we need to write jpk? (we decide max = 9) 
     283         clk = REPEAT(' ', ifmtk)                               ! create the equivalent in blank string 
     284         WRITE(cdline,9300) TRIM(ADJUSTL(clkt)), TRIM(ADJUSTL(cdprefix)), pval, TRIM(cli), TRIM(clj), clk(1:ifmtk), TRIM(clsuff) 
    108285      ELSE 
    109          zmax(1) = MINVAL( e3t(:,:,1,Kmm)  )                                         ! ssh min 
    110       ENDIF 
    111       zmax(2) = MAXVAL(  ABS( uu(:,:,:,Kmm) )  )                                  ! velocity max (zonal only) 
    112       zmax(3) = REAL( nstop , wp )                                            ! stop indicator 
    113       ! 
    114       IF( ll_colruns ) THEN 
    115          CALL mpp_max( "stpctl", zmax )          ! max over the global domain 
    116          nstop = NINT( zmax(3) )                 ! nstop indicator sheared among all local domains 
    117       ENDIF 
    118       !                                   !==  run statistics  ==!   ("run.stat" files) 
    119       IF( ll_wrtruns ) THEN 
    120          WRITE(numrun,9500) kt, zmax(1), zmax(2) 
    121          istatus = NF90_PUT_VAR( idrun, idssh, (/ zmax(1)/), (/kt/), (/1/) ) 
    122          istatus = NF90_PUT_VAR( idrun,   idu, (/ zmax(2)/), (/kt/), (/1/) ) 
    123          IF( MOD( kt , 100 ) == 0 ) istatus = NF90_SYNC(idrun) 
    124          IF( kt == nitend         ) istatus = NF90_CLOSE(idrun) 
    125       END IF 
    126       !                                   !==  error handling  ==! 
    127       IF( ( sn_cfctl%l_glochk .OR. lsomeoce ) .AND. (   &  ! domain contains some ocean points, check for sensible ranges 
    128          &  zmax(1) <    0._wp .OR.   &                    ! negative sea surface height  
    129          &  zmax(2) >   10._wp .OR.   &                    ! too large velocity ( > 10 m/s) 
    130          &  ISNAN( zmax(1) + zmax(2) ) ) ) THEN            ! NaN encounter in the tests 
    131          IF( lk_mpp .AND. sn_cfctl%l_glochk ) THEN 
    132             ! have use mpp_max (because sn_cfctl%l_glochk=.T. and distributed) 
    133             CALL mpp_maxloc( 'stpctl', ABS(ssh(:,:,Kmm))        , ssmask(:,:)  , zzz, ih  ) 
    134             CALL mpp_maxloc( 'stpctl', ABS(uu(:,:,:,Kmm))          , umask (:,:,:), zzz, iu  ) 
    135          ELSE 
    136             ! find local min and max locations 
    137             ih(:)  = MAXLOC( ABS( ssh(:,:,Kmm)   )                              ) + (/ nimpp - 1, njmpp - 1    /) 
    138             iu(:)  = MAXLOC( ABS( uu  (:,:,:,Kmm) )                              ) + (/ nimpp - 1, njmpp - 1, 0 /) 
    139          ENDIF 
    140           
    141          WRITE(ctmp1,*) ' stp_ctl: (e3t0) ssh < 0 m  or  |U| > 10 m/s  or  NaN encounter in the tests' 
    142          WRITE(ctmp2,9100) kt,   zmax(1), ih(1) , ih(2) 
    143          WRITE(ctmp3,9200) kt,   zmax(2), iu(1) , iu(2) , iu(3) 
    144          WRITE(ctmp4,*) '      ===> output of last computed fields in output.abort.nc file' 
    145           
    146          CALL dia_wri_state( Kmm, 'output.abort' )     ! create an output.abort file 
    147           
    148          IF( .NOT. sn_cfctl%l_glochk ) THEN 
    149             WRITE(ctmp8,*) 'E R R O R message from sub-domain: ', narea 
    150             CALL ctl_stop( 'STOP', ctmp1, ' ', ctmp8, ' ', ctmp2, ctmp3, ctmp4 ) 
    151          ELSE 
    152             CALL ctl_stop( ctmp1, ' ', ctmp2, ctmp3, ctmp4 ) 
    153          ENDIF 
    154  
    155          kindic = -3 
    156          ! 
    157       ENDIF 
    158       ! 
    159 9100  FORMAT (' kt=',i8,'   |ssh| min: ',1pg11.4,', at  i j  : ',2i5) 
    160 9200  FORMAT (' kt=',i8,'   |U|   max: ',1pg11.4,', at  i j k: ',3i5) 
    161 9500  FORMAT(' it :', i8, '    |ssh|_max: ', D23.16, ' |U|_max: ', D23.16) 
    162       ! 
    163    END SUBROUTINE stp_ctl 
     286         WRITE(clfmt, '(i1)') INT(LOG10(REAL(jpk,wp))) + 1      ! how many digits to we need to write jpk? (we decide max = 9) 
     287         !!! WRITE(clk, '(i'//clfmt//')') kloc(3)               ! this is creating a compilation error with AGRIF 
     288         cl4 = '(i'//clfmt//')'   ;   WRITE(clk, cl4) kloc(3)   ! this is ok with AGRIF 
     289         WRITE(cdline,9400) TRIM(ADJUSTL(clkt)), TRIM(ADJUSTL(cdprefix)), pval, TRIM(cli), TRIM(clj),    TRIM(clk), TRIM(clsuff) 
     290      ENDIF 
     291      ! 
     2929100  FORMAT('MPI rank ', a) 
     2939200  FORMAT('found in ', a, ' MPI tasks, spread out among ranks ', a, ' to ', a) 
     2949300  FORMAT('kt ', a, ' ', a, ' ', 1pg11.4, ' at i j   ', a, ' ', a, ' ', a, ' ', a) 
     2959400  FORMAT('kt ', a, ' ', a, ' ', 1pg11.4, ' at i j k ', a, ' ', a, ' ', a, ' ', a) 
     296      ! 
     297   END SUBROUTINE wrt_line 
     298 
    164299 
    165300   !!====================================================================== 
Note: See TracChangeset for help on using the changeset viewer.