New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Changeset 15574 for NEMO/branches/2021/dev_r14318_RK3_stage1/src/OCE/stpmlf.F90 – NEMO

Ignore:
Timestamp:
2021-12-03T20:32:50+01:00 (3 years ago)
Author:
techene
Message:

#2605 #2715 trunk merged into dev_r14318_RK3_stage1

Location:
NEMO/branches/2021/dev_r14318_RK3_stage1
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2021/dev_r14318_RK3_stage1

    • Property svn:externals
      •  

        old new  
        99 
        1010# SETTE 
        11 ^/utils/CI/sette@14244        sette 
         11^/utils/CI/sette@HEAD        sette 
         12 
  • NEMO/branches/2021/dev_r14318_RK3_stage1/src/OCE/stpmlf.F90

    r15193 r15574  
    9595      INTEGER ::   ji, jj, jk, jn, jtile   ! dummy loop indice 
    9696      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:)       ::   zgdept 
    97       REAL(wp), TARGET     , DIMENSION(jpi,jpj,jpk) ::   zau, zav, zaw   ! advective velocity 
    9897      !! --------------------------------------------------------------------- 
    9998#if defined key_agrif 
     
    142141#endif 
    143142      ENDIF 
     143      IF( kstp + nn_fsbc - 1 == nitrst .AND. lwxios ) THEN 
    144144#if defined key_si3 
    145       IF( kstp + nn_fsbc - 1 == nitrst .AND. lwxios ) THEN 
    146145                             CALL iom_swap(                     cw_icerst_cxt ) 
    147146                             CALL iom_init_closedef(            cw_icerst_cxt ) 
    148147                             CALL iom_setkt( kstp - nit000 + 1, cw_icerst_cxt ) 
    149       ENDIF 
    150 #endif 
     148#endif 
     149         IF( ln_abl      ) THEN 
     150                             CALL iom_swap(                     cw_ablrst_cxt ) 
     151                             CALL iom_init_closedef(            cw_ablrst_cxt ) 
     152                             CALL iom_setkt( kstp - nit000 + 1, cw_ablrst_cxt ) 
     153         ENDIF 
     154      ENDIF 
    151155      IF( kstp /= nit000 )   CALL day( kstp )         ! Calendar (day was already called at nit000 in day_init) 
    152156                             CALL iom_setkt( kstp - nit000 + 1,      cxios_context          )   ! tell IOM we are at time step kstp 
     
    178182 
    179183      !  VERTICAL PHYSICS 
     184      IF( ln_tile ) CALL dom_tile_start         ! [tiling] ZDF tiling loop 
     185      DO jtile = 1, nijtile 
     186         IF( ln_tile ) CALL dom_tile( ntsi, ntsj, ntei, ntej, ktile = jtile ) 
    180187                         CALL zdf_phy( kstp, Nbb, Nnn, Nrhs )   ! vertical physics update (top/bot drag, avt, avs, avm + MLD) 
     188      END DO 
     189      IF( ln_tile ) CALL dom_tile_stop 
    181190 
    182191      !  LATERAL  PHYSICS 
    183192      ! 
    184       IF( l_ldfslp ) THEN                             ! slope of lateral mixing 
    185                          CALL eos( ts(:,:,:,:,Nbb), rhd, gdept_0(:,:,:) )               ! before in situ density 
    186  
    187          IF( ln_zps .AND. .NOT. ln_isfcav)                                    & 
     193      IF( ln_zps .OR. l_ldfslp ) CALL eos( ts(:,:,:,:,Nbb), rhd, gdept_0(:,:,:) )               ! before in situ density 
     194 
     195      IF( ln_zps .AND. .NOT. ln_isfcav)                                    & 
    188196            &            CALL zps_hde    ( kstp, Nnn, jpts, ts(:,:,:,:,Nbb), gtsu, gtsv,  &  ! Partial steps: before horizontal gradient 
    189197            &                                          rhd, gru , grv    )       ! of t, s, rd at the last ocean level 
    190198 
    191          IF( ln_zps .AND.       ln_isfcav)                                                & 
     199      IF( ln_zps .AND.       ln_isfcav)                                                & 
    192200            &            CALL zps_hde_isf( kstp, Nnn, jpts, ts(:,:,:,:,Nbb), gtsu, gtsv, gtui, gtvi,  &  ! Partial steps for top cell (ISF) 
    193201            &                                          rhd, gru , grv , grui, grvi   )       ! of t, s, rd at the first ocean level 
     202 
     203      IF( l_ldfslp ) THEN                             ! slope of lateral mixing 
    194204         IF( ln_traldf_triad ) THEN 
    195205                         CALL ldf_slp_triad( kstp, Nbb, Nnn )             ! before slope for triad operator 
     
    223233                         uu(:,:,:,Nrhs) = 0._wp            ! set dynamics trends to zero 
    224234                         vv(:,:,:,Nrhs) = 0._wp 
    225 !!st                          
    226                          zau(:,:,:) = uu(:,:,:,Nnn)        !!st   patch for MLF will be computed in RK3 
    227                          zav(:,:,:) = vv(:,:,:,Nnn) 
    228                          zaw(:,:,:) = ww(:,:,:    ) 
    229 !!st 
    230       IF(  lk_asminc .AND. ln_asmiau .AND. ln_dyninc )   & 
    231                &         CALL dyn_asm_inc   ( kstp, Nbb, Nnn, uu, vv, Nrhs )  ! apply dynamics assimilation increment 
    232       IF( ln_bdy     )   CALL bdy_dyn3d_dmp ( kstp, Nbb,      uu, vv, Nrhs )  ! bdy damping trends 
    233 #if defined key_agrif 
     235 
     236      IF( ln_dyndmp .AND. ln_c1d )  CALL dyn_dmp( kstp, Nbb, Nnn, uu(:,:,:,Nrhs), vv(:,:,:,Nrhs), Nrhs )   ! internal damping trends- momentum 
     237 
     238      IF( ln_tile ) CALL dom_tile_start         ! [tiling] DYN tiling loop (1) 
     239      DO jtile = 1, nijtile 
     240         IF( ln_tile ) CALL dom_tile( ntsi, ntsj, ntei, ntej, ktile = jtile ) 
     241 
     242         IF(  lk_asminc .AND. ln_asmiau .AND. ln_dyninc )   & 
     243                  &         CALL dyn_asm_inc   ( kstp, Nbb, Nnn, uu, vv, Nrhs )  ! apply dynamics assimilation increment 
     244         IF( ln_bkgwri )    CALL asm_bkg_wri( kstp, Nnn )     ! output background fields 
     245         IF( ln_bdy     )   CALL bdy_dyn3d_dmp ( kstp, Nbb,      uu, vv, Nrhs )  ! bdy damping trends 
     246#if defined key_agrif 
     247      END DO 
     248      IF( ln_tile ) CALL dom_tile_stop 
     249 
    234250      IF(.NOT. Agrif_Root())  & 
    235251               &         CALL Agrif_Sponge_dyn        ! momentum sponge 
    236 #endif 
    237                          CALL dyn_adv( kstp, Nbb, Nnn      , uu, vv, Nrhs, zau, zav, zaw )  ! advection (VF or FF)   ==> RHS 
    238                          CALL dyn_vor( kstp,      Nnn      , uu, vv, Nrhs )  ! vorticity              ==> RHS 
    239                          CALL dyn_ldf( kstp, Nbb, Nnn      , uu, vv, Nrhs )  ! lateral mixing 
    240       IF( ln_zdfosm  )   CALL dyn_osm( kstp,      Nnn      , uu, vv, Nrhs )  ! OSMOSIS non-local velocity fluxes ==> RHS 
    241                          CALL dyn_hpg( kstp,      Nnn      , uu, vv, Nrhs )  ! horizontal gradient of Hydrostatic pressure 
    242                          CALL dyn_spg( kstp, Nbb, Nnn, Nrhs, uu, vv, ssh, uu_b, vv_b, Naa )  ! surface pressure gradient 
    243                           
    244       IF( ln_dynspg_ts ) THEN      ! With split-explicit free surface, since now transports have been updated and ssh(:,:,Nrhs) 
    245                                    ! as well as vertical scale factors and vertical velocity need to be updated 
    246                             CALL div_hor    ( kstp, Nbb, Nnn )                ! Horizontal divergence  (2nd call in time-split case) 
    247          IF(.NOT.lk_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  
    248       ENDIF 
     252 
     253      IF( ln_tile ) CALL dom_tile_start         ! [tiling] DYN tiling loop (1, continued) 
     254      DO jtile = 1, nijtile 
     255         IF( ln_tile ) CALL dom_tile( ntsi, ntsj, ntei, ntej, ktile = jtile ) 
     256#endif 
     257                            CALL dyn_adv( kstp, Nbb, Nnn      , uu, vv, Nrhs )  ! advection (VF or FF)   ==> RHS 
     258                            CALL dyn_vor( kstp,      Nnn      , uu, vv, Nrhs )  ! vorticity              ==> RHS 
     259                            CALL dyn_ldf( kstp, Nbb, Nnn      , uu, vv, Nrhs )  ! lateral mixing 
     260         IF( ln_zdfosm  )   CALL dyn_osm( kstp,      Nnn      , uu, vv, Nrhs )  ! OSMOSIS non-local velocity fluxes ==> RHS 
     261                            CALL dyn_hpg( kstp,      Nnn      , uu, vv, Nrhs )  ! horizontal gradient of Hydrostatic pressure 
     262      END DO 
     263      IF( ln_tile ) CALL dom_tile_stop 
     264 
     265                            CALL dyn_spg( kstp, Nbb, Nnn, Nrhs, uu, vv, ssh, uu_b, vv_b, Naa )  ! surface pressure gradient 
     266 
     267      IF( ln_tile ) CALL dom_tile_start         ! [tiling] DYN tiling loop (2) 
     268      DO jtile = 1, nijtile 
     269         IF( ln_tile ) CALL dom_tile( ntsi, ntsj, ntei, ntej, ktile = jtile ) 
     270 
     271         IF( ln_dynspg_ts ) THEN      ! With split-explicit free surface, since now transports have been updated and ssh(:,:,Nrhs) 
     272                                      ! as well as vertical scale factors and vertical velocity need to be updated 
     273                            CALL div_hor    ( kstp, Nbb, Nnn )                  ! Horizontal divergence  (2nd call in time-split case) 
     274            IF(.NOT.lk_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 
     275         ENDIF 
    249276                            CALL dyn_zdf    ( kstp, Nbb, Nnn, Nrhs, uu, vv, Naa  )  ! vertical diffusion 
     277      END DO 
     278      IF( ln_tile ) CALL dom_tile_stop 
     279 
    250280      IF( ln_dynspg_ts ) THEN                                                       ! vertical scale factors and vertical velocity need to be updated 
    251281                            CALL wzv        ( kstp, Nbb, Nnn, Naa, ww )             ! Nnn cross-level velocity 
    252282         IF( ln_zad_Aimp )  CALL wAimp      ( kstp,      Nnn )                      ! Adaptive-implicit vertical advection partitioning 
    253283      ENDIF 
    254 !!st 
    255                             zau(:,:,:) = uu(:,:,:,Nnn)        !!st   patch for MLF will be computed in RK3 
    256                             zav(:,:,:) = vv(:,:,:,Nnn) 
    257                             zaw(:,:,:) = ww(:,:,:    ) 
    258 !!st 
     284 
    259285      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    260286      ! cool skin 
     
    291317      ! Active tracers 
    292318      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    293       ! Loop over tile domains 
     319 
     320      IF( ln_tile )   CALL dom_tile_start         ! [tiling] TRA tiling loop (1) 
     321 
    294322      DO jtile = 1, nijtile 
    295          IF( ln_tile    )   CALL dom_tile( ntsi, ntsj, ntei, ntej, ktile = jtile ) 
     323         IF( ln_tile )   CALL dom_tile( ntsi, ntsj, ntei, ntej, ktile = jtile ) 
    296324 
    297325         DO jn = 1, jpts 
     
    311339         IF( ln_bdy     )   CALL bdy_tra_dmp( kstp, Nbb,      ts, Nrhs )  ! bdy damping trends 
    312340      END DO 
     341      IF( ln_tile ) CALL dom_tile_stop 
    313342 
    314343#if defined key_agrif 
    315344      IF(.NOT. Agrif_Root() ) THEN 
    316          IF( ln_tile    )   CALL dom_tile( ntsi, ntsj, ntei, ntej, ktile = 0 ) 
    317345                            CALL Agrif_Sponge_tra        ! tracers sponge 
    318346      ENDIF 
     
    320348 
    321349      ! TEMP: [tiling] Separate loop over tile domains (due to tra_adv workarounds for tiling) 
     350      IF( ln_tile ) CALL dom_tile_start         ! [tiling] TRA tiling loop (2) 
    322351      DO jtile = 1, nijtile 
    323          IF( ln_tile    )  CALL dom_tile( ntsi, ntsj, ntei, ntej, ktile = jtile ) 
    324  
    325                             CALL tra_adv    ( kstp, Nbb, Nnn, ts, Nrhs, zau, zav, zaw )  ! hor. + vert. advection ==> RHS 
     352         IF( ln_tile ) CALL dom_tile( ntsi, ntsj, ntei, ntej, ktile = jtile ) 
     353 
     354                            CALL tra_adv    ( kstp, Nbb, Nnn, ts, Nrhs )  ! hor. + vert. advection ==> RHS 
    326355         IF( ln_zdfmfc  )   CALL tra_mfc    ( kstp, Nbb,      ts, Nrhs )  ! Mass Flux Convection 
    327356         IF( ln_zdfosm  ) THEN 
     
    334363         IF( ln_zdfnpc  )   CALL tra_npc    ( kstp,      Nnn, Nrhs, ts, Naa  )  ! update after fields by non-penetrative convection 
    335364      END DO 
    336  
    337       IF( ln_tile ) CALL dom_tile( ntsi, ntsj, ntei, ntej, ktile = 0 ) ! Revert to tile over full domain 
     365      IF( ln_tile ) CALL dom_tile_stop 
     366 
    338367      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    339368      ! Set boundary conditions, time filter and swap time levels 
     
    456485      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) ::   puu, pvv   ! velocities 
    457486      ! 
    458       INTEGER  ::   jk   ! dummy loop indices 
     487      INTEGER  ::   ji,jj, jk   ! dummy loop indices 
    459488      REAL(wp), DIMENSION(jpi,jpj) ::   zue, zve 
    460489      !!---------------------------------------------------------------------- 
     
    462491      ! Ensure below that barotropic velocities match time splitting estimate 
    463492      ! Compute actual transport and replace it with ts estimate at "after" time step 
    464       zue(:,:) = e3u(:,:,1,Kaa) * puu(:,:,1,Kaa) * umask(:,:,1) 
    465       zve(:,:) = e3v(:,:,1,Kaa) * pvv(:,:,1,Kaa) * vmask(:,:,1) 
     493      DO_2D( 0, 0, 0, 0 ) 
     494         zue(ji,jj) = e3u(ji,jj,1,Kaa) * puu(ji,jj,1,Kaa) * umask(ji,jj,1) 
     495         zve(ji,jj) = e3v(ji,jj,1,Kaa) * pvv(ji,jj,1,Kaa) * vmask(ji,jj,1) 
     496      END_2D 
    466497      DO jk = 2, jpkm1 
    467          zue(:,:) = zue(:,:) + e3u(:,:,jk,Kaa) * puu(:,:,jk,Kaa) * umask(:,:,jk) 
    468          zve(:,:) = zve(:,:) + e3v(:,:,jk,Kaa) * pvv(:,:,jk,Kaa) * vmask(:,:,jk) 
     498         DO_2D( 0, 0, 0, 0 ) 
     499            zue(ji,jj) = zue(ji,jj) + e3u(ji,jj,jk,Kaa) * puu(ji,jj,jk,Kaa) * umask(ji,jj,jk) 
     500            zve(ji,jj) = zve(ji,jj) + e3v(ji,jj,jk,Kaa) * pvv(ji,jj,jk,Kaa) * vmask(ji,jj,jk) 
     501         END_2D 
    469502      END DO 
    470503      DO jk = 1, jpkm1 
    471          puu(:,:,jk,Kaa) = ( puu(:,:,jk,Kaa) - zue(:,:) * r1_hu(:,:,Kaa) + uu_b(:,:,Kaa) ) * umask(:,:,jk) 
    472          pvv(:,:,jk,Kaa) = ( pvv(:,:,jk,Kaa) - zve(:,:) * r1_hv(:,:,Kaa) + vv_b(:,:,Kaa) ) * vmask(:,:,jk) 
     504         DO_2D( 0, 0, 0, 0 ) 
     505            puu(ji,jj,jk,Kaa) = ( puu(ji,jj,jk,Kaa) - zue(ji,jj) * r1_hu(ji,jj,Kaa) + uu_b(ji,jj,Kaa) ) * umask(ji,jj,jk) 
     506            pvv(ji,jj,jk,Kaa) = ( pvv(ji,jj,jk,Kaa) - zve(ji,jj) * r1_hv(ji,jj,Kaa) + vv_b(ji,jj,Kaa) ) * vmask(ji,jj,jk) 
     507         END_2D 
    473508      END DO 
    474509      ! 
     
    518553# endif 
    519554      !                                        ! local domain boundaries  (T-point, unchanged sign) 
    520       CALL lbc_lnk_multi( 'finalize_lbc', puu(:,:,:,       Kaa), 'U', -1., pvv(:,:,:       ,Kaa), 'V', -1.   & 
    521                        &                , pts(:,:,:,jp_tem,Kaa), 'T',  1., pts(:,:,:,jp_sal,Kaa), 'T',  1. ) 
    522       ! 
     555      CALL lbc_lnk( 'finalize_lbc', puu(:,:,:,       Kaa), 'U', -1., pvv(:,:,:       ,Kaa), 'V', -1.   & 
     556                       &          , pts(:,:,:,jp_tem,Kaa), 'T',  1., pts(:,:,:,jp_sal,Kaa), 'T',  1. ) 
     557      ! 
     558      ! lbc_lnk needed for zdf_sh2 when using nn_hls = 2, moved here to allow tiling in zdf_phy 
     559      IF( nn_hls == 2 .AND. l_zdfsh2 ) CALL lbc_lnk( 'stp', avm_k, 'W', 1.0_wp ) 
     560 
     561      ! dom_qco_r3c defines over [nn_hls, nn_hls-1, nn_hls, nn_hls-1] 
     562      IF( nn_hls == 2 .AND. .NOT. lk_linssh ) THEN 
     563         CALL lbc_lnk( 'finalize_lbc', r3u(:,:,Kaa), 'U', 1._wp, r3v(:,:,Kaa), 'V', 1._wp, & 
     564            &                          r3u_f(:,:),   'U', 1._wp, r3v_f(:,:),   'V', 1._wp ) 
     565      ENDIF 
    523566      !                                        !* BDY open boundaries 
    524567      IF( ln_bdy )   THEN 
Note: See TracChangeset for help on using the changeset viewer.