- Timestamp:
- 2020-09-14T17:40:34+02:00 (4 years ago)
- Location:
- NEMO/branches/2019/dev_r11351_fldread_with_XIOS
- Files:
-
- 11 edited
- 1 copied
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_r11351_fldread_with_XIOS
- Property svn:externals
-
old new 3 3 ^/utils/build/mk@HEAD mk 4 4 ^/utils/tools@HEAD tools 5 ^/vendors/AGRIF/dev @HEADext/AGRIF5 ^/vendors/AGRIF/dev_r12970_AGRIF_CMEMS ext/AGRIF 6 6 ^/vendors/FCM@HEAD ext/FCM 7 7 ^/vendors/IOIPSL@HEAD ext/IOIPSL 8 9 # SETTE 10 ^/utils/CI/sette@13382 sette
-
- Property svn:externals
-
NEMO/branches/2019/dev_r11351_fldread_with_XIOS/tests/VORTEX/EXPREF/1_context_nemo.xml
r9930 r13463 5 5 --> 6 6 <context id="1_nemo"> 7 <!-- $id$ --> 7 <!-- $id$ --> 8 <variable_definition> 9 <!-- Year/Month/Day of time origin for NetCDF files; defaults to 1800-01-01 --> 10 <variable id="ref_year" type="int"> 1900 </variable> 11 <variable id="ref_month" type="int"> 01 </variable> 12 <variable id="ref_day" type="int"> 01 </variable> 13 <variable id="rau0" type="float" > 1026.0 </variable> 14 <variable id="cpocean" type="float" > 3991.86795711963 </variable> 15 <variable id="convSpsu" type="float" > 0.99530670233846 </variable> 16 <variable id="rhoic" type="float" > 917.0 </variable> 17 <variable id="rhosn" type="float" > 330.0 </variable> 18 <variable id="missval" type="float" > 1.e20 </variable> 19 </variable_definition> 20 8 21 <!-- Fields definition --> 9 22 <field_definition src="./field_def_nemo-oce.xml"/> <!-- NEMO ocean dynamics --> … … 11 24 <!-- Files definition --> 12 25 <file_definition src="./file_def_nemo-oce.xml"/> <!-- NEMO ocean dynamics --> 13 <!-- 14 ============================================================================================================ 15 = grid definition = = DO NOT CHANGE = 16 ============================================================================================================ 17 --> 18 19 <axis_definition> 20 <axis id="deptht" long_name="Vertical T levels" unit="m" positive="down" /> 21 <axis id="depthu" long_name="Vertical U levels" unit="m" positive="down" /> 22 <axis id="depthv" long_name="Vertical V levels" unit="m" positive="down" /> 23 <axis id="depthw" long_name="Vertical W levels" unit="m" positive="down" /> 24 <axis id="nfloat" long_name="Float number" unit="-" /> 25 <axis id="icbcla" long_name="Iceberg class" unit="1" /> 26 <axis id="ncatice" long_name="Ice category" unit="1" /> 27 <axis id="iax_20C" long_name="20 degC isotherm" unit="degC" /> 28 <axis id="iax_28C" long_name="28 degC isotherm" unit="degC" /> 29 </axis_definition> 26 27 <!-- Axis definition --> 28 <axis_definition src="./axis_def_nemo.xml"/> 30 29 30 <!-- Domain definition --> 31 31 <domain_definition src="./domain_def_nemo.xml"/> 32 33 <!-- Grids definition --> 34 <grid_definition src="./grid_def_nemo.xml"/> 32 35 33 <grid_definition src="./grid_def_nemo.xml"/> 34 36 35 37 </context> -
NEMO/branches/2019/dev_r11351_fldread_with_XIOS/tests/VORTEX/EXPREF/1_namelist_cfg
r10075 r13463 40 40 !----------------------------------------------------------------------- 41 41 ln_linssh = .false. ! =T linear free surface ==>> model level are fixed in time 42 rn_ rdt = 480. ! time step for the dynamics (and tracer if nn_acc=0)42 rn_Dt = 480. ! time step for the dynamics (and tracer if nn_acc=0) 43 43 rn_atfp = 0.05 ! asselin time filter parameter 44 44 / … … 98 98 &namagrif ! AGRIF zoom ("key_agrif") 99 99 !----------------------------------------------------------------------- 100 ln_spc_dyn = .true. ! use 0 as special value for dynamics 101 rn_sponge_tra = 800. ! coefficient for tracer sponge layer [m2/s] 102 rn_sponge_dyn = 800. ! coefficient for dynamics sponge layer [m2/s] 103 ln_chk_bathy = .FALSE. ! 100 rn_sponge_tra = 0.00768 ! coefficient for tracer sponge layer [] 101 rn_sponge_dyn = 0.00768 ! coefficient for dynamics sponge layer [] 104 102 / 105 103 !!====================================================================== … … 107 105 !! !! 108 106 !! namdrg top/bottom drag coefficient (default: NO selection) 109 !! namdrg_top top friction (ln_ OFF=F & ln_isfcav=T)110 !! namdrg_bot bottom friction (ln_ OFF=F)107 !! namdrg_top top friction (ln_drg_OFF=F & ln_isfcav=T) 108 !! namdrg_bot bottom friction (ln_drg_OFF=F) 111 109 !! nambbc bottom temperature boundary condition (default: OFF) 112 110 !! nambbl bottom boundary layer scheme (default: OFF) … … 116 114 &namdrg ! top/bottom drag coefficient (default: NO selection) 117 115 !----------------------------------------------------------------------- 118 ln_ OFF = .true. ! free-slip : Cd = 0116 ln_drg_OFF = .true. ! free-slip : Cd = 0 (F => fill namdrg_bot 119 117 / 120 118 !!====================================================================== … … 213 211 ln_bt_av = .true. ! Time filtering of barotropic variables 214 212 nn_bt_flt = 1 ! Time filter choice = 0 None 215 ! ! = 1 Boxcar over nn_ barosub-steps216 ! ! = 2 Boxcar over 2*nn_ baro" "213 ! ! = 1 Boxcar over nn_e sub-steps 214 ! ! = 2 Boxcar over 2*nn_e " " 217 215 ln_bt_auto = .false. ! Number of sub-step defined from: 218 nn_ baro = 24 ! =F : the number of sub-step in rn_rdt seconds216 nn_e = 24 ! =F : the number of sub-step in rn_Dt seconds 219 217 / 220 218 !----------------------------------------------------------------------- … … 272 270 !! namdiu Cool skin and warm layer models (default: OFF) 273 271 !! namdiu Cool skin and warm layer models (default: OFF) 274 !! namflo float parameters ("key_float") 275 !! nam_diaharm Harmonic analysis of tidal constituents ("key_diaharm") 276 !! namdct transports through some sections ("key_diadct") 277 !! nam_diatmb Top Middle Bottom Output (default: OFF) 272 !! namflo float parameters (default: OFF) 273 !! nam_diadct transports through some sections (default: OFF) 278 274 !! nam_dia25h 25h Mean Output (default: OFF) 279 275 !! namnc4 netcdf4 chunking and compression settings ("key_netcdf4") -
NEMO/branches/2019/dev_r11351_fldread_with_XIOS/tests/VORTEX/EXPREF/context_nemo.xml
r9930 r13463 5 5 --> 6 6 <context id="nemo"> 7 <!-- $id$ --> 7 <!-- $id$ --> 8 <variable_definition> 9 <!-- Year/Month/Day of time origin for NetCDF files; defaults to 1800-01-01 --> 10 <variable id="ref_year" type="int"> 1900 </variable> 11 <variable id="ref_month" type="int"> 01 </variable> 12 <variable id="ref_day" type="int"> 01 </variable> 13 <variable id="rau0" type="float" > 1026.0 </variable> 14 <variable id="cpocean" type="float" > 3991.86795711963 </variable> 15 <variable id="convSpsu" type="float" > 0.99530670233846 </variable> 16 <variable id="rhoic" type="float" > 917.0 </variable> 17 <variable id="rhosn" type="float" > 330.0 </variable> 18 <variable id="missval" type="float" > 1.e20 </variable> 19 </variable_definition> 20 8 21 <!-- Fields definition --> 9 22 <field_definition src="./field_def_nemo-oce.xml"/> <!-- NEMO ocean dynamics --> … … 11 24 <!-- Files definition --> 12 25 <file_definition src="./file_def_nemo-oce.xml"/> <!-- NEMO ocean dynamics --> 13 <!-- 14 ============================================================================================================ 15 = grid definition = = DO NOT CHANGE = 16 ============================================================================================================ 17 --> 18 19 <axis_definition> 20 <axis id="deptht" long_name="Vertical T levels" unit="m" positive="down" /> 21 <axis id="depthu" long_name="Vertical U levels" unit="m" positive="down" /> 22 <axis id="depthv" long_name="Vertical V levels" unit="m" positive="down" /> 23 <axis id="depthw" long_name="Vertical W levels" unit="m" positive="down" /> 24 <axis id="nfloat" long_name="Float number" unit="-" /> 25 <axis id="icbcla" long_name="Iceberg class" unit="1" /> 26 <axis id="ncatice" long_name="Ice category" unit="1" /> 27 <axis id="iax_20C" long_name="20 degC isotherm" unit="degC" /> 28 <axis id="iax_28C" long_name="28 degC isotherm" unit="degC" /> 29 </axis_definition> 26 27 <!-- Axis definition --> 28 <axis_definition src="./axis_def_nemo.xml"/> 30 29 30 <!-- Domain definition --> 31 31 <domain_definition src="./domain_def_nemo.xml"/> 32 33 <!-- Grids definition --> 34 <grid_definition src="./grid_def_nemo.xml"/> 32 35 33 <grid_definition src="./grid_def_nemo.xml"/>34 36 35 37 </context> -
NEMO/branches/2019/dev_r11351_fldread_with_XIOS/tests/VORTEX/EXPREF/namelist_cfg
r10075 r13463 40 40 !----------------------------------------------------------------------- 41 41 ln_linssh = .false. ! =T linear free surface ==>> model level are fixed in time 42 rn_ rdt = 1440. ! time step for the dynamics (and tracer if nn_acc=0)42 rn_Dt = 1440. ! time step for the dynamics (and tracer if nn_acc=0) 43 43 rn_atfp = 0.05 ! asselin time filter parameter 44 44 / … … 99 99 !! !! 100 100 !! namdrg top/bottom drag coefficient (default: NO selection) 101 !! namdrg_top top friction (ln_ OFF=F & ln_isfcav=T)102 !! namdrg_bot bottom friction (ln_ OFF=F)101 !! namdrg_top top friction (ln_drg_OFF=F & ln_isfcav=T) 102 !! namdrg_bot bottom friction (ln_drg_OFF=F) 103 103 !! nambbc bottom temperature boundary condition (default: OFF) 104 104 !! nambbl bottom boundary layer scheme (default: OFF) … … 108 108 &namdrg ! top/bottom drag coefficient (default: NO selection) 109 109 !----------------------------------------------------------------------- 110 ln_ OFF = .true. ! free-slip : Cd = 0110 ln_drg_OFF = .true. ! free-slip : Cd = 0 (F => fill namdrg_bot 111 111 / 112 112 !!====================================================================== … … 204 204 ln_bt_av = .true. ! Time filtering of barotropic variables 205 205 nn_bt_flt = 1 ! Time filter choice = 0 None 206 ! ! = 1 Boxcar over nn_ barosub-steps207 ! ! = 2 Boxcar over 2*nn_ baro" "206 ! ! = 1 Boxcar over nn_e sub-steps 207 ! ! = 2 Boxcar over 2*nn_e " " 208 208 ln_bt_auto = .false. ! Number of sub-step defined from: 209 nn_ baro = 24 ! =F : the number of sub-step in rn_rdt seconds209 nn_e = 24 ! =F : the number of sub-step in rn_Dt seconds 210 210 / 211 211 !----------------------------------------------------------------------- … … 263 263 !! namdiu Cool skin and warm layer models (default: OFF) 264 264 !! namdiu Cool skin and warm layer models (default: OFF) 265 !! namflo float parameters ("key_float") 266 !! nam_diaharm Harmonic analysis of tidal constituents ("key_diaharm") 267 !! namdct transports through some sections ("key_diadct") 268 !! nam_diatmb Top Middle Bottom Output (default: OFF) 265 !! namflo float parameters (default: OFF) 266 !! nam_diadct transports through some sections (default: OFF) 269 267 !! nam_dia25h 25h Mean Output (default: OFF) 270 268 !! namnc4 netcdf4 chunking and compression settings ("key_netcdf4") -
NEMO/branches/2019/dev_r11351_fldread_with_XIOS/tests/VORTEX/MY_SRC/domvvl.F90
r10572 r13463 8 8 !! 3.3 ! 2011-10 (M. Leclair) totally rewrote domvvl: vvl option includes z_star and z_tilde coordinates 9 9 !! 3.6 ! 2014-11 (P. Mathiot) add ice shelf capability 10 !! 4.1 ! 2019-08 (A. Coward, D. Storkey) rename dom_vvl_sf_swp -> dom_vvl_sf_update for new timestepping 11 !! 4.x ! 2020-02 (G. Madec, S. Techene) introduce ssh to h0 ratio 10 12 !!---------------------------------------------------------------------- 11 13 12 !!----------------------------------------------------------------------13 !! dom_vvl_init : define initial vertical scale factors, depths and column thickness14 !! dom_vvl_sf_nxt : Compute next vertical scale factors15 !! dom_vvl_sf_swp : Swap vertical scale factors and update the vertical grid16 !! dom_vvl_interpol : Interpolate vertical scale factors from one grid point to another17 !! dom_vvl_rst : read/write restart file18 !! dom_vvl_ctl : Check the vvl options19 !!----------------------------------------------------------------------20 14 USE oce ! ocean dynamics and tracers 21 15 USE phycst ! physical constant … … 35 29 PRIVATE 36 30 37 PUBLIC dom_vvl_init ! called by domain.F9038 PUBLIC dom_vvl_sf_nxt ! called by step.F9039 PUBLIC dom_vvl_sf_swp ! called by step.F9040 PUBLIC dom_vvl_interpol ! called by dynnxt.F9041 42 31 ! !!* Namelist nam_vvl 43 32 LOGICAL , PUBLIC :: ln_vvl_zstar = .FALSE. ! zstar vertical coordinate … … 61 50 REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:) :: frq_rst_hdv ! retoring period for low freq. divergence 62 51 52 #if defined key_qco 53 !!---------------------------------------------------------------------- 54 !! 'key_qco' EMPTY MODULE Quasi-Eulerian vertical coordonate 55 !!---------------------------------------------------------------------- 56 #else 57 !!---------------------------------------------------------------------- 58 !! Default key Old management of time varying vertical coordinate 59 !!---------------------------------------------------------------------- 60 61 !!---------------------------------------------------------------------- 62 !! dom_vvl_init : define initial vertical scale factors, depths and column thickness 63 !! dom_vvl_sf_nxt : Compute next vertical scale factors 64 !! dom_vvl_sf_update : Swap vertical scale factors and update the vertical grid 65 !! dom_vvl_interpol : Interpolate vertical scale factors from one grid point to another 66 !! dom_vvl_rst : read/write restart file 67 !! dom_vvl_ctl : Check the vvl options 68 !!---------------------------------------------------------------------- 69 70 PUBLIC dom_vvl_init ! called by domain.F90 71 PUBLIC dom_vvl_zgr ! called by isfcpl.F90 72 PUBLIC dom_vvl_sf_nxt ! called by step.F90 73 PUBLIC dom_vvl_sf_update ! called by step.F90 74 PUBLIC dom_vvl_interpol ! called by dynnxt.F90 75 63 76 !! * Substitutions 64 # include " vectopt_loop_substitute.h90"77 # include "do_loop_substitute.h90" 65 78 !!---------------------------------------------------------------------- 66 79 !! NEMO/OCE 4.0 , NEMO Consortium (2018) … … 93 106 94 107 95 SUBROUTINE dom_vvl_init 108 SUBROUTINE dom_vvl_init( Kbb, Kmm, Kaa ) 96 109 !!---------------------------------------------------------------------- 97 110 !! *** ROUTINE dom_vvl_init *** … … 102 115 !! ** Method : - use restart file and/or initialize 103 116 !! - interpolate scale factors 117 !! 118 !! ** Action : - e3t_(n/b) and tilde_e3t_(n/b) 119 !! - Regrid: e3[u/v](:,:,:,Kmm) 120 !! e3[u/v](:,:,:,Kmm) 121 !! e3w(:,:,:,Kmm) 122 !! e3[u/v]w_b 123 !! e3[u/v]w_n 124 !! gdept(:,:,:,Kmm), gdepw(:,:,:,Kmm) and gde3w 125 !! - h(t/u/v)_0 126 !! - frq_rst_e3t and frq_rst_hdv 127 !! 128 !! Reference : Leclair, M., and G. Madec, 2011, Ocean Modelling. 129 !!---------------------------------------------------------------------- 130 INTEGER, INTENT(in) :: Kbb, Kmm, Kaa 131 ! 132 IF(lwp) WRITE(numout,*) 133 IF(lwp) WRITE(numout,*) 'dom_vvl_init : Variable volume activated' 134 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~' 135 ! 136 CALL dom_vvl_ctl ! choose vertical coordinate (z_star, z_tilde or layer) 137 ! 138 ! ! Allocate module arrays 139 IF( dom_vvl_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'dom_vvl_init : unable to allocate arrays' ) 140 ! 141 ! ! Read or initialize e3t_(b/n), tilde_e3t_(b/n) and hdiv_lf 142 CALL dom_vvl_rst( nit000, Kbb, Kmm, 'READ' ) 143 e3t(:,:,jpk,Kaa) = e3t_0(:,:,jpk) ! last level always inside the sea floor set one for all 144 ! 145 CALL dom_vvl_zgr(Kbb, Kmm, Kaa) ! interpolation scale factor, depth and water column 146 ! 147 END SUBROUTINE dom_vvl_init 148 149 150 SUBROUTINE dom_vvl_zgr(Kbb, Kmm, Kaa) 151 !!---------------------------------------------------------------------- 152 !! *** ROUTINE dom_vvl_init *** 153 !! 154 !! ** Purpose : Interpolation of all scale factors, 155 !! depths and water column heights 156 !! 157 !! ** Method : - interpolate scale factors 104 158 !! 105 159 !! ** Action : - e3t_(n/b) and tilde_e3t_(n/b) … … 115 169 !! Reference : Leclair, M., and G. Madec, 2011, Ocean Modelling. 116 170 !!---------------------------------------------------------------------- 171 INTEGER, INTENT(in) :: Kbb, Kmm, Kaa 172 !!---------------------------------------------------------------------- 117 173 INTEGER :: ji, jj, jk 118 174 INTEGER :: ii0, ii1, ij0, ij1 … … 120 176 !!---------------------------------------------------------------------- 121 177 ! 122 IF(lwp) WRITE(numout,*)123 IF(lwp) WRITE(numout,*) 'dom_vvl_init : Variable volume activated'124 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'125 !126 CALL dom_vvl_ctl ! choose vertical coordinate (z_star, z_tilde or layer)127 !128 ! ! Allocate module arrays129 IF( dom_vvl_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'dom_vvl_init : unable to allocate arrays' )130 !131 ! ! Read or initialize e3t_(b/n), tilde_e3t_(b/n) and hdiv_lf132 CALL dom_vvl_rst( nit000, 'READ' )133 e3t_a(:,:,jpk) = e3t_0(:,:,jpk) ! last level always inside the sea floor set one for all134 !135 178 ! !== Set of all other vertical scale factors ==! (now and before) 136 179 ! ! Horizontal interpolation of e3t 137 CALL dom_vvl_interpol( e3t _b(:,:,:), e3u_b(:,:,:), 'U' ) ! from T to U138 CALL dom_vvl_interpol( e3t _n(:,:,:), e3u_n(:,:,:), 'U' )139 CALL dom_vvl_interpol( e3t _b(:,:,:), e3v_b(:,:,:), 'V' ) ! from T to V140 CALL dom_vvl_interpol( e3t _n(:,:,:), e3v_n(:,:,:), 'V' )141 CALL dom_vvl_interpol( e3u _n(:,:,:), e3f_n(:,:,:), 'F' ) ! from U to F180 CALL dom_vvl_interpol( e3t(:,:,:,Kbb), e3u(:,:,:,Kbb), 'U' ) ! from T to U 181 CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3u(:,:,:,Kmm), 'U' ) 182 CALL dom_vvl_interpol( e3t(:,:,:,Kbb), e3v(:,:,:,Kbb), 'V' ) ! from T to V 183 CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3v(:,:,:,Kmm), 'V' ) 184 CALL dom_vvl_interpol( e3u(:,:,:,Kmm), e3f(:,:,:), 'F' ) ! from U to F 142 185 ! ! Vertical interpolation of e3t,u,v 143 CALL dom_vvl_interpol( e3t _n(:,:,:), e3w_n (:,:,:), 'W' ) ! from T to W144 CALL dom_vvl_interpol( e3t _b(:,:,:), e3w_b (:,:,:), 'W' )145 CALL dom_vvl_interpol( e3u _n(:,:,:), e3uw_n(:,:,:), 'UW' ) ! from U to UW146 CALL dom_vvl_interpol( e3u _b(:,:,:), e3uw_b(:,:,:), 'UW' )147 CALL dom_vvl_interpol( e3v _n(:,:,:), e3vw_n(:,:,:), 'VW' ) ! from V to UW148 CALL dom_vvl_interpol( e3v _b(:,:,:), e3vw_b(:,:,:), 'VW' )186 CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3w (:,:,:,Kmm), 'W' ) ! from T to W 187 CALL dom_vvl_interpol( e3t(:,:,:,Kbb), e3w (:,:,:,Kbb), 'W' ) 188 CALL dom_vvl_interpol( e3u(:,:,:,Kmm), e3uw(:,:,:,Kmm), 'UW' ) ! from U to UW 189 CALL dom_vvl_interpol( e3u(:,:,:,Kbb), e3uw(:,:,:,Kbb), 'UW' ) 190 CALL dom_vvl_interpol( e3v(:,:,:,Kmm), e3vw(:,:,:,Kmm), 'VW' ) ! from V to UW 191 CALL dom_vvl_interpol( e3v(:,:,:,Kbb), e3vw(:,:,:,Kbb), 'VW' ) 149 192 150 193 ! We need to define e3[tuv]_a for AGRIF initialisation (should not be a problem for the restartability...) 151 e3t _a(:,:,:) = e3t_n(:,:,:)152 e3u _a(:,:,:) = e3u_n(:,:,:)153 e3v _a(:,:,:) = e3v_n(:,:,:)194 e3t(:,:,:,Kaa) = e3t(:,:,:,Kmm) 195 e3u(:,:,:,Kaa) = e3u(:,:,:,Kmm) 196 e3v(:,:,:,Kaa) = e3v(:,:,:,Kmm) 154 197 ! 155 198 ! !== depth of t and w-point ==! (set the isf depth as it is in the initial timestep) 156 gdept_n(:,:,1) = 0.5_wp * e3w_n(:,:,1) ! reference to the ocean surface (used for MLD and light penetration) 157 gdepw_n(:,:,1) = 0.0_wp 158 gde3w_n(:,:,1) = gdept_n(:,:,1) - sshn(:,:) ! reference to a common level z=0 for hpg 159 gdept_b(:,:,1) = 0.5_wp * e3w_b(:,:,1) 160 gdepw_b(:,:,1) = 0.0_wp 161 DO jk = 2, jpk ! vertical sum 162 DO jj = 1,jpj 163 DO ji = 1,jpi 164 ! zcoef = tmask - wmask ! 0 everywhere tmask = wmask, ie everywhere expect at jk = mikt 165 ! ! 1 everywhere from mbkt to mikt + 1 or 1 (if no isf) 166 ! ! 0.5 where jk = mikt 167 !!gm ??????? BUG ? gdept_n as well as gde3w_n does not include the thickness of ISF ?? 168 zcoef = ( tmask(ji,jj,jk) - wmask(ji,jj,jk) ) 169 gdepw_n(ji,jj,jk) = gdepw_n(ji,jj,jk-1) + e3t_n(ji,jj,jk-1) 170 gdept_n(ji,jj,jk) = zcoef * ( gdepw_n(ji,jj,jk ) + 0.5 * e3w_n(ji,jj,jk)) & 171 & + (1-zcoef) * ( gdept_n(ji,jj,jk-1) + e3w_n(ji,jj,jk)) 172 gde3w_n(ji,jj,jk) = gdept_n(ji,jj,jk) - sshn(ji,jj) 173 gdepw_b(ji,jj,jk) = gdepw_b(ji,jj,jk-1) + e3t_b(ji,jj,jk-1) 174 gdept_b(ji,jj,jk) = zcoef * ( gdepw_b(ji,jj,jk ) + 0.5 * e3w_b(ji,jj,jk)) & 175 & + (1-zcoef) * ( gdept_b(ji,jj,jk-1) + e3w_b(ji,jj,jk)) 176 END DO 177 END DO 199 gdept(:,:,1,Kmm) = 0.5_wp * e3w(:,:,1,Kmm) ! reference to the ocean surface (used for MLD and light penetration) 200 gdepw(:,:,1,Kmm) = 0.0_wp 201 gde3w(:,:,1) = gdept(:,:,1,Kmm) - ssh(:,:,Kmm) ! reference to a common level z=0 for hpg 202 gdept(:,:,1,Kbb) = 0.5_wp * e3w(:,:,1,Kbb) 203 gdepw(:,:,1,Kbb) = 0.0_wp 204 DO_3D( 1, 1, 1, 1, 2, jpk ) 205 ! zcoef = tmask - wmask ! 0 everywhere tmask = wmask, ie everywhere expect at jk = mikt 206 ! ! 1 everywhere from mbkt to mikt + 1 or 1 (if no isf) 207 ! ! 0.5 where jk = mikt 208 !!gm ??????? BUG ? gdept(:,:,:,Kmm) as well as gde3w does not include the thickness of ISF ?? 209 zcoef = ( tmask(ji,jj,jk) - wmask(ji,jj,jk) ) 210 gdepw(ji,jj,jk,Kmm) = gdepw(ji,jj,jk-1,Kmm) + e3t(ji,jj,jk-1,Kmm) 211 gdept(ji,jj,jk,Kmm) = zcoef * ( gdepw(ji,jj,jk ,Kmm) + 0.5 * e3w(ji,jj,jk,Kmm)) & 212 & + (1-zcoef) * ( gdept(ji,jj,jk-1,Kmm) + e3w(ji,jj,jk,Kmm)) 213 gde3w(ji,jj,jk) = gdept(ji,jj,jk,Kmm) - ssh(ji,jj,Kmm) 214 gdepw(ji,jj,jk,Kbb) = gdepw(ji,jj,jk-1,Kbb) + e3t(ji,jj,jk-1,Kbb) 215 gdept(ji,jj,jk,Kbb) = zcoef * ( gdepw(ji,jj,jk ,Kbb) + 0.5 * e3w(ji,jj,jk,Kbb)) & 216 & + (1-zcoef) * ( gdept(ji,jj,jk-1,Kbb) + e3w(ji,jj,jk,Kbb)) 217 END_3D 218 ! 219 ! !== thickness of the water column !! (ocean portion only) 220 ht(:,:) = e3t(:,:,1,Kmm) * tmask(:,:,1) !!gm BUG : this should be 1/2 * e3w(k=1) .... 221 hu(:,:,Kbb) = e3u(:,:,1,Kbb) * umask(:,:,1) 222 hu(:,:,Kmm) = e3u(:,:,1,Kmm) * umask(:,:,1) 223 hv(:,:,Kbb) = e3v(:,:,1,Kbb) * vmask(:,:,1) 224 hv(:,:,Kmm) = e3v(:,:,1,Kmm) * vmask(:,:,1) 225 DO jk = 2, jpkm1 226 ht(:,:) = ht(:,:) + e3t(:,:,jk,Kmm) * tmask(:,:,jk) 227 hu(:,:,Kbb) = hu(:,:,Kbb) + e3u(:,:,jk,Kbb) * umask(:,:,jk) 228 hu(:,:,Kmm) = hu(:,:,Kmm) + e3u(:,:,jk,Kmm) * umask(:,:,jk) 229 hv(:,:,Kbb) = hv(:,:,Kbb) + e3v(:,:,jk,Kbb) * vmask(:,:,jk) 230 hv(:,:,Kmm) = hv(:,:,Kmm) + e3v(:,:,jk,Kmm) * vmask(:,:,jk) 178 231 END DO 179 232 ! 180 ! !== thickness of the water column !! (ocean portion only)181 ht_n(:,:) = e3t_n(:,:,1) * tmask(:,:,1) !!gm BUG : this should be 1/2 * e3w(k=1) ....182 hu_b(:,:) = e3u_b(:,:,1) * umask(:,:,1)183 hu_n(:,:) = e3u_n(:,:,1) * umask(:,:,1)184 hv_b(:,:) = e3v_b(:,:,1) * vmask(:,:,1)185 hv_n(:,:) = e3v_n(:,:,1) * vmask(:,:,1)186 DO jk = 2, jpkm1187 ht_n(:,:) = ht_n(:,:) + e3t_n(:,:,jk) * tmask(:,:,jk)188 hu_b(:,:) = hu_b(:,:) + e3u_b(:,:,jk) * umask(:,:,jk)189 hu_n(:,:) = hu_n(:,:) + e3u_n(:,:,jk) * umask(:,:,jk)190 hv_b(:,:) = hv_b(:,:) + e3v_b(:,:,jk) * vmask(:,:,jk)191 hv_n(:,:) = hv_n(:,:) + e3v_n(:,:,jk) * vmask(:,:,jk)192 END DO193 !194 233 ! !== inverse of water column thickness ==! (u- and v- points) 195 r1_hu _b(:,:) = ssumask(:,:) / ( hu_b(:,:) + 1._wp - ssumask(:,:) ) ! _i mask due to ISF196 r1_hu _n(:,:) = ssumask(:,:) / ( hu_n(:,:) + 1._wp - ssumask(:,:) )197 r1_hv _b(:,:) = ssvmask(:,:) / ( hv_b(:,:) + 1._wp - ssvmask(:,:) )198 r1_hv _n(:,:) = ssvmask(:,:) / ( hv_n(:,:) + 1._wp - ssvmask(:,:) )234 r1_hu(:,:,Kbb) = ssumask(:,:) / ( hu(:,:,Kbb) + 1._wp - ssumask(:,:) ) ! _i mask due to ISF 235 r1_hu(:,:,Kmm) = ssumask(:,:) / ( hu(:,:,Kmm) + 1._wp - ssumask(:,:) ) 236 r1_hv(:,:,Kbb) = ssvmask(:,:) / ( hv(:,:,Kbb) + 1._wp - ssvmask(:,:) ) 237 r1_hv(:,:,Kmm) = ssvmask(:,:) / ( hv(:,:,Kmm) + 1._wp - ssvmask(:,:) ) 199 238 200 239 ! !== z_tilde coordinate case ==! (Restoring frequencies) … … 208 247 IF( ln_vvl_ztilde_as_zstar ) THEN ! z-star emulation using z-tile 209 248 frq_rst_e3t(:,:) = 0._wp !Ignore namelist settings 210 frq_rst_hdv(:,:) = 1._wp / r dt249 frq_rst_hdv(:,:) = 1._wp / rn_Dt 211 250 ENDIF 212 251 IF ( ln_vvl_zstar_at_eqtor ) THEN ! use z-star in vicinity of the Equator 213 DO jj = 1, jpj 214 DO ji = 1, jpi 252 DO_2D( 1, 1, 1, 1 ) 215 253 !!gm case |gphi| >= 6 degrees is useless initialized just above by default 216 IF( ABS(gphit(ji,jj)) >= 6.) THEN 217 ! values outside the equatorial band and transition zone (ztilde) 218 frq_rst_e3t(ji,jj) = 2.0_wp * rpi / ( MAX( rn_rst_e3t , rsmall ) * 86400.e0_wp ) 219 frq_rst_hdv(ji,jj) = 2.0_wp * rpi / ( MAX( rn_lf_cutoff, rsmall ) * 86400.e0_wp ) 220 ELSEIF( ABS(gphit(ji,jj)) <= 2.5) THEN ! Equator strip ==> z-star 221 ! values inside the equatorial band (ztilde as zstar) 222 frq_rst_e3t(ji,jj) = 0.0_wp 223 frq_rst_hdv(ji,jj) = 1.0_wp / rdt 224 ELSE ! transition band (2.5 to 6 degrees N/S) 225 ! ! (linearly transition from z-tilde to z-star) 226 frq_rst_e3t(ji,jj) = 0.0_wp + (frq_rst_e3t(ji,jj)-0.0_wp)*0.5_wp & 227 & * ( 1.0_wp - COS( rad*(ABS(gphit(ji,jj))-2.5_wp) & 228 & * 180._wp / 3.5_wp ) ) 229 frq_rst_hdv(ji,jj) = (1.0_wp / rdt) & 230 & + ( frq_rst_hdv(ji,jj)-(1.e0_wp / rdt) )*0.5_wp & 231 & * ( 1._wp - COS( rad*(ABS(gphit(ji,jj))-2.5_wp) & 232 & * 180._wp / 3.5_wp ) ) 233 ENDIF 234 END DO 235 END DO 254 IF( ABS(gphit(ji,jj)) >= 6.) THEN 255 ! values outside the equatorial band and transition zone (ztilde) 256 frq_rst_e3t(ji,jj) = 2.0_wp * rpi / ( MAX( rn_rst_e3t , rsmall ) * 86400.e0_wp ) 257 frq_rst_hdv(ji,jj) = 2.0_wp * rpi / ( MAX( rn_lf_cutoff, rsmall ) * 86400.e0_wp ) 258 ELSEIF( ABS(gphit(ji,jj)) <= 2.5) THEN ! Equator strip ==> z-star 259 ! values inside the equatorial band (ztilde as zstar) 260 frq_rst_e3t(ji,jj) = 0.0_wp 261 frq_rst_hdv(ji,jj) = 1.0_wp / rn_Dt 262 ELSE ! transition band (2.5 to 6 degrees N/S) 263 ! ! (linearly transition from z-tilde to z-star) 264 frq_rst_e3t(ji,jj) = 0.0_wp + (frq_rst_e3t(ji,jj)-0.0_wp)*0.5_wp & 265 & * ( 1.0_wp - COS( rad*(ABS(gphit(ji,jj))-2.5_wp) & 266 & * 180._wp / 3.5_wp ) ) 267 frq_rst_hdv(ji,jj) = (1.0_wp / rn_Dt) & 268 & + ( frq_rst_hdv(ji,jj)-(1.e0_wp / rn_Dt) )*0.5_wp & 269 & * ( 1._wp - COS( rad*(ABS(gphit(ji,jj))-2.5_wp) & 270 & * 180._wp / 3.5_wp ) ) 271 ENDIF 272 END_2D 236 273 IF( cn_cfg == "orca" .OR. cn_cfg == "ORCA" ) THEN 237 274 IF( nn_cfg == 3 ) THEN ! ORCA2: Suppress ztilde in the Foxe Basin for ORCA2 238 ii0 = 103 ; ii1 = 111239 ij0 = 128 ; ij1 = 135 ;275 ii0 = 103 + nn_hls - 1 ; ii1 = 111 + nn_hls - 1 276 ij0 = 128 + nn_hls ; ij1 = 135 + nn_hls 240 277 frq_rst_e3t( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.0_wp 241 frq_rst_hdv( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1.e0_wp / r dt278 frq_rst_hdv( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1.e0_wp / rn_Dt 242 279 ENDIF 243 280 ENDIF … … 263 300 ENDIF 264 301 ! 265 END SUBROUTINE dom_vvl_ init266 267 268 SUBROUTINE dom_vvl_sf_nxt( kt, kcall )302 END SUBROUTINE dom_vvl_zgr 303 304 305 SUBROUTINE dom_vvl_sf_nxt( kt, Kbb, Kmm, Kaa, kcall ) 269 306 !!---------------------------------------------------------------------- 270 307 !! *** ROUTINE dom_vvl_sf_nxt *** … … 288 325 !! Reference : Leclair, M., and Madec, G. 2011, Ocean Modelling. 289 326 !!---------------------------------------------------------------------- 290 INTEGER, INTENT( in ) :: kt ! time step 291 INTEGER, INTENT( in ), OPTIONAL :: kcall ! optional argument indicating call sequence 327 INTEGER, INTENT( in ) :: kt ! time step 328 INTEGER, INTENT( in ) :: Kbb, Kmm, Kaa ! time step 329 INTEGER, INTENT( in ), OPTIONAL :: kcall ! optional argument indicating call sequence 292 330 ! 293 331 INTEGER :: ji, jj, jk ! dummy loop indices 294 332 INTEGER , DIMENSION(3) :: ijk_max, ijk_min ! temporary integers 295 REAL(wp) :: z 2dt, z_tmin, z_tmax! local scalars333 REAL(wp) :: z_tmin, z_tmax ! local scalars 296 334 LOGICAL :: ll_do_bclinic ! local logical 297 335 REAL(wp), DIMENSION(jpi,jpj) :: zht, z_scale, zwu, zwv, zhdiv 298 REAL(wp), DIMENSION(jpi,jpj,jpk) :: ze3t 336 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ze3t 337 LOGICAL , DIMENSION(:,:,:), ALLOCATABLE :: llmsk 299 338 !!---------------------------------------------------------------------- 300 339 ! … … 321 360 ! ! --------------------------------------------- ! 322 361 ! 323 z_scale(:,:) = ( ssh a(:,:) - sshb(:,:) ) * ssmask(:,:) / ( ht_0(:,:) + sshn(:,:) + 1. - ssmask(:,:) )362 z_scale(:,:) = ( ssh(:,:,Kaa) - ssh(:,:,Kbb) ) * ssmask(:,:) / ( ht_0(:,:) + ssh(:,:,Kmm) + 1. - ssmask(:,:) ) 324 363 DO jk = 1, jpkm1 325 ! formally this is the same as e3t _a= e3t_0*(1+ssha/ht_0)326 e3t _a(:,:,jk) = e3t_b(:,:,jk) + e3t_n(:,:,jk) * z_scale(:,:) * tmask(:,:,jk)364 ! formally this is the same as e3t(:,:,:,Kaa) = e3t_0*(1+ssha/ht_0) 365 e3t(:,:,jk,Kaa) = e3t(:,:,jk,Kbb) + e3t(:,:,jk,Kmm) * z_scale(:,:) * tmask(:,:,jk) 327 366 END DO 328 367 ! 329 IF( ln_vvl_ztilde .OR. ln_vvl_layer.AND. ll_do_bclinic ) THEN ! z_tilde or layer coordinate !330 ! ! ------baroclinic part------ !368 IF( (ln_vvl_ztilde .OR. ln_vvl_layer) .AND. ll_do_bclinic ) THEN ! z_tilde or layer coordinate ! 369 ! ! ------baroclinic part------ ! 331 370 ! I - initialization 332 371 ! ================== … … 337 376 zht(:,:) = 0._wp 338 377 DO jk = 1, jpkm1 339 zhdiv(:,:) = zhdiv(:,:) + e3t _n(:,:,jk) * hdivn(:,:,jk)340 zht (:,:) = zht (:,:) + e3t _n(:,:,jk) * tmask(:,:,jk)378 zhdiv(:,:) = zhdiv(:,:) + e3t(:,:,jk,Kmm) * hdiv(:,:,jk) 379 zht (:,:) = zht (:,:) + e3t(:,:,jk,Kmm) * tmask(:,:,jk) 341 380 END DO 342 381 zhdiv(:,:) = zhdiv(:,:) / ( zht(:,:) + 1. - tmask_i(:,:) ) … … 347 386 IF( kt > nit000 ) THEN 348 387 DO jk = 1, jpkm1 349 hdiv_lf(:,:,jk) = hdiv_lf(:,:,jk) - r dt * frq_rst_hdv(:,:) &350 & * ( hdiv_lf(:,:,jk) - e3t _n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) )388 hdiv_lf(:,:,jk) = hdiv_lf(:,:,jk) - rn_Dt * frq_rst_hdv(:,:) & 389 & * ( hdiv_lf(:,:,jk) - e3t(:,:,jk,Kmm) * ( hdiv(:,:,jk) - zhdiv(:,:) ) ) 351 390 END DO 352 391 ENDIF … … 361 400 IF( ln_vvl_ztilde ) THEN ! z_tilde case 362 401 DO jk = 1, jpkm1 363 tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - ( e3t _n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) - hdiv_lf(:,:,jk) )402 tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - ( e3t(:,:,jk,Kmm) * ( hdiv(:,:,jk) - zhdiv(:,:) ) - hdiv_lf(:,:,jk) ) 364 403 END DO 365 404 ELSE ! layer case 366 405 DO jk = 1, jpkm1 367 tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - e3t _n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) * tmask(:,:,jk)406 tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - e3t(:,:,jk,Kmm) * ( hdiv(:,:,jk) - zhdiv(:,:) ) * tmask(:,:,jk) 368 407 END DO 369 408 ENDIF … … 381 420 zwu(:,:) = 0._wp 382 421 zwv(:,:) = 0._wp 383 DO jk = 1, jpkm1 ! a - first derivative: diffusive fluxes 384 DO jj = 1, jpjm1 385 DO ji = 1, fs_jpim1 ! vector opt. 386 un_td(ji,jj,jk) = rn_ahe3 * umask(ji,jj,jk) * e2_e1u(ji,jj) & 387 & * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji+1,jj ,jk) ) 388 vn_td(ji,jj,jk) = rn_ahe3 * vmask(ji,jj,jk) * e1_e2v(ji,jj) & 389 & * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji ,jj+1,jk) ) 390 zwu(ji,jj) = zwu(ji,jj) + un_td(ji,jj,jk) 391 zwv(ji,jj) = zwv(ji,jj) + vn_td(ji,jj,jk) 392 END DO 393 END DO 394 END DO 395 DO jj = 1, jpj ! b - correction for last oceanic u-v points 396 DO ji = 1, jpi 397 un_td(ji,jj,mbku(ji,jj)) = un_td(ji,jj,mbku(ji,jj)) - zwu(ji,jj) 398 vn_td(ji,jj,mbkv(ji,jj)) = vn_td(ji,jj,mbkv(ji,jj)) - zwv(ji,jj) 399 END DO 400 END DO 401 DO jk = 1, jpkm1 ! c - second derivative: divergence of diffusive fluxes 402 DO jj = 2, jpjm1 403 DO ji = fs_2, fs_jpim1 ! vector opt. 404 tilde_e3t_a(ji,jj,jk) = tilde_e3t_a(ji,jj,jk) + ( un_td(ji-1,jj ,jk) - un_td(ji,jj,jk) & 405 & + vn_td(ji ,jj-1,jk) - vn_td(ji,jj,jk) & 406 & ) * r1_e1e2t(ji,jj) 407 END DO 408 END DO 409 END DO 422 DO_3D( 1, 0, 1, 0, 1, jpkm1 ) 423 un_td(ji,jj,jk) = rn_ahe3 * umask(ji,jj,jk) * e2_e1u(ji,jj) & 424 & * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji+1,jj ,jk) ) 425 vn_td(ji,jj,jk) = rn_ahe3 * vmask(ji,jj,jk) * e1_e2v(ji,jj) & 426 & * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji ,jj+1,jk) ) 427 zwu(ji,jj) = zwu(ji,jj) + un_td(ji,jj,jk) 428 zwv(ji,jj) = zwv(ji,jj) + vn_td(ji,jj,jk) 429 END_3D 430 DO_2D( 1, 1, 1, 1 ) 431 un_td(ji,jj,mbku(ji,jj)) = un_td(ji,jj,mbku(ji,jj)) - zwu(ji,jj) 432 vn_td(ji,jj,mbkv(ji,jj)) = vn_td(ji,jj,mbkv(ji,jj)) - zwv(ji,jj) 433 END_2D 434 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 435 tilde_e3t_a(ji,jj,jk) = tilde_e3t_a(ji,jj,jk) + ( un_td(ji-1,jj ,jk) - un_td(ji,jj,jk) & 436 & + vn_td(ji ,jj-1,jk) - vn_td(ji,jj,jk) & 437 & ) * r1_e1e2t(ji,jj) 438 END_3D 410 439 ! ! d - thickness diffusion transport: boundary conditions 411 440 ! (stored for tracer advction and continuity equation) … … 414 443 ! 4 - Time stepping of baroclinic scale factors 415 444 ! --------------------------------------------- 416 ! Leapfrog time stepping417 ! ~~~~~~~~~~~~~~~~~~~~~~418 IF( neuler == 0 .AND. kt == nit000 ) THEN419 z2dt = rdt420 ELSE421 z2dt = 2.0_wp * rdt422 ENDIF423 445 CALL lbc_lnk( 'domvvl', tilde_e3t_a(:,:,:), 'T', 1._wp ) 424 tilde_e3t_a(:,:,:) = tilde_e3t_b(:,:,:) + z2dt * tmask(:,:,:) * tilde_e3t_a(:,:,:)446 tilde_e3t_a(:,:,:) = tilde_e3t_b(:,:,:) + rDt * tmask(:,:,:) * tilde_e3t_a(:,:,:) 425 447 426 448 ! Maximum deformation control 427 449 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~ 428 ze3t(:,:,jpk) = 0._wp 429 DO jk = 1, jpkm1 430 ze3t(:,:,jk) = tilde_e3t_a(:,:,jk) / e3t_0(:,:,jk) * tmask(:,:,jk) * tmask_i(:,:) 431 END DO 432 z_tmax = MAXVAL( ze3t(:,:,:) ) 433 CALL mpp_max( 'domvvl', z_tmax ) ! max over the global domain 434 z_tmin = MINVAL( ze3t(:,:,:) ) 435 CALL mpp_min( 'domvvl', z_tmin ) ! min over the global domain 450 ALLOCATE( ze3t(jpi,jpj,jpk), llmsk(jpi,jpj,jpk) ) 451 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 452 ze3t(ji,jj,jk) = tilde_e3t_a(ji,jj,jk) / e3t_0(ji,jj,jk) * tmask(ji,jj,jk) * tmask_i(ji,jj) 453 END_3D 454 ! 455 llmsk( 1:Nis1,:,:) = .FALSE. ! exclude halos from the checked region 456 llmsk(Nie1: jpi,:,:) = .FALSE. 457 llmsk(:, 1:Njs1,:) = .FALSE. 458 llmsk(:,Nje1: jpj,:) = .FALSE. 459 ! 460 llmsk(Nis0:Nie0,Njs0:Nje0,:) = tmask(Nis0:Nie0,Njs0:Nje0,:) == 1._wp ! define only the inner domain 461 z_tmax = MAXVAL( ze3t(:,:,:), mask = llmsk ) ; CALL mpp_max( 'domvvl', z_tmax ) ! max over the global domain 462 z_tmin = MINVAL( ze3t(:,:,:), mask = llmsk ) ; CALL mpp_min( 'domvvl', z_tmin ) ! min over the global domain 436 463 ! - ML - test: for the moment, stop simulation for too large e3_t variations 437 464 IF( ( z_tmax > rn_zdef_max ) .OR. ( z_tmin < - rn_zdef_max ) ) THEN 438 IF( lk_mpp ) THEN 439 CALL mpp_maxloc( 'domvvl', ze3t, tmask, z_tmax, ijk_max ) 440 CALL mpp_minloc( 'domvvl', ze3t, tmask, z_tmin, ijk_min ) 441 ELSE 442 ijk_max = MAXLOC( ze3t(:,:,:) ) 443 ijk_max(1) = ijk_max(1) + nimpp - 1 444 ijk_max(2) = ijk_max(2) + njmpp - 1 445 ijk_min = MINLOC( ze3t(:,:,:) ) 446 ijk_min(1) = ijk_min(1) + nimpp - 1 447 ijk_min(2) = ijk_min(2) + njmpp - 1 448 ENDIF 465 CALL mpp_maxloc( 'domvvl', ze3t, llmsk, z_tmax, ijk_max ) 466 CALL mpp_minloc( 'domvvl', ze3t, llmsk, z_tmin, ijk_min ) 449 467 IF (lwp) THEN 450 468 WRITE(numout, *) 'MAX( tilde_e3t_a(:,:,:) / e3t_0(:,:,:) ) =', z_tmax … … 455 473 ENDIF 456 474 ENDIF 475 DEALLOCATE( ze3t, llmsk ) 457 476 ! - ML - end test 458 477 ! - ML - Imposing these limits will cause a baroclinicity error which is corrected for below … … 476 495 zht(:,:) = zht(:,:) + tilde_e3t_a(:,:,jk) * tmask(:,:,jk) 477 496 END DO 478 z_scale(:,:) = - zht(:,:) / ( ht_0(:,:) + ssh n(:,:) + 1. - ssmask(:,:) )497 z_scale(:,:) = - zht(:,:) / ( ht_0(:,:) + ssh(:,:,Kmm) + 1. - ssmask(:,:) ) 479 498 DO jk = 1, jpkm1 480 dtilde_e3t_a(:,:,jk) = dtilde_e3t_a(:,:,jk) + e3t _n(:,:,jk) * z_scale(:,:) * tmask(:,:,jk)499 dtilde_e3t_a(:,:,jk) = dtilde_e3t_a(:,:,jk) + e3t(:,:,jk,Kmm) * z_scale(:,:) * tmask(:,:,jk) 481 500 END DO 482 501 … … 486 505 ! ! ---baroclinic part--------- ! 487 506 DO jk = 1, jpkm1 488 e3t _a(:,:,jk) = e3t_a(:,:,jk) + dtilde_e3t_a(:,:,jk) * tmask(:,:,jk)507 e3t(:,:,jk,Kaa) = e3t(:,:,jk,Kaa) + dtilde_e3t_a(:,:,jk) * tmask(:,:,jk) 489 508 END DO 490 509 ENDIF … … 501 520 zht(:,:) = 0.0_wp 502 521 DO jk = 1, jpkm1 503 zht(:,:) = zht(:,:) + e3t _n(:,:,jk) * tmask(:,:,jk)522 zht(:,:) = zht(:,:) + e3t(:,:,jk,Kmm) * tmask(:,:,jk) 504 523 END DO 505 z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + ssh n(:,:) - zht(:,:) ) )524 z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + ssh(:,:,Kmm) - zht(:,:) ) ) 506 525 CALL mpp_max( 'domvvl', z_tmax ) ! max over the global domain 507 IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+sshn-SUM(e3t _n))) =', z_tmax526 IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+sshn-SUM(e3t(:,:,:,Kmm)))) =', z_tmax 508 527 ! 509 528 zht(:,:) = 0.0_wp 510 529 DO jk = 1, jpkm1 511 zht(:,:) = zht(:,:) + e3t _a(:,:,jk) * tmask(:,:,jk)530 zht(:,:) = zht(:,:) + e3t(:,:,jk,Kaa) * tmask(:,:,jk) 512 531 END DO 513 z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + ssh a(:,:) - zht(:,:) ) )532 z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + ssh(:,:,Kaa) - zht(:,:) ) ) 514 533 CALL mpp_max( 'domvvl', z_tmax ) ! max over the global domain 515 IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+ssha-SUM(e3t _a))) =', z_tmax534 IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+ssha-SUM(e3t(:,:,:,Kaa)))) =', z_tmax 516 535 ! 517 536 zht(:,:) = 0.0_wp 518 537 DO jk = 1, jpkm1 519 zht(:,:) = zht(:,:) + e3t _b(:,:,jk) * tmask(:,:,jk)538 zht(:,:) = zht(:,:) + e3t(:,:,jk,Kbb) * tmask(:,:,jk) 520 539 END DO 521 z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + ssh b(:,:) - zht(:,:) ) )540 z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + ssh(:,:,Kbb) - zht(:,:) ) ) 522 541 CALL mpp_max( 'domvvl', z_tmax ) ! max over the global domain 523 IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+sshb-SUM(e3t _b))) =', z_tmax524 ! 525 z_tmax = MAXVAL( tmask(:,:,1) * ABS( ssh b(:,:) ) )542 IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+sshb-SUM(e3t(:,:,:,Kbb)))) =', z_tmax 543 ! 544 z_tmax = MAXVAL( tmask(:,:,1) * ABS( ssh(:,:,Kbb) ) ) 526 545 CALL mpp_max( 'domvvl', z_tmax ) ! max over the global domain 527 IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(ssh b))) =', z_tmax528 ! 529 z_tmax = MAXVAL( tmask(:,:,1) * ABS( ssh n(:,:) ) )546 IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(ssh(:,:,Kbb)))) =', z_tmax 547 ! 548 z_tmax = MAXVAL( tmask(:,:,1) * ABS( ssh(:,:,Kmm) ) ) 530 549 CALL mpp_max( 'domvvl', z_tmax ) ! max over the global domain 531 IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(ssh n))) =', z_tmax532 ! 533 z_tmax = MAXVAL( tmask(:,:,1) * ABS( ssh a(:,:) ) )550 IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(ssh(:,:,Kmm)))) =', z_tmax 551 ! 552 z_tmax = MAXVAL( tmask(:,:,1) * ABS( ssh(:,:,Kaa) ) ) 534 553 CALL mpp_max( 'domvvl', z_tmax ) ! max over the global domain 535 IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(ssh a))) =', z_tmax554 IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(ssh(:,:,Kaa)))) =', z_tmax 536 555 END IF 537 556 … … 540 559 ! *********************************** ! 541 560 542 CALL dom_vvl_interpol( e3t _a(:,:,:), e3u_a(:,:,:), 'U' )543 CALL dom_vvl_interpol( e3t _a(:,:,:), e3v_a(:,:,:), 'V' )561 CALL dom_vvl_interpol( e3t(:,:,:,Kaa), e3u(:,:,:,Kaa), 'U' ) 562 CALL dom_vvl_interpol( e3t(:,:,:,Kaa), e3v(:,:,:,Kaa), 'V' ) 544 563 545 564 ! *********************************** ! … … 547 566 ! *********************************** ! 548 567 549 hu _a(:,:) = e3u_a(:,:,1) * umask(:,:,1)550 hv _a(:,:) = e3v_a(:,:,1) * vmask(:,:,1)568 hu(:,:,Kaa) = e3u(:,:,1,Kaa) * umask(:,:,1) 569 hv(:,:,Kaa) = e3v(:,:,1,Kaa) * vmask(:,:,1) 551 570 DO jk = 2, jpkm1 552 hu _a(:,:) = hu_a(:,:) + e3u_a(:,:,jk) * umask(:,:,jk)553 hv _a(:,:) = hv_a(:,:) + e3v_a(:,:,jk) * vmask(:,:,jk)571 hu(:,:,Kaa) = hu(:,:,Kaa) + e3u(:,:,jk,Kaa) * umask(:,:,jk) 572 hv(:,:,Kaa) = hv(:,:,Kaa) + e3v(:,:,jk,Kaa) * vmask(:,:,jk) 554 573 END DO 555 574 ! ! Inverse of the local depth 556 575 !!gm BUG ? don't understand the use of umask_i here ..... 557 r1_hu _a(:,:) = ssumask(:,:) / ( hu_a(:,:) + 1._wp - ssumask(:,:) )558 r1_hv _a(:,:) = ssvmask(:,:) / ( hv_a(:,:) + 1._wp - ssvmask(:,:) )576 r1_hu(:,:,Kaa) = ssumask(:,:) / ( hu(:,:,Kaa) + 1._wp - ssumask(:,:) ) 577 r1_hv(:,:,Kaa) = ssvmask(:,:) / ( hv(:,:,Kaa) + 1._wp - ssvmask(:,:) ) 559 578 ! 560 579 IF( ln_timing ) CALL timing_stop('dom_vvl_sf_nxt') … … 563 582 564 583 565 SUBROUTINE dom_vvl_sf_ swp( kt)566 !!---------------------------------------------------------------------- 567 !! *** ROUTINE dom_vvl_sf_ swp***584 SUBROUTINE dom_vvl_sf_update( kt, Kbb, Kmm, Kaa ) 585 !!---------------------------------------------------------------------- 586 !! *** ROUTINE dom_vvl_sf_update *** 568 587 !! 569 !! ** Purpose : compute time filter and swap of scale factors588 !! ** Purpose : for z tilde case: compute time filter and swap of scale factors 570 589 !! compute all depths and related variables for next time step 571 590 !! write outputs and restart file 572 591 !! 573 !! ** Method : - swap of e3t with trick for volume/tracer conservation 592 !! ** Method : - swap of e3t with trick for volume/tracer conservation (ONLY FOR Z TILDE CASE) 574 593 !! - reconstruct scale factor at other grid points (interpolate) 575 594 !! - recompute depths and water height fields 576 595 !! 577 !! ** Action : - e3t_(b/n), tilde_e3t_(b/n) and e3(u/v)_nready for next time step596 !! ** Action : - tilde_e3t_(b/n) ready for next time step 578 597 !! - Recompute: 579 598 !! e3(u/v)_b 580 !! e3w _n599 !! e3w(:,:,:,Kmm) 581 600 !! e3(u/v)w_b 582 601 !! e3(u/v)w_n 583 !! gdept _n, gdepw_n and gde3w_n602 !! gdept(:,:,:,Kmm), gdepw(:,:,:,Kmm) and gde3w 584 603 !! h(u/v) and h(u/v)r 585 604 !! … … 587 606 !! Leclair, M., and G. Madec, 2011, Ocean Modelling. 588 607 !!---------------------------------------------------------------------- 589 INTEGER, INTENT( in ) :: kt ! time step 608 INTEGER, INTENT( in ) :: kt ! time step 609 INTEGER, INTENT( in ) :: Kbb, Kmm, Kaa ! time level indices 590 610 ! 591 611 INTEGER :: ji, jj, jk ! dummy loop indices … … 595 615 IF( ln_linssh ) RETURN ! No calculation in linear free surface 596 616 ! 597 IF( ln_timing ) CALL timing_start('dom_vvl_sf_ swp')617 IF( ln_timing ) CALL timing_start('dom_vvl_sf_update') 598 618 ! 599 619 IF( kt == nit000 ) THEN 600 620 IF(lwp) WRITE(numout,*) 601 IF(lwp) WRITE(numout,*) 'dom_vvl_sf_ swp : - time filter and swap of scale factors'602 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~ - interpolate scale factors and compute depths for next time step'621 IF(lwp) WRITE(numout,*) 'dom_vvl_sf_update : - interpolate scale factors and compute depths for next time step' 622 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~~~' 603 623 ENDIF 604 624 ! … … 607 627 ! - ML - e3(t/u/v)_b are allready computed in dynnxt. 608 628 IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN 609 IF( neuler == 0 .AND. kt == nit000) THEN629 IF( l_1st_euler ) THEN 610 630 tilde_e3t_b(:,:,:) = tilde_e3t_n(:,:,:) 611 631 ELSE 612 632 tilde_e3t_b(:,:,:) = tilde_e3t_n(:,:,:) & 613 & + atfp * ( tilde_e3t_b(:,:,:) - 2.0_wp * tilde_e3t_n(:,:,:) + tilde_e3t_a(:,:,:) )633 & + rn_atfp * ( tilde_e3t_b(:,:,:) - 2.0_wp * tilde_e3t_n(:,:,:) + tilde_e3t_a(:,:,:) ) 614 634 ENDIF 615 635 tilde_e3t_n(:,:,:) = tilde_e3t_a(:,:,:) 616 636 ENDIF 617 gdept_b(:,:,:) = gdept_n(:,:,:)618 gdepw_b(:,:,:) = gdepw_n(:,:,:)619 620 e3t_n(:,:,:) = e3t_a(:,:,:)621 e3u_n(:,:,:) = e3u_a(:,:,:)622 e3v_n(:,:,:) = e3v_a(:,:,:)623 637 624 638 ! Compute all missing vertical scale factor and depths … … 626 640 ! Horizontal scale factor interpolations 627 641 ! -------------------------------------- 628 ! - ML - e3u _b and e3v_b are allready computed in dynnxt629 ! - JC - hu _b, hv_b, hur_b, hvr_b also642 ! - ML - e3u(:,:,:,Kbb) and e3v(:,:,:,Kbb) are already computed in dynnxt 643 ! - JC - hu(:,:,:,Kbb), hv(:,:,:,:,Kbb), hur_b, hvr_b also 630 644 631 CALL dom_vvl_interpol( e3u _n(:,:,:), e3f_n(:,:,:), 'F' )645 CALL dom_vvl_interpol( e3u(:,:,:,Kmm), e3f(:,:,:), 'F' ) 632 646 633 647 ! Vertical scale factor interpolations 634 CALL dom_vvl_interpol( e3t _n(:,:,:), e3w_n(:,:,:), 'W' )635 CALL dom_vvl_interpol( e3u _n(:,:,:), e3uw_n(:,:,:), 'UW' )636 CALL dom_vvl_interpol( e3v _n(:,:,:), e3vw_n(:,:,:), 'VW' )637 CALL dom_vvl_interpol( e3t _b(:,:,:), e3w_b(:,:,:), 'W' )638 CALL dom_vvl_interpol( e3u _b(:,:,:), e3uw_b(:,:,:), 'UW' )639 CALL dom_vvl_interpol( e3v _b(:,:,:), e3vw_b(:,:,:), 'VW' )648 CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3w(:,:,:,Kmm), 'W' ) 649 CALL dom_vvl_interpol( e3u(:,:,:,Kmm), e3uw(:,:,:,Kmm), 'UW' ) 650 CALL dom_vvl_interpol( e3v(:,:,:,Kmm), e3vw(:,:,:,Kmm), 'VW' ) 651 CALL dom_vvl_interpol( e3t(:,:,:,Kbb), e3w(:,:,:,Kbb), 'W' ) 652 CALL dom_vvl_interpol( e3u(:,:,:,Kbb), e3uw(:,:,:,Kbb), 'UW' ) 653 CALL dom_vvl_interpol( e3v(:,:,:,Kbb), e3vw(:,:,:,Kbb), 'VW' ) 640 654 641 655 ! t- and w- points depth (set the isf depth as it is in the initial step) 642 gdept_n(:,:,1) = 0.5_wp * e3w_n(:,:,1) 643 gdepw_n(:,:,1) = 0.0_wp 644 gde3w_n(:,:,1) = gdept_n(:,:,1) - sshn(:,:) 645 DO jk = 2, jpk 646 DO jj = 1,jpj 647 DO ji = 1,jpi 648 ! zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk)) ! 0 everywhere tmask = wmask, ie everywhere expect at jk = mikt 649 ! 1 for jk = mikt 650 zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk)) 651 gdepw_n(ji,jj,jk) = gdepw_n(ji,jj,jk-1) + e3t_n(ji,jj,jk-1) 652 gdept_n(ji,jj,jk) = zcoef * ( gdepw_n(ji,jj,jk ) + 0.5 * e3w_n(ji,jj,jk) ) & 653 & + (1-zcoef) * ( gdept_n(ji,jj,jk-1) + e3w_n(ji,jj,jk) ) 654 gde3w_n(ji,jj,jk) = gdept_n(ji,jj,jk) - sshn(ji,jj) 655 END DO 656 END DO 657 END DO 656 gdept(:,:,1,Kmm) = 0.5_wp * e3w(:,:,1,Kmm) 657 gdepw(:,:,1,Kmm) = 0.0_wp 658 gde3w(:,:,1) = gdept(:,:,1,Kmm) - ssh(:,:,Kmm) 659 DO_3D( 1, 1, 1, 1, 2, jpk ) 660 ! zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk)) ! 0 everywhere tmask = wmask, ie everywhere expect at jk = mikt 661 ! 1 for jk = mikt 662 zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk)) 663 gdepw(ji,jj,jk,Kmm) = gdepw(ji,jj,jk-1,Kmm) + e3t(ji,jj,jk-1,Kmm) 664 gdept(ji,jj,jk,Kmm) = zcoef * ( gdepw(ji,jj,jk ,Kmm) + 0.5 * e3w(ji,jj,jk,Kmm) ) & 665 & + (1-zcoef) * ( gdept(ji,jj,jk-1,Kmm) + e3w(ji,jj,jk,Kmm) ) 666 gde3w(ji,jj,jk) = gdept(ji,jj,jk,Kmm) - ssh(ji,jj,Kmm) 667 END_3D 658 668 659 669 ! Local depth and Inverse of the local depth of the water 660 670 ! ------------------------------------------------------- 661 hu_n(:,:) = hu_a(:,:) ; r1_hu_n(:,:) = r1_hu_a(:,:) 662 hv_n(:,:) = hv_a(:,:) ; r1_hv_n(:,:) = r1_hv_a(:,:) 663 ! 664 ht_n(:,:) = e3t_n(:,:,1) * tmask(:,:,1) 671 ! 672 ht(:,:) = e3t(:,:,1,Kmm) * tmask(:,:,1) 665 673 DO jk = 2, jpkm1 666 ht _n(:,:) = ht_n(:,:) + e3t_n(:,:,jk) * tmask(:,:,jk)674 ht(:,:) = ht(:,:) + e3t(:,:,jk,Kmm) * tmask(:,:,jk) 667 675 END DO 668 676 669 677 ! write restart file 670 678 ! ================== 671 IF( lrst_oce ) CALL dom_vvl_rst( kt, 'WRITE' )672 ! 673 IF( ln_timing ) CALL timing_stop('dom_vvl_sf_ swp')674 ! 675 END SUBROUTINE dom_vvl_sf_ swp679 IF( lrst_oce ) CALL dom_vvl_rst( kt, Kbb, Kmm, 'WRITE' ) 680 ! 681 IF( ln_timing ) CALL timing_stop('dom_vvl_sf_update') 682 ! 683 END SUBROUTINE dom_vvl_sf_update 676 684 677 685 … … 704 712 ! 705 713 CASE( 'U' ) !* from T- to U-point : hor. surface weighted mean 706 DO jk = 1, jpk 707 DO jj = 1, jpjm1 708 DO ji = 1, fs_jpim1 ! vector opt. 709 pe3_out(ji,jj,jk) = 0.5_wp * ( umask(ji,jj,jk) * (1.0_wp - zlnwd) + zlnwd ) * r1_e1e2u(ji,jj) & 710 & * ( e1e2t(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3t_0(ji ,jj,jk) ) & 711 & + e1e2t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) ) 712 END DO 713 END DO 714 END DO 714 DO_3D( 1, 0, 1, 0, 1, jpk ) 715 pe3_out(ji,jj,jk) = 0.5_wp * ( umask(ji,jj,jk) * (1.0_wp - zlnwd) + zlnwd ) * r1_e1e2u(ji,jj) & 716 & * ( e1e2t(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3t_0(ji ,jj,jk) ) & 717 & + e1e2t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) ) 718 END_3D 715 719 CALL lbc_lnk( 'domvvl', pe3_out(:,:,:), 'U', 1._wp ) 716 720 pe3_out(:,:,:) = pe3_out(:,:,:) + e3u_0(:,:,:) 717 721 ! 718 722 CASE( 'V' ) !* from T- to V-point : hor. surface weighted mean 719 DO jk = 1, jpk 720 DO jj = 1, jpjm1 721 DO ji = 1, fs_jpim1 ! vector opt. 722 pe3_out(ji,jj,jk) = 0.5_wp * ( vmask(ji,jj,jk) * (1.0_wp - zlnwd) + zlnwd ) * r1_e1e2v(ji,jj) & 723 & * ( e1e2t(ji,jj ) * ( pe3_in(ji,jj ,jk) - e3t_0(ji,jj ,jk) ) & 724 & + e1e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) ) 725 END DO 726 END DO 727 END DO 723 DO_3D( 1, 0, 1, 0, 1, jpk ) 724 pe3_out(ji,jj,jk) = 0.5_wp * ( vmask(ji,jj,jk) * (1.0_wp - zlnwd) + zlnwd ) * r1_e1e2v(ji,jj) & 725 & * ( e1e2t(ji,jj ) * ( pe3_in(ji,jj ,jk) - e3t_0(ji,jj ,jk) ) & 726 & + e1e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) ) 727 END_3D 728 728 CALL lbc_lnk( 'domvvl', pe3_out(:,:,:), 'V', 1._wp ) 729 729 pe3_out(:,:,:) = pe3_out(:,:,:) + e3v_0(:,:,:) 730 730 ! 731 731 CASE( 'F' ) !* from U-point to F-point : hor. surface weighted mean 732 DO jk = 1, jpk 733 DO jj = 1, jpjm1 734 DO ji = 1, fs_jpim1 ! vector opt. 735 pe3_out(ji,jj,jk) = 0.5_wp * ( umask(ji,jj,jk) * umask(ji,jj+1,jk) * (1.0_wp - zlnwd) + zlnwd ) & 736 & * r1_e1e2f(ji,jj) & 737 & * ( e1e2u(ji,jj ) * ( pe3_in(ji,jj ,jk) - e3u_0(ji,jj ,jk) ) & 738 & + e1e2u(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3u_0(ji,jj+1,jk) ) ) 739 END DO 740 END DO 741 END DO 732 DO_3D( 1, 0, 1, 0, 1, jpk ) 733 pe3_out(ji,jj,jk) = 0.5_wp * ( umask(ji,jj,jk) * umask(ji,jj+1,jk) * (1.0_wp - zlnwd) + zlnwd ) & 734 & * r1_e1e2f(ji,jj) & 735 & * ( e1e2u(ji,jj ) * ( pe3_in(ji,jj ,jk) - e3u_0(ji,jj ,jk) ) & 736 & + e1e2u(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3u_0(ji,jj+1,jk) ) ) 737 END_3D 742 738 CALL lbc_lnk( 'domvvl', pe3_out(:,:,:), 'F', 1._wp ) 743 739 pe3_out(:,:,:) = pe3_out(:,:,:) + e3f_0(:,:,:) … … 783 779 784 780 785 SUBROUTINE dom_vvl_rst( kt, cdrw )781 SUBROUTINE dom_vvl_rst( kt, Kbb, Kmm, cdrw ) 786 782 !!--------------------------------------------------------------------- 787 783 !! *** ROUTINE dom_vvl_rst *** … … 795 791 !! they are set to 0. 796 792 !!---------------------------------------------------------------------- 797 INTEGER , INTENT(in) :: kt ! ocean time-step 798 CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag 793 INTEGER , INTENT(in) :: kt ! ocean time-step 794 INTEGER , INTENT(in) :: Kbb, Kmm ! ocean time level indices 795 CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag 799 796 ! 800 797 INTEGER :: ji, jj, jk … … 806 803 IF( ln_rstart ) THEN !* Read the restart file 807 804 CALL rst_read_open ! open the restart file if necessary 808 CALL iom_get( numror, jpdom_auto glo, 'sshn' , sshn, ldxios = lrxios )805 CALL iom_get( numror, jpdom_auto, 'sshn' , ssh(:,:,Kmm), ldxios = lrxios ) 809 806 ! 810 807 id1 = iom_varid( numror, 'e3t_b', ldstop = .FALSE. ) … … 813 810 id4 = iom_varid( numror, 'tilde_e3t_n', ldstop = .FALSE. ) 814 811 id5 = iom_varid( numror, 'hdiv_lf', ldstop = .FALSE. ) 812 ! 815 813 ! ! --------- ! 816 814 ! ! all cases ! 817 815 ! ! --------- ! 816 ! 818 817 IF( MIN( id1, id2 ) > 0 ) THEN ! all required arrays exist 819 CALL iom_get( numror, jpdom_auto glo, 'e3t_b', e3t_b(:,:,:), ldxios = lrxios )820 CALL iom_get( numror, jpdom_auto glo, 'e3t_n', e3t_n(:,:,:), ldxios = lrxios )818 CALL iom_get( numror, jpdom_auto, 'e3t_b', e3t(:,:,:,Kbb), ldxios = lrxios ) 819 CALL iom_get( numror, jpdom_auto, 'e3t_n', e3t(:,:,:,Kmm), ldxios = lrxios ) 821 820 ! needed to restart if land processor not computed 822 IF(lwp) write(numout,*) 'dom_vvl_rst : e3t _b and e3t_nfound in restart files'821 IF(lwp) write(numout,*) 'dom_vvl_rst : e3t(:,:,:,Kbb) and e3t(:,:,:,Kmm) found in restart files' 823 822 WHERE ( tmask(:,:,:) == 0.0_wp ) 824 e3t _n(:,:,:) = e3t_0(:,:,:)825 e3t _b(:,:,:) = e3t_0(:,:,:)823 e3t(:,:,:,Kmm) = e3t_0(:,:,:) 824 e3t(:,:,:,Kbb) = e3t_0(:,:,:) 826 825 END WHERE 827 IF( neuler == 0) THEN828 e3t _b(:,:,:) = e3t_n(:,:,:)826 IF( l_1st_euler ) THEN 827 e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm) 829 828 ENDIF 830 829 ELSE IF( id1 > 0 ) THEN 831 IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t _nnot found in restart files'830 IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t(:,:,:,Kmm) not found in restart files' 832 831 IF(lwp) write(numout,*) 'e3t_n set equal to e3t_b.' 833 IF(lwp) write(numout,*) ' neuler is forced to 0'834 CALL iom_get( numror, jpdom_auto glo, 'e3t_b', e3t_b(:,:,:), ldxios = lrxios )835 e3t _n(:,:,:) = e3t_b(:,:,:)836 neuler = 0832 IF(lwp) write(numout,*) 'l_1st_euler is forced to true' 833 CALL iom_get( numror, jpdom_auto, 'e3t_b', e3t(:,:,:,Kbb), ldxios = lrxios ) 834 e3t(:,:,:,Kmm) = e3t(:,:,:,Kbb) 835 l_1st_euler = .true. 837 836 ELSE IF( id2 > 0 ) THEN 838 IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t _bnot found in restart files'837 IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t(:,:,:,Kbb) not found in restart files' 839 838 IF(lwp) write(numout,*) 'e3t_b set equal to e3t_n.' 840 IF(lwp) write(numout,*) ' neuler is forced to 0'841 CALL iom_get( numror, jpdom_auto glo, 'e3t_n', e3t_n(:,:,:), ldxios = lrxios )842 e3t _b(:,:,:) = e3t_n(:,:,:)843 neuler = 0839 IF(lwp) write(numout,*) 'l_1st_euler is forced to true' 840 CALL iom_get( numror, jpdom_auto, 'e3t_n', e3t(:,:,:,Kmm), ldxios = lrxios ) 841 e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm) 842 l_1st_euler = .true. 844 843 ELSE 845 IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t _nnot found in restart file'844 IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t(:,:,:,Kmm) not found in restart file' 846 845 IF(lwp) write(numout,*) 'Compute scale factor from sshn' 847 IF(lwp) write(numout,*) ' neuler is forced to 0'846 IF(lwp) write(numout,*) 'l_1st_euler is forced to true' 848 847 DO jk = 1, jpk 849 e3t _n(:,:,jk) = e3t_0(:,:,jk) * ( ht_0(:,:) + sshn(:,:) ) &848 e3t(:,:,jk,Kmm) = e3t_0(:,:,jk) * ( ht_0(:,:) + ssh(:,:,Kmm) ) & 850 849 & / ( ht_0(:,:) + 1._wp - ssmask(:,:) ) * tmask(:,:,jk) & 851 850 & + e3t_0(:,:,jk) * (1._wp -tmask(:,:,jk)) 852 851 END DO 853 e3t _b(:,:,:) = e3t_n(:,:,:)854 neuler = 0852 e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm) 853 l_1st_euler = .true. 855 854 ENDIF 856 855 ! ! ----------- ! … … 864 863 ! ! ----------------------- ! 865 864 IF( MIN( id3, id4 ) > 0 ) THEN ! all required arrays exist 866 CALL iom_get( numror, jpdom_auto glo, 'tilde_e3t_b', tilde_e3t_b(:,:,:), ldxios = lrxios )867 CALL iom_get( numror, jpdom_auto glo, 'tilde_e3t_n', tilde_e3t_n(:,:,:), ldxios = lrxios )865 CALL iom_get( numror, jpdom_auto, 'tilde_e3t_b', tilde_e3t_b(:,:,:), ldxios = lrxios ) 866 CALL iom_get( numror, jpdom_auto, 'tilde_e3t_n', tilde_e3t_n(:,:,:), ldxios = lrxios ) 868 867 ELSE ! one at least array is missing 869 868 tilde_e3t_b(:,:,:) = 0.0_wp … … 874 873 ! ! ------------ ! 875 874 IF( id5 > 0 ) THEN ! required array exists 876 CALL iom_get( numror, jpdom_auto glo, 'hdiv_lf', hdiv_lf(:,:,:), ldxios = lrxios )875 CALL iom_get( numror, jpdom_auto, 'hdiv_lf', hdiv_lf(:,:,:), ldxios = lrxios ) 877 876 ELSE ! array is missing 878 877 hdiv_lf(:,:,:) = 0.0_wp … … 888 887 IF( cn_cfg == 'wad' ) THEN 889 888 ! Wetting and drying test case 890 CALL usr_def_istate( gdept _b, tmask, tsb, ub, vb, sshb)891 ts n (:,:,:,:) = tsb (:,:,:,:) ! set now values from to before ones892 ssh n (:,:) = sshb(:,:)893 u n (:,:,:) = ub (:,:,:)894 v n (:,:,:) = vb (:,:,:)889 CALL usr_def_istate( gdept(:,:,:,Kbb), tmask, ts(:,:,:,:,Kbb), uu(:,:,:,Kbb), vv(:,:,:,Kbb), ssh(:,:,Kbb) ) 890 ts (:,:,:,:,Kmm) = ts (:,:,:,:,Kbb) ! set now values from to before ones 891 ssh (:,:,Kmm) = ssh(:,:,Kbb) 892 uu (:,:,:,Kmm) = uu (:,:,:,Kbb) 893 vv (:,:,:,Kmm) = vv (:,:,:,Kbb) 895 894 ELSE 896 895 ! if not test case 897 sshn(:,:) = -ssh_ref 898 sshb(:,:) = -ssh_ref 899 900 DO jj = 1, jpj 901 DO ji = 1, jpi 902 IF( ht_0(ji,jj)-ssh_ref < rn_wdmin1 ) THEN ! if total depth is less than min depth 903 904 sshb(ji,jj) = rn_wdmin1 - (ht_0(ji,jj) ) 905 sshn(ji,jj) = rn_wdmin1 - (ht_0(ji,jj) ) 906 ssha(ji,jj) = rn_wdmin1 - (ht_0(ji,jj) ) 907 ENDIF 908 ENDDO 909 ENDDO 896 ssh(:,:,Kmm) = -ssh_ref 897 ssh(:,:,Kbb) = -ssh_ref 898 899 DO_2D( 1, 1, 1, 1 ) 900 IF( ht_0(ji,jj)-ssh_ref < rn_wdmin1 ) THEN ! if total depth is less than min depth 901 ssh(ji,jj,Kbb) = rn_wdmin1 - (ht_0(ji,jj) ) 902 ssh(ji,jj,Kmm) = rn_wdmin1 - (ht_0(ji,jj) ) 903 ENDIF 904 END_2D 910 905 ENDIF !If test case else 911 906 912 907 ! Adjust vertical metrics for all wad 913 908 DO jk = 1, jpk 914 e3t _n(:,:,jk) = e3t_0(:,:,jk) * ( ht_0(:,:) + sshn(:,:) ) &909 e3t(:,:,jk,Kmm) = e3t_0(:,:,jk) * ( ht_0(:,:) + ssh(:,:,Kmm) ) & 915 910 & / ( ht_0(:,:) + 1._wp - ssmask(:,:) ) * tmask(:,:,jk) & 916 911 & + e3t_0(:,:,jk) * ( 1._wp - tmask(:,:,jk) ) 917 912 END DO 918 e3t_b(:,:,:) = e3t_n(:,:,:) 919 920 DO ji = 1, jpi 921 DO jj = 1, jpj 922 IF ( ht_0(ji,jj) .LE. 0.0 .AND. NINT( ssmask(ji,jj) ) .EQ. 1) THEN 923 CALL ctl_stop( 'dom_vvl_rst: ht_0 must be positive at potentially wet points' ) 924 ENDIF 925 END DO 926 END DO 913 e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm) 914 915 DO_2D( 1, 1, 1, 1 ) 916 IF ( ht_0(ji,jj) .LE. 0.0 .AND. NINT( ssmask(ji,jj) ) .EQ. 1) THEN 917 CALL ctl_stop( 'dom_vvl_rst: ht_0 must be positive at potentially wet points' ) 918 ENDIF 919 END_2D 927 920 ! 928 921 ELSE 929 922 ! 930 ! usr_def_istate called here only to get sshb, that is needed to initialize e3t_b and e3t_n 931 CALL usr_def_istate( gdept_0, tmask, tsb, ub, vb, sshb ) 932 ! usr_def_istate will be called again in istate_init to initialize ts(bn), ssh(bn), u(bn) and v(bn) 923 ! usr_def_istate called here only to get ssh(Kbb) needed to initialize e3t(Kbb) and e3t(Kmm) 924 ! 925 CALL usr_def_istate( gdept_0, tmask, ts(:,:,:,:,Kbb), uu(:,:,:,Kbb), vv(:,:,:,Kbb), ssh(:,:,Kbb) ) 926 ! 927 ! usr_def_istate will be called again in istate_init to initialize ts, ssh, u and v 933 928 ! 934 929 DO jk=1,jpk 935 e3t _b(:,:,jk) = e3t_0(:,:,jk) * ( ht_0(:,:) + sshb(:,:) ) &930 e3t(:,:,jk,Kbb) = e3t_0(:,:,jk) * ( ht_0(:,:) + ssh(:,:,Kbb) ) & 936 931 & / ( ht_0(:,:) + 1._wp - ssmask(:,:) ) * tmask(:,:,jk) & 937 & + e3t_0(:,:,jk) * ( 1._wp - tmask(:,:,jk) ) ! make sure e3t _b!= 0 on land points932 & + e3t_0(:,:,jk) * ( 1._wp - tmask(:,:,jk) ) ! make sure e3t(:,:,:,Kbb) != 0 on land points 938 933 END DO 939 e3t_n(:,:,:) = e3t_b(:,:,:) 940 sshn(:,:) = sshb(:,:) ! needed later for gde3w 941 !!$ e3t_n(:,:,:)=e3t_0(:,:,:) 942 !!$ e3t_b(:,:,:)=e3t_0(:,:,:) 934 e3t(:,:,:,Kmm) = e3t(:,:,:,Kbb) 935 ssh(:,:,Kmm) = ssh(:,:,Kbb) ! needed later for gde3w 943 936 ! 944 937 END IF ! end of ll_wd edits … … 958 951 ! ! all cases ! 959 952 ! ! --------- ! 960 CALL iom_rstput( kt, nitrst, numrow, 'e3t_b', e3t _b(:,:,:), ldxios = lwxios )961 CALL iom_rstput( kt, nitrst, numrow, 'e3t_n', e3t _n(:,:,:), ldxios = lwxios )953 CALL iom_rstput( kt, nitrst, numrow, 'e3t_b', e3t(:,:,:,Kbb), ldxios = lwxios ) 954 CALL iom_rstput( kt, nitrst, numrow, 'e3t_n', e3t(:,:,:,Kmm), ldxios = lwxios ) 962 955 ! ! ----------------------- ! 963 956 IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN ! z_tilde and layer cases ! … … 992 985 !!---------------------------------------------------------------------- 993 986 ! 994 REWIND( numnam_ref ) ! Namelist nam_vvl in reference namelist :995 987 READ ( numnam_ref, nam_vvl, IOSTAT = ios, ERR = 901) 996 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_vvl in reference namelist', lwp ) 997 REWIND( numnam_cfg ) ! Namelist nam_vvl in configuration namelist : Parameters of the run 988 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_vvl in reference namelist' ) 998 989 READ ( numnam_cfg, nam_vvl, IOSTAT = ios, ERR = 902 ) 999 902 IF( ios > 0 ) CALL ctl_nam ( ios , 'nam_vvl in configuration namelist' , lwp)990 902 IF( ios > 0 ) CALL ctl_nam ( ios , 'nam_vvl in configuration namelist' ) 1000 991 IF(lwm) WRITE ( numond, nam_vvl ) 1001 992 ! … … 1019 1010 WRITE(numout,*) ' rn_rst_e3t = 0.e0' 1020 1011 WRITE(numout,*) ' hard-wired : z-tilde cutoff frequency of low-pass filter (days)' 1021 WRITE(numout,*) ' rn_lf_cutoff = 1.0/r dt'1012 WRITE(numout,*) ' rn_lf_cutoff = 1.0/rn_Dt' 1022 1013 ELSE 1023 1014 WRITE(numout,*) ' z-tilde to zstar restoration timescale (days) rn_rst_e3t = ', rn_rst_e3t … … 1034 1025 ! 1035 1026 IF( ioptio /= 1 ) CALL ctl_stop( 'Choose ONE vertical coordinate in namelist nam_vvl' ) 1036 IF( .NOT. ln_vvl_zstar .AND. ln_isf ) CALL ctl_stop( 'Only vvl_zstar has been tested with ice shelf cavity' )1037 1027 ! 1038 1028 IF(lwp) THEN ! Print the choice … … 1050 1040 END SUBROUTINE dom_vvl_ctl 1051 1041 1042 #endif 1043 1052 1044 !!====================================================================== 1053 1045 END MODULE domvvl -
NEMO/branches/2019/dev_r11351_fldread_with_XIOS/tests/VORTEX/MY_SRC/usrdef_hgr.F90
r10074 r13463 26 26 PUBLIC usr_def_hgr ! called by domhgr.F90 27 27 28 !! * Substitutions 29 # include "do_loop_substitute.h90" 28 30 !!---------------------------------------------------------------------- 29 31 !! NEMO/OCE 4.0 , NEMO Consortium (2018) … … 61 63 REAL(wp), DIMENSION(:,:), INTENT(out) :: pe1e2u, pe1e2v ! u- & v-surfaces (if reduction in strait) [m2] 62 64 ! 63 INTEGER :: ji, jj ! dummy loop indices65 INTEGER :: ji, jj ! dummy loop indices 64 66 REAL(wp) :: zphi0, zlam0, zbeta, zf0 65 REAL(wp) :: zti, z ui, ztj, zvj ! local scalars67 REAL(wp) :: zti, ztj ! local scalars 66 68 !!------------------------------------------------------------------------------- 67 69 ! … … 75 77 ! Position coordinates (in kilometers) 76 78 ! ========== 77 zlam0 = -(jpiglo-1)/2 * 1.e-3 * rn_dx 78 zphi0 = -(jpjglo-1)/2 * 1.e-3 * rn_dy 79 79 #if defined key_agrif 80 IF( Agrif_Root() ) THEN 81 #endif 82 ! Compatibility WITH old version: 83 ! jperio = 0 => Ni0glo = jpigo_old_version 84 ! => jpiglo-1 replaced by Ni0glo-1 85 zlam0 = -REAL( (Ni0glo-1)/2, wp) * 1.e-3 * rn_dx 86 zphi0 = -REAL( (Nj0glo-1)/2, wp) * 1.e-3 * rn_dy 80 87 #if defined key_agrif 81 ! ! let lower left longitude and latitude from parent 82 IF (.NOT.Agrif_root()) THEN 83 zlam0 = (0.5_wp-(Agrif_parent(jpiglo)-1)/2)*1.e-3*Agrif_irhox()*rn_dx & 84 &+(Agrif_Ix()+nbghostcells-1)*Agrif_irhox()*rn_dx*1.e-3-(0.5_wp+nbghostcells)*rn_dx*1.e-3 85 zphi0 = (0.5_wp-(Agrif_parent(jpjglo)-1)/2)*1.e-3*Agrif_irhoy()*rn_dy & 86 &+(Agrif_Iy()+nbghostcells-1)*Agrif_irhoy()*rn_dy*1.e-3-(0.5_wp+nbghostcells)*rn_dy*1.e-3 88 ELSE 89 ! ! let lower left longitude and latitude from parent 90 ! Compatibility WITH old version: 91 ! jperio = 0 => Ni0glo = jpigo_old_version 92 ! => Agrif_parent(jpiglo)-1 replaced by Agrif_parent(Ni0glo)-1 93 zlam0 = ( 0.5_wp - REAL( ( Agrif_parent(Ni0glo)-1 ) / 2, wp) ) * 1.e-3 * Agrif_irhox() * rn_dx & 94 & + ( Agrif_Ix() + nbghostcells - 1 ) * Agrif_irhox() * rn_dx * 1.e-3 - ( 0.5_wp + nbghostcells ) * rn_dx * 1.e-3 95 zphi0 = ( 0.5_wp - REAL( ( Agrif_parent(Nj0glo)-1 ) / 2, wp) ) * 1.e-3 * Agrif_irhoy() * rn_dy & 96 & + ( Agrif_Iy() + nbghostcells - 1 ) * Agrif_irhoy() * rn_dy * 1.e-3 - ( 0.5_wp + nbghostcells ) * rn_dy * 1.e-3 87 97 ENDIF 88 98 #endif 89 99 90 DO jj = 1, jpj 91 DO ji = 1, jpi 92 zti = FLOAT( ji - 1 + nimpp - 1 ) ; ztj = FLOAT( jj - 1 + njmpp - 1 ) 93 zui = FLOAT( ji - 1 + nimpp - 1 ) + 0.5_wp ; zvj = FLOAT( jj - 1 + njmpp - 1 ) + 0.5_wp 94 95 plamt(ji,jj) = zlam0 + rn_dx * 1.e-3 * zti 96 plamu(ji,jj) = zlam0 + rn_dx * 1.e-3 * zui 97 plamv(ji,jj) = plamt(ji,jj) 98 plamf(ji,jj) = plamu(ji,jj) 99 100 pphit(ji,jj) = zphi0 + rn_dy * 1.e-3 * ztj 101 pphiv(ji,jj) = zphi0 + rn_dy * 1.e-3 * zvj 102 pphiu(ji,jj) = pphit(ji,jj) 103 pphif(ji,jj) = pphiv(ji,jj) 104 END DO 105 END DO 100 DO_2D( 1, 1, 1, 1 ) 101 zti = REAL( mig0_oldcmp(ji) - 1, wp ) ! start at i=0 in the global grid without halos 102 ztj = REAL( mjg0_oldcmp(jj) - 1, wp ) ! start at j=0 in the global grid without halos 103 104 plamt(ji,jj) = zlam0 + rn_dx * 1.e-3 * zti 105 plamu(ji,jj) = zlam0 + rn_dx * 1.e-3 * ( zti + 0.5_wp ) 106 plamv(ji,jj) = plamt(ji,jj) 107 plamf(ji,jj) = plamu(ji,jj) 108 109 pphit(ji,jj) = zphi0 + rn_dy * 1.e-3 * ztj 110 pphiv(ji,jj) = zphi0 + rn_dy * 1.e-3 * ( ztj + 0.5_wp ) 111 pphiu(ji,jj) = pphit(ji,jj) 112 pphif(ji,jj) = pphiv(ji,jj) 113 END_2D 106 114 ! 107 115 ! Horizontal scale factors (in meters) -
NEMO/branches/2019/dev_r11351_fldread_with_XIOS/tests/VORTEX/MY_SRC/usrdef_istate.F90
r10425 r13463 28 28 PUBLIC usr_def_istate ! called by istate.F90 29 29 30 !! * Substitutions 31 # include "do_loop_substitute.h90" 30 32 !!---------------------------------------------------------------------- 31 33 !! NEMO/OCE 4.0 , NEMO Consortium (2018) … … 69 71 zH = 0.5_wp * 5000._wp 70 72 ! 71 zP0 = r au0 * zf0 * zumax * zlambda * SQRT(EXP(1._wp)/2._wp)73 zP0 = rho0 * zf0 * zumax * zlambda * SQRT(EXP(1._wp)/2._wp) 72 74 ! 73 75 ! Sea level: 74 76 za = -zP0 * (1._wp-EXP(-zH)) / (grav*(zH-1._wp + EXP(-zH))) 75 DO ji=1, jpi 76 DO jj=1, jpj 77 zx = glamt(ji,jj) * 1.e3 78 zy = gphit(ji,jj) * 1.e3 79 zrho1 = rau0 + za * EXP(-(zx**2+zy**2)/zlambda**2) 80 pssh(ji,jj) = zP0 * EXP(-(zx**2+zy**2)/zlambda**2)/(zrho1*grav) * ptmask(ji,jj,1) 81 END DO 82 END DO 77 DO_2D( 1, 1, 1, 1 ) 78 zx = glamt(ji,jj) * 1.e3 79 zy = gphit(ji,jj) * 1.e3 80 zrho1 = rho0 + za * EXP(-(zx**2+zy**2)/zlambda**2) 81 pssh(ji,jj) = zP0 * EXP(-(zx**2+zy**2)/zlambda**2)/(zrho1*grav) * ptmask(ji,jj,1) 82 END_2D 83 83 ! 84 84 ! temperature: 85 DO ji=1, jpi 86 DO jj=1, jpj 87 zx = glamt(ji,jj) * 1.e3 88 zy = gphit(ji,jj) * 1.e3 89 DO jk=1,jpk 90 zdt = pdept(ji,jj,jk) 91 zrho1 = rau0 * (1._wp + zn2*zdt/grav) 92 IF (zdt < zH) THEN 93 zrho1 = zrho1 - zP0 * (1._wp-EXP(zdt-zH)) & 94 & * EXP(-(zx**2+zy**2)/zlambda**2) / (grav*(zH -1._wp + exp(-zH))); 95 ENDIF 96 pts(ji,jj,jk,jp_tem) = (20._wp + (rau0-zrho1) / 0.28_wp) * ptmask(ji,jj,jk) 97 END DO 85 DO_2D( 1, 1, 1, 1 ) 86 zx = glamt(ji,jj) * 1.e3 87 zy = gphit(ji,jj) * 1.e3 88 DO jk=1,jpk 89 zdt = pdept(ji,jj,jk) 90 zrho1 = rho0 * (1._wp + zn2*zdt/grav) 91 IF (zdt < zH) THEN 92 zrho1 = zrho1 - zP0 * (1._wp-EXP(zdt-zH)) & 93 & * EXP(-(zx**2+zy**2)/zlambda**2) / (grav*(zH -1._wp + EXP(-zH))); 94 ENDIF 95 pts(ji,jj,jk,jp_tem) = (20._wp + (rho0-zrho1) / 0.28_wp) * ptmask(ji,jj,jk) 98 96 END DO 99 END DO97 END_2D 100 98 ! 101 99 ! salinity: … … 103 101 ! 104 102 ! velocities: 105 za = 2._wp * zP0 / (zf0 * rau0 * zlambda**2) 106 DO ji=1, jpim1 107 DO jj=1, jpj 108 zx = glamu(ji,jj) * 1.e3 109 zy = gphiu(ji,jj) * 1.e3 110 DO jk=1, jpk 111 zdu = 0.5_wp * (pdept(ji ,jj,jk) + pdept(ji+1,jj,jk)) 112 IF (zdu < zH) THEN 113 zf = (zH-1._wp-zdu+EXP(zdu-zH)) / (zH-1._wp+EXP(-zH)) 114 pu(ji,jj,jk) = (za * zf * zy * EXP(-(zx**2+zy**2)/zlambda**2)) * ptmask(ji,jj,jk) * ptmask(ji+1,jj,jk) 115 ELSE 116 pu(ji,jj,jk) = 0._wp 117 ENDIF 118 END DO 103 za = 2._wp * zP0 / (zf0 * rho0 * zlambda**2) 104 DO_2D( 0, 0, 0, 0 ) 105 zx = glamu(ji,jj) * 1.e3 106 zy = gphiu(ji,jj) * 1.e3 107 DO jk=1, jpk 108 zdu = 0.5_wp * (pdept(ji ,jj,jk) + pdept(ji+1,jj,jk)) 109 IF (zdu < zH) THEN 110 zf = (zH-1._wp-zdu+EXP(zdu-zH)) / (zH-1._wp+EXP(-zH)) 111 pu(ji,jj,jk) = (za * zf * zy * EXP(-(zx**2+zy**2)/zlambda**2)) * ptmask(ji,jj,jk) * ptmask(ji+1,jj,jk) 112 ELSE 113 pu(ji,jj,jk) = 0._wp 114 ENDIF 119 115 END DO 120 END DO116 END_2D 121 117 ! 122 DO ji=1, jpi 123 DO jj=1, jpjm1 124 zx = glamv(ji,jj) * 1.e3 125 zy = gphiv(ji,jj) * 1.e3 126 DO jk=1, jpk 127 zdv = 0.5_wp * (pdept(ji ,jj,jk) + pdept(ji,jj+1,jk)) 128 IF (zdv < zH) THEN 129 zf = (zH-1._wp-zdv+EXP(zdv-zH)) / (zH-1._wp+EXP(-zH)) 130 pv(ji,jj,jk) = -(za * zf * zx * EXP(-(zx**2+zy**2)/zlambda**2)) * ptmask(ji,jj,jk) * ptmask(ji,jj+1,jk) 131 ELSE 132 pv(ji,jj,jk) = 0._wp 133 ENDIF 134 END DO 118 DO_2D( 0, 0, 0, 0 ) 119 zx = glamv(ji,jj) * 1.e3 120 zy = gphiv(ji,jj) * 1.e3 121 DO jk=1, jpk 122 zdv = 0.5_wp * (pdept(ji ,jj,jk) + pdept(ji,jj+1,jk)) 123 IF (zdv < zH) THEN 124 zf = (zH-1._wp-zdv+EXP(zdv-zH)) / (zH-1._wp+EXP(-zH)) 125 pv(ji,jj,jk) = -(za * zf * zx * EXP(-(zx**2+zy**2)/zlambda**2)) * ptmask(ji,jj,jk) * ptmask(ji,jj+1,jk) 126 ELSE 127 pv(ji,jj,jk) = 0._wp 128 ENDIF 135 129 END DO 136 END DO 137 138 CALL lbc_lnk( 'usrdef_istate', pu, 'U', -1. ) 139 CALL lbc_lnk( 'usrdef_istate', pv, 'V', -1. ) 130 END_2D 131 ! 132 CALL lbc_lnk_multi( 'usrdef_istate', pu, 'U', -1., pv, 'V', -1. ) 140 133 ! 141 134 END SUBROUTINE usr_def_istate -
NEMO/branches/2019/dev_r11351_fldread_with_XIOS/tests/VORTEX/MY_SRC/usrdef_nam.F90
r10074 r13463 14 14 !! usr_def_hgr : initialize the horizontal mesh 15 15 !!---------------------------------------------------------------------- 16 USE dom_oce , ONLY: nimpp , njmpp ! i- & j-indices of the local domain16 USE dom_oce 17 17 USE par_oce ! ocean space and time domain 18 18 USE phycst ! physical constants … … 40 40 CONTAINS 41 41 42 SUBROUTINE usr_def_nam( ldtxt, ldnam,cd_cfg, kk_cfg, kpi, kpj, kpk, kperio )42 SUBROUTINE usr_def_nam( cd_cfg, kk_cfg, kpi, kpj, kpk, kperio ) 43 43 !!---------------------------------------------------------------------- 44 44 !! *** ROUTINE dom_nam *** … … 52 52 !! ** input : - namusr_def namelist found in namelist_cfg 53 53 !!---------------------------------------------------------------------- 54 CHARACTER(len=*), DIMENSION(:), INTENT(out) :: ldtxt, ldnam ! stored print information55 54 CHARACTER(len=*) , INTENT(out) :: cd_cfg ! configuration name 56 55 INTEGER , INTENT(out) :: kk_cfg ! configuration resolution … … 58 57 INTEGER , INTENT(out) :: kperio ! lateral global domain b.c. 59 58 ! 60 INTEGER :: ios , ii! Local integer59 INTEGER :: ios ! Local integer 61 60 REAL(wp):: zlx, zly, zh ! Local scalars 62 61 !! … … 64 63 !!---------------------------------------------------------------------- 65 64 ! 66 ii = 167 !68 REWIND( numnam_cfg ) ! Namelist namusr_def (exist in namelist_cfg only)69 65 READ ( numnam_cfg, namusr_def, IOSTAT = ios, ERR = 902 ) 70 902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namusr_def in configuration namelist' , .TRUE.)66 902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namusr_def in configuration namelist' ) 71 67 ! 72 68 #if defined key_agrif … … 80 76 #endif 81 77 ! 82 WRITE( ldnam(:), namusr_def )78 IF(lwm) WRITE( numond, namusr_def ) 83 79 ! 84 80 cd_cfg = 'VORTEX' ! name & resolution (not used) 85 kk_cfg = INT( rn_dx )81 kk_cfg = nINT( rn_dx ) 86 82 ! 87 ! Global Domain size: VORTEX global domain is 1800 km x 1800 Km x 5000 m88 kpi =INT( 1800.e3 / rn_dx ) + 389 kpj =INT( 1800.e3 / rn_dy ) + 390 kpk = INT( 5000._wp / rn_dz ) + 191 #if defined key_agrif 92 IF( .NOT. Agrif_Root() ) THEN93 kpi = nbcellsx + 2 + 2*nbghostcells 94 kpj = nbcellsy + 2 + 2*nbghostcells 83 IF( Agrif_Root() ) THEN ! Global Domain size: VORTEX global domain is 1800 km x 1800 Km x 5000 m 84 kpi = NINT( 1800.e3 / rn_dx ) + 3 85 kpj = NINT( 1800.e3 / rn_dy ) + 3 86 ELSE ! Global Domain size: add nbghostcells + 1 "land" point on each side 87 kpi = nbcellsx + 2 * ( nbghostcells + 1 ) 88 kpj = nbcellsy + 2 * ( nbghostcells + 1 ) 89 !!$ kpi = nbcellsx + nbghostcells_x + nbghostcells_x + 2 90 !!$ kpj = nbcellsy + nbghostcells_y_s + nbghostcells_y_n + 2 95 91 ENDIF 96 #endif 92 kpk = NINT( 5000._wp / rn_dz ) + 1 97 93 ! 98 94 zlx = (kpi-2)*rn_dx*1.e-3 99 95 zly = (kpj-2)*rn_dy*1.e-3 100 96 zh = (kpk-1)*rn_dz 101 ! ! control print102 WRITE(ldtxt(ii),*) ' ' ; ii = ii + 1103 WRITE(ldtxt(ii),*) 'usr_def_nam : read the user defined namelist (namusr_def) in namelist_cfg' ; ii = ii + 1104 WRITE(ldtxt(ii),*) '~~~~~~~~~~~ ' ; ii = ii + 1105 WRITE(ldtxt(ii),*) ' Namelist namusr_def : VORTEX test case' ; ii = ii + 1106 WRITE(ldtxt(ii),*) ' horizontal resolution rn_dx = ', rn_dx, ' m' ; ii = ii + 1107 WRITE(ldtxt(ii),*) ' horizontal resolution rn_dy = ', rn_dy, ' m' ; ii = ii + 1108 WRITE(ldtxt(ii),*) ' vertical resolution rn_dz = ', rn_dz, ' m' ; ii = ii + 1109 WRITE(ldtxt(ii),*) ' VORTEX domain: ' ; ii = ii + 1110 WRITE(ldtxt(ii),*) ' LX [km]: ', zlx ; ii = ii + 1111 WRITE(ldtxt(ii),*) ' LY [km]: ', zly ; ii = ii + 1112 WRITE(ldtxt(ii),*) ' H [m] : ', zh ; ii = ii + 1113 WRITE(ldtxt(ii),*) ' Reference latitude rn_ppgphi0 = ', rn_ppgphi0 ; ii = ii + 1114 !115 97 ! ! Set the lateral boundary condition of the global domain 116 98 kperio = 0 ! VORTEX configuration : closed basin 117 ! 118 WRITE(ldtxt(ii),*) ' ' ; ii = ii + 1 119 WRITE(ldtxt(ii),*) ' Lateral boundary condition of the global domain' ; ii = ii + 1 120 WRITE(ldtxt(ii),*) ' VORTEX : closed basin jperio = ', kperio ; ii = ii + 1 99 ! ! control print 100 IF(lwp) THEN 101 WRITE(numout,*) ' ' 102 WRITE(numout,*) 'usr_def_nam : read the user defined namelist (namusr_def) in namelist_cfg' 103 WRITE(numout,*) '~~~~~~~~~~~ ' 104 WRITE(numout,*) ' Namelist namusr_def : VORTEX test case' 105 WRITE(numout,*) ' horizontal resolution rn_dx = ', rn_dx, ' m' 106 WRITE(numout,*) ' horizontal resolution rn_dy = ', rn_dy, ' m' 107 WRITE(numout,*) ' vertical resolution rn_dz = ', rn_dz, ' m' 108 WRITE(numout,*) ' resulting global domain size : Ni0glo = ', kpi 109 WRITE(numout,*) ' Nj0glo = ', kpj 110 WRITE(numout,*) ' jpkglo = ', kpk 111 WRITE(numout,*) ' VORTEX domain: ' 112 WRITE(numout,*) ' LX [km]: ', zlx 113 WRITE(numout,*) ' LY [km]: ', zly 114 WRITE(numout,*) ' H [m] : ', zh 115 WRITE(numout,*) ' Reference latitude rn_ppgphi0 = ', rn_ppgphi0 116 WRITE(numout,*) ' ' 117 WRITE(numout,*) ' Lateral boundary condition of the global domain' 118 WRITE(numout,*) ' VORTEX : closed basin jperio = ', kperio 119 ENDIF 121 120 ! 122 121 END SUBROUTINE usr_def_nam -
NEMO/branches/2019/dev_r11351_fldread_with_XIOS/tests/VORTEX/MY_SRC/usrdef_sbc.F90
r10074 r13463 30 30 PUBLIC usrdef_sbc_ice_flx ! routine called by icestp.F90 for ice thermo 31 31 32 !! * Substitutions33 # include "vectopt_loop_substitute.h90"34 32 !!---------------------------------------------------------------------- 35 33 !! NEMO/OCE 4.0 , NEMO Consortium (2018) … … 39 37 CONTAINS 40 38 41 SUBROUTINE usrdef_sbc_oce( kt )39 SUBROUTINE usrdef_sbc_oce( kt, Kbb ) 42 40 !!--------------------------------------------------------------------- 43 41 !! *** ROUTINE usr_def_sbc *** … … 54 52 !!---------------------------------------------------------------------- 55 53 INTEGER, INTENT(in) :: kt ! ocean time step 54 INTEGER, INTENT(in) :: Kbb ! ocean time index 56 55 !!--------------------------------------------------------------------- 57 56 ! … … 79 78 END SUBROUTINE usrdef_sbc_ice_tau 80 79 81 SUBROUTINE usrdef_sbc_ice_flx( kt ) 80 81 SUBROUTINE usrdef_sbc_ice_flx( kt, phs, phi ) 82 82 INTEGER, INTENT(in) :: kt ! ocean time step 83 REAL(wp), DIMENSION(:,:,:), INTENT(in) :: phs ! snow thickness 84 REAL(wp), DIMENSION(:,:,:), INTENT(in) :: phi ! ice thickness 83 85 END SUBROUTINE usrdef_sbc_ice_flx 84 86 -
NEMO/branches/2019/dev_r11351_fldread_with_XIOS/tests/VORTEX/MY_SRC/usrdef_zgr.F90
r10425 r13463 29 29 PUBLIC usr_def_zgr ! called by domzgr.F90 30 30 31 !! * Substitutions32 # include "vectopt_loop_substitute.h90"33 31 !!---------------------------------------------------------------------- 34 32 !! NEMO/OCE 4.0 , NEMO Consortium (2018) … … 194 192 CALL lbc_lnk( 'usrdef_zgr', z2d, 'T', 1. ) ! set surrounding land to zero (here jperio=0 ==>> closed) 195 193 ! 196 k_bot(:,:) = INT( z2d(:,:) )! =jpkm1 over the ocean point, =0 elsewhere194 k_bot(:,:) = NINT( z2d(:,:) ) ! =jpkm1 over the ocean point, =0 elsewhere 197 195 ! 198 196 k_top(:,:) = MIN( 1 , k_bot(:,:) ) ! = 1 over the ocean point, =0 elsewhere
Note: See TracChangeset
for help on using the changeset viewer.