Changeset 13005 for NEMO/branches/2020
- Timestamp:
- 2020-06-02T17:46:26+02:00 (4 years ago)
- Location:
- NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater
- Files:
-
- 12 edited
-
cfgs/AM98/EXPREF/namelist_cfg (modified) (6 diffs)
-
cfgs/AM98/MY_SRC/usrdef_hgr.F90 (modified) (3 diffs)
-
cfgs/AM98/MY_SRC/usrdef_nam.F90 (modified) (5 diffs)
-
cfgs/AM98/MY_SRC/usrdef_sbc.F90 (modified) (1 diff)
-
cfgs/AM98/MY_SRC/usrdef_zgr.F90 (modified) (9 diffs)
-
src/SWE/domvvl.F90 (modified) (1 diff)
-
src/SWE/dynadv.F90 (modified) (1 diff)
-
src/SWE/dynkeg.F90 (modified) (2 diffs)
-
src/SWE/dynldf_lap_blp.F90 (modified) (6 diffs)
-
src/SWE/dynvor.F90 (modified) (17 diffs)
-
src/SWE/ldfdyn.F90 (modified) (3 diffs)
-
src/SWE/step.F90 (modified) (3 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/cfgs/AM98/EXPREF/namelist_cfg
r12667 r13005 23 23 nn_it000 = 1 ! first time step 24 24 nn_itend = 1 ! 10 ans - 2 min - last time step 25 ! nn_itend = 2592000 ! 10 ans - 2 min - last time step26 25 ! nn_itend = 1036800 ! 10 ans - 5 min 27 ! nn_itend = 34560 ! 1 an - 15 min28 ! nn_itend = 172800 ! 10 ans - 30 min 26 ! nn_itend = 345600 ! 10 ans - 15 min - AM98 1/4deg + P.Marchand 27 ! nn_itend = 172800 ! 10 ans - 30 min - AM98 1/4deg 29 28 ! nn_itend = 34560 ! 10 ans - 2h 30 29 ! nn_itend = 86400 ! 10 ans - 1h … … 43 42 &namusr_def ! AM98 user defined namelist 44 43 !----------------------------------------------------------------------- 45 nn_AM98 = 4 ! AM98 resolution [1/ degrees]44 nn_AM98 = 4 ! AM98 resolution [1/nn_AM98] 46 45 jpkglo = 2 ! number of model levels 47 46 rn_theta = 0. ! rotation angle fo the grid [deg] 47 ! 48 nn_gc = 2 ! number of ghostcells 49 rn_domsiz = 2000000 ! size of the domain (default 2000km) [m] 50 rn_dx = 100000 ! gridspacing (default 100km) [m] 48 51 ! 49 52 rn_tau = 0.20 ! wind stress on the surface [N/m2] … … 53 56 rn_beta = 2.e-11 ! beta-plan coriolis parameter [1/m.s] 54 57 rn_f0 = 0.7e-4 ! f-plan coriolis parameter [1/s] 58 ! 59 nn_dynldf_lap_typ = 1 ! choose type of laplacian (ideally from namelist) 60 ! ! = 1 divrot laplacian 61 ! ! = 2 symmetric laplacian (Griffies&Hallberg 2000) 62 ! ! = 3 symmetric laplacian (cartesian) 63 ln_dynldf_lap_PM =.false. ! if True - apply the P.Marchand boundary condition on the laplacian viscosity 64 ! 55 65 / 56 66 !----------------------------------------------------------------------- … … 60 70 ! 61 71 ! rn_Dt = 7200. ! 2h - 1 deg cfg - time step for the dynamics [s] 62 ! rn_Dt = 3600. ! 1h - 1 deg cfg63 rn_Dt = 1800. ! 30min - 1/4 deg cfg64 ! rn_Dt = 900. ! 15min - 1/4 deg cfg72 ! rn_Dt = 3600. ! 1h - AM98 1 deg 73 rn_Dt = 1800. ! 30min - AM98 1/4deg 74 ! rn_Dt = 900. ! 15min - AM98 1/4deg + P.Marchand 65 75 ! rn_Dt = 600. ! 10min - 1/4 deg cfg 66 76 ! rn_Dt = 300. ! 5min - 1/4 deg cfg … … 109 119 &namlbc ! lateral momentum boundary condition (default: NO selection) 110 120 !----------------------------------------------------------------------- 111 rn_shlat = 0. ! free slip 121 ! ! free slip ! partial slip ! no slip ! strong slip 122 rn_shlat = 0. ! shlat = 0 ! 0 < shlat < 2 ! shlat = 2 ! 2 < shlat 112 123 / 113 124 !!====================================================================== … … 141 152 &namdyn_adv ! formulation of the momentum advection (default: NO selection) 142 153 !----------------------------------------------------------------------- 143 ln_dynadv_vec = .true. ! vector form - 2nd centered scheme 154 ln_dynadv_OFF = .false. ! linear dynamics (no momentum advection) 155 ln_dynadv_vec = .false. ! vector form - 2nd centered scheme 144 156 nn_dynkeg = 0 ! grad(KE) scheme: =0 C2 ; =1 Hollingsworth correction 157 ln_dynadv_cen2 = .false. ! flux form - 2nd order centered scheme 158 ln_dynadv_ubs = .false. ! flux form - 3rd order UBS scheme 145 159 / 146 160 !----------------------------------------------------------------------- 147 161 &namdyn_vor ! Vorticity / Coriolis scheme (default: NO selection) 148 162 !----------------------------------------------------------------------- 149 ln_dynvor_ene = .true. ! energy conserving scheme 150 ln_dynvor_ens = .false. ! enstrophy conserving scheme 151 ln_dynvor_een = .false. ! energy & enstrophy scheme 163 ln_dynvor_ene = .false. ! energy conserving scheme 164 ln_dynvor_ens = .false. ! enstrophy conserving scheme 165 ln_dynvor_enT = .false. ! energy conserving scheme (T-point) 166 ln_dynvor_eeT = .false. ! energy conserving scheme (een using e3t) 167 ln_dynvor_een = .false. ! energy & enstrophy scheme 152 168 nn_een_e3f = 0 ! =0 e3f = mi(mj(e3t))/4 153 169 ! ! =1 e3f = mi(mj(e3t))/mi(mj( tmask)) 154 ln_dynvor_msk = .false. ! vorticity multiplied by fmask (=T) 155 ! ! (f-point vorticity schemes only) 170 ln_dynvor_msk = .false. ! vorticity multiplied by fmask (=T) 171 ! ! (f-point vorticity schemes only) 172 ! 173 ln_dynvor_ens_adVO = .false. ! Alternative Direction 174 ln_dynvor_ens_adKE = .false. 175 ln_dynvor_ens_adKEVO = .false. 176 ln_dynvor_ene_adVO = .false. 177 ln_dynvor_ene_adKE = .false. 178 ln_dynvor_ene_adKEVO = .false. 156 179 / 157 180 !----------------------------------------------------------------------- -
NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/cfgs/AM98/MY_SRC/usrdef_hgr.F90
r12614 r13005 68 68 INTEGER :: ji, jj ! dummy loop indices 69 69 REAL(wp) :: zlam1, zlam0, zcos_theta, zim1 , zjm1 , ze1 , ze1deg ! local scalars 70 REAL(wp) :: zphi1, zphi0, zsin_theta, zim05, zjm05, znorme ! - - 71 REAL(wp) :: zl, zgl, zbl ! - - 70 REAL(wp) :: zphi1, zphi0, zsin_theta, zim05, zjm05, znorme ! - - 71 REAL(wp) :: zgl, zbl ! - - 72 72 73 !!------------------------------------------------------------------------------- 73 74 ! … … 81 82 ! !== grid point position ==! 82 83 ! 83 ze1 = 100000._wp / REAL(nn_AM98, wp) ! (100km) gridspacing [m] 84 zl = 2000000._wp ! lentgh side square domain [m] 85 zgl = 2000000._wp + 4._wp*ze1 ! length side square + ghost cells [m] 86 ! 2000km x 2000km + 2 cells around sides 87 ! 88 ! biggest L = 3500 km (worst case) 89 !zbl = 3500000._wp ! length side bigger domain [m] 90 ! non rotated (0deg) - worst case (centered in the 3500kmx3500km) 91 !zlam0 = 875000._wp ! m 92 !zphi0 = 875000._wp ! m 93 ! rotated case (45eg) 94 ! origin (in meters) of the 2000km x 2000km domain 95 !zlam0 = - 1000000._wp ! - 1000km 96 !zphi0 = 1474873.7341529163_wp ! 3500*sin(45deg) - 1000km 97 98 84 ze1 = rn_dx / REAL(nn_AM98, wp) ! [m] gridspacing used 85 zgl = rn_domsiz + 2._wp * REAL(nn_gc, wp) * ze1 ! [m] length of the square with ghostcells 99 86 ! fit the best square around the square + ghost cells 100 87 zbl = zgl * ( COS( rn_theta * rad ) + SIN( rn_theta * rad ) ) ! length side bigger domain [m] … … 111 98 112 99 ! exact origin in meters 113 zlam1 = zbl * COS((rn_theta + 45 )* rad ) / SQRT( 2._wp ) - zl/2._wp114 zphi1 = zbl * SIN((rn_theta + 45 )* rad ) / SQRT( 2._wp ) - zl/2._wp100 zlam1 = zbl * COS((rn_theta + 45 )* rad ) / SQRT( 2._wp ) - rn_domsiz/2._wp 101 zphi1 = zbl * SIN((rn_theta + 45 )* rad ) / SQRT( 2._wp ) - rn_domsiz/2._wp 115 102 ! origin put in the true corner of a cell so there will be no cropping 116 103 ! of the edge cells 117 zlam0 = REAL( anint( zlam1 / ze1 ) + 1, wp ) * ze1118 zphi0 = REAl( anint( zphi1 / ze1 ) + 1, wp ) * ze1104 zlam0 = REAL( anint( zlam1 / ze1 ), wp ) * ze1 105 zphi0 = REAl( anint( zphi1 / ze1 ), wp ) * ze1 119 106 120 107 IF(lwp) WRITE(numout,*) ' origin position zlam0 = ', zlam0/1000, ' km' -
NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/cfgs/AM98/MY_SRC/usrdef_nam.F90
r12614 r13005 28 28 INTEGER , PUBLIC :: nn_AM98 ! 1/nn_AM98 = the resolution chosen in degrees and thus defining the horizontal domain size 29 29 REAL(wp), PUBLIC :: rn_theta ! rotation angle (in degree) of the grid 30 REAL(wp), PUBLIC :: rn_tau ! wind stress on the surface [N/m2] 31 REAL(wp), PUBLIC :: rn_f0 ! f-plan coriolis parameter [1/s] 32 REAL(wp), PUBLIC :: rn_beta ! beta-plan coriolis parameter [1/m.s] 33 REAL(wp), PUBLIC :: rn_modified_grav ! modified gravity [m/s2] 34 REAL(wp), PUBLIC :: rn_rfr ! layer friction [1/s] 30 INTEGER , PUBLIC :: nn_gc ! number of ghostcells 31 REAL(wp), PUBLIC :: rn_domsiz ! size of the domain (default 2000km) [m] 32 REAL(wp), PUBLIC :: rn_dx ! gridspacing (default 100km) [m] 33 REAL(wp), PUBLIC :: rn_tau ! wind stress on the surface [N/m2] 34 REAL(wp), PUBLIC :: rn_f0 ! f-plan coriolis parameter [1/s] 35 REAL(wp), PUBLIC :: rn_beta ! beta-plan coriolis parameter [1/m.s] 36 REAL(wp), PUBLIC :: rn_modified_grav ! modified gravity [m/s2] 37 REAL(wp), PUBLIC :: rn_rfr ! layer friction [1/s] 38 ! !!* temporary *!! 39 INTEGER , PUBLIC :: nn_dynldf_lap_typ ! choose type of laplacian (ideally from namelist) 40 ! ! = 1 divrot laplacian 41 ! ! = 2 symmetric laplacian (Griffies&Hallberg 2000) 42 ! ! = 3 symmetric laplacian (cartesian) 43 LOGICAL , PUBLIC :: ln_dynldf_lap_PM ! if True - apply the P.Marchand boundary condition on the laplacian 35 44 ! 36 45 !!---------------------------------------------------------------------- … … 61 70 REAL(wp) :: ze1, zgl, zbl ! gridspacing, length of the biggest square 62 71 !! 63 NAMELIST/namusr_def/ nn_AM98, rn_theta, jpkglo, rn_tau, & ! 72 NAMELIST/namusr_def/ nn_AM98, rn_theta, jpkglo, & ! 73 & nn_gc ,rn_domsiz, rn_dx, & ! domain parameters 64 74 & rn_f0 ,rn_beta, & ! coriolis parameter 65 & rn_modified_grav, rn_rfr ! 75 & rn_modified_grav, rn_rfr , rn_tau, & ! reduced gravity, friction, wind 76 & nn_dynldf_lap_typ, ln_dynldf_lap_PM ! temporary parameter 66 77 !!---------------------------------------------------------------------- 67 78 ! … … 78 89 79 90 kk_cfg = nn_AM98 80 81 ! square domain of 2000km x 2000km with dx = dy = 100km 82 ! kpi = 20 * nn_AM98 + 2 ! Global Domain size 83 ! kpj = 20 * nn_AM98 + 2 84 85 ! Number of cells 86 ! lenght of the biggest square + 2 ghost cells on each side 87 ze1 = 100000._wp / REAL(nn_AM98, wp) ! (100km) gridspacing [m] 88 zgl = 2000000._wp + 4._wp*ze1 ! length side square + ghostcells [m] 89 !zbl = zgl / COS( rn_theta * rad ) ! length side bigger domain [m] 91 ! 92 ze1 = rn_dx / REAL(nn_AM98, wp) ! [m] gridspacing used 93 zgl = rn_domsiz + 2._wp * REAL(nn_gc, wp) * ze1 ! [m] length of the square with ghostcells 94 ! rotation 90 95 zbl = zgl * ( COS( rn_theta * rad ) + SIN( rn_theta * rad ) ) ! length side bigger domain [m] 91 96 ! 92 97 kpi = ceiling(zbl / ze1 ) ! Global Domain size + ghost cells 93 98 kpj = ceiling(zbl / ze1 ) ! Global Domain size + ghost cells … … 104 109 IF( rn_modified_grav /= 0._wp) grav = rn_modified_grav ! update the gravity 105 110 ! 111 ! !== Temporary namelist parameter ==! 112 ! 113 ! ll_dynldf_lap_PM = ln_dynldf_lap_PM 114 ! 106 115 kpk = jpkglo 107 116 ! ! Set the lateral boundary condition of the global domain … … 114 123 WRITE(numout,*) '~~~~~~~~~~~ ' 115 124 WRITE(numout,*) ' Namelist namusr_def : AM98 case' 125 WRITE(numout,*) ' domain size rn_domsiz = ', rn_domsiz, 'm' 126 WRITE(numout,*) ' gridspacing rn_dx = ', rn_dx, 'm' 116 127 WRITE(numout,*) ' inverse resolution & implied domain size nn_AM98 = ', nn_AM98 128 WRITE(numout,*) ' implied gridspacing rn_dx = ', rn_dx, 'm' 129 WRITE(numout,*) ' number of ghostcells nn_gc = ', nn_gc 130 WRITE(numout,*) ' ' 131 WRITE(numout,*) ' rotation angle chosen rn_theta = ', rn_theta, 'deg' 132 WRITE(numout,*) ' modified gravity rn_modified_grav = ', rn_modified_grav, 'm2/s' 117 133 #if defined key_agrif 118 134 IF( Agrif_Root() ) THEN 119 135 #endif 120 WRITE(numout,*) ' jpiglo = Nx*nn_AM98 +4jpiglo = ', kpi121 WRITE(numout,*) ' jpjglo = N y*nn_AM98+4jpjglo = ', kpj136 WRITE(numout,*) ' jpiglo = Nx*nn_AM98 + 2*n_ghost jpiglo = ', kpi 137 WRITE(numout,*) ' jpjglo = Nx*nn_AM98 + 2*n_ghost jpjglo = ', kpj 122 138 #if defined key_agrif 123 139 ENDIF -
NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/cfgs/AM98/MY_SRC/usrdef_sbc.F90
r12614 r13005 92 92 END DO 93 93 94 95 94 ! module of wind stress and wind speed at T-point 96 95 zcoef = 1. / ( zrhoa * zcdrag ) -
NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/cfgs/AM98/MY_SRC/usrdef_zgr.F90
r12614 r13005 3 3 !! *** MODULE usrdef_zgr *** 4 4 !! 5 !! === GYREconfiguration ===5 !! === AM98 configuration === 6 6 !! 7 7 !! User defined : vertical coordinate system of a user configuration … … 64 64 ! 65 65 IF(lwp) WRITE(numout,*) 66 IF(lwp) WRITE(numout,*) 'usr_def_zgr : GYREconfiguration (z-coordinate closed flat box ocean without cavities)'66 IF(lwp) WRITE(numout,*) 'usr_def_zgr : AM98 configuration (z-coordinate closed flat box ocean without cavities)' 67 67 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 68 68 ! … … 70 70 ! type of vertical coordinate 71 71 ! --------------------------- 72 ld_zco = .FALSE. ! GYREcase: z-coordinate without ocean cavities72 ld_zco = .FALSE. ! AM98 case: z-coordinate without ocean cavities 73 73 ld_zps = .FALSE. 74 74 ld_sco = .TRUE. … … 170 170 !! ** Purpose : set the masked top and bottom ocean t-levels 171 171 !! 172 !! ** Method : GYREcase = closed flat box ocean without ocean cavities172 !! ** Method : AM98 case = closed flat box ocean without ocean cavities 173 173 !! k_top = 1 except along north, south, east and west boundaries 174 174 !! k_bot = jpk-1 except along north, south, east and west boundaries … … 182 182 INTEGER :: ji, jj ! dummy loop indices 183 183 REAL(wp) :: zylim0, zylim1, zxlim0, zxlim1 ! limit of the domain [m] 184 REAL(WP) :: zcoeff ! local scalar 184 185 !!---------------------------------------------------------------------- 185 186 ! … … 187 188 IF(lwp) WRITE(numout,*) ' zgr_top_bot : defines the top and bottom wet ocean levels.' 188 189 IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~' 189 IF(lwp) WRITE(numout,*) ' GYREcase : closed flat box ocean without ocean cavities'190 IF(lwp) WRITE(numout,*) ' AM98 case : closed flat box ocean without ocean cavities' 190 191 ! 191 192 z2d(:,:) = REAL( jpkm1 , wp ) ! flat bottom … … 193 194 CALL lbc_lnk( 'usrdef_zgr', z2d, 'T', 1. ) ! set surrounding land to zero (here jperio=0 ==>> closed) 194 195 ! 195 196 196 ! rotated domain (45deg) 197 197 !zylim0 = 100000._wp * SQRT( 2._wp ) / 3 ! dy*sqrt(2)/2 * 2/3 … … 214 214 DO jj = 1, jpj 215 215 DO ji = 1, jpi 216 217 216 ! if T point in the 2000 km x 2000 km domain 218 217 ! IF ( gphit(ji,jj) > zylim0 .AND. gphit(ji,jj) < zylim1 .AND. & … … 227 226 k_bot(ji,jj) = 0 228 227 END IF 229 230 228 END DO 231 229 END DO 232 230 ! mask the lonely corners 231 DO jj = 1, jpj 232 DO ji = 1, jpi 233 zcoeff = k_top(ji+1,jj) + k_top(ji,jj+1) & 234 + k_top(ji-1,jj) + k_top(ji,jj-1) 235 IF ( zcoeff <= 1._wp ) THEN 236 k_top(ji,jj) = 0 ! = land 237 k_bot(ji,jj) = 0 238 END IF 239 END DO 240 END DO 233 241 ! k_bot(:,:) = NINT( z2d(:,:) ) ! =jpkm1 over the ocean point, =0 elsewhere 234 242 ! -
NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/domvvl.F90
r12614 r13005 723 723 & + e1e2t(ji+1,jj+1) * pssh(ji+1,jj+1) ) * r1_hf_0(ji,jj) * r1_e1e2f(ji,jj) 724 724 END_2D 725 !!an dans le cas tourné, hf augmente et trend VOR diminue 726 ! DO_2D_10_10 727 ! zc3(ji,jj) = ( e1e2t(ji ,jj ) * pssh(ji ,jj ) & 728 ! & + e1e2t(ji+1,jj ) * pssh(ji+1,jj ) & 729 ! & + e1e2t(ji ,jj+1) * pssh(ji ,jj+1) & 730 ! & + e1e2t(ji+1,jj+1) * pssh(ji+1,jj+1) ) * r1_hf_0(ji,jj) * r1_e1e2f(ji,jj) & 731 ! & / MAX( tmask(ji,jj) + tmask(ji+1,jj) + tmask(ji,jj+1) + tmask(ji+1,jj+1), 1._wp ) 732 ! END_2D 733 725 734 CALL lbc_lnk( 'domvvl', zc3(:,:), 'F', 1._wp ) 726 735 ! -
NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/dynadv.F90
r12822 r13005 74 74 CASE( np_VEC_c2 ) 75 75 CALL dyn_keg ( kt, nn_dynkeg, Kmm, puu, pvv, Krhs ) ! vector form : horizontal gradient of kinetic energy 76 !!an SWE : w = 0 76 77 CALL dyn_zad ( kt, Kmm, puu, pvv, Krhs ) ! vector form : vertical advection 77 78 CASE( np_FLX_c2 ) -
NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/dynkeg.F90
r12822 r13005 29 29 30 30 PUBLIC dyn_keg ! routine called by step module 31 PUBLIC dyn_kegAD ! routine called by step module 31 32 32 33 INTEGER, PARAMETER, PUBLIC :: nkeg_C2 = 0 !: 2nd order centered scheme (standard scheme) … … 155 156 ! 156 157 END SUBROUTINE dyn_keg 157 158 159 160 SUBROUTINE dyn_kegAD( kt, kscheme, puu, pvv, pu_rhs, pv_rhs ) 161 !!---------------------------------------------------------------------- 162 !! *** ROUTINE dyn_kegAD *** 163 !! 164 !! ** Purpose : Compute the now momentum trend due to the horizontal 165 !! gradient of the horizontal kinetic energy and add it to the 166 !! general momentum trend. 167 !! 168 !! ** Method : * kscheme = nkeg_C2 : 2nd order centered scheme that 169 !! conserve kinetic energy. Compute the now horizontal kinetic energy 170 !! zhke = 1/2 [ mi-1( un^2 ) + mj-1( vn^2 ) ] 171 !! * kscheme = nkeg_HW : Hollingsworth correction following 172 !! Arakawa (2001). The now horizontal kinetic energy is given by: 173 !! zhke = 1/6 [ mi-1( 2 * un^2 + ((u(j+1)+u(j-1))/2)^2 ) 174 !! + mj-1( 2 * vn^2 + ((v(i+1)+v(i-1))/2)^2 ) ] 175 !! 176 !! Take its horizontal gradient and add it to the general momentum 177 !! trend. 178 !! u(rhs) = u(rhs) - 1/e1u di[ zhke ] 179 !! v(rhs) = v(rhs) - 1/e2v dj[ zhke ] 180 !! 181 !! ** Action : - Update the (puu(:,:,:,Krhs), pvv(:,:,:,Krhs)) with the hor. ke gradient trend 182 !! - send this trends to trd_dyn (l_trddyn=T) for post-processing 183 !! 184 !! ** References : Arakawa, A., International Geophysics 2001. 185 !! Hollingsworth et al., Quart. J. Roy. Meteor. Soc., 1983. 186 !!---------------------------------------------------------------------- 187 ! 188 INTEGER , INTENT( in ) :: kt ! ocean time-step index 189 INTEGER , INTENT( in ) :: kscheme ! =0/1/2 type of KEG scheme 190 REAL(wp), DIMENSION(jpi,jpj,jpk) , INTENT(inout) :: puu, pvv ! ocean velocities at Kmm 191 REAL(wp), DIMENSION(jpi,jpj,jpk),OPTIONAL, INTENT(inout) :: pu_rhs, pv_rhs ! RHS 192 ! 193 INTEGER :: ji, jj, jk ! dummy loop indices 194 REAL(wp) :: zu, zv ! local scalars 195 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zhke 196 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: ztrdu, ztrdv 197 !!---------------------------------------------------------------------- 198 ! 199 IF( ln_timing ) CALL timing_start('dyn_keg') 200 ! 201 IF( kt == nit000 ) THEN 202 IF(lwp) WRITE(numout,*) 203 IF(lwp) WRITE(numout,*) 'dyn_kegAD : kinetic energy gradient trend, scheme number=', kscheme 204 IF(lwp) WRITE(numout,*) '~~~~~~~~~' 205 ENDIF 206 207 zhke(:,:,jpk) = 0._wp 208 209 SELECT CASE ( kscheme ) !== Horizontal kinetic energy at T-point ==! 210 !!an45 to be ADDED : que cas C2 - "wet points only" il suffit de x2 le terme quadratic a la coast (nn_dynkeg_adv = 2) 211 CASE ( nkeg_C2_wpo ) !-- Standard scheme --! 212 DO_3D_01_01( 1, jpkm1 ) 213 zu = ( puu(ji-1,jj ,jk) * puu(ji-1,jj ,jk) & 214 & + puu(ji ,jj ,jk) * puu(ji ,jj ,jk) ) * ( 2._wp - umask(ji-1,jj,jk) * umask(ji,jj,jk) ) 215 zv = ( pvv(ji ,jj-1,jk) * pvv(ji ,jj-1,jk) & 216 & + pvv(ji ,jj ,jk) * pvv(ji ,jj ,jk) ) * ( 2._wp - vmask(ji,jj-1,jk) * vmask(ji,jj,jk) ) 217 zhke(ji,jj,jk) = 0.25_wp * ( zv + zu ) 218 END_3D 219 !!an45 220 ! 221 CASE ( nkeg_C2 ) !-- Standard scheme --! 222 DO_3D_01_01( 1, jpkm1 ) 223 zu = puu(ji-1,jj ,jk) * puu(ji-1,jj ,jk) & 224 & + puu(ji ,jj ,jk) * puu(ji ,jj ,jk) 225 zv = pvv(ji ,jj-1,jk) * pvv(ji ,jj-1,jk) & 226 & + pvv(ji ,jj ,jk) * pvv(ji ,jj ,jk) 227 zhke(ji,jj,jk) = 0.25_wp * ( zv + zu ) 228 END_3D 229 !!an 00_00 ? 230 CASE ( nkeg_HW ) !-- Hollingsworth scheme --! 231 DO_3D_00_00( 1, jpkm1 ) 232 zu = 8._wp * ( puu(ji-1,jj ,jk) * puu(ji-1,jj ,jk) & 233 & + puu(ji ,jj ,jk) * puu(ji ,jj ,jk) ) & 234 & + ( puu(ji-1,jj-1,jk) + puu(ji-1,jj+1,jk) ) * ( puu(ji-1,jj-1,jk) + puu(ji-1,jj+1,jk) ) & 235 & + ( puu(ji ,jj-1,jk) + puu(ji ,jj+1,jk) ) * ( puu(ji ,jj-1,jk) + puu(ji ,jj+1,jk) ) 236 ! 237 zv = 8._wp * ( pvv(ji ,jj-1,jk) * pvv(ji ,jj-1,jk) & 238 & + pvv(ji ,jj ,jk) * pvv(ji ,jj ,jk) ) & 239 & + ( pvv(ji-1,jj-1,jk) + pvv(ji+1,jj-1,jk) ) * ( pvv(ji-1,jj-1,jk) + pvv(ji+1,jj-1,jk) ) & 240 & + ( pvv(ji-1,jj ,jk) + pvv(ji+1,jj ,jk) ) * ( pvv(ji-1,jj ,jk) + pvv(ji+1,jj ,jk) ) 241 zhke(ji,jj,jk) = r1_48 * ( zv + zu ) 242 END_3D 243 CALL lbc_lnk( 'dynkeg', zhke, 'T', 1. ) 244 ! 245 END SELECT 246 ! 247 IF( PRESENT( pu_rhs ) .AND. PRESENT( pv_rhs ) ) THEN !*** NO alternating direction ***! 248 ! 249 DO_3D_00_00( 1, jpkm1 ) 250 pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) - ( zhke(ji+1,jj ,jk) - zhke(ji,jj,jk) ) / e1u(ji,jj) 251 pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) - ( zhke(ji ,jj+1,jk) - zhke(ji,jj,jk) ) / e2v(ji,jj) 252 END_3D 253 ! 254 ELSEIF( PRESENT( pu_rhs ) .AND. .NOT. PRESENT( pv_rhs ) ) THEN !*** Alternating direction : i-component ***! 255 ! 256 DO_3D_00_00( 1, jpkm1 ) 257 pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) - ( zhke(ji+1,jj ,jk) - zhke(ji,jj,jk) ) / e1u(ji,jj) 258 END_3D 259 ! 260 ELSEIF( .NOT. PRESENT( pu_rhs ) .AND. PRESENT( pv_rhs ) ) THEN !*** Alternating direction : j-component ***! 261 ! 262 DO_3D_00_00( 1, jpkm1 ) 263 pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) - ( zhke(ji ,jj+1,jk) - zhke(ji,jj,jk) ) / e2v(ji,jj) 264 END_3D 265 ! 266 ENDIF 267 IF( ln_timing ) CALL timing_stop('dyn_kegAD') 268 ! 269 END SUBROUTINE dyn_kegAD 158 270 !!====================================================================== 159 271 END MODULE dynkeg -
NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/dynldf_lap_blp.F90
r12822 r13005 20 20 USE in_out_manager ! I/O manager 21 21 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 22 22 USE lib_mpp ! MPP library 23 ! 24 USE usrdef_nam , ONLY : nn_dynldf_lap_typ ! use laplacian parameter 25 ! 23 26 IMPLICIT NONE 24 27 PRIVATE … … 31 34 INTEGER, PUBLIC, PARAMETER :: np_dynldf_lap_symN = 3 ! symmetric laplacian (cartesian) 32 35 33 INTEGER, PUBLIC, PARAMETER :: ln_dynldf_lap_typ = 1 ! choose type of laplacian (ideally from namelist)36 !INTEGER, PUBLIC, PARAMETER :: nn_dynldf_lap_typ = 1 ! choose type of laplacian (ideally from namelist) 34 37 !!anSYM 35 38 !! * Substitutions … … 77 80 WRITE(numout,*) 'dyn_ldf : iso-level harmonic (laplacian) operator, pass=', kpass 78 81 WRITE(numout,*) '~~~~~~~ ' 79 WRITE(numout,*) ' ln_dynldf_lap_typ = ', ln_dynldf_lap_typ80 SELECT CASE( ln_dynldf_lap_typ ) ! print the choice of operator82 WRITE(numout,*) ' nn_dynldf_lap_typ = ', nn_dynldf_lap_typ 83 SELECT CASE( nn_dynldf_lap_typ ) ! print the choice of operator 81 84 CASE( np_dynldf_lap_rot ) ; WRITE(numout,*) ' ==>>> div-rot laplacian' 82 85 CASE( np_dynldf_lap_sym ) ; WRITE(numout,*) ' ==>>> symmetric laplacian (covariant form)' 83 CASE( np_dynldf_lap_symN) ; WRITE(numout,*) ' ==>>> symmetric laplacian ( simpleform)'86 CASE( np_dynldf_lap_symN) ; WRITE(numout,*) ' ==>>> symmetric laplacian (cartesian form)' 84 87 END SELECT 85 88 ENDIF … … 89 92 ENDIF 90 93 ! 91 SELECT CASE( ln_dynldf_lap_typ )94 SELECT CASE( nn_dynldf_lap_typ ) 92 95 ! 93 96 CASE ( np_dynldf_lap_rot ) !== Vorticity-Divergence form ==! … … 157 160 END DO ! End of slab 158 161 ! 159 CASE ( np_dynldf_lap_symN ) !== Symmetric form ==! ( naiveway)162 CASE ( np_dynldf_lap_symN ) !== Symmetric form ==! (cartesian way) 160 163 ! 161 164 DO jk = 1, jpkm1 ! Horizontal slab … … 190 193 ! 191 194 CASE DEFAULT ! error 192 CALL ctl_stop('STOP','dyn_ldf_lap: wrong value for ln_dynldf_lap_typ' )195 CALL ctl_stop('STOP','dyn_ldf_lap: wrong value for nn_dynldf_lap_typ' ) 193 196 END SELECT 194 197 ! -
NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/dynvor.F90
r12667 r13005 22 22 !! - ! 2018-04 (G. Madec) add pre-computed gradient for metric term calculation 23 23 !! 4.x ! 2020-03 (G. Madec, A. Nasser) make ln_dynvor_msk truly efficient on relative vorticity 24 !! 4.x ! 2020-03 (G. Madec, A. Nasser) alternate direction computation of vorticity tendancy 25 !! ! for ENS, ENE 24 26 !!---------------------------------------------------------------------- 25 27 … … 45 47 USE lib_mpp ! MPP library 46 48 USE timing ! Timing 49 !!anAD only 50 USE dynkeg, ONLY : dyn_kegAD 51 USE dynadv, ONLY : nn_dynkeg 47 52 48 53 IMPLICIT NONE … … 53 58 54 59 ! !!* Namelist namdyn_vor: vorticity term 55 LOGICAL, PUBLIC :: ln_dynvor_ens !: enstrophy conserving scheme (ENS) 56 LOGICAL, PUBLIC :: ln_dynvor_ene !: f-point energy conserving scheme (ENE) 57 LOGICAL, PUBLIC :: ln_dynvor_enT !: t-point energy conserving scheme (ENT) 58 LOGICAL, PUBLIC :: ln_dynvor_eeT !: t-point energy conserving scheme (EET) 59 LOGICAL, PUBLIC :: ln_dynvor_een !: energy & enstrophy conserving scheme (EEN) 60 INTEGER, PUBLIC :: nn_een_e3f !: e3f=masked averaging of e3t divided by 4 (=0) or by the sum of mask (=1) 61 LOGICAL, PUBLIC :: ln_dynvor_mix !: mixed scheme (MIX) 62 LOGICAL, PUBLIC :: ln_dynvor_msk !: vorticity multiplied by fmask (=T) or not (=F) (all vorticity schemes) 60 LOGICAL, PUBLIC :: ln_dynvor_ens !: enstrophy conserving scheme (ENS) 61 LOGICAL, PUBLIC :: ln_dynvor_ens_adVO = .FALSE. !: AD enstrophy conserving scheme (ENS_ad) 62 LOGICAL, PUBLIC :: ln_dynvor_ens_adKE = .FALSE. !: AD enstrophy conserving scheme (ENS_ad) 63 LOGICAL, PUBLIC :: ln_dynvor_ens_adKEVO = .FALSE. !: AD enstrophy conserving scheme (ENS_ad) 64 LOGICAL, PUBLIC :: ln_dynvor_ene !: f-point energy conserving scheme (ENE) 65 LOGICAL, PUBLIC :: ln_dynvor_ene_adVO = .FALSE. !: f-point AD energy conserving scheme (ENE_ad) 66 LOGICAL, PUBLIC :: ln_dynvor_ene_adKE = .FALSE. !: f-point AD energy conserving scheme (ENE_ad) 67 LOGICAL, PUBLIC :: ln_dynvor_ene_adKEVO = .FALSE. !: f-point AD energy conserving scheme (ENE_ad) 68 LOGICAL, PUBLIC :: ln_dynvor_enT !: t-point energy conserving scheme (ENT) 69 LOGICAL, PUBLIC :: ln_dynvor_eeT !: t-point energy conserving scheme (EET) 70 LOGICAL, PUBLIC :: ln_dynvor_een !: energy & enstrophy conserving scheme (EEN) 71 INTEGER, PUBLIC :: nn_een_e3f !: e3f=masked averaging of e3t divided by 4 (=0) or by the sum of mask (=1) 72 LOGICAL, PUBLIC :: ln_dynvor_mix !: mixed scheme (MIX) 73 LOGICAL, PUBLIC :: ln_dynvor_msk !: vorticity multiplied by fmask (=T) or not (=F) (all vorticity schemes) 63 74 64 75 INTEGER, PUBLIC :: nvor_scheme !: choice of the type of advection scheme … … 70 81 INTEGER, PUBLIC, PARAMETER :: np_EEN = 4 ! EEN scheme 71 82 INTEGER, PUBLIC, PARAMETER :: np_MIX = 5 ! MIX scheme 72 73 INTEGER :: ncor, nrvm, ntot ! choice of calculated vorticity 83 !!an 84 INTEGER, PUBLIC, PARAMETER :: np_ENS_adKE = 11 ! ENS scheme - AD scheme (KE only) 85 INTEGER, PUBLIC, PARAMETER :: np_ENS_adVO = 12 ! ENS scheme - AD scheme (VOR only) 86 INTEGER, PUBLIC, PARAMETER :: np_ENS_adKEVO = 13 ! ENS scheme - AD scheme (KE+VOR) 87 INTEGER, PUBLIC, PARAMETER :: np_ENE_adKE = 21 ! ENE scheme - AD scheme (KE only) 88 INTEGER, PUBLIC, PARAMETER :: np_ENE_adVO = 22 ! ENE scheme - AD scheme (VOR only) 89 INTEGER, PUBLIC, PARAMETER :: np_ENE_adKEVO = 23 ! ENE scheme - AD scheme (KE+VOR) 90 !!an 91 92 !!an ds step on pourra spécifier la valeur de ntot = np_COR ou np_COR + np_RVO 93 INTEGER, PUBLIC :: ncor, nrvm, ntot ! choice of calculated vorticity 74 94 ! ! associated indices: 75 95 INTEGER, PUBLIC, PARAMETER :: np_COR = 1 ! Coriolis (planetary) … … 97 117 CONTAINS 98 118 99 SUBROUTINE dyn_vor( kt, K mm, puu, pvv, Krhs)119 SUBROUTINE dyn_vor( kt, Kbb, Kmm, puu, pvv, Krhs) 100 120 !!---------------------------------------------------------------------- 101 121 !! … … 107 127 !! for futher diagnostics (l_trddyn=T) 108 128 !!---------------------------------------------------------------------- 109 INTEGER , INTENT( in ) :: kt ! ocean time-step index 110 INTEGER , INTENT( in ) :: Kmm, Krhs ! ocean time level indices 111 REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: puu, pvv ! ocean velocity field and RHS of momentum equation 129 INTEGER :: ji, jj, jk ! dummy loop indice 130 INTEGER , INTENT( in ) :: kt ! ocean time-step index 131 INTEGER , INTENT( in ) :: Kmm, Krhs, Kbb ! ocean time level indices 132 REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: puu, pvv ! ocean velocity field and RHS of momentum equation 112 133 ! 113 134 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: ztrdu, ztrdv 135 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zuu, zvv 114 136 !!---------------------------------------------------------------------- 115 137 ! … … 165 187 IF( ln_stcor ) CALL vor_eeT( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend 166 188 CASE( np_ENE ) !* energy conserving scheme 189 CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pu_rhs=puu(:,:,:,Krhs), pv_rhs=pvv(:,:,:,Krhs) ) 190 CALL vor_ene( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! total vorticity trend 191 IF( ln_stcor ) CALL vor_ene( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend 192 193 CASE( np_ENE_adVO ) !* energy conserving scheme with alternating direction 194 IF( MOD( kt , 2 ) == 1 ) THEN ! even time step: u-vor then v-vor components 195 196 !== Alternative direction - VOR only ==! 197 198 CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pu_rhs=puu(:,:,:,Krhs), pv_rhs=pvv(:,:,:,Krhs) ) ! compute and add uu-vorticity trend 199 200 ALLOCATE( zuu(jpi,jpj,jpk) ) 201 CALL vor_ene( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pu_rhs=puu(:,:,:,Krhs) ) ! compute and add uu-vorticity trend 202 zuu(:,:,:) = puu(:,:,:,Kbb) + rDt * puu(:,:,:,Krhs) * umask(:,:,:) 203 CALL lbc_lnk( 'dynvor', zuu(:,:,:) , 'U', -1._wp ) 204 CALL vor_ene( kt, Kmm, ntot, zuu(:,:,:) , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) ) ! compute and add vv-vorticity trend 205 DEALLOCATE( zuu ) 206 ELSE 207 CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pu_rhs=puu(:,:,:,Krhs), pv_rhs=pvv(:,:,:,Krhs) ) ! compute and add uu-vorticity trend 208 209 ALLOCATE( zvv(jpi,jpj,jpk) ) 210 CALL vor_ene( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) ) ! compute and add vv-vorticity trend 211 zvv(:,:,:) = pvv(:,:,:,Kbb) + rDt * pvv(:,:,:,Krhs) * vmask(:,:,:) 212 CALL lbc_lnk( 'dynvor', zvv(:,:,:) , 'V', -1._wp ) 213 CALL vor_ene( kt, Kmm, ntot, puu(:,:,:,Kmm), zvv(:,:,:) , pu_rhs=puu(:,:,:,Krhs) ) ! compute and add uu-vorticity trend 214 DEALLOCATE( zvv ) 215 ENDIF 216 IF( ln_stcor ) CALL vor_ene( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend 217 CASE( np_ENE_adKE ) !* energy conserving scheme with alternating direction 218 IF( MOD( kt , 2 ) == 1 ) THEN ! even time step: u-vor then v-vor components 219 220 !== Alternative direction - KEG only ==! 221 167 222 CALL vor_ene( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! total vorticity trend 223 224 ALLOCATE( zuu(jpi,jpj,jpk) ) 225 CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pu_rhs=puu(:,:,:,Krhs) ) ! compute and add uu-vorticity trend 226 zuu(:,:,:) = puu(:,:,:,Kbb) + rDt * puu(:,:,:,Krhs) * umask(:,:,:) 227 CALL lbc_lnk( 'dynvor', zuu(:,:,:) , 'U', -1._wp ) 228 CALL dyn_kegAD( kt, nn_dynkeg, zuu(:,:,:) , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) ) ! compute and add vv-vorticity trend 229 DEALLOCATE( zuu ) 230 ELSE 231 232 CALL vor_ene( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! total vorticity trend 233 234 ALLOCATE( zvv(jpi,jpj,jpk) ) 235 CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) ) ! compute and add vv-vorticity trend 236 zvv(:,:,:) = pvv(:,:,:,Kbb) + rDt * pvv(:,:,:,Krhs) * vmask(:,:,:) 237 CALL lbc_lnk( 'dynvor', zvv(:,:,:) , 'V', -1._wp ) 238 CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , zvv(:,:,:) , pu_rhs=puu(:,:,:,Krhs) ) ! compute and add uu-vorticity trend 239 DEALLOCATE( zvv ) 240 ENDIF 168 241 IF( ln_stcor ) CALL vor_ene( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend 242 243 CASE( np_ENE_adKEVO ) !* energy conserving scheme with alternating direction 244 IF( MOD( kt , 2 ) == 1 ) THEN ! even time step: u-vor then v-vor components 245 246 !== Alternative direction - KE + VOR ==! 247 248 ALLOCATE( zuu(jpi,jpj,jpk) ) 249 CALL vor_ene( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pu_rhs=puu(:,:,:,Krhs) ) ! compute and add uu-vorticity trend 250 CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pu_rhs=puu(:,:,:,Krhs) ) ! 251 zuu(:,:,:) = puu(:,:,:,Kbb) + rDt * puu(:,:,:,Krhs) * umask(:,:,:) 252 CALL lbc_lnk( 'dynvor', zuu(:,:,:) , 'U', -1._wp ) 253 CALL vor_ene( kt, Kmm, ntot, zuu(:,:,:) , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) ) ! compute and add vv-vorticity trend 254 CALL dyn_kegAD( kt, nn_dynkeg, zuu(:,:,:) , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) ) 255 DEALLOCATE( zuu ) 256 ELSE 257 258 ALLOCATE( zvv(jpi,jpj,jpk) ) 259 CALL vor_ene( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) ) ! compute and add vv-vorticity trend 260 CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) ) 261 zvv(:,:,:) = pvv(:,:,:,Kbb) + rDt * pvv(:,:,:,Krhs) * vmask(:,:,:) 262 CALL lbc_lnk( 'dynvor', zvv(:,:,:) , 'V', -1._wp ) 263 CALL vor_ene( kt, Kmm, ntot, puu(:,:,:,Kmm), zvv(:,:,:) , pu_rhs=puu(:,:,:,Krhs) ) ! compute and add uu-vorticity trend 264 CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm), zvv(:,:,:) , pu_rhs=puu(:,:,:,Krhs) ) 265 DEALLOCATE( zvv ) 266 ENDIF 169 267 CASE( np_ENS ) !* enstrophy conserving scheme 268 CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pu_rhs=puu(:,:,:,Krhs), pv_rhs=pvv(:,:,:,Krhs) ) 269 CALL vor_ens( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! total vorticity trend 270 IF( ln_stcor ) CALL vor_ens( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend 271 CASE( np_ENS_adVO ) !* enstrophy conserving scheme with alternating direction 272 IF( MOD( kt , 2 ) == 1 ) THEN ! even time step: u-vor then v-vor components 273 274 !== Alternative direction - VOR only ==! 275 276 CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pu_rhs=puu(:,:,:,Krhs), pv_rhs=pvv(:,:,:,Krhs) ) 277 278 ALLOCATE( zuu(jpi,jpj,jpk) ) 279 CALL vor_ens( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pu_rhs=puu(:,:,:,Krhs) ) ! compute and add uu-vorticity trend 280 zuu(:,:,:) = puu(:,:,:,Kbb) + rDt * puu(:,:,:,Krhs) * umask(:,:,:) 281 CALL lbc_lnk( 'dynvor', zuu(:,:,:) , 'U', -1._wp ) 282 CALL vor_ens( kt, Kmm, ntot, zuu(:,:,:) , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) ) ! compute and add vv-vorticity trend 283 DEALLOCATE( zuu ) 284 ELSE 285 CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pu_rhs=puu(:,:,:,Krhs), pv_rhs=pvv(:,:,:,Krhs) ) 286 287 ALLOCATE( zvv(jpi,jpj,jpk) ) 288 CALL vor_ens( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) ) ! compute and add vv-vorticity trend 289 zvv(:,:,:) = pvv(:,:,:,Kbb) + rDt * pvv(:,:,:,Krhs) * vmask(:,:,:) 290 CALL lbc_lnk( 'dynvor', zvv(:,:,:) , 'V', -1._wp ) 291 CALL vor_ens( kt, Kmm, ntot, puu(:,:,:,Kmm) , zvv(:,:,:), pu_rhs=puu(:,:,:,Krhs) ) ! compute and add uu-vorticity trend 292 DEALLOCATE( zvv ) 293 ENDIF 294 IF( ln_stcor ) CALL vor_ens( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend 295 CASE( np_ENS_adKE ) !* enstrophy conserving scheme with alternating direction 296 IF( MOD( kt , 2 ) == 1 ) THEN ! even time step: u-vor then v-vor components 297 298 !== Alternative direction - KEG only ==! 299 170 300 CALL vor_ens( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! total vorticity trend 301 302 ALLOCATE( zuu(jpi,jpj,jpk) ) 303 CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pu_rhs=puu(:,:,:,Krhs) ) ! compute and add uu-vorticity trend 304 zuu(:,:,:) = puu(:,:,:,Kbb) + rDt * puu(:,:,:,Krhs) * umask(:,:,:) 305 CALL lbc_lnk( 'dynvor', zuu(:,:,:) , 'U', -1._wp ) 306 CALL dyn_kegAD( kt, nn_dynkeg, zuu(:,:,:) , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) ) ! compute and add vv-vorticity trend 307 DEALLOCATE( zuu ) 308 ELSE 309 310 CALL vor_ens( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! total vorticity trend 311 312 ALLOCATE( zvv(jpi,jpj,jpk) ) 313 CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) ) ! compute and add vv-vorticity trend 314 zvv(:,:,:) = pvv(:,:,:,Kbb) + rDt * pvv(:,:,:,Krhs) * vmask(:,:,:) 315 CALL lbc_lnk( 'dynvor', zvv(:,:,:) , 'V', -1._wp ) 316 CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , zvv(:,:,:) , pu_rhs=puu(:,:,:,Krhs) ) ! compute and add uu-vorticity trend 317 DEALLOCATE( zvv ) 318 ENDIF 319 CASE( np_ENS_adKEVO ) !* enstrophy conserving scheme with alternating direction 320 IF( MOD( kt , 2 ) == 1 ) THEN ! even time step: u-vor then v-vor components 321 322 !== Alternative direction - KE + VOR ==! 323 324 ALLOCATE( zuu(jpi,jpj,jpk) ) 325 CALL vor_ens( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pu_rhs=puu(:,:,:,Krhs) ) ! compute and add uu-vorticity trend 326 CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pu_rhs=puu(:,:,:,Krhs) ) ! 327 zuu(:,:,:) = puu(:,:,:,Kbb) + rDt * puu(:,:,:,Krhs) * umask(:,:,:) 328 CALL lbc_lnk( 'dynvor', zuu(:,:,:) , 'U', -1._wp ) 329 CALL vor_ens( kt, Kmm, ntot, zuu(:,:,:) , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) ) ! compute and add vv-vorticity trend 330 CALL dyn_kegAD( kt, nn_dynkeg, zuu(:,:,:) , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) ) 331 DEALLOCATE( zuu ) 332 ELSE 333 334 ALLOCATE( zvv(jpi,jpj,jpk) ) 335 CALL vor_ens( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) ) ! compute and add vv-vorticity trend 336 CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) ) 337 zvv(:,:,:) = pvv(:,:,:,Kbb) + rDt * pvv(:,:,:,Krhs) * vmask(:,:,:) 338 CALL lbc_lnk( 'dynvor', zvv(:,:,:) , 'V', -1._wp ) 339 CALL vor_ens( kt, Kmm, ntot, puu(:,:,:,Kmm), zvv(:,:,:) , pu_rhs=puu(:,:,:,Krhs) ) ! compute and add uu-vorticity trend 340 CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm), zvv(:,:,:) , pu_rhs=puu(:,:,:,Krhs) ) 341 DEALLOCATE( zvv ) 342 ENDIF 171 343 IF( ln_stcor ) CALL vor_ens( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend 172 344 CASE( np_MIX ) !* mixed ene-ens scheme … … 313 485 !! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689. 314 486 !!---------------------------------------------------------------------- 315 INTEGER , INTENT(in ) :: kt ! ocean time-step index316 INTEGER , INTENT(in ) :: Kmm ! ocean time level index317 INTEGER , INTENT(in ) :: kvor ! total, planetary, relative, or metric318 REAL(wp), DIMENSION(jpi,jpj,jpk) , INTENT(inout) :: pu, pv ! now velocities319 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pu_rhs, pv_rhs ! total v-trend487 INTEGER , INTENT(in ) :: kt ! ocean time-step index 488 INTEGER , INTENT(in ) :: Kmm ! ocean time level index 489 INTEGER , INTENT(in ) :: kvor ! total, planetary, relative, or metric 490 REAL(wp), DIMENSION(jpi,jpj,jpk) , INTENT(inout) :: pu, pv ! now velocities 491 REAL(wp), DIMENSION(jpi,jpj,jpk),OPTIONAL, INTENT(inout) :: pu_rhs, pv_rhs ! total v-trend 320 492 ! 321 493 INTEGER :: ji, jj, jk ! dummy loop indices … … 371 543 END SELECT 372 544 ! 373 ! !== horizontal fluxes and potential vorticity ==! 374 zwx(:,:) = e2u(:,:) * e3u(:,:,jk,Kmm) * pu(:,:,jk) 375 zwy(:,:) = e1v(:,:) * e3v(:,:,jk,Kmm) * pv(:,:,jk) 376 zwz(:,:) = zwz(:,:) / e3f(:,:,jk) 377 ! 378 ! !== compute and add the vorticity term trend =! 379 DO_2D_00_00 380 zy1 = zwy(ji,jj-1) + zwy(ji+1,jj-1) 381 zy2 = zwy(ji,jj ) + zwy(ji+1,jj ) 382 zx1 = zwx(ji-1,jj) + zwx(ji-1,jj+1) 383 zx2 = zwx(ji ,jj) + zwx(ji ,jj+1) 384 pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + r1_4 * r1_e1u(ji,jj) * ( zwz(ji ,jj-1) * zy1 + zwz(ji,jj) * zy2 ) 385 pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) - r1_4 * r1_e2v(ji,jj) * ( zwz(ji-1,jj ) * zx1 + zwz(ji,jj) * zx2 ) 386 END_2D 545 IF( PRESENT( pu_rhs ) .AND. PRESENT( pv_rhs ) ) THEN !*** NO alternating direction ***! 546 ! 547 ! !== horizontal fluxes and potential vorticity ==! 548 zwx(:,:) = e2u(:,:) * e3u(:,:,jk,Kmm) * pu(:,:,jk) 549 zwy(:,:) = e1v(:,:) * e3v(:,:,jk,Kmm) * pv(:,:,jk) 550 zwz(:,:) = zwz(:,:) / e3f(:,:,jk) 551 ! 552 ! !== compute and add the vorticity term trend =! 553 DO_2D_00_00 554 zy1 = zwy(ji,jj-1) + zwy(ji+1,jj-1) 555 zy2 = zwy(ji,jj ) + zwy(ji+1,jj ) 556 zx1 = zwx(ji-1,jj) + zwx(ji-1,jj+1) 557 zx2 = zwx(ji ,jj) + zwx(ji ,jj+1) 558 pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + r1_4 * r1_e1u(ji,jj) * ( zwz(ji ,jj-1) * zy1 + zwz(ji,jj) * zy2 ) 559 pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) - r1_4 * r1_e2v(ji,jj) * ( zwz(ji-1,jj ) * zx1 + zwz(ji,jj) * zx2 ) 560 END_2D 561 ! 562 ! 563 ELSEIF( PRESENT( pu_rhs ) .AND. .NOT. PRESENT( pv_rhs ) ) THEN !*** Alternating direction : i-component ***! 564 ! 565 ! 566 ! !== horizontal fluxes and potential vorticity ==! 567 zwy(:,:) = e1v(:,:) * e3v(:,:,jk,Kmm) * pv(:,:,jk) 568 zwz(:,:) = zwz(:,:) / e3f(:,:,jk) 569 ! 570 ! !== compute and add the vorticity term trend =! 571 DO_2D_00_00 572 zy1 = zwy(ji,jj-1) + zwy(ji+1,jj-1) 573 zy2 = zwy(ji,jj ) + zwy(ji+1,jj ) 574 pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + r1_4 * r1_e1u(ji,jj) * ( zwz(ji ,jj-1) * zy1 + zwz(ji,jj) * zy2 ) 575 END_2D 576 ! 577 ELSEIF( .NOT. PRESENT( pu_rhs ) .AND. PRESENT( pv_rhs ) ) THEN !*** Alternating direction : j-component ***! 578 ! 579 ! 580 ! !== horizontal fluxes and potential vorticity ==! 581 zwx(:,:) = e2u(:,:) * e3u(:,:,jk,Kmm) * pu(:,:,jk) 582 zwz(:,:) = zwz(:,:) / e3f(:,:,jk) 583 ! 584 ! !== compute and add the vorticity term trend =! 585 DO_2D_00_00 586 zx1 = zwx(ji-1,jj) + zwx(ji-1,jj+1) 587 zx2 = zwx(ji ,jj) + zwx(ji ,jj+1) 588 pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) - r1_4 * r1_e2v(ji,jj) * ( zwz(ji-1,jj ) * zx1 + zwz(ji,jj) * zx2 ) 589 END_2D 590 ! 591 ENDIF 387 592 ! ! =============== 388 593 END DO ! End of slab … … 411 616 !! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689. 412 617 !!---------------------------------------------------------------------- 413 INTEGER , INTENT(in ) :: kt ! ocean time-step index414 INTEGER , INTENT(in ) :: Kmm ! ocean time level index415 INTEGER , INTENT(in ) :: kvor ! total, planetary, relative, or metric416 REAL(wp), DIMENSION(jpi,jpj,jpk) , INTENT(inout) :: pu, pv ! now velocities417 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pu_rhs, pv_rhs ! total v-trend618 INTEGER , INTENT(in ) :: kt ! ocean time-step index 619 INTEGER , INTENT(in ) :: Kmm ! ocean time level index 620 INTEGER , INTENT(in ) :: kvor ! total, planetary, relative, or metric 621 REAL(wp), DIMENSION(jpi,jpj,jpk) , INTENT(inout) :: pu, pv ! now velocities 622 REAL(wp), DIMENSION(jpi,jpj,jpk),OPTIONAL, INTENT(inout) :: pu_rhs, pv_rhs ! total v-trend 418 623 ! 419 624 INTEGER :: ji, jj, jk ! dummy loop indices … … 469 674 ! 470 675 ! 471 ! !== horizontal fluxes and potential vorticity ==! 472 zwx(:,:) = e2u(:,:) * e3u(:,:,jk,Kmm) * pu(:,:,jk) 473 zwy(:,:) = e1v(:,:) * e3v(:,:,jk,Kmm) * pv(:,:,jk) 474 zwz(:,:) = zwz(:,:) / e3f(:,:,jk) 475 ! 476 ! !== compute and add the vorticity term trend =! 477 DO_2D_00_00 478 zuav = r1_8 * r1_e1u(ji,jj) * ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) & 479 & + zwy(ji ,jj ) + zwy(ji+1,jj ) ) 480 zvau =-r1_8 * r1_e2v(ji,jj) * ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) & 481 & + zwx(ji ,jj ) + zwx(ji ,jj+1) ) 482 pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + zuav * ( zwz(ji ,jj-1) + zwz(ji,jj) ) 483 pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) + zvau * ( zwz(ji-1,jj ) + zwz(ji,jj) ) 484 END_2D 676 !!an wut ? v et u 677 IF( PRESENT( pu_rhs ) .AND. PRESENT( pv_rhs ) ) THEN !*** NO alternating direction ***! 678 ! 679 ! !== horizontal fluxes and potential vorticity ==! 680 zwx(:,:) = e2u(:,:) * e3u(:,:,jk,Kmm) * pu(:,:,jk) 681 zwy(:,:) = e1v(:,:) * e3v(:,:,jk,Kmm) * pv(:,:,jk) 682 zwz(:,:) = zwz(:,:) / e3f(:,:,jk) 683 ! 684 ! !== compute and add the vorticity term trend =! 685 DO_2D_00_00 686 zuav = r1_8 * r1_e1u(ji,jj) * ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) & 687 & + zwy(ji ,jj ) + zwy(ji+1,jj ) ) 688 zvau =-r1_8 * r1_e2v(ji,jj) * ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) & 689 & + zwx(ji ,jj ) + zwx(ji ,jj+1) ) 690 pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + zuav * ( zwz(ji ,jj-1) + zwz(ji,jj) ) 691 pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) + zvau * ( zwz(ji-1,jj ) + zwz(ji,jj) ) 692 END_2D 693 ! 694 ELSEIF( PRESENT( pu_rhs ) .AND. .NOT. PRESENT( pv_rhs ) ) THEN !*** Alternating direction : i-component ***! 695 ! 696 ! 697 ! !== horizontal fluxes and potential vorticity ==! 698 zwy(:,:) = e1v(:,:) * e3v(:,:,jk,Kmm) * pv(:,:,jk) 699 zwz(:,:) = zwz(:,:) / e3f(:,:,jk) 700 ! 701 ! !== compute and add the vorticity term trend =! 702 DO_2D_00_00 703 zuav = r1_8 * r1_e1u(ji,jj) * ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) & 704 & + zwy(ji ,jj ) + zwy(ji+1,jj ) ) 705 pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + zuav * ( zwz(ji ,jj-1) + zwz(ji,jj) ) 706 END_2D 707 ! 708 ELSEIF( .NOT. PRESENT( pu_rhs ) .AND. PRESENT( pv_rhs ) ) THEN !*** Alternating direction : j-component ***! 709 ! 710 ! 711 ! !== horizontal fluxes and potential vorticity ==! 712 zwx(:,:) = e2u(:,:) * e3u(:,:,jk,Kmm) * pu(:,:,jk) 713 zwz(:,:) = zwz(:,:) / e3f(:,:,jk) 714 ! 715 ! !== compute and add the vorticity term trend =! 716 DO_2D_00_00 717 zvau =-r1_8 * r1_e2v(ji,jj) * ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) & 718 & + zwx(ji ,jj ) + zwx(ji ,jj+1) ) 719 pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) + zvau * ( zwz(ji-1,jj ) + zwz(ji,jj) ) 720 END_2D 721 ! 722 ENDIF 485 723 ! ! =============== 486 724 END DO ! End of slab … … 762 1000 !! 763 1001 NAMELIST/namdyn_vor/ ln_dynvor_ens, ln_dynvor_ene, ln_dynvor_enT, ln_dynvor_eeT, & 764 & ln_dynvor_een, nn_een_e3f , ln_dynvor_mix, ln_dynvor_msk 1002 & ln_dynvor_een, nn_een_e3f , ln_dynvor_mix, ln_dynvor_msk, & 1003 & ln_dynvor_ens_adVO, ln_dynvor_ens_adKE, ln_dynvor_ens_adKEVO, & ! Alternative direction parameters 1004 & ln_dynvor_ene_adVO, ln_dynvor_ene_adKE, ln_dynvor_ene_adKEVO 765 1005 !!---------------------------------------------------------------------- 766 1006 ! … … 779 1019 IF(lwp) THEN ! Namelist print 780 1020 WRITE(numout,*) ' Namelist namdyn_vor : choice of the vorticity term scheme' 781 WRITE(numout,*) ' enstrophy conserving scheme ln_dynvor_ens = ', ln_dynvor_ens782 WRITE(numout,*) ' f-point energy conserving scheme ln_dynvor_ene = ', ln_dynvor_ene783 WRITE(numout,*) ' t-point energy conserving scheme ln_dynvor_enT = ', ln_dynvor_enT784 WRITE(numout,*) ' energy conserving scheme (een using e3t) ln_dynvor_eeT = ', ln_dynvor_eeT785 WRITE(numout,*) ' enstrophy and energy conserving scheme ln_dynvor_een = ', ln_dynvor_een786 WRITE(numout,*) ' e3f = averaging /4 (=0) or /sum(tmask) (=1) nn_een_e3f = ', nn_een_e3f787 WRITE(numout,*) ' mixed enstrophy/energy conserving scheme ln_dynvor_mix = ', ln_dynvor_mix788 WRITE(numout,*) ' masked (=T) or unmasked(=F) vorticity ln_dynvor_msk = ', ln_dynvor_msk1021 WRITE(numout,*) ' enstrophy conserving scheme ln_dynvor_ens = ', ln_dynvor_ens 1022 WRITE(numout,*) ' f-point energy conserving scheme ln_dynvor_ene = ', ln_dynvor_ene 1023 WRITE(numout,*) ' t-point energy conserving scheme ln_dynvor_enT = ', ln_dynvor_enT 1024 WRITE(numout,*) ' energy conserving scheme (een using e3t) ln_dynvor_eeT = ', ln_dynvor_eeT 1025 WRITE(numout,*) ' enstrophy and energy conserving scheme ln_dynvor_een = ', ln_dynvor_een 1026 WRITE(numout,*) ' e3f = averaging /4 (=0) or /sum(tmask) (=1) nn_een_e3f = ', nn_een_e3f 1027 WRITE(numout,*) ' mixed enstrophy/energy conserving scheme ln_dynvor_mix = ', ln_dynvor_mix 1028 WRITE(numout,*) ' masked (=T) or unmasked(=F) vorticity ln_dynvor_msk = ', ln_dynvor_msk 789 1029 ENDIF 790 1030 … … 799 1039 DO_3D_10_10( 1, jpk ) 800 1040 IF( tmask(ji,jj+1,jk) + tmask(ji+1,jj+1,jk) & 801 & + tmask(ji,jj ,jk) + tmask(ji+1,jj +1,jk) == 3._wp ) fmask(ji,jj,jk) = 1._wp1041 & + tmask(ji,jj ,jk) + tmask(ji+1,jj ,jk) == 3._wp ) fmask(ji,jj,jk) = 1._wp 802 1042 END_3D 803 1043 ! … … 809 1049 ioptio = 0 ! type of scheme for vorticity (set nvor_scheme) 810 1050 IF( ln_dynvor_ens ) THEN ; ioptio = ioptio + 1 ; nvor_scheme = np_ENS ; ENDIF 1051 IF( ln_dynvor_ens_adVO ) THEN ; ioptio = ioptio + 1 ; nvor_scheme = np_ENS_adVO ; ENDIF 1052 IF( ln_dynvor_ens_adKE ) THEN ; ioptio = ioptio + 1 ; nvor_scheme = np_ENS_adKE ; ENDIF 1053 IF( ln_dynvor_ens_adKEVO ) THEN ; ioptio = ioptio + 1 ; nvor_scheme = np_ENS_adKEVO ; ENDIF 811 1054 IF( ln_dynvor_ene ) THEN ; ioptio = ioptio + 1 ; nvor_scheme = np_ENE ; ENDIF 1055 IF( ln_dynvor_ene_adVO ) THEN ; ioptio = ioptio + 1 ; nvor_scheme = np_ENE_adVO ; ENDIF 1056 IF( ln_dynvor_ene_adKE ) THEN ; ioptio = ioptio + 1 ; nvor_scheme = np_ENE_adKE ; ENDIF 1057 IF( ln_dynvor_ene_adKEVO ) THEN ; ioptio = ioptio + 1 ; nvor_scheme = np_ENE_adKEVO ; ENDIF 812 1058 IF( ln_dynvor_enT ) THEN ; ioptio = ioptio + 1 ; nvor_scheme = np_ENT ; ENDIF 813 1059 IF( ln_dynvor_eeT ) THEN ; ioptio = ioptio + 1 ; nvor_scheme = np_EET ; ENDIF … … 826 1072 CASE( np_VEC_c2 ) 827 1073 IF(lwp) WRITE(numout,*) ' ==>>> vector form dynamics : total vorticity = Coriolis + relative vorticity' 828 nrvm = np_RVO ! relative vorticity 1074 nrvm = np_RVO ! relative vorticity 829 1075 ntot = np_CRV ! relative + planetary vorticity 830 1076 CASE( np_FLX_c2 , np_FLX_ubs ) … … 856 1102 WRITE(numout,*) 857 1103 SELECT CASE( nvor_scheme ) 858 CASE( np_ENS ) ; WRITE(numout,*) ' ==>>> enstrophy conserving scheme (ENS)' 859 CASE( np_ENE ) ; WRITE(numout,*) ' ==>>> energy conserving scheme (Coriolis at F-points) (ENE)' 860 CASE( np_ENT ) ; WRITE(numout,*) ' ==>>> energy conserving scheme (Coriolis at T-points) (ENT)' 861 CASE( np_EET ) ; WRITE(numout,*) ' ==>>> energy conserving scheme (EEN scheme using e3t) (EET)' 862 CASE( np_EEN ) ; WRITE(numout,*) ' ==>>> energy and enstrophy conserving scheme (EEN)' 863 CASE( np_MIX ) ; WRITE(numout,*) ' ==>>> mixed enstrophy/energy conserving scheme (MIX)' 1104 1105 CASE( np_ENS ) ; WRITE(numout,*) ' ==>>> enstrophy conserving scheme (ENS)' 1106 CASE( np_ENS_adVO ) ; WRITE(numout,*) ' ==>>> AD enstrophy conserving scheme (ENS_adVO) on vorticity only' 1107 CASE( np_ENS_adKE ) ; WRITE(numout,*) ' ==>>> AD enstrophy conserving scheme (ENS_adKE) on kinetic energy only' 1108 CASE( np_ENS_adKEVO ) ; WRITE(numout,*) ' ==>>> AD enstrophy conserving scheme (ENS_adKEVO) on kinetic energy and vorticity' 1109 1110 CASE( np_ENE ) ; WRITE(numout,*) ' ==>>> energy conserving scheme (Coriolis at F-points) (ENE)' 1111 CASE( np_ENE_adVO ) ; WRITE(numout,*) ' ==>>> AD energy conserving scheme (Coriolis at F-points) (ENE_adVO) on vorticity only' 1112 CASE( np_ENE_adKE ) ; WRITE(numout,*) ' ==>>> AD energy conserving scheme (Coriolis at F-points) (ENE_adKE) on kinetic energy only' 1113 CASE( np_ENE_adKEVO ) ; WRITE(numout,*) ' ==>>> AD energy conserving scheme (Coriolis at F-points) (ENE_adKEVO) on kinetic energy and vorticity' 1114 1115 CASE( np_ENT ) ; WRITE(numout,*) ' ==>>> energy conserving scheme (Coriolis at T-points) (ENT)' 1116 CASE( np_EET ) ; WRITE(numout,*) ' ==>>> energy conserving scheme (EEN scheme using e3t) (EET)' 1117 CASE( np_EEN ) ; WRITE(numout,*) ' ==>>> energy and enstrophy conserving scheme (EEN)' 1118 CASE( np_MIX ) ; WRITE(numout,*) ' ==>>> mixed enstrophy/energy conserving scheme (MIX)' 864 1119 END SELECT 865 1120 ENDIF -
NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/ldfdyn.F90
r12822 r13005 25 25 USE lib_mpp ! distribued memory computing library 26 26 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 27 27 ! 28 USE usrdef_nam , ONLY : ln_dynldf_lap_PM 29 ! 28 30 IMPLICIT NONE 29 31 PRIVATE … … 60 62 INTEGER , PUBLIC :: nldf_dyn !: type of lateral diffusion used defined from ln_dynldf_... (namlist logicals) 61 63 LOGICAL , PUBLIC :: l_ldfdyn_time !: flag for time variation of the lateral eddy viscosity coef. 62 64 !!an 65 !LOGICAL , PUBLIC :: ll_dynldf_lap_PM !: flag for P.Marchand modification on viscosity 66 !!an 63 67 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: ahmt, ahmf !: eddy viscosity coef. at T- and F-points [m2/s or m4/s] 64 68 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: dtensq !: horizontal tension squared (Smagorinsky only) … … 323 327 IF( .NOT.l_ldfdyn_time ) THEN !* No time variation 324 328 IF( ln_dynldf_lap ) THEN ! laplacian operator (mask only) 325 ahmt(:,:,1:jpkm1) = ahmt(:,:,1:jpkm1) * tmask(:,:,1:jpkm1) 326 WRITE(numout,*) ' ahmt tmask ' 329 !!an ! 330 WRITE(numout,*) ' ln_dynldf_lap_PM = ',ln_dynldf_lap_PM 331 IF( ln_dynldf_lap_PM ) THEN ! laplacian operator (mask only) 327 332 !! mask tension at the coast (equivalent of ghostpoints at T) 328 ! DO jk = 1, jpk 329 ! DO jj = 1, jpjm1 330 ! DO ji = 1, jpim1 ! NO vector opt. 331 ! ! si sum(fmask)==3 = mouillé (on touche pas) 332 ! ! si sum = 2 alors on met a 0 zsum = fmask + fmask + fmask,.. et si zsum < 2 -> 0 sinon = 1 333 ! zsum = fmask(ji,jj ,jk) + fmask(ji+1,jj ,jk) & 334 ! & + fmask(ji,jj+1,jk) + fmask(ji+1,jj+1,jk) 335 ! IF ( zsum < 2._wp ) ahmt(ji,jj,jk) = 0 336 ! ! 337 ! !ahmt(ji,jj,jk) = ahmt(ji,jj,jk) * fmask(ji,jj ,jk) * fmask(ji+1,jj ,jk) & 338 ! ! & * fmask(ji,jj+1,jk) * fmask(ji+1,jj+1,jk) 339 ! END DO 340 ! END DO 341 ! END DO 342 ! ahmt(jpi,:,1:jpkm1) = 0._wp 343 ! ahmt(:,jpj,1:jpkm1) = 0._wp 344 ! WRITE(numout,*) ' an45 ahmt x0' 345 346 ahmf(:,:,1:jpkm1) = ahmf(:,:,1:jpkm1) * fmask(:,:,1:jpkm1) 347 WRITE(numout,*) ' ahmf fmask ' 348 !!an apply no slip at the coast (ssfmask = 1 within the domain and at the coast contrary to fmask in free slip) 349 ! DO jk = 1, jpkm1 350 ! ahmf(:,:,jk) = ahmf(:,:,jk) * ( 2._wp * ssfmask(:,:) - fmask(:,:,jk) ) 351 ! END DO 352 ! WRITE(numout,*) ' an45 ahmf x2' 353 333 DO jk = 1, jpk 334 DO jj = 1, jpjm1 335 DO ji = 1, jpim1 ! NO vector opt. 336 ! si sum(fmask)==3 = mouillé (on touche pas) 337 ! si sum = 2 alors on met a 0 zsum = fmask + fmask + fmask,.. et si zsum < 2 -> 0 sinon = 1 338 zsum = fmask(ji,jj ,jk) + fmask(ji+1,jj ,jk) & 339 & + fmask(ji,jj+1,jk) + fmask(ji+1,jj+1,jk) 340 IF ( zsum < 2._wp ) ahmt(ji,jj,jk) = 0 341 ! 342 !ahmt(ji,jj,jk) = ahmt(ji,jj,jk) * fmask(ji,jj ,jk) * fmask(ji+1,jj ,jk) & 343 ! & * fmask(ji,jj+1,jk) * fmask(ji+1,jj+1,jk) 344 END DO 345 END DO 346 END DO 347 ahmt(jpi,:,1:jpkm1) = 0._wp 348 ahmt(:,jpj,1:jpkm1) = 0._wp 349 WRITE(numout,*) ' ahmt x0' 350 !! apply no slip at the coast (ssfmask = 1 within the domain and at the coast contrary to fmask in free slip) 351 DO jk = 1, jpkm1 352 ahmf(:,:,jk) = ahmf(:,:,jk) * ( 2._wp * ssfmask(:,:) - fmask(:,:,jk) ) 353 END DO 354 WRITE(numout,*) ' ahmf x2' 355 ELSE 356 ! classic boundary condition on the viscosity coefficient 357 ahmt(:,:,1:jpkm1) = ahmt(:,:,1:jpkm1) * tmask(:,:,1:jpkm1) 358 WRITE(numout,*) ' ahmt tmasked ' 359 ahmf(:,:,1:jpkm1) = ahmf(:,:,1:jpkm1) * fmask(:,:,1:jpkm1) 360 WRITE(numout,*) ' ahmf fmasked ' 361 ENDIF 362 !!an ! 354 363 ELSEIF( ln_dynldf_blp ) THEN ! bilaplacian operator (square root + mask) 355 364 ahmt(:,:,1:jpkm1) = SQRT( ahmt(:,:,1:jpkm1) ) * tmask(:,:,1:jpkm1) -
NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/step.F90
r12614 r13005 13 13 USE phycst ! physical constants 14 14 USE usrdef_nam 15 USE lib_mpp ! MPP library 16 USE dynvor , ONLY : ln_dynvor_ens_adVO, ln_dynvor_ens_adKE, ln_dynvor_ens_adKEVO, & 17 & ln_dynvor_ene_adVO, ln_dynvor_ene_adKE, ln_dynvor_ene_adKEVO 15 18 ! 16 19 USE iom ! xIOs server … … 138 141 & CALL Agrif_Sponge_dyn ! momentum sponge 139 142 #endif 140 CALL dyn_adv( kstp, Nbb, Nnn , uu, vv, Nrhs ) ! advection (VF or FF) ==> RHS141 142 CALL dyn_vor( kstp, Nnn , uu, vv, Nrhs ) ! vorticity ==> RHS143 144 CALL dyn_ldf( kstp, Nbb, Nnn , uu, vv, Nrhs ) ! lateral mixing145 143 146 144 !!an - calcul du gradient de pression horizontal (explicit) … … 150 148 END_3D 151 149 ! 150 151 ! IF( kstp == nit000 .AND. lwp ) THEN 152 ! WRITE(numout,*) 153 ! WRITE(numout,*) 'step.F90 : classic script used' 154 ! WRITE(numout,*) '~~~~~~~' 155 ! IF( ln_dynvor_ens_adVO .OR. ln_dynvor_ens_adKE .OR. ln_dynvor_ens_adKEVO & 156 ! & .OR. ln_dynvor_ene_adVO .OR. ln_dynvor_ene_adKE .OR. ln_dynvor_ene_adKEVO ) THEN 157 ! CALL ctl_stop('STOP','step : alternative direction asked but classis step' ) 158 ! ENDIF 159 ! ENDIF 160 !!an 161 ! CALL dyn_adv( kstp, Nbb, Nnn , uu, vv, Nrhs ) ! advection (VF or FF) ==> RHS 162 ! 163 ! CALL dyn_vor( kstp, Nbb, Nnn , uu, vv, Nrhs ) ! vorticity ==> RHS 164 ! 165 !!an In dynvor, dynkegAD is called even if not AD, so we keep the same step.F90 166 167 CALL dyn_vor( kstp, Nbb, Nnn , uu, vv, Nrhs) ! vorticity ==> RHS 168 169 CALL dyn_ldf( kstp, Nbb, Nnn , uu, vv, Nrhs ) ! lateral mixing 170 152 171 ! add wind stress forcing and layer linear friction to the RHS 153 172 z1_2rho0 = 0.5_wp * r1_rho0
Note: See TracChangeset
for help on using the changeset viewer.
