Changeset 13604 for NEMO/branches/2020
- Timestamp:
- 2020-10-14T18:07:31+02:00 (4 years ago)
- Location:
- NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/src/SWE
- Files:
-
- 1 added
- 14 deleted
- 3 edited
- 1 copied
-
asminc.F90 (deleted)
-
diawri.F90 (deleted)
-
dom_oce.F90 (deleted)
-
domain.F90 (deleted)
-
dommsk.F90 (deleted)
-
domvvl.F90 (deleted)
-
dynadv.F90 (deleted)
-
dynatf.F90 (deleted)
-
dynkeg.F90 (deleted)
-
dynldf_lap_blp.F90 (deleted)
-
dynvor.F90 (modified) (1 diff)
-
ldfdyn.F90 (deleted)
-
phycst.F90 (deleted)
-
restart.F90 (added)
-
sbcice_cice.F90 (deleted)
-
step.F90 (deleted)
-
stpRK3.F90 (modified) (15 diffs)
-
stp_oce.F90 (copied) (copied from NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/src/OCE/step_oce.F90) (8 diffs)
-
stpctl.F90 (modified) (5 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/src/SWE/dynvor.F90
r13295 r13604 447 447 IF( ln_dynvor_msk ) THEN ! mask the relative vorticity 448 448 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) 450 450 END_2D 451 451 ENDIF -
NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/src/SWE/stpRK3.F90
r13295 r13604 1 1 MODULE stpRK3 2 2 !!====================================================================== 3 !! *** MODULE st ep***3 !! *** MODULE stpRK3 *** 4 4 !! Time-stepping : manager of the shallow water equation time stepping 5 !! 3rd order Runge-Kutta scheme 5 6 !!====================================================================== 6 7 !! 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 15 15 ! 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 18 19 19 20 IMPLICIT NONE … … 21 22 22 23 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 28 27 29 28 !! * Substitutions … … 37 36 CONTAINS 38 37 39 #if defined key_agrif40 RECURSIVE SUBROUTINE stp_RK3( )41 INTEGER :: kstp ! ocean time-step index42 #else43 38 SUBROUTINE stp_RK3( kstp ) 44 INTEGER, INTENT(in) :: kstp ! ocean time-step index45 #endif46 39 !!---------------------------------------------------------------------- 47 40 !! *** ROUTINE stp_RK3 *** 48 41 !! 49 !! ** Purpose : - Time stepping of shallow water (SHW) (momentum and ssh eqs.)42 !! ** Purpose : - RK3 Time stepping scheme for shallow water Eq. 50 43 !! 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 57 65 !!---------------------------------------------------------------------- 66 INTEGER, INTENT(in ) :: kstp ! ocean time-step index 67 ! 58 68 INTEGER :: ji, jj, jk ! dummy loop indice 59 69 INTEGER :: indic ! error indicator if < 0 60 !!gm kcall can be removed, I guess61 INTEGER :: kcall ! optional integer argument (dom_vvl_sf_nxt)62 70 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 ! - - 67 73 !! --------------------------------------------------------------------- 68 #if defined key_agrif69 kstp = nit000 + Agrif_Nb_Step()70 Kbb_a = Nbb; Kmm_a = Nnn; Krhs_a = Nrhs ! agrif_oce module copies of time level indices71 IF( lk_agrif_debug ) THEN72 IF( Agrif_Root() .and. lwp) WRITE(*,*) '---'73 IF(lwp) WRITE(*,*) 'Grid Number', Agrif_Fixed(),' time step ', kstp, 'int tstep', Agrif_NbStepint()74 ENDIF75 IF( kstp == nit000 + 1 ) lk_agrif_fstep = .FALSE.76 # if defined key_iomput77 IF( Agrif_Nbstepint() == 0 ) CALL iom_swap( cxios_context )78 # endif79 #endif80 74 ! 81 75 IF( ln_timing ) CALL timing_start('stp_RK3') … … 93 87 indic = 0 ! reset to no error condition 94 88 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) 98 91 CALL iom_init_closedef 99 IF( ln_crs ) CALL iom_init( TRIM(cxios_context)//"_crs" ) ! for coarse grid100 92 ENDIF 101 93 IF( kstp /= nit000 ) CALL day( kstp ) ! Calendar (day was already called at nit000 in day_init) 102 94 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 118 106 IF( l_ldfdyn_time ) CALL ldf_dyn( kstp, Nbb ) ! eddy viscosity coeff. 119 120 107 121 108 !====================================================================== … … 127 114 128 115 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 129 ! RK3 1st stage Ocean dynamics : hdiv, ssh, e3, u, v, w116 ! RK3 1st stage Ocean dynamics : u, v, ssh 130 117 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 131 118 rDt = rn_Dt / 3._wp 132 119 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 142 128 #if defined key_RK3all 143 CALL dyn_ldf( kstp, Nbb, Nbb, uu, vv, Nrhs ) ! lateral mixing144 #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 148 134 uu(ji,jj,jk,Nrhs) = uu(ji,jj,jk,Nrhs) - grav * ( ssh(ji+1,jj,Nbb) - ssh(ji,jj,Nbb) ) * r1_e1u(ji,jj) 149 135 vv(ji,jj,jk,Nrhs) = vv(ji,jj,jk,Nrhs) - grav * ( ssh(ji,jj+1,Nbb) - ssh(ji,jj,Nbb) ) * r1_e2v(ji,jj) … … 151 137 ! 152 138 #if defined key_RK3all 153 ! add wind stress forcing and layer linear friction to the RHS139 ! ! wind stress and layer friction 154 140 z5_6 = 5._wp/6._wp 155 141 DO_3D( 0, 0, 0, 0,1,jpkm1) … … 160 146 END_3D 161 147 #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 165 177 DO_3D( 0, 0, 0, 0,1,jpkm1) 166 178 uu(ji,jj,jk,Naa) = uu(ji,jj,jk,Nbb) + rDt * uu(ji,jj,jk,Nrhs) * umask(ji,jj,jk) … … 168 180 END_3D 169 181 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) 177 190 END_3D 178 191 ENDIF 179 ! Swap time levels 192 ! 193 ! !== Swap time levels ==! 194 ! 180 195 Nrhs= Nnn 181 196 Nnn = Naa 182 197 Naa = Nrhs 183 198 184 185 199 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 186 200 ! RK3 2nd stage Ocean dynamics : hdiv, ssh, e3, u, v, w … … 188 202 rDt = rn_Dt / 2._wp 189 203 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 199 211 #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 ! 205 215 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 216 ! ! horizontal pressure gradient 206 217 uu(ji,jj,jk,Nrhs) = uu(ji,jj,jk,Nrhs) - grav * ( ssh(ji+1,jj,Nnn) - ssh(ji,jj,Nnn) ) * r1_e1u(ji,jj) 207 218 vv(ji,jj,jk,Nrhs) = vv(ji,jj,jk,Nrhs) - grav * ( ssh(ji,jj+1,Nnn) - ssh(ji,jj,Nnn) ) * r1_e2v(ji,jj) 208 219 END_3D 209 220 ! 210 ! add wind stress forcing and layer linear friction to the RHS211 221 #if defined key_RK3all 222 ! ! wind stress and layer friction 212 223 z3_4 = 3._wp/4._wp 213 224 DO_3D( 0, 0, 0, 0,1,jpkm1) … … 218 229 END_3D 219 230 #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 223 241 DO_3D( 0, 0, 0, 0,1,jpkm1) 224 242 uu(ji,jj,jk,Naa) = uu(ji,jj,jk,Nbb) + rDt * uu(ji,jj,jk,Nrhs) * umask(ji,jj,jk) … … 226 244 END_3D 227 245 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) 235 254 END_3D 236 255 ENDIF 237 ! Swap time levels 256 ! 257 ! !== Swap time levels ==! 258 ! 238 259 Nrhs= Nnn 239 260 Nnn = Naa … … 245 266 rDt = rn_Dt 246 267 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 266 278 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 279 ! ! horizontal pressure gradient 267 280 uu(ji,jj,jk,Nrhs) = uu(ji,jj,jk,Nrhs) - grav * ( ssh(ji+1,jj,Nnn) - ssh(ji,jj,Nnn) ) * r1_e1u(ji,jj) 268 281 vv(ji,jj,jk,Nrhs) = vv(ji,jj,jk,Nrhs) - grav * ( ssh(ji,jj+1,Nnn) - ssh(ji,jj,Nnn) ) * r1_e2v(ji,jj) 269 282 END_3D 270 ! 271 ! add wind stress forcing and layer linear friction to the RHS 283 ! ! wind stress and layer friction 272 284 z1_2rho0 = 0.5_wp * r1_rho0 273 285 DO_3D( 0, 0, 0, 0,1,jpkm1) … … 276 288 vv(ji,jj,jk,Nrhs) = vv(ji,jj,jk,Nrhs) + z1_2rho0 * ( vtau_b(ji,jj) + vtau(ji,jj) ) / e3v(ji,jj,jk,Nnn) & 277 289 & - 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 282 301 DO_3D( 1, 1, 1, 1,1,jpkm1) 283 302 zua = uu(ji,jj,jk,Nbb) + rDt * uu(ji,jj,jk,Nrhs) * umask(ji,jj,jk) … … 291 310 END_3D 292 311 ! 293 ELSE ! flux form : applied on thickness weighted velocity312 ELSE ! flux form : applied on thickness weighted velocity 294 313 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)297 314 zue3b = e3u(ji,jj,jk,Nbb) * uu(ji,jj,jk,Nbb) 298 315 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) 302 318 ! 303 uu(ji,jj,jk,Naa) = zue3a / e3 t(ji,jj,jk,Naa)304 vv(ji,jj,jk,Naa) = zve3a / e3 t(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) 305 321 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 309 323 310 324 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 ! 317 329 Nrhs = Nbb 318 330 Nbb = Naa 319 331 Naa = Nrhs 320 ! 321 ! CALL dom_vvl_sf_update_st( kstp, Nbb, Nnn, Naa ) ! recompute vertical scale factors 332 322 333 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 323 334 ! diagnostics and outputs 324 335 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 325 IF( ln_floats ) CALL flo_stp ( kstp, Nbb, Nnn ) ! drifting Floats336 326 337 IF( ln_diacfl ) CALL dia_cfl ( kstp, Nnn ) ! Courant number diagnostics 327 328 338 CALL dia_wri ( kstp, Nnn ) ! ocean model: outputs 329 330 339 ! 331 340 IF( lrst_oce ) CALL rst_write ( kstp, Nbb, Nnn ) ! write output ocean restart file 332 333 334 #if defined key_agrif335 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>336 ! AGRIF337 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<338 Kbb_a = Nbb; Kmm_a = Nnn; Krhs_a = Nrhs ! agrif_oce module copies of time level indices339 CALL Agrif_Integrate_ChildGrids( stp_RK3 ) ! allows to finish all the Child Grids before updating340 341 IF( Agrif_NbStepint() == 0 ) THEN342 CALL Agrif_update_all( ) ! Update all components343 ENDIF344 #endif345 IF( ln_diaobs ) CALL dia_obs ( kstp, Nnn ) ! obs-minus-model (assimilation) diagnostics (call after dynamics update)346 341 347 342 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 348 343 ! Control 349 344 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 350 CALL stp_ctl ( kstp, Nbb, Nnn, indic ) 351 352 345 CALL stp_ctl ( kstp, Nnn ) 346 353 347 IF( kstp == nit000 ) THEN ! 1st time step only 354 348 CALL iom_close( numror ) ! close input ocean restart file … … 360 354 #if defined key_iomput 361 355 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 372 359 ! 373 360 IF( ln_timing ) CALL timing_stop('stp_RK3') 374 361 ! 375 362 END SUBROUTINE stp_RK3 376 ! 363 377 364 !!====================================================================== 378 365 END MODULE stpRK3 -
NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/src/SWE/stp_oce.F90
r13426 r13604 1 MODULE st ep_oce1 MODULE stp_oce 2 2 !!====================================================================== 3 !! *** MODULE st ep_oce ***3 !! *** MODULE stp_oce *** 4 4 !! Ocean time-stepping : module used in both initialisation phase and time stepping 5 !! (i.e. nemo_init and stp or stp_MLF routines) 5 6 !!====================================================================== 6 7 !! History : 3.3 ! 2010-08 (C. Ethe) Original code - reorganisation of the initial phase … … 9 10 USE oce ! ocean dynamics and tracers variables 10 11 USE dom_oce ! ocean space and time domain variables 11 USE zdf_oce ! ocean vertical physics variables12 USE zdfdrg , ONLY : ln_drgimp ! implicit top/bottom friction13 12 14 13 USE daymod ! calendar (day routine) … … 19 18 USE sbccpl ! surface boundary condition: coupled formulation (call send at end of step) 20 19 USE sbcapr ! surface boundary condition: atmospheric pressure 21 USE tide_mod, ONLY : ln_tide, tide_update22 20 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) 23 27 24 28 USE isf_oce ! ice shelf boundary condition 25 29 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) 26 45 27 46 USE traqsr ! solar radiation penetration (tra_qsr routine) … … 39 58 USE eosbn2 ! equation of state (eos_bn2 routine) 40 59 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 51 60 USE stopar ! Stochastic parametrization (sto_par routine) 52 61 USE stopts 53 54 USE bdy_oce , ONLY : ln_bdy55 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)64 62 65 63 USE ldfslp ! iso-neutral slopes (ldf_slp routine) … … 67 65 USE ldftra ! lateral eddy diffusive coef. (ldf_tra routine) 68 66 67 USE zdf_oce ! ocean vertical physics variables 69 68 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 71 71 72 72 USE diu_layers ! diurnal SST bulk and coolskin routines … … 81 81 USE diahth ! thermocline depth (dia_hth routine) 82 82 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) 85 85 USE diadetide ! Weights computation for daily detiding of model diagnostics 86 86 USE diamlr ! IOM context management for multiple-linear-regression analysis … … 92 92 USE asminc ! assimilation increments (tra_asm_inc routine) 93 93 ! (dyn_asm_inc routine) 94 USE asmbkg 94 USE asmbkg ! writing out state trajectory 95 95 USE stpctl ! time stepping control (stp_ctl routine) 96 96 USE restart ! ocean restart (rst_wri routine) … … 117 117 !! Software governed by the CeCILL license (see ./LICENSE) 118 118 !!====================================================================== 119 END MODULE st ep_oce119 END MODULE stp_oce -
NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/src/SWE/stpctl.F90
r12983 r13604 3 3 !! *** MODULE stpctl *** 4 4 !! Ocean run control : gross check of the ocean time stepping 5 !! *** Shallow Water Equation (SWE) case *** 6 !! ( No test on temperature and salinity ) 5 7 !!====================================================================== 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 13 9 !!---------------------------------------------------------------------- 14 10 … … 19 15 USE dom_oce ! ocean space and time domain variables 20 16 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 ! 21 20 USE diawri ! Standard run outputs (dia_wri_state routine) 22 !23 21 USE in_out_manager ! I/O manager 24 22 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 25 23 USE lib_mpp ! distributed memory computing 26 USE zdf_oce , ONLY : ln_zad_Aimp ! ocean vertical physics variables27 USE wet_dry, ONLY : ll_wd, ssh_ref ! reference depth for negative bathy28 29 24 USE netcdf ! NetCDF library 25 30 26 IMPLICIT NONE 31 27 PRIVATE … … 33 29 PUBLIC stp_ctl ! routine called by step.F90 34 30 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 38 35 # include "domzgr_substitute.h90" 39 36 !!---------------------------------------------------------------------- … … 44 41 CONTAINS 45 42 46 SUBROUTINE stp_ctl( kt, K bb, Kmm, kindic)43 SUBROUTINE stp_ctl( kt, Kmm ) 47 44 !!---------------------------------------------------------------------- 48 45 !! *** ROUTINE stp_ctl *** … … 52 49 !! ** Method : - Save the time step in numstp 53 50 !! - Print it each 50 time steps 54 !! - Stop the run IF problem encountered by setting indic=-355 !! Problems checked: |ssh| maximum larger than 10 m51 !! - Stop the run IF problem encountered by setting nstop > 0 52 !! Problems checked: e3t0+ssh minimum smaller that 0 56 53 !! |U| maximum larger than 10 m/s 57 !! negative sea surface salinity54 !! ( not for SWE : negative sea surface salinity ) 58 55 !! 59 56 !! ** Actions : "time.step" file = last ocean time-step 60 57 !! "run.stat" file = run statistics 61 !! nstop indicator sheared among all local domain (lk_mpp=T)58 !! nstop indicator sheared among all local domain 62 59 !!---------------------------------------------------------------------- 63 60 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 88 94 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 89 96 clname = 'run.stat.nc' 90 97 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 101 118 WRITE ( numstp, '(1x, i8)' ) kt 102 119 REWIND( numstp ) 103 120 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) 239 9500 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) 108 285 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 ! 292 9100 FORMAT('MPI rank ', a) 293 9200 FORMAT('found in ', a, ' MPI tasks, spread out among ranks ', a, ' to ', a) 294 9300 FORMAT('kt ', a, ' ', a, ' ', 1pg11.4, ' at i j ', a, ' ', a, ' ', a, ' ', a) 295 9400 FORMAT('kt ', a, ' ', a, ' ', 1pg11.4, ' at i j k ', a, ' ', a, ' ', a, ' ', a) 296 ! 297 END SUBROUTINE wrt_line 298 164 299 165 300 !!======================================================================
Note: See TracChangeset
for help on using the changeset viewer.
