Changeset 12614 for NEMO/branches
- Timestamp:
- 2020-03-26T15:59:52+01:00 (4 years ago)
- Location:
- NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater
- Files:
-
- 22 added
- 13 copied
-
cfgs/AM98 (added)
-
cfgs/AM98/EXPREF (added)
-
cfgs/AM98/EXPREF/axis_def_nemo.xml (added)
-
cfgs/AM98/EXPREF/communication_report.txt (added)
-
cfgs/AM98/EXPREF/context_nemo.xml (added)
-
cfgs/AM98/EXPREF/domain_def_nemo.xml (added)
-
cfgs/AM98/EXPREF/field_def_nemo-oce.xml (added)
-
cfgs/AM98/EXPREF/file_def_nemo-oce.xml (added)
-
cfgs/AM98/EXPREF/grid_def_nemo.xml (added)
-
cfgs/AM98/EXPREF/iodef.xml (added)
-
cfgs/AM98/EXPREF/layout.dat (added)
-
cfgs/AM98/EXPREF/namelist_cfg (added)
-
cfgs/AM98/EXPREF/namelist_ref (added)
-
cfgs/AM98/MY_SRC (added)
-
cfgs/AM98/MY_SRC/usrdef_fmask.F90 (added)
-
cfgs/AM98/MY_SRC/usrdef_hgr.F90 (added)
-
cfgs/AM98/MY_SRC/usrdef_istate.F90 (added)
-
cfgs/AM98/MY_SRC/usrdef_nam.F90 (added)
-
cfgs/AM98/MY_SRC/usrdef_sbc.F90 (added)
-
cfgs/AM98/MY_SRC/usrdef_zgr.F90 (added)
-
cfgs/AM98/cpp_AM98.fcm (added)
-
src/SWE (added)
-
src/SWE/asminc.F90 (copied) (copied from NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/OCE/ASM/asminc.F90) (1 diff)
-
src/SWE/diawri.F90 (copied) (copied from NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/OCE/DIA/diawri.F90) (3 diffs)
-
src/SWE/dom_oce.F90 (copied) (copied from NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/OCE/DOM/dom_oce.F90) (4 diffs)
-
src/SWE/domain.F90 (copied) (copied from NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/OCE/DOM/domain.F90) (3 diffs)
-
src/SWE/dommsk.F90 (copied) (copied from NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/OCE/DOM/dommsk.F90) (1 diff)
-
src/SWE/domvvl.F90 (copied) (copied from NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/OCE/DOM/domvvl.F90) (7 diffs)
-
src/SWE/dynatf.F90 (copied) (copied from NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/OCE/DYN/dynatf.F90) (2 diffs)
-
src/SWE/dynvor.F90 (copied) (copied from NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/OCE/DYN/dynvor.F90) (1 diff)
-
src/SWE/nemogcm.F90 (copied) (copied from NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/OCE/nemogcm.F90) (5 diffs)
-
src/SWE/phycst.F90 (copied) (copied from NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/OCE/DOM/phycst.F90) (3 diffs)
-
src/SWE/sbcice_cice.F90 (copied) (copied from NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/OCE/SBC/sbcice_cice.F90) (1 diff)
-
src/SWE/step.F90 (copied) (copied from NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/OCE/step.F90) (12 diffs)
-
src/SWE/stpctl.F90 (copied) (copied from NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/OCE/stpctl.F90) (6 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/asminc.F90
r12529 r12614 896 896 IF ( kt == nitdin_r ) THEN 897 897 ! 898 l_1st_euler = 0! Force Euler forward step898 l_1st_euler = .TRUE. ! Force Euler forward step 899 899 ! 900 900 ! Sea-ice : SI3 case -
NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/diawri.F90
r12529 r12614 119 119 REAL(wp):: zztmp , zztmpx ! local scalar 120 120 REAL(wp):: zztmp2, zztmpy ! - - 121 REAL(wp):: ze3 ! - - 121 122 REAL(wp), DIMENSION(jpi,jpj) :: z2d ! 2D workspace 122 123 REAL(wp), DIMENSION(jpi,jpj,jpk) :: z3d ! 3D workspace … … 146 147 CALL iom_put( "ssh" , (ssh(:,:,Kmm)+ssh_ref)*tmask(:,:,1) ) ! sea surface height (brought back to the reference used for wetting and drying) 147 148 ELSE 148 CALL iom_put( "ssh" , ssh(:,:,Kmm) ) ! sea surface height 149 ENDIF 149 CALL iom_put( "ssh" , ssh(:,:,Kmm) ) ! sea surface height 150 ENDIF 151 152 !!an 153 IF( iom_use("ht") ) & ! water column at t-point 154 CALL iom_put( "ht" , ht_0(:,:) + ssh(:,:,Kmm) ) 155 ! 156 IF( iom_use("hu") ) THEN ! water column at u-point 157 z2d(:,:) = 0._wp 158 DO_2D_10_10 159 z2d(ji,jj) = 0.5_wp * ( e3t(ji ,jj,1,Kmm) * e1e2t(ji ,jj) & 160 & + e3t(ji+1,jj,1,Kmm) * e1e2t(ji+1,jj) ) * r1_e1e2u(ji,jj) 161 END_2D 162 CALL lbc_lnk( 'diawri', z2d, 'U', 1._wp ) 163 CALL iom_put( "hu", z2d ) 164 ENDIF 165 ! 166 IF( iom_use("hv") ) THEN ! water column at v-point 167 z2d(:,:) = 0._wp 168 DO_2D_10_10 169 z2d(ji,jj) = 0.5_wp * ( e3t(ji,jj+1,1,Kmm) * e1e2t(ji,jj+1) & 170 & + e3t(ji,jj ,1,Kmm) * e1e2t(ji,jj ) ) * r1_e1e2v(ji,jj) 171 END_2D 172 CALL lbc_lnk( 'diawri', z2d, 'V', 1._wp ) 173 CALL iom_put( "hv", z2d ) 174 ENDIF 175 ! 176 IF( iom_use("hf") ) THEN ! water column at f-point 177 z2d(:,:) = 0._wp 178 DO_2D_10_10 179 z2d(ji,jj) = 0.25_wp * ( e3t(ji,jj+1,1,Kmm) * e1e2t(ji,jj+1) + e3t(ji+1,jj+1,1,Kmm) * e1e2t(ji+1,jj+1) & 180 & + e3t(ji,jj ,1,Kmm) * e1e2t(ji,jj ) + e3t(ji+1,jj ,1,Kmm) * e1e2t(ji+1,jj ) ) * r1_e1e2f(ji,jj) 181 END_2D 182 CALL lbc_lnk( 'diawri', z2d, 'F', 1._wp ) 183 CALL iom_put( "hf", z2d ) 184 ENDIF 185 !!an 150 186 151 187 IF( iom_use("wetdep") ) & ! wet depth 152 188 CALL iom_put( "wetdep" , ht_0(:,:) + ssh(:,:,Kmm) ) 153 154 CALL iom_put( "toce", ts(:,:,:,jp_tem,Kmm) ) ! 3D temperature155 CALL iom_put( "sst", ts(:,:,1,jp_tem,Kmm) ) ! surface temperature156 IF ( iom_use("sbt") ) THEN157 DO_2D_11_11158 ikbot = mbkt(ji,jj)159 z2d(ji,jj) = ts(ji,jj,ikbot,jp_tem,Kmm)160 END_2D161 CALL iom_put( "sbt", z2d ) ! bottom temperature162 ENDIF163 164 CALL iom_put( "soce", ts(:,:,:,jp_sal,Kmm) ) ! 3D salinity165 CALL iom_put( "sss", ts(:,:,1,jp_sal,Kmm) ) ! surface salinity166 IF ( iom_use("sbs") ) THEN167 DO_2D_11_11168 ikbot = mbkt(ji,jj)169 z2d(ji,jj) = ts(ji,jj,ikbot,jp_sal,Kmm)170 END_2D171 CALL iom_put( "sbs", z2d ) ! bottom salinity172 ENDIF173 189 174 190 IF ( iom_use("taubot") ) THEN ! bottom stress … … 186 202 CALL iom_put( "taubot", z2d ) 187 203 ENDIF 188 189 CALL iom_put( "uoce", uu(:,:,:,Kmm) ) ! 3D i-current 190 CALL iom_put( "ssu", uu(:,:,1,Kmm) ) ! surface i-current 191 IF ( iom_use("sbu") ) THEN 192 DO_2D_11_11 193 ikbot = mbku(ji,jj) 194 z2d(ji,jj) = uu(ji,jj,ikbot,Kmm) 204 ! 205 CALL iom_put( "ssu" , uu(:,:,1,Kmm) ) ! surface i-current 206 CALL iom_put( "ssv" , vv(:,:,1,Kmm) ) ! surface j-current 207 CALL iom_put( "woce", ww ) ! vertical velocity 208 ! 209 210 IF ( iom_use("sKE") ) THEN ! surface kinetic energy at T point 211 z2d(:,:) = 0._wp 212 DO_2D_00_00 213 z2d(ji,jj) = 0.25_wp * ( uu(ji ,jj,1,Kmm) * uu(ji ,jj,1,Kmm) * e1e2u(ji ,jj) * e3u(ji ,jj,1,Kmm) & 214 & + uu(ji-1,jj,1,Kmm) * uu(ji-1,jj,1,Kmm) * e1e2u(ji-1,jj) * e3u(ji-1,jj,1,Kmm) & 215 & + vv(ji,jj ,1,Kmm) * vv(ji,jj ,1,Kmm) * e1e2v(ji,jj ) * e3v(ji,jj ,1,Kmm) & 216 & + vv(ji,jj-1,1,Kmm) * vv(ji,jj-1,1,Kmm) * e1e2v(ji,jj-1) * e3v(ji,jj-1,1,Kmm) ) & 217 & * r1_e1e2t(ji,jj) / e3t(ji,jj,1,Kmm) * tmask(ji,jj,1) 195 218 END_2D 196 CALL iom_put( "sbu", z2d ) ! bottom i-current 219 ! 220 CALL lbc_lnk( 'diawri', z2d, 'T', 1. ) 221 IF ( iom_use("sKE" ) ) CALL iom_put( "sKE" , z2d ) 222 223 ENDIF 224 225 !!an sKEf 226 IF ( iom_use("sKEf") ) THEN ! surface kinetic energy at F point 227 z2d(:,:) = 0._wp ! CAUTION : only valid in SWE, not with bathymetry 228 DO_2D_00_00 229 z2d(ji,jj) = 0.25_wp * ( uu(ji,jj ,1,Kmm) * uu(ji,jj ,1,Kmm) * e1e2u(ji,jj ) * e3u(ji,jj ,1,Kmm) & 230 & + uu(ji,jj+1,1,Kmm) * uu(ji,jj+1,1,Kmm) * e1e2u(ji,jj+1) * e3u(ji,jj+1,1,Kmm) & 231 & + vv(ji ,jj,1,Kmm) * vv(ji,jj ,1,Kmm) * e1e2v(ji ,jj) * e3v(ji ,jj,1,Kmm) & 232 & + vv(ji+1,jj,1,Kmm) * vv(ji+1,jj,1,Kmm) * e1e2v(ji+1,jj) * e3v(ji+1,jj,1,Kmm) ) & 233 & * r1_e1e2f(ji,jj) / e3f(ji,jj,1) * ssfmask(ji,jj) 234 END_2D 235 ! 236 CALL lbc_lnk( 'diawri', z2d, 'F', 1. ) 237 CALL iom_put( "sKEf", z2d ) 238 ENDIF 239 !!an 240 ! 241 CALL iom_put( "hdiv", hdiv ) ! Horizontal divergence 242 ! 243 ! Output of vorticity terms 244 IF ( iom_use("relvor") .OR. iom_use("plavor") .OR. & 245 & iom_use("relpotvor") .OR. iom_use("abspotvor") .OR. & 246 & iom_use("Ens") ) THEN 247 ! 248 z2d(:,:) = 0._wp 249 ze3 = 0._wp 250 DO_2D_10_10 251 z2d(ji,jj) = ( e2v(ji+1,jj ) * vv(ji+1,jj ,1,Kmm) - e2v(ji,jj) * vv(ji,jj,1,Kmm) & 252 & - e1u(ji ,jj+1) * uu(ji ,jj+1,1,Kmm) + e1u(ji,jj) * uu(ji,jj,1,Kmm) ) * r1_e1e2f(ji,jj) 253 END_2D 254 CALL lbc_lnk( 'diawri', z2d, 'F', 1. ) 255 CALL iom_put( "relvor", z2d ) ! relative vorticity ( zeta ) 256 ! 257 CALL iom_put( "plavor", ff_f ) ! planetary vorticity ( f ) 258 ! 259 DO_2D_10_10 260 ze3 = ( e3t(ji,jj+1,1,Kmm) * e1e2t(ji,jj+1) + e3t(ji+1,jj+1,1,Kmm) * e1e2t(ji+1,jj+1) & 261 & + e3t(ji,jj ,1,Kmm) * e1e2t(ji,jj ) + e3t(ji+1,jj ,1,Kmm) * e1e2t(ji+1,jj ) ) * r1_e1e2f(ji,jj) 262 IF( ze3 /= 0._wp ) THEN ; ze3 = 4._wp / ze3 263 ELSE ; ze3 = 0._wp 264 ENDIF 265 z2d(ji,jj) = ze3 * z2d(ji,jj) 266 END_2D 267 CALL lbc_lnk( 'diawri', z2d, 'F', 1. ) 268 CALL iom_put( "relpotvor", z2d ) ! relative potential vorticity (zeta/h) 269 ! 270 DO_2D_10_10 271 ze3 = ( e3t(ji,jj+1,1,Kmm) * e1e2t(ji,jj+1) + e3t(ji+1,jj+1,1,Kmm) * e1e2t(ji+1,jj+1) & 272 & + e3t(ji,jj ,1,Kmm) * e1e2t(ji,jj ) + e3t(ji+1,jj ,1,Kmm) * e1e2t(ji+1,jj ) ) * r1_e1e2f(ji,jj) 273 IF( ze3 /= 0._wp ) THEN ; ze3 = 4._wp / ze3 274 ELSE ; ze3 = 0._wp 275 ENDIF 276 z2d(ji,jj) = ze3 * ff_f(ji,jj) + z2d(ji,jj) 277 END_2D 278 CALL lbc_lnk( 'diawri', z2d, 'F', 1. ) 279 CALL iom_put( "abspotvor", z2d ) ! absolute potential vorticity ( q ) 280 ! 281 DO_2D_10_10 282 z2d(ji,jj) = 0.5_wp * z2d(ji,jj) * z2d(ji,jj) 283 END_2D 284 CALL lbc_lnk( 'diawri', z2d, 'F', 1. ) 285 CALL iom_put( "Ens", z2d ) ! potential enstrophy ( 1/2*q2 ) 286 ! 197 287 ENDIF 198 288 199 CALL iom_put( "voce", vv(:,:,:,Kmm) ) ! 3D j-current 200 CALL iom_put( "ssv", vv(:,:,1,Kmm) ) ! surface j-current 201 IF ( iom_use("sbv") ) THEN 202 DO_2D_11_11 203 ikbot = mbkv(ji,jj) 204 z2d(ji,jj) = vv(ji,jj,ikbot,Kmm) 205 END_2D 206 CALL iom_put( "sbv", z2d ) ! bottom j-current 207 ENDIF 208 209 IF( ln_zad_Aimp ) ww = ww + wi ! Recombine explicit and implicit parts of vertical velocity for diagnostic output 210 ! 211 CALL iom_put( "woce", ww ) ! vertical velocity 212 IF( iom_use('w_masstr') .OR. iom_use('w_masstr2') ) THEN ! vertical mass transport & its square value 213 ! Caution: in the VVL case, it only correponds to the baroclinic mass transport. 214 z2d(:,:) = rho0 * e1e2t(:,:) 215 DO jk = 1, jpk 216 z3d(:,:,jk) = ww(:,:,jk) * z2d(:,:) 217 END DO 218 CALL iom_put( "w_masstr" , z3d ) 219 IF( iom_use('w_masstr2') ) CALL iom_put( "w_masstr2", z3d(:,:,:) * z3d(:,:,:) ) 220 ENDIF 221 ! 222 IF( ln_zad_Aimp ) ww = ww - wi ! Remove implicit part of vertical velocity that was added for diagnostic output 223 224 CALL iom_put( "avt" , avt ) ! T vert. eddy diff. coef. 225 CALL iom_put( "avs" , avs ) ! S vert. eddy diff. coef. 226 CALL iom_put( "avm" , avm ) ! T vert. eddy visc. coef. 227 228 IF( iom_use('logavt') ) CALL iom_put( "logavt", LOG( MAX( 1.e-20_wp, avt(:,:,:) ) ) ) 229 IF( iom_use('logavs') ) CALL iom_put( "logavs", LOG( MAX( 1.e-20_wp, avs(:,:,:) ) ) ) 230 231 IF ( iom_use("sstgrad") .OR. iom_use("sstgrad2") ) THEN 232 DO_2D_00_00 233 zztmp = ts(ji,jj,1,jp_tem,Kmm) 234 zztmpx = ( ts(ji+1,jj,1,jp_tem,Kmm) - zztmp ) * r1_e1u(ji,jj) + ( zztmp - ts(ji-1,jj ,1,jp_tem,Kmm) ) * r1_e1u(ji-1,jj) 235 zztmpy = ( ts(ji,jj+1,1,jp_tem,Kmm) - zztmp ) * r1_e2v(ji,jj) + ( zztmp - ts(ji ,jj-1,1,jp_tem,Kmm) ) * r1_e2v(ji,jj-1) 236 z2d(ji,jj) = 0.25 * ( zztmpx * zztmpx + zztmpy * zztmpy ) & 237 & * umask(ji,jj,1) * umask(ji-1,jj,1) * vmask(ji,jj,1) * umask(ji,jj-1,1) 238 END_2D 239 CALL lbc_lnk( 'diawri', z2d, 'T', 1. ) 240 CALL iom_put( "sstgrad2", z2d ) ! square of module of sst gradient 241 z2d(:,:) = SQRT( z2d(:,:) ) 242 CALL iom_put( "sstgrad" , z2d ) ! module of sst gradient 243 ENDIF 244 245 ! heat and salt contents 246 IF( iom_use("heatc") ) THEN 247 z2d(:,:) = 0._wp 248 DO_3D_11_11( 1, jpkm1 ) 249 z2d(ji,jj) = z2d(ji,jj) + e3t(ji,jj,jk,Kmm) * ts(ji,jj,jk,jp_tem,Kmm) * tmask(ji,jj,jk) 250 END_3D 251 CALL iom_put( "heatc", rho0_rcp * z2d ) ! vertically integrated heat content (J/m2) 252 ENDIF 253 254 IF( iom_use("saltc") ) THEN 255 z2d(:,:) = 0._wp 256 DO_3D_11_11( 1, jpkm1 ) 257 z2d(ji,jj) = z2d(ji,jj) + e3t(ji,jj,jk,Kmm) * ts(ji,jj,jk,jp_sal,Kmm) * tmask(ji,jj,jk) 258 END_3D 259 CALL iom_put( "saltc", rho0 * z2d ) ! vertically integrated salt content (PSU*kg/m2) 260 ENDIF 261 ! 262 IF ( iom_use("eken") ) THEN 263 z3d(:,:,jpk) = 0._wp 264 DO_3D_00_00( 1, jpkm1 ) 265 zztmp = 0.25_wp * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm) 266 z3d(ji,jj,jk) = zztmp * ( uu(ji-1,jj,jk,Kmm)**2 * e2u(ji-1,jj) * e3u(ji-1,jj,jk,Kmm) & 267 & + uu(ji ,jj,jk,Kmm)**2 * e2u(ji ,jj) * e3u(ji ,jj,jk,Kmm) & 268 & + vv(ji,jj-1,jk,Kmm)**2 * e1v(ji,jj-1) * e3v(ji,jj-1,jk,Kmm) & 269 & + vv(ji,jj ,jk,Kmm)**2 * e1v(ji,jj ) * e3v(ji,jj ,jk,Kmm) ) 270 END_3D 271 CALL lbc_lnk( 'diawri', z3d, 'T', 1. ) 272 CALL iom_put( "eken", z3d ) ! kinetic energy 273 ENDIF 274 ! 275 CALL iom_put( "hdiv", hdiv ) ! Horizontal divergence 276 ! 277 IF( iom_use("u_masstr") .OR. iom_use("u_masstr_vint") .OR. iom_use("u_heattr") .OR. iom_use("u_salttr") ) THEN 278 z3d(:,:,jpk) = 0.e0 279 z2d(:,:) = 0.e0 280 DO jk = 1, jpkm1 281 z3d(:,:,jk) = rho0 * uu(:,:,jk,Kmm) * e2u(:,:) * e3u(:,:,jk,Kmm) * umask(:,:,jk) 282 z2d(:,:) = z2d(:,:) + z3d(:,:,jk) 283 END DO 284 CALL iom_put( "u_masstr" , z3d ) ! mass transport in i-direction 285 CALL iom_put( "u_masstr_vint", z2d ) ! mass transport in i-direction vertical sum 286 ENDIF 287 288 IF( iom_use("u_heattr") ) THEN 289 z2d(:,:) = 0._wp 290 DO_3D_00_00( 1, jpkm1 ) 291 z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( ts(ji,jj,jk,jp_tem,Kmm) + ts(ji+1,jj,jk,jp_tem,Kmm) ) 292 END_3D 293 CALL lbc_lnk( 'diawri', z2d, 'U', -1. ) 294 CALL iom_put( "u_heattr", 0.5*rcp * z2d ) ! heat transport in i-direction 295 ENDIF 296 297 IF( iom_use("u_salttr") ) THEN 298 z2d(:,:) = 0.e0 299 DO_3D_00_00( 1, jpkm1 ) 300 z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( ts(ji,jj,jk,jp_sal,Kmm) + ts(ji+1,jj,jk,jp_sal,Kmm) ) 301 END_3D 302 CALL lbc_lnk( 'diawri', z2d, 'U', -1. ) 303 CALL iom_put( "u_salttr", 0.5 * z2d ) ! heat transport in i-direction 304 ENDIF 305 306 307 IF( iom_use("v_masstr") .OR. iom_use("v_heattr") .OR. iom_use("v_salttr") ) THEN 308 z3d(:,:,jpk) = 0.e0 309 DO jk = 1, jpkm1 310 z3d(:,:,jk) = rho0 * vv(:,:,jk,Kmm) * e1v(:,:) * e3v(:,:,jk,Kmm) * vmask(:,:,jk) 311 END DO 312 CALL iom_put( "v_masstr", z3d ) ! mass transport in j-direction 313 ENDIF 314 315 IF( iom_use("v_heattr") ) THEN 316 z2d(:,:) = 0.e0 317 DO_3D_00_00( 1, jpkm1 ) 318 z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( ts(ji,jj,jk,jp_tem,Kmm) + ts(ji,jj+1,jk,jp_tem,Kmm) ) 319 END_3D 320 CALL lbc_lnk( 'diawri', z2d, 'V', -1. ) 321 CALL iom_put( "v_heattr", 0.5*rcp * z2d ) ! heat transport in j-direction 322 ENDIF 323 324 IF( iom_use("v_salttr") ) THEN 325 z2d(:,:) = 0._wp 326 DO_3D_00_00( 1, jpkm1 ) 327 z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( ts(ji,jj,jk,jp_sal,Kmm) + ts(ji,jj+1,jk,jp_sal,Kmm) ) 328 END_3D 329 CALL lbc_lnk( 'diawri', z2d, 'V', -1. ) 330 CALL iom_put( "v_salttr", 0.5 * z2d ) ! heat transport in j-direction 331 ENDIF 332 333 IF( iom_use("tosmint") ) THEN 334 z2d(:,:) = 0._wp 335 DO_3D_00_00( 1, jpkm1 ) 336 z2d(ji,jj) = z2d(ji,jj) + e3t(ji,jj,jk,Kmm) * ts(ji,jj,jk,jp_tem,Kmm) 337 END_3D 338 CALL lbc_lnk( 'diawri', z2d, 'T', -1. ) 339 CALL iom_put( "tosmint", rho0 * z2d ) ! Vertical integral of temperature 340 ENDIF 341 IF( iom_use("somint") ) THEN 342 z2d(:,:)=0._wp 343 DO_3D_00_00( 1, jpkm1 ) 344 z2d(ji,jj) = z2d(ji,jj) + e3t(ji,jj,jk,Kmm) * ts(ji,jj,jk,jp_sal,Kmm) 345 END_3D 346 CALL lbc_lnk( 'diawri', z2d, 'T', -1. ) 347 CALL iom_put( "somint", rho0 * z2d ) ! Vertical integral of salinity 348 ENDIF 349 350 CALL iom_put( "bn2", rn2 ) ! Brunt-Vaisala buoyancy frequency (N^2) 351 ! 352 353 IF (ln_dia25h) CALL dia_25h( kt, Kmm ) ! 25h averaging 354 289 ! 355 290 IF( ln_timing ) CALL timing_stop('dia_wri') 356 291 ! -
NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/dom_oce.F90
r12529 r12614 151 151 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hu_0 !: u-depth [m] 152 152 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hv_0 !: v-depth [m] 153 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hf_0 !: f-depth [m] 154 ! ! inverse of reference heights of water column 155 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: r1_ht_0 !: t-depth [1/m] 156 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: r1_hu_0 !: u-depth [1/m] 157 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: r1_hv_0 !: v-depth [1/m] 158 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: r1_hf_0 !: f-depth [1/m] 159 153 160 ! time-dependent heights of water column 154 161 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ht !: height of water column at T-points [m] … … 178 185 INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: mikt, miku, mikv, mikf !: top first wet T-, U-, V-, F-level (ISF) 179 186 180 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ssmask, ssumask, ssvmask 187 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ssmask, ssumask, ssvmask, ssfmask !: surface mask at T-,U-, V- and F-pts 181 188 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:), TARGET :: tmask, umask, vmask, fmask !: land/ocean mask at T-, U-, V- and F-pts 182 189 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:), TARGET :: wmask, wumask, wvmask !: land/ocean mask at WT-, WU- and WV-pts … … 268 275 & e3uw_0(jpi,jpj,jpk) , e3vw_0(jpi,jpj,jpk) , & 269 276 & e3uw (jpi,jpj,jpk,jpt) , e3vw (jpi,jpj,jpk,jpt) , STAT=ierr(5) ) 270 ! 271 ALLOCATE( ht_0(jpi,jpj) , hu_0(jpi,jpj) , hv_0(jpi,jpj) , & 272 & ht (jpi,jpj) , hu( jpi,jpj,jpt), hv( jpi,jpj,jpt) , r1_hu(jpi,jpj,jpt) , r1_hv(jpi,jpj,jpt) , & 273 & STAT=ierr(6) ) 277 ! 278 ALLOCATE( ht_0(jpi,jpj) , hu_0(jpi,jpj) , hv_0(jpi,jpj) , hf_0(jpi,jpj) , & 279 & r1_ht_0(jpi,jpj) , r1_hu_0(jpi,jpj) , r1_hv_0(jpi,jpj) , r1_hf_0(jpi,jpj) , & 280 & ht (jpi,jpj) , hu (jpi,jpj,jpt) , hv (jpi,jpj,jpt) , & 281 & r1_hu (jpi,jpj,jpt) , r1_hv (jpi,jpj,jpt) , STAT=ierr(6) ) 274 282 ! 275 283 ALLOCATE( risfdep(jpi,jpj) , bathy(jpi,jpj) , STAT=ierr(7) ) … … 277 285 ALLOCATE( gdept_1d(jpk) , gdepw_1d(jpk) , e3t_1d(jpk) , e3w_1d(jpk) , STAT=ierr(8) ) 278 286 ! 279 ALLOCATE( tmask_i(jpi,jpj) , tmask_h(jpi,jpj) , &280 & ssmask (jpi,jpj) , ssumask(jpi,jpj) , ssvmask(jpi,jpj) , &281 & mbkt (jpi,jpj) , mbku (jpi,jpj) , mbkv (jpi,jpj) , STAT=ierr(9) )287 ALLOCATE( tmask_i(jpi,jpj) , tmask_h(jpi,jpj) , & 288 & ssmask (jpi,jpj) , ssumask(jpi,jpj) , ssvmask(jpi,jpj) , ssfmask(jpi,jpj) , & 289 & mbkt (jpi,jpj) , mbku (jpi,jpj) , mbkv (jpi,jpj) , STAT=ierr(9) ) 282 290 ! 283 291 ALLOCATE( mikt(jpi,jpj), miku(jpi,jpj), mikv(jpi,jpj), mikf(jpi,jpj), STAT=ierr(10) ) -
NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/domain.F90
r12529 r12614 143 143 144 144 CALL dom_msk( ik_top, ik_bot ) ! Masks 145 ! 146 ht_0(:,:) = 0._wp ! Reference ocean thickness145 146 ht_0(:,:) = 0._wp ! Reference ocean thickness 147 147 hu_0(:,:) = 0._wp 148 148 hv_0(:,:) = 0._wp 149 DO jk = 1, jpk 149 hf_0(:,:) = 0._wp 150 DO jk = 1, jpkm1 150 151 ht_0(:,:) = ht_0(:,:) + e3t_0(:,:,jk) * tmask(:,:,jk) 151 152 hu_0(:,:) = hu_0(:,:) + e3u_0(:,:,jk) * umask(:,:,jk) 152 153 hv_0(:,:) = hv_0(:,:) + e3v_0(:,:,jk) * vmask(:,:,jk) 154 hf_0(:,:) = hf_0(:,:) + e3f_0(:,:,jk) * ssfmask(:,:) ! CAUTION : only valid in SWE, not with bathymetry 153 155 END DO 154 ! 156 ! ! Inverse of reference ocean thickness 157 r1_ht_0(:,:) = ssmask(:,:) / ( ht_0(:,:) + 1._wp - ssmask(:,:) ) 158 r1_hu_0(:,:) = ssumask(:,:) / ( hu_0(:,:) + 1._wp - ssumask(:,:) ) 159 r1_hv_0(:,:) = ssvmask(:,:) / ( hv_0(:,:) + 1._wp - ssvmask(:,:) ) 160 r1_hf_0(:,:) = ssfmask(:,:) / ( hf_0(:,:) + 1._wp - ssfmask(:,:) ) 161 155 162 ! !== time varying part of coordinate system ==! 156 163 ! … … 190 197 ! 191 198 IF( ln_meshmask ) CALL dom_wri ! Create a domain file 199 192 200 IF( .NOT.ln_rstart ) CALL dom_ctl ! Domain control 193 201 ! … … 200 208 WRITE(numout,*) 201 209 ENDIF 210 202 211 ! 203 212 END SUBROUTINE dom_init -
NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/dommsk.F90
r12529 r12614 192 192 ssumask(:,:) = MAXVAL( umask(:,:,:), DIM=3 ) 193 193 ssvmask(:,:) = MAXVAL( vmask(:,:,:), DIM=3 ) 194 194 !!an 195 ! ssfmask(:,:) = MAXVAL( fmask(:,:,:), DIM=3 ) 196 DO_2D_10_10 197 ssfmask(ji,jj) = MAX( tmask(ji,jj+1,1), tmask(ji+1,jj+1,1), & 198 & tmask(ji,jj ,1), tmask(ji+1,jj ,1) ) 199 END_2D 200 CALL lbc_lnk( 'dommsk', ssfmask, 'F', 1._wp ) 201 202 !!an 195 203 196 204 ! Interior domain mask (used for global sum) -
NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/domvvl.F90
r12529 r12614 166 166 ! !== Set of all other vertical scale factors ==! (now and before) 167 167 ! ! Horizontal interpolation of e3t 168 CALL dom_vvl_interpol( e3t(:,:,:,Kbb), e3u(:,:,:,Kbb), 'U' ) ! from T to U169 CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3u(:,:,:,Kmm), 'U' )170 CALL dom_vvl_interpol( e3t(:,:,:,Kbb), e3v(:,:,:,Kbb), 'V' ) ! from T to V171 CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3v(:,:,:,Kmm), 'V' )172 CALL dom_vvl_interpol( e3u(:,:,:,Kmm), e3f(:,:,:), 'F' ) ! from U to F168 CALL dom_vvl_interpol( ssh(:,:,Kbb), e3u(:,:,:,Kbb), 'U' ) ! from T to U 169 CALL dom_vvl_interpol( ssh(:,:,Kmm), e3u(:,:,:,Kmm), 'U' ) 170 CALL dom_vvl_interpol( ssh(:,:,Kbb), e3v(:,:,:,Kbb), 'V' ) ! from T to V 171 CALL dom_vvl_interpol( ssh(:,:,Kmm), e3v(:,:,:,Kmm), 'V' ) 172 CALL dom_vvl_interpol( ssh(:,:,Kmm), e3f(:,:,:), 'F' ) ! from U to F 173 173 ! ! Vertical interpolation of e3t,u,v 174 CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3w (:,:,:,Kmm), 'W' ) ! from T to W175 CALL dom_vvl_interpol( e3t(:,:,:,Kbb), e3w (:,:,:,Kbb), 'W' )176 CALL dom_vvl_interpol( e3u(:,:,:,Kmm), e3uw(:,:,:,Kmm), 'UW' ) ! from U to UW177 CALL dom_vvl_interpol( e3u(:,:,:,Kbb), e3uw(:,:,:,Kbb), 'UW' )178 CALL dom_vvl_interpol( e3v(:,:,:,Kmm), e3vw(:,:,:,Kmm), 'VW' ) ! from V to UW179 CALL dom_vvl_interpol( e3v(:,:,:,Kbb), e3vw(:,:,:,Kbb), 'VW' )174 CALL dom_vvl_interpol( ssh(:,:,Kmm), e3w (:,:,:,Kmm), 'W' ) ! from T to W 175 CALL dom_vvl_interpol( ssh(:,:,Kbb), e3w (:,:,:,Kbb), 'W' ) 176 CALL dom_vvl_interpol( ssh(:,:,Kmm), e3uw(:,:,:,Kmm), 'UW' ) ! from U to UW 177 CALL dom_vvl_interpol( ssh(:,:,Kbb), e3uw(:,:,:,Kbb), 'UW' ) 178 CALL dom_vvl_interpol( ssh(:,:,Kmm), e3vw(:,:,:,Kmm), 'VW' ) ! from V to UW 179 CALL dom_vvl_interpol( ssh(:,:,Kbb), e3vw(:,:,:,Kbb), 'VW' ) 180 180 181 181 ! We need to define e3[tuv]_a for AGRIF initialisation (should not be a problem for the restartability...) … … 549 549 ! *********************************** ! 550 550 551 CALL dom_vvl_interpol( e3t(:,:,:,Kaa), e3u(:,:,:,Kaa), 'U' )552 CALL dom_vvl_interpol( e3t(:,:,:,Kaa), e3v(:,:,:,Kaa), 'V' )551 CALL dom_vvl_interpol( ssh(:,:,Kaa), e3u(:,:,:,Kaa), 'U' ) 552 CALL dom_vvl_interpol( ssh(:,:,Kaa), e3v(:,:,:,Kaa), 'V' ) 553 553 554 554 ! *********************************** ! … … 633 633 ! - JC - hu(:,:,:,Kbb), hv(:,:,:,:,Kbb), hur_b, hvr_b also 634 634 635 CALL dom_vvl_interpol( e3u(:,:,:,Kmm), e3f(:,:,:), 'F' )635 CALL dom_vvl_interpol( ssh(:,:,Kmm), e3f(:,:,:), 'F' ) 636 636 637 637 ! Vertical scale factor interpolations 638 CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3w(:,:,:,Kmm), 'W' )639 CALL dom_vvl_interpol( e3u(:,:,:,Kmm), e3uw(:,:,:,Kmm), 'UW' )640 CALL dom_vvl_interpol( e3v(:,:,:,Kmm), e3vw(:,:,:,Kmm), 'VW' )641 CALL dom_vvl_interpol( e3t(:,:,:,Kbb), e3w(:,:,:,Kbb), 'W' )642 CALL dom_vvl_interpol( e3u(:,:,:,Kbb), e3uw(:,:,:,Kbb), 'UW' )643 CALL dom_vvl_interpol( e3v(:,:,:,Kbb), e3vw(:,:,:,Kbb), 'VW' )638 CALL dom_vvl_interpol( ssh(:,:,Kmm), e3w(:,:,:,Kmm), 'W' ) 639 CALL dom_vvl_interpol( ssh(:,:,Kmm), e3uw(:,:,:,Kmm), 'UW' ) 640 CALL dom_vvl_interpol( ssh(:,:,Kmm), e3vw(:,:,:,Kmm), 'VW' ) 641 CALL dom_vvl_interpol( ssh(:,:,Kbb), e3w(:,:,:,Kbb), 'W' ) 642 CALL dom_vvl_interpol( ssh(:,:,Kbb), e3uw(:,:,:,Kbb), 'UW' ) 643 CALL dom_vvl_interpol( ssh(:,:,Kbb), e3vw(:,:,:,Kbb), 'VW' ) 644 644 645 645 ! t- and w- points depth (set the isf depth as it is in the initial step) … … 674 674 675 675 676 SUBROUTINE dom_vvl_interpol( p e3_in, pe3_out, pout)676 SUBROUTINE dom_vvl_interpol( pssh, pe3, cdp ) 677 677 !!--------------------------------------------------------------------- 678 678 !! *** ROUTINE dom_vvl__interpol *** … … 684 684 !! - vertical interpolation: simple averaging 685 685 !!---------------------------------------------------------------------- 686 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in ) :: pe3_in ! input e3 to be interpolated 687 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pe3_out ! output interpolated e3 688 CHARACTER(LEN=*) , INTENT(in ) :: pout ! grid point of out scale factors 689 ! ! = 'U', 'V', 'W, 'F', 'UW' or 'VW' 690 ! 691 INTEGER :: ji, jj, jk ! dummy loop indices 692 REAL(wp) :: zlnwd ! =1./0. when ln_wd_il = T/F 693 !!---------------------------------------------------------------------- 694 ! 695 IF(ln_wd_il) THEN 696 zlnwd = 1.0_wp 697 ELSE 698 zlnwd = 0.0_wp 699 END IF 700 ! 701 SELECT CASE ( pout ) !== type of interpolation ==! 686 REAL(wp), DIMENSION(:,:) , INTENT(in ) :: pssh ! input e3 NOT used here (ssh is used instead) 687 REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: pe3 ! scale factor e3 to be updated [m] 688 CHARACTER(LEN=*) , INTENT(in ) :: cdp ! grid point of the scale factor ( 'U', 'V', 'W, 'F', 'UW' or 'VW' ) 689 ! 690 INTEGER :: ji, jj, jk ! dummy loop indices 691 REAL(wp), DIMENSION(jpi,jpj) :: zc3 ! 2D workspace 692 !!---------------------------------------------------------------------- 693 ! 694 SELECT CASE ( cdp ) !== type of interpolation ==! 702 695 ! 703 696 CASE( 'U' ) !* from T- to U-point : hor. surface weighted mean 704 DO_3D_10_10( 1, jpk ) 705 pe3_out(ji,jj,jk) = 0.5_wp * ( umask(ji,jj,jk) * (1.0_wp - zlnwd) + zlnwd ) * r1_e1e2u(ji,jj) & 706 & * ( e1e2t(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3t_0(ji ,jj,jk) ) & 707 & + e1e2t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) ) 708 END_3D 709 CALL lbc_lnk( 'domvvl', pe3_out(:,:,:), 'U', 1._wp ) 710 pe3_out(:,:,:) = pe3_out(:,:,:) + e3u_0(:,:,:) 697 DO_2D_00_00 698 zc3(ji,jj) = 0.5_wp * ( e1e2t(ji ,jj) * pssh(ji ,jj) & 699 & + e1e2t(ji+1,jj) * pssh(ji+1,jj) ) * r1_hu_0(ji,jj) * r1_e1e2u(ji,jj) 700 END_2D 701 CALL lbc_lnk( 'domvvl', zc3(:,:), 'U', 1._wp ) 702 ! 703 DO jk = 1, jpkm1 704 pe3(:,:,jk) = e3u_0(:,:,jk) * ( 1.0_wp + zc3(:,:) ) 705 END DO 711 706 ! 712 707 CASE( 'V' ) !* from T- to V-point : hor. surface weighted mean 713 DO_3D_10_10( 1, jpk ) 714 pe3_out(ji,jj,jk) = 0.5_wp * ( vmask(ji,jj,jk) * (1.0_wp - zlnwd) + zlnwd ) * r1_e1e2v(ji,jj) & 715 & * ( e1e2t(ji,jj ) * ( pe3_in(ji,jj ,jk) - e3t_0(ji,jj ,jk) ) & 716 & + e1e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) ) 717 END_3D 718 CALL lbc_lnk( 'domvvl', pe3_out(:,:,:), 'V', 1._wp ) 719 pe3_out(:,:,:) = pe3_out(:,:,:) + e3v_0(:,:,:) 708 DO_2D_00_00 709 zc3(ji,jj) = 0.5_wp * ( e1e2t(ji,jj ) * pssh(ji,jj ) & 710 & + e1e2t(ji,jj+1) * pssh(ji,jj+1) ) * r1_hv_0(ji,jj) * r1_e1e2v(ji,jj) 711 END_2D 712 CALL lbc_lnk( 'domvvl', zc3(:,:), 'V', 1._wp ) 713 ! 714 DO jk = 1, jpkm1 715 pe3(:,:,jk) = e3v_0(:,:,jk) * ( 1.0_wp + zc3(:,:) ) 716 END DO 720 717 ! 721 718 CASE( 'F' ) !* from U-point to F-point : hor. surface weighted mean 722 DO_3D_10_10( 1, jpk ) 723 pe3_out(ji,jj,jk) = 0.5_wp * ( umask(ji,jj,jk) * umask(ji,jj+1,jk) * (1.0_wp - zlnwd) + zlnwd ) & 724 & * r1_e1e2f(ji,jj) & 725 & * ( e1e2u(ji,jj ) * ( pe3_in(ji,jj ,jk) - e3u_0(ji,jj ,jk) ) & 726 & + e1e2u(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3u_0(ji,jj+1,jk) ) ) 727 END_3D 728 CALL lbc_lnk( 'domvvl', pe3_out(:,:,:), 'F', 1._wp ) 729 pe3_out(:,:,:) = pe3_out(:,:,:) + e3f_0(:,:,:) 719 DO_2D_10_10 720 zc3(ji,jj) = 0.25_wp * ( e1e2t(ji ,jj ) * pssh(ji ,jj ) & 721 & + e1e2t(ji+1,jj ) * pssh(ji+1,jj ) & 722 & + e1e2t(ji ,jj+1) * pssh(ji ,jj+1) & 723 & + e1e2t(ji+1,jj+1) * pssh(ji+1,jj+1) ) * r1_hf_0(ji,jj) * r1_e1e2f(ji,jj) 724 END_2D 725 CALL lbc_lnk( 'domvvl', zc3(:,:), 'F', 1._wp ) 726 ! 727 DO jk = 1, jpkm1 ! Horizontal interpolation of e3f from ssh 728 e3f(:,:,jk) = e3f_0(:,:,jk) * ( 1._wp + zc3(:,:) ) 729 END DO 730 730 ! 731 731 CASE( 'W' ) !* from T- to W-point : vertical simple mean 732 ! 733 pe3_out(:,:,1) = e3w_0(:,:,1) + pe3_in(:,:,1) - e3t_0(:,:,1) 734 ! - ML - The use of mask in this formulea enables the special treatment of the last w-point without indirect adressing 735 !!gm BUG? use here wmask in case of ISF ? to be checked 736 DO jk = 2, jpk 737 pe3_out(:,:,jk) = e3w_0(:,:,jk) + ( 1.0_wp - 0.5_wp * ( tmask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd ) ) & 738 & * ( pe3_in(:,:,jk-1) - e3t_0(:,:,jk-1) ) & 739 & + 0.5_wp * ( tmask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd ) & 740 & * ( pe3_in(:,:,jk ) - e3t_0(:,:,jk ) ) 741 END DO 742 ! 743 CASE( 'UW' ) !* from U- to UW-point : vertical simple mean 744 ! 745 pe3_out(:,:,1) = e3uw_0(:,:,1) + pe3_in(:,:,1) - e3u_0(:,:,1) 746 ! - ML - The use of mask in this formaula enables the special treatment of the last w- point without indirect adressing 747 !!gm BUG? use here wumask in case of ISF ? to be checked 748 DO jk = 2, jpk 749 pe3_out(:,:,jk) = e3uw_0(:,:,jk) + ( 1.0_wp - 0.5_wp * ( umask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd ) ) & 750 & * ( pe3_in(:,:,jk-1) - e3u_0(:,:,jk-1) ) & 751 & + 0.5_wp * ( umask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd ) & 752 & * ( pe3_in(:,:,jk ) - e3u_0(:,:,jk ) ) 753 END DO 754 ! 755 CASE( 'VW' ) !* from V- to VW-point : vertical simple mean 756 ! 757 pe3_out(:,:,1) = e3vw_0(:,:,1) + pe3_in(:,:,1) - e3v_0(:,:,1) 758 ! - ML - The use of mask in this formaula enables the special treatment of the last w- point without indirect adressing 759 !!gm BUG? use here wvmask in case of ISF ? to be checked 760 DO jk = 2, jpk 761 pe3_out(:,:,jk) = e3vw_0(:,:,jk) + ( 1.0_wp - 0.5_wp * ( vmask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd ) ) & 762 & * ( pe3_in(:,:,jk-1) - e3v_0(:,:,jk-1) ) & 763 & + 0.5_wp * ( vmask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd ) & 764 & * ( pe3_in(:,:,jk ) - e3v_0(:,:,jk ) ) 765 END DO 732 zc3(:,:) = pssh(:,:) * r1_ht_0(:,:) 733 ! 734 DO jk = 1, jpk 735 pe3(:,:,jk) = e3w_0(:,:,jk) * ( 1.0_wp + zc3(:,:) ) 736 END DO 737 ! 738 CASE( 'UW' ) !* from U- to UW-point 739 ! 740 DO_2D_00_00 741 zc3(ji,jj) = 0.5_wp * ( e1e2t(ji ,jj) * pssh(ji ,jj) & 742 & + e1e2t(ji+1,jj) * pssh(ji+1,jj) ) * r1_hu_0(ji,jj) * r1_e1e2u(ji,jj) 743 END_2D 744 CALL lbc_lnk( 'domvvl', zc3(:,:), 'U', 1._wp ) 745 ! 746 DO jk = 1, jpk 747 pe3(:,:,jk) = e3uw_0(:,:,jk) * ( 1.0_wp + zc3(:,:) ) 748 END DO 749 CASE( 'VW' ) !* from U- to UW-point : vertical simple mean 750 ! 751 DO_2D_00_00 752 zc3(ji,jj) = 0.5_wp * ( e1e2t(ji,jj ) * pssh(ji,jj ) & 753 & + e1e2t(ji,jj+1) * pssh(ji,jj+1) ) * r1_hv_0(ji,jj) * r1_e1e2v(ji,jj) 754 END_2D 755 CALL lbc_lnk( 'domvvl', zc3(:,:), 'V', 1._wp ) 756 ! 757 DO jk = 1, jpk 758 pe3(:,:,jk) = e3vw_0(:,:,jk) * ( 1.0_wp + zc3(:,:) ) 759 END DO 760 ! 766 761 END SELECT 767 762 ! … … 878 873 ! Wetting and drying test case 879 874 CALL usr_def_istate( gdept(:,:,:,Kbb), tmask, ts(:,:,:,:,Kbb), uu(:,:,:,Kbb), vv(:,:,:,Kbb), ssh(:,:,Kbb) ) 880 ts (:,:,:,:,Kmm) = ts (:,:,:,:,Kbb) ! set now values from to before ones875 !!an ts (:,:,:,:,Kmm) = ts (:,:,:,:,Kbb) ! set now values from to before ones 881 876 ssh (:,:,Kmm) = ssh(:,:,Kbb) 882 877 uu (:,:,:,Kmm) = uu (:,:,:,Kbb) … … 923 918 ! e3t(:,:,:,Kmm) = e3t(:,:,:,Kbb) 924 919 ssh(:,:,Kmm)=0._wp 920 ssh(:,:,Kbb)=0._wp 925 921 e3t(:,:,:,Kmm)=e3t_0(:,:,:) 926 922 e3t(:,:,:,Kbb)=e3t_0(:,:,:) -
NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/dynatf.F90
r12529 r12614 213 213 IF( ln_dynadv_vec ) THEN ! Asselin filter applied on velocity 214 214 ! Before filtered scale factor at (u/v)-points 215 CALL dom_vvl_interpol( pe3t(:,:,:,Kmm), pe3u(:,:,:,Kmm), 'U' )216 CALL dom_vvl_interpol( pe3t(:,:,:,Kmm), pe3v(:,:,:,Kmm), 'V' )215 CALL dom_vvl_interpol( ssh(:,:,Kmm), pe3u(:,:,:,Kmm), 'U' ) 216 CALL dom_vvl_interpol( ssh(:,:,Kmm), pe3v(:,:,:,Kmm), 'V' ) 217 217 DO_3D_11_11( 1, jpkm1 ) 218 218 puu(ji,jj,jk,Kmm) = puu(ji,jj,jk,Kmm) + rn_atfp * ( puu(ji,jj,jk,Kbb) - 2._wp * puu(ji,jj,jk,Kmm) + puu(ji,jj,jk,Kaa) ) … … 224 224 ALLOCATE( ze3u_f(jpi,jpj,jpk) , ze3v_f(jpi,jpj,jpk) ) 225 225 ! Now filtered scale factor at (u/v)-points stored in ze3u_f, ze3v_f 226 CALL dom_vvl_interpol( pe3t(:,:,:,Kmm), ze3u_f, 'U' )227 CALL dom_vvl_interpol( pe3t(:,:,:,Kmm), ze3v_f, 'V' )226 CALL dom_vvl_interpol( ssh(:,:,Kmm), ze3u_f, 'U' ) 227 CALL dom_vvl_interpol( ssh(:,:,Kmm), ze3v_f, 'V' ) 228 228 DO_3D_11_11( 1, jpkm1 ) 229 229 zue3a = pe3u(ji,jj,jk,Kaa) * puu(ji,jj,jk,Kaa) -
NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/dynvor.F90
r12529 r12614 461 461 zwz(ji,jj) = ff_f(ji,jj) + ( e2v(ji+1,jj ) * pv(ji+1,jj ,jk) - e2v(ji,jj) * pv(ji,jj,jk) & 462 462 & - e1u(ji ,jj+1) * pu(ji ,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk) ) * r1_e1e2f(ji,jj) 463 !!an & * fmask(ji,jj,jk) 463 464 END_2D 464 465 CASE ( np_CME ) !* Coriolis + metric -
NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/nemogcm.F90
r12529 r12614 71 71 USE stopts ! Stochastic param.: ??? 72 72 USE diu_layers ! diurnal bulk SST and coolskin 73 USE step_diu ! diurnal bulk SST timestepping (called from here if run offline)74 73 USE crsini ! initialise grid coarsening utility 75 74 USE dia25h ! 25h mean output … … 140 139 CALL nemo_init !== Initialisations ==! 141 140 ! !-----------------------! 141 142 142 #if defined key_agrif 143 143 Kbb_a = Nbb; Kmm_a = Nnn; Krhs_a = Nrhs ! agrif_oce module copies of time level indices … … 165 165 istp = nit000 166 166 ! 167 #if defined key_c1d 168 DO WHILE ( istp <= nitend .AND. nstop == 0 ) !== C1D time-stepping ==! 169 CALL stp_c1d( istp ) 167 ! !== Standard time-stepping ==! 168 ! 169 DO WHILE( istp <= nitend .AND. nstop == 0 ) 170 171 ncom_stp = istp 172 IF( ln_timing ) THEN 173 zstptiming = MPI_Wtime() 174 IF ( istp == ( nit000 + 1 ) ) elapsed_time = zstptiming 175 IF ( istp == nitend ) elapsed_time = zstptiming - elapsed_time 176 ENDIF 177 178 CALL stp ( istp ) 170 179 istp = istp + 1 180 181 IF( lwp .AND. ln_timing ) WRITE(numtime,*) 'timing step ', istp-1, ' : ', MPI_Wtime() - zstptiming 182 171 183 END DO 172 #else 173 ! 174 # if defined key_agrif 175 ! !== AGRIF time-stepping ==! 176 CALL Agrif_Regrid() 177 ! 178 ! Recursive update from highest nested level to lowest: 179 Kbb_a = Nbb; Kmm_a = Nnn; Krhs_a = Nrhs ! agrif_oce module copies of time level indices 180 CALL Agrif_step_child_adj(Agrif_Update_All) 181 ! 182 DO WHILE( istp <= nitend .AND. nstop == 0 ) 183 CALL stp 184 istp = istp + 1 185 END DO 186 ! 187 IF( .NOT. Agrif_Root() ) THEN 188 CALL Agrif_ParentGrid_To_ChildGrid() 189 IF( ln_diaobs ) CALL dia_obs_wri 190 IF( ln_timing ) CALL timing_finalize 191 CALL Agrif_ChildGrid_To_ParentGrid() 192 ENDIF 193 ! 194 # else 195 ! 196 IF( .NOT.ln_diurnal_only ) THEN !== Standard time-stepping ==! 197 ! 198 DO WHILE( istp <= nitend .AND. nstop == 0 ) 199 200 ncom_stp = istp 201 IF( ln_timing ) THEN 202 zstptiming = MPI_Wtime() 203 IF ( istp == ( nit000 + 1 ) ) elapsed_time = zstptiming 204 IF ( istp == nitend ) elapsed_time = zstptiming - elapsed_time 205 ENDIF 206 207 CALL stp ( istp ) 208 istp = istp + 1 209 210 IF( lwp .AND. ln_timing ) WRITE(numtime,*) 'timing step ', istp-1, ' : ', MPI_Wtime() - zstptiming 211 212 END DO 213 ! 214 ELSE !== diurnal SST time-steeping only ==! 215 ! 216 DO WHILE( istp <= nitend .AND. nstop == 0 ) 217 CALL stp_diurnal( istp ) ! time step only the diurnal SST 218 istp = istp + 1 219 END DO 220 ! 221 ENDIF 222 ! 223 # endif 224 ! 225 #endif 226 ! 227 IF( ln_diaobs ) CALL dia_obs_wri 228 ! 229 IF( ln_icebergs ) CALL icb_end( nitend ) 230 184 ! 185 ! 231 186 ! !------------------------! 232 187 ! !== finalize the run ==! … … 244 199 ! 245 200 #if defined key_iomput 246 CALL xios_finalize ! end mpp communications with xios 247 IF( lk_oasis ) CALL cpl_finalize ! end coupling and mpp communications with OASIS 201 CALL xios_finalize ! end mpp communications with xios 248 202 #else 249 IF ( lk_oasis ) THEN ; CALL cpl_finalize ! end coupling and mpp communications with OASIS 250 ELSEIF( lk_mpp ) THEN ; CALL mppstop ! end mpp communications 251 ENDIF 203 IF( lk_mpp ) CALL mppstop ! end mpp communications 252 204 #endif 253 205 ! … … 419 371 ! 420 372 CALL phy_cst ! Physical constants 421 CALL eos_init ! Equation of state 422 IF( lk_c1d ) CALL c1d_init ! 1D column configuration 423 CALL wad_init ! Wetting and drying options 373 424 374 CALL dom_init( Nbb, Nnn, Naa, "OPA") ! Domain 425 IF( ln_crs ) CALL crs_init( Nnn ) ! coarsened grid: domain initialization 375 426 376 IF( sn_cfctl%l_prtctl ) & 427 377 & CALL prt_ctl_init ! Print control 428 429 CALL diurnal_sst_bulk_init ! diurnal sst 430 IF( ln_diurnal ) CALL diurnal_sst_coolskin_init ! cool skin 431 ! 432 IF( ln_diurnal_only ) THEN ! diurnal only: a subset of the initialisation routines 433 CALL istate_init( Nbb, Nnn, Naa ) ! ocean initial state (Dynamics and tracers) 434 CALL sbc_init( Nbb, Nnn, Naa ) ! Forcings : surface module 435 CALL tra_qsr_init ! penetrative solar radiation qsr 436 IF( ln_diaobs ) THEN ! Observation & model comparison 437 CALL dia_obs_init( Nnn ) ! Initialize observational data 438 CALL dia_obs( nit000 - 1, Nnn ) ! Observation operator for restart 439 ENDIF 440 IF( lk_asminc ) CALL asm_inc_init( Nbb, Nnn, Nrhs ) ! Assimilation increments 441 ! 442 RETURN ! end of initialization 443 ENDIF 444 ! 445 378 446 379 CALL istate_init( Nbb, Nnn, Naa ) ! ocean initial state (Dynamics and tracers) 447 380 448 381 ! ! external forcing 449 382 CALL tide_init ! tidal harmonics 383 450 384 CALL sbc_init( Nbb, Nnn, Naa ) ! surface boundary conditions (including sea-ice) 451 CALL bdy_init ! Open boundaries initialisation 452 453 ! ! Ocean physics 454 CALL zdf_phy_init( Nnn ) ! Vertical physics 455 385 386 387 ! ! Ocean physics 456 388 ! ! Lateral physics 457 CALL ldf_tra_init ! Lateral ocean tracer physics458 CALL ldf_eiv_init ! eddy induced velocity param.459 389 CALL ldf_dyn_init ! Lateral ocean momentum physics 460 461 ! ! Active tracers 462 IF( ln_traqsr ) CALL tra_qsr_init ! penetrative solar radiation qsr 463 CALL tra_bbc_init ! bottom heat flux 464 CALL tra_bbl_init ! advective (and/or diffusive) bottom boundary layer scheme 465 CALL tra_dmp_init ! internal tracer damping 466 CALL tra_adv_init ! horizontal & vertical advection 467 CALL tra_ldf_init ! lateral mixing 468 390 391 469 392 ! ! Dynamics 470 IF( lk_c1d ) CALL dyn_dmp_init ! internal momentum damping471 393 CALL dyn_adv_init ! advection (vector or flux form) 394 472 395 CALL dyn_vor_init ! vorticity term including Coriolis 396 473 397 CALL dyn_ldf_init ! lateral mixing 474 CALL dyn_hpg_init( Nnn ) ! horizontal gradient of Hydrostatic pressure 398 475 399 CALL dyn_spg_init ! surface pressure gradient 476 400 477 #if defined key_top478 ! ! Passive tracers479 CALL trc_init( Nbb, Nnn, Naa )480 #endif481 IF( l_ldfslp ) CALL ldf_slp_init ! slope of lateral mixing482 483 ! ! Icebergs484 CALL icb_init( rn_Dt, nit000) ! initialise icebergs instance485 486 ! ice shelf487 CALL isf_init( Nbb, Nnn, Naa )488 489 ! ! Misc. options490 CALL sto_par_init ! Stochastic parametrization491 IF( ln_sto_eos ) CALL sto_pts_init ! RRandom T/S fluctuations492 493 401 ! ! Diagnostics 494 402 CALL flo_init( Nnn ) ! drifting Floats 403 495 404 IF( ln_diacfl ) CALL dia_cfl_init ! Initialise CFL diagnostics 496 ! CALL dia_ptr_init ! Poleward TRansports initialization 497 CALL dia_dct_init ! Sections tranports 498 CALL dia_hsb_init( Nnn ) ! heat content, salt content and volume budgets 405 499 406 CALL trd_init( Nnn ) ! Mixed-layer/Vorticity/Integral constraints trends 500 CALL dia_obs_init( Nnn ) ! Initialize observational data 501 CALL dia_25h_init( Nbb ) ! 25h mean outputs 502 CALL dia_detide_init ! Weights computation for daily detiding of model diagnostics 503 IF( ln_diaobs ) CALL dia_obs( nit000-1, Nnn ) ! Observation operator for restart 504 CALL dia_mlr_init ! Initialisation of IOM context management for multiple-linear-regression analysis 505 506 ! ! Assimilation increments 507 IF( lk_asminc ) CALL asm_inc_init( Nbb, Nnn, Nrhs ) ! Initialize assimilation increments 508 ! 407 408 509 409 IF(lwp) WRITE(numout,cform_aaa) ! Flag AAAAAAA 510 410 ! 511 411 IF( ln_timing ) CALL timing_stop( 'nemo_init') 512 412 ! 413 513 414 END SUBROUTINE nemo_init 514 415 -
NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/phycst.F90
r12529 r12614 16 16 !!---------------------------------------------------------------------- 17 17 USE par_oce ! ocean parameters 18 USE in_out_manager ! I/O manager 18 USE in_out_manager ! I/O manager 19 19 20 20 IMPLICIT NONE … … 66 66 REAL(wp), PUBLIC :: r1_rhos !: 1 / rhos 67 67 REAL(wp), PUBLIC :: r1_rcpi !: 1 / rcpi 68 68 69 !!---------------------------------------------------------------------- 69 70 !! NEMO/OCE 4.0 , NEMO Consortium (2018) … … 92 93 r1_rhos = 1._wp / rhos 93 94 r1_rcpi = 1._wp / rcpi 94 95 ! 96 rho0 = 1026._wp !: volumic mass of reference [kg/m3] 97 rcp = 3991.86795711963_wp !: heat capacity [J/K] 98 ! 99 rho0_rcp = rho0 * rcp 100 r1_rho0 = 1._wp / rho0 101 r1_rcp = 1._wp / rcp 102 r1_rho0_rcp = 1._wp / rho0_rcp 103 ! 95 104 IF(lwp) THEN 96 105 WRITE(numout,*) -
NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/sbcice_cice.F90
r12529 r12614 244 244 ! Horizontal scale factor interpolations 245 245 ! -------------------------------------- 246 CALL dom_vvl_interpol( e3t(:,:,:,Kbb), e3u(:,:,:,Kbb), 'U' )247 CALL dom_vvl_interpol( e3t(:,:,:,Kbb), e3v(:,:,:,Kbb), 'V' )248 CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3u(:,:,:,Kmm), 'U' )249 CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3v(:,:,:,Kmm), 'V' )250 CALL dom_vvl_interpol( e3u(:,:,:,Kmm), e3f(:,:,:), 'F' )246 CALL dom_vvl_interpol( ssh(:,:,Kbb), e3u(:,:,:,Kbb), 'U' ) 247 CALL dom_vvl_interpol( ssh(:,:,Kbb), e3v(:,:,:,Kbb), 'V' ) 248 CALL dom_vvl_interpol( ssh(:,:,Kmm), e3u(:,:,:,Kmm), 'U' ) 249 CALL dom_vvl_interpol( ssh(:,:,Kmm), e3v(:,:,:,Kmm), 'V' ) 250 CALL dom_vvl_interpol( ssh(:,:,Kmm), e3f(:,:,:), 'F' ) 251 251 ! Vertical scale factor interpolations 252 252 ! ------------------------------------ 253 CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3w (:,:,:,Kmm), 'W' )254 CALL dom_vvl_interpol( e3u(:,:,:,Kmm), e3uw(:,:,:,Kmm), 'UW' )255 CALL dom_vvl_interpol( e3v(:,:,:,Kmm), e3vw(:,:,:,Kmm), 'VW' )256 CALL dom_vvl_interpol( e3u(:,:,:,Kbb), e3uw(:,:,:,Kbb), 'UW' )257 CALL dom_vvl_interpol( e3v(:,:,:,Kbb), e3vw(:,:,:,Kbb), 'VW' )253 CALL dom_vvl_interpol( ssh(:,:,Kmm), e3w (:,:,:,Kmm), 'W' ) 254 CALL dom_vvl_interpol( ssh(:,:,Kmm), e3uw(:,:,:,Kmm), 'UW' ) 255 CALL dom_vvl_interpol( ssh(:,:,Kmm), e3vw(:,:,:,Kmm), 'VW' ) 256 CALL dom_vvl_interpol( ssh(:,:,Kbb), e3uw(:,:,:,Kbb), 'UW' ) 257 CALL dom_vvl_interpol( ssh(:,:,Kbb), e3vw(:,:,:,Kbb), 'VW' ) 258 258 ! t- and w- points depth 259 259 ! ---------------------- -
NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/step.F90
r12529 r12614 2 2 !!====================================================================== 3 3 !! *** MODULE step *** 4 !! Time-stepping : manager of the ocean, tracer and icetime stepping4 !! Time-stepping : manager of the shallow water equation time stepping 5 5 !!====================================================================== 6 !! History : OPA ! 1991-03 (G. Madec) Original code 7 !! - ! 1991-11 (G. Madec) 8 !! - ! 1992-06 (M. Imbard) add a first output record 9 !! - ! 1996-04 (G. Madec) introduction of dynspg 10 !! - ! 1996-04 (M.A. Foujols) introduction of passive tracer 11 !! 8.0 ! 1997-06 (G. Madec) new architecture of call 12 !! 8.2 ! 1997-06 (G. Madec, M. Imbard, G. Roullet) free surface 13 !! - ! 1999-02 (G. Madec, N. Grima) hpg implicit 14 !! - ! 2000-07 (J-M Molines, M. Imbard) Open Bondary Conditions 15 !! NEMO 1.0 ! 2002-06 (G. Madec) free form, suppress macro-tasking 16 !! - ! 2004-08 (C. Talandier) New trends organization 17 !! - ! 2005-01 (C. Ethe) Add the KPP closure scheme 18 !! - ! 2005-11 (G. Madec) Reorganisation of tra and dyn calls 19 !! - ! 2006-01 (L. Debreu, C. Mazauric) Agrif implementation 20 !! - ! 2006-07 (S. Masson) restart using iom 21 !! 3.2 ! 2009-02 (G. Madec, R. Benshila) reintroduicing z*-coordinate 22 !! - ! 2009-06 (S. Masson, G. Madec) TKE restart compatible with key_cpl 23 !! 3.3 ! 2010-05 (K. Mogensen, A. Weaver, M. Martin, D. Lea) Assimilation interface 24 !! - ! 2010-10 (C. Ethe, G. Madec) reorganisation of initialisation phase + merge TRC-TRA 25 !! 3.4 ! 2011-04 (G. Madec, C. Ethe) Merge of dtatem and dtasal 26 !! 3.6 ! 2012-07 (J. Simeon, G. Madec. C. Ethe) Online coarsening of outputs 27 !! 3.6 ! 2014-04 (F. Roquet, G. Madec) New equations of state 28 !! 3.6 ! 2014-10 (E. Clementi, P. Oddo) Add Qiao vertical mixing in case of waves 29 !! 3.7 ! 2014-10 (G. Madec) LDF simplication 30 !! - ! 2014-12 (G. Madec) remove KPP scheme 31 !! - ! 2015-11 (J. Chanut) free surface simplification (remove filtered free surface) 32 !! 4.0 ! 2017-05 (G. Madec) introduction of the vertical physics manager (zdfphy) 33 !! 4.1 ! 2019-08 (A. Coward, D. Storkey) rewrite in preparation for new timestepping scheme 34 !!---------------------------------------------------------------------- 35 36 !!---------------------------------------------------------------------- 37 !! stp : OPA system time-stepping 6 !! History : NEMO ! 2020-03 (A. Nasser, G. Madec) Original code from 4.0.2 7 !!---------------------------------------------------------------------- 8 9 !!---------------------------------------------------------------------- 10 !! stp : Shallow Water time-stepping 38 11 !!---------------------------------------------------------------------- 39 12 USE step_oce ! time stepping definition modules 13 USE phycst ! physical constants 14 USE usrdef_nam 40 15 ! 41 USE iom ! xIOs server 16 USE iom ! xIOs server 42 17 43 18 IMPLICIT NONE … … 45 20 46 21 PUBLIC stp ! called by nemogcm.F90 47 22 48 23 !!---------------------------------------------------------------------- 49 24 !! time level indices 50 25 !!---------------------------------------------------------------------- 51 INTEGER, PUBLIC :: Nbb, Nnn, Naa, Nrhs !! used by nemo_init 52 26 INTEGER, PUBLIC :: Nbb, Nnn, Naa, Nrhs !! used by nemo_init 27 28 !! * Substitutions 29 # include "do_loop_substitute.h90" 53 30 !!---------------------------------------------------------------------- 54 31 !! NEMO/OCE 4.0 , NEMO Consortium (2018) … … 68 45 !! *** ROUTINE stp *** 69 46 !! 70 !! ** Purpose : - Time stepping of OPA (momentum and active tracer eqs.) 71 !! - Time stepping of SI3 (dynamic and thermodynamic eqs.) 72 !! - Time stepping of TRC (passive tracer eqs.) 47 !! ** Purpose : - Time stepping of shallow water (SHW) (momentum and ssh eqs.) 73 48 !! 74 !! ** Method : -1- Update forcings and data 75 !! -2- Update ocean physics 76 !! -3- Compute the t and s trends 77 !! -4- Update t and s 78 !! -5- Compute the momentum trends 79 !! -6- Update the horizontal velocity 80 !! -7- Compute the diagnostics variables (rd,N2, hdiv,w) 81 !! -8- Outputs and diagnostics 49 !! ** Method : -1- Update forcings 50 !! -2- Update the ssh at Naa 51 !! -3- Compute the momentum trends (Nrhs) 52 !! -4- Update the horizontal velocity 53 !! -5- Apply Asselin time filter to uu,vv,ssh 54 !! -6- Outputs and diagnostics 82 55 !!---------------------------------------------------------------------- 83 56 INTEGER :: ji, jj, jk ! dummy loop indice … … 85 58 !!gm kcall can be removed, I guess 86 59 INTEGER :: kcall ! optional integer argument (dom_vvl_sf_nxt) 60 REAL(wp):: z1_2rho0 ! local scalars 61 62 REAL(wp) :: zue3a, zue3n, zue3b ! local scalars 63 REAL(wp) :: zve3a, zve3n, zve3b ! - - 64 REAL(wp) :: ze3t_tf, ze3u_tf, ze3v_tf, zua, zva 87 65 !! --------------------------------------------------------------------- 88 66 #if defined key_agrif … … 105 83 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 106 84 ! 107 IF( l_1st_euler ) THEN 108 ! start or restart with Euler 1st time-step 109 rDt = rn_Dt 85 IF( l_1st_euler ) THEN ! start or restart with Euler 1st time-step 86 rDt = rn_Dt 110 87 r1_Dt = 1._wp / rDt 111 88 ENDIF 89 90 IF ( kstp == nit000 ) ww(:,:,:) = 0._wp ! initialize vertical velocity one for all to zero 91 112 92 ! 113 93 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> … … 132 112 IF( ln_apr_dyn ) CALL sbc_apr ( kstp ) ! atmospheric pressure (NB: call before bdy_dta which needs ssh_ib) 133 113 IF( ln_bdy ) CALL bdy_dta ( kstp, Nnn ) ! update dynamic & tracer data at open boundaries 134 IF( ln_isf ) CALL isf_stp ( kstp, Nnn ) 135 CALL sbc ( kstp, Nbb, Nnn ) ! Sea Boundary Condition (including sea-ice) 136 137 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 138 ! Update stochastic parameters and random T/S fluctuations 139 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 140 IF( ln_sto_eos ) CALL sto_par( kstp ) ! Stochastic parameters 141 IF( ln_sto_eos ) CALL sto_pts( ts(:,:,:,:,Nnn) ) ! Random T/S fluctuations 114 CALL sbc ( kstp, Nbb, Nnn ) ! Sea Boundary Condition 142 115 143 116 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 144 117 ! Ocean physics update 145 118 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 146 ! THERMODYNAMICS147 CALL eos_rab( ts(:,:,:,:,Nbb), rab_b, Nnn ) ! before local thermal/haline expension ratio at T-points148 CALL eos_rab( ts(:,:,:,:,Nnn), rab_n, Nnn ) ! now local thermal/haline expension ratio at T-points149 CALL bn2 ( ts(:,:,:,:,Nbb), rab_b, rn2b, Nnn ) ! before Brunt-Vaisala frequency150 CALL bn2 ( ts(:,:,:,:,Nnn), rab_n, rn2, Nnn ) ! now Brunt-Vaisala frequency151 152 ! VERTICAL PHYSICS153 CALL zdf_phy( kstp, Nbb, Nnn, Nrhs ) ! vertical physics update (top/bot drag, avt, avs, avm + MLD)154 155 119 ! LATERAL PHYSICS 156 !157 IF( l_ldfslp ) THEN ! slope of lateral mixing158 CALL eos( ts(:,:,:,:,Nbb), rhd, gdept_0(:,:,:) ) ! before in situ density159 160 IF( ln_zps .AND. .NOT. ln_isfcav) &161 & CALL zps_hde ( kstp, Nnn, jpts, ts(:,:,:,:,Nbb), gtsu, gtsv, & ! Partial steps: before horizontal gradient162 & rhd, gru , grv ) ! of t, s, rd at the last ocean level163 164 IF( ln_zps .AND. ln_isfcav) &165 & CALL zps_hde_isf( kstp, Nnn, jpts, ts(:,:,:,:,Nbb), gtsu, gtsv, gtui, gtvi, & ! Partial steps for top cell (ISF)166 & rhd, gru , grv , grui, grvi ) ! of t, s, rd at the first ocean level167 IF( ln_traldf_triad ) THEN168 CALL ldf_slp_triad( kstp, Nbb, Nnn ) ! before slope for triad operator169 ELSE170 CALL ldf_slp ( kstp, rhd, rn2b, Nbb, Nnn ) ! before slope for standard operator171 ENDIF172 ENDIF173 120 ! ! eddy diffusivity coeff. 174 IF( l_ldftra_time .OR. l_ldfeiv_time ) CALL ldf_tra( kstp, Nbb, Nnn ) ! and/or eiv coeff.175 121 IF( l_ldfdyn_time ) CALL ldf_dyn( kstp, Nbb ) ! eddy viscosity coeff. 176 122 … … 180 126 181 127 CALL ssh_nxt ( kstp, Nbb, Nnn, ssh, Naa ) ! after ssh (includes call to div_hor) 128 182 129 IF( .NOT.ln_linssh ) CALL dom_vvl_sf_nxt( kstp, Nbb, Nnn, Naa ) ! after vertical scale factors 183 CALL wzv ( kstp, Nbb, Nnn, ww, Naa ) ! now cross-level velocity 184 IF( ln_zad_Aimp ) CALL wAimp ( kstp, Nnn ) ! Adaptive-implicit vertical advection partitioning 185 CALL eos ( ts(:,:,:,:,Nnn), rhd, rhop, gdept(:,:,:,Nnn) ) ! now in situ density for hpg computation 186 187 130 188 131 uu(:,:,:,Nrhs) = 0._wp ! set dynamics trends to zero 189 132 vv(:,:,:,Nrhs) = 0._wp 190 133 191 IF( lk_asminc .AND. ln_asmiau .AND. ln_dyninc ) &192 & CALL dyn_asm_inc ( kstp, Nbb, Nnn, uu, vv, Nrhs ) ! apply dynamics assimilation increment193 134 IF( ln_bdy ) CALL bdy_dyn3d_dmp ( kstp, Nbb, uu, vv, Nrhs ) ! bdy damping trends 135 194 136 #if defined key_agrif 195 137 IF(.NOT. Agrif_Root()) & … … 197 139 #endif 198 140 CALL dyn_adv( kstp, Nbb, Nnn , uu, vv, Nrhs ) ! advection (VF or FF) ==> RHS 141 199 142 CALL dyn_vor( kstp, Nnn , uu, vv, Nrhs ) ! vorticity ==> RHS 143 200 144 CALL dyn_ldf( kstp, Nbb, Nnn , uu, vv, Nrhs ) ! lateral mixing 201 IF( ln_zdfosm ) CALL dyn_osm( kstp, Nnn , uu, vv, Nrhs ) ! OSMOSIS non-local velocity fluxes ==> RHS 202 CALL dyn_hpg( kstp, Nnn , uu, vv, Nrhs ) ! horizontal gradient of Hydrostatic pressure 203 CALL dyn_spg( kstp, Nbb, Nnn, Nrhs, uu, vv, ssh, uu_b, vv_b, Naa ) ! surface pressure gradient 204 205 ! With split-explicit free surface, since now transports have been updated and ssh(:,:,Nrhs) as well 206 IF( ln_dynspg_ts ) THEN ! vertical scale factors and vertical velocity need to be updated 207 CALL div_hor ( kstp, Nbb, Nnn ) ! Horizontal divergence (2nd call in time-split case) 208 IF(.NOT.ln_linssh) CALL dom_vvl_sf_nxt( kstp, Nbb, Nnn, Naa, kcall=2 ) ! after vertical scale factors (update depth average component) 209 ENDIF 210 CALL dyn_zdf ( kstp, Nbb, Nnn, Nrhs, uu, vv, Naa ) ! vertical diffusion 211 IF( ln_dynspg_ts ) THEN ! vertical scale factors and vertical velocity need to be updated 212 CALL wzv ( kstp, Nbb, Nnn, ww, Naa ) ! now cross-level velocity 213 IF( ln_zad_Aimp ) CALL wAimp ( kstp, Nnn ) ! Adaptive-implicit vertical advection partitioning 214 ENDIF 215 216 217 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 218 ! cool skin 219 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 220 IF ( ln_diurnal ) CALL diurnal_layers( kstp ) 221 222 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 223 ! diagnostics and outputs 224 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 225 IF( ln_floats ) CALL flo_stp ( kstp, Nbb, Nnn ) ! drifting Floats 226 IF( ln_diacfl ) CALL dia_cfl ( kstp, Nnn ) ! Courant number diagnostics 227 CALL dia_hth ( kstp, Nnn ) ! Thermocline depth (20 degres isotherm depth) 228 IF( ln_diadct ) CALL dia_dct ( kstp, Nnn ) ! Transports 229 CALL dia_ar5 ( kstp, Nnn ) ! ar5 diag 230 CALL dia_ptr ( kstp, Nnn ) ! Poleward adv/ldf TRansports diagnostics 231 CALL dia_wri ( kstp, Nnn ) ! ocean model: outputs 232 IF( ln_crs ) CALL crs_fld ( kstp, Nnn ) ! ocean model: online field coarsening & output 233 IF( lk_diadetide ) CALL dia_detide( kstp ) ! Weights computation for daily detiding of model diagnostics 234 IF( lk_diamlr ) CALL dia_mlr ! Update time used in multiple-linear-regression analysis 235 236 #if defined key_top 237 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 238 ! Passive Tracer Model 239 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 240 CALL trc_stp ( kstp, Nbb, Nnn, Nrhs, Naa ) ! time-stepping 241 #endif 242 243 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 244 ! Active tracers 245 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 246 ts(:,:,:,:,Nrhs) = 0._wp ! set tracer trends to zero 247 248 IF( lk_asminc .AND. ln_asmiau .AND. & 249 & ln_trainc ) CALL tra_asm_inc( kstp, Nbb, Nnn, ts, Nrhs ) ! apply tracer assimilation increment 250 CALL tra_sbc ( kstp, Nnn, ts, Nrhs ) ! surface boundary condition 251 IF( ln_traqsr ) CALL tra_qsr ( kstp, Nnn, ts, Nrhs ) ! penetrative solar radiation qsr 252 IF( ln_isf ) CALL tra_isf ( kstp, Nnn, ts, Nrhs ) ! ice shelf heat flux 253 IF( ln_trabbc ) CALL tra_bbc ( kstp, Nnn, ts, Nrhs ) ! bottom heat flux 254 IF( ln_trabbl ) CALL tra_bbl ( kstp, Nbb, Nnn, ts, Nrhs ) ! advective (and/or diffusive) bottom boundary layer scheme 255 IF( ln_tradmp ) CALL tra_dmp ( kstp, Nbb, Nnn, ts, Nrhs ) ! internal damping trends 256 IF( ln_bdy ) CALL bdy_tra_dmp( kstp, Nbb, ts, Nrhs ) ! bdy damping trends 257 #if defined key_agrif 258 IF(.NOT. Agrif_Root()) & 259 & CALL Agrif_Sponge_tra ! tracers sponge 260 #endif 261 CALL tra_adv ( kstp, Nbb, Nnn, ts, Nrhs ) ! hor. + vert. advection ==> RHS 262 IF( ln_zdfosm ) CALL tra_osm ( kstp, Nnn, ts, Nrhs ) ! OSMOSIS non-local tracer fluxes ==> RHS 263 IF( lrst_oce .AND. ln_zdfosm ) & 264 & CALL osm_rst ( kstp, Nnn, 'WRITE' ) ! write OSMOSIS outputs + ww (so must do here) to restarts 265 CALL tra_ldf ( kstp, Nbb, Nnn, ts, Nrhs ) ! lateral mixing 266 267 CALL tra_zdf ( kstp, Nbb, Nnn, Nrhs, ts, Naa ) ! vertical mixing and after tracer fields 268 IF( ln_zdfnpc ) CALL tra_npc ( kstp, Nnn, Nrhs, ts, Naa ) ! update after fields by non-penetrative convection 145 146 !!an - calcul du gradient de pression horizontal (explicit) 147 DO_3D_00_00( 1, jpkm1 ) 148 uu(ji,jj,jk,Nrhs) = uu(ji,jj,jk,Nrhs) - grav * ( ssh(ji+1,jj,Nnn) - ssh(ji,jj,Nnn) ) * r1_e1u(ji,jj) 149 vv(ji,jj,jk,Nrhs) = vv(ji,jj,jk,Nrhs) - grav * ( ssh(ji,jj+1,Nnn) - ssh(ji,jj,Nnn) ) * r1_e2v(ji,jj) 150 END_3D 151 ! 152 ! add wind stress forcing and layer linear friction to the RHS 153 z1_2rho0 = 0.5_wp * r1_rho0 154 DO_3D_00_00(1,jpkm1) 155 uu(ji,jj,jk,Nrhs) = uu(ji,jj,jk,Nrhs) + z1_2rho0 * ( utau_b(ji,jj) + utau(ji,jj) ) / e3u(ji,jj,jk,Nnn) & 156 & - rn_rfr * uu(ji,jj,jk,Nbb) 157 vv(ji,jj,jk,Nrhs) = vv(ji,jj,jk,Nrhs) + z1_2rho0 * ( vtau_b(ji,jj) + vtau(ji,jj) ) / e3v(ji,jj,jk,Nnn) & 158 & - rn_rfr * vv(ji,jj,jk,Nbb) 159 END_3D 160 !!an 161 162 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 163 ! Leap-Frog time splitting + Robert-Asselin time filter on u,v,e3 164 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 165 166 !!an futur module dyn_nxt (a la place de dyn_atf) 167 168 IF( ln_dynadv_vec ) THEN ! vector invariant form : applied on velocity 169 IF( l_1st_euler ) THEN ! Euler time stepping (no Asselin filter) 170 DO_3D_00_00(1,jpkm1) 171 uu(ji,jj,jk,Naa) = uu(ji,jj,jk,Nbb) + rDt * uu(ji,jj,jk,Nrhs) * umask(ji,jj,jk) 172 vv(ji,jj,jk,Naa) = vv(ji,jj,jk,Nbb) + rDt * vv(ji,jj,jk,Nrhs) * vmask(ji,jj,jk) 173 END_3D 174 ELSE ! Leap Frog time stepping + Asselin filter 175 DO_3D_11_11(1,jpkm1) 176 zua = uu(ji,jj,jk,Nbb) + rDt * uu(ji,jj,jk,Nrhs) * umask(ji,jj,jk) 177 zva = vv(ji,jj,jk,Nbb) + rDt * vv(ji,jj,jk,Nrhs) * vmask(ji,jj,jk) 178 ! ! Asselin time filter on u,v (Nnn) 179 uu(ji,jj,jk,Nnn) = uu(ji,jj,jk,Nnn) + rn_atfp * (uu(ji,jj,jk,Nbb) - 2._wp * uu(ji,jj,jk,Nnn) + zua) 180 vv(ji,jj,jk,Nnn) = vv(ji,jj,jk,Nnn) + rn_atfp * (vv(ji,jj,jk,Nbb) - 2._wp * vv(ji,jj,jk,Nnn) + zva) 181 ! 182 ze3u_tf = e3u(ji,jj,jk,Nnn) + rn_atfp * ( e3u(ji,jj,jk,Nbb) - 2._wp * e3u(ji,jj,jk,Nnn) + e3u(ji,jj,jk,Naa) ) 183 ze3v_tf = e3v(ji,jj,jk,Nnn) + rn_atfp * ( e3v(ji,jj,jk,Nbb) - 2._wp * e3v(ji,jj,jk,Nnn) + e3v(ji,jj,jk,Naa) ) 184 ze3t_tf = e3t(ji,jj,jk,Nnn) + rn_atfp * ( e3t(ji,jj,jk,Nbb) - 2._wp * e3t(ji,jj,jk,Nnn) + e3t(ji,jj,jk,Naa) ) 185 ! 186 e3u(ji,jj,jk,Nnn) = ze3u_tf 187 e3v(ji,jj,jk,Nnn) = ze3v_tf 188 e3t(ji,jj,jk,Nnn) = ze3t_tf 189 ! 190 uu(ji,jj,jk,Naa) = zua 191 vv(ji,jj,jk,Naa) = zva 192 END_3D 193 ENDIF 194 ! 195 ELSE ! flux form : applied on thickness weighted velocity 196 IF( l_1st_euler ) THEN ! Euler time stepping (no Asselin filter) 197 DO_3D_00_00(1,jpkm1) 198 zue3b = e3u(ji,jj,jk,Nbb) * uu(ji,jj,jk,Nbb) 199 zve3b = e3v(ji,jj,jk,Nbb) * vv(ji,jj,jk,Nbb) 200 ! ! LF time stepping 201 zue3a = zue3b + rDt * e3u(ji,jj,jk,Nrhs) * uu(ji,jj,jk,Nrhs) * umask(ji,jj,jk) 202 zve3a = zve3b + rDt * e3v(ji,jj,jk,Nrhs) * vv(ji,jj,jk,Nrhs) * vmask(ji,jj,jk) 203 ! 204 uu(ji,jj,jk,Naa) = zue3a / e3u(ji,jj,jk,Naa) 205 vv(ji,jj,jk,Naa) = zve3a / e3v(ji,jj,jk,Naa) 206 END_3D 207 ELSE ! Leap Frog time stepping + Asselin filter 208 DO_3D_11_11(1,jpkm1) 209 zue3n = e3u(ji,jj,jk,Nnn) * uu(ji,jj,jk,Nnn) 210 zve3n = e3v(ji,jj,jk,Nnn) * vv(ji,jj,jk,Nnn) 211 zue3b = e3u(ji,jj,jk,Nbb) * uu(ji,jj,jk,Nbb) 212 zve3b = e3v(ji,jj,jk,Nbb) * vv(ji,jj,jk,Nbb) 213 ! ! LF time stepping 214 zue3a = zue3b + rDt * e3u(ji,jj,jk,Nrhs) * uu(ji,jj,jk,Nrhs) * umask(ji,jj,jk) 215 zve3a = zve3b + rDt * e3v(ji,jj,jk,Nrhs) * vv(ji,jj,jk,Nrhs) * vmask(ji,jj,jk) 216 ! ! Asselin time filter on e3u/v/t 217 ze3u_tf = e3u(ji,jj,jk,Nnn) + rn_atfp * ( e3u(ji,jj,jk,Nbb) - 2._wp * e3u(ji,jj,jk,Nnn) + e3u(ji,jj,jk,Naa) ) 218 ze3v_tf = e3v(ji,jj,jk,Nnn) + rn_atfp * ( e3v(ji,jj,jk,Nbb) - 2._wp * e3v(ji,jj,jk,Nnn) + e3v(ji,jj,jk,Naa) ) 219 ze3t_tf = e3t(ji,jj,jk,Nnn) + rn_atfp * ( e3t(ji,jj,jk,Nbb) - 2._wp * e3t(ji,jj,jk,Nnn) + e3t(ji,jj,jk,Naa) ) 220 ! ! Asselin time filter on u,v (Nnn) 221 uu(ji,jj,jk,Nnn) = ( zue3n + rn_atfp * ( zue3b - 2._wp * zue3n + zue3a ) ) / ze3u_tf 222 vv(ji,jj,jk,Nnn) = ( zve3n + rn_atfp * ( zve3b - 2._wp * zve3n + zve3a ) ) / ze3v_tf 223 ! 224 e3u(ji,jj,jk,Nnn) = ze3u_tf 225 e3v(ji,jj,jk,Nnn) = ze3v_tf 226 e3t(ji,jj,jk,Nnn) = ze3t_tf 227 ! 228 uu(ji,jj,jk,Naa) = zue3a / e3u(ji,jj,jk,Naa) 229 vv(ji,jj,jk,Naa) = zve3a / e3v(ji,jj,jk,Naa) 230 END_3D 231 ENDIF 232 ENDIF 233 234 CALL lbc_lnk_multi( 'stp', uu(:,:,:,Nnn), 'U', -1., vv(:,:,:,Nnn), 'V', -1., & !* local domain boundaries 235 & uu(:,:,:,Naa), 'U', -1., vv(:,:,:,Naa), 'V', -1. ) 236 237 !!an 269 238 270 239 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 271 240 ! Set boundary conditions, time filter and swap time levels 272 241 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 273 !!jc1: For agrif, it would be much better to finalize tracers/momentum here (e.g. bdy conditions) and move the swap 274 !! (and time filtering) after Agrif update. Then restart would be done after and would contain updated fields. 275 !! If so: 276 !! (i) no need to call agrif update at initialization time 277 !! (ii) no need to update "before" fields 278 !! 279 !! Apart from creating new tra_swp/dyn_swp routines, this however: 280 !! (i) makes boundary conditions at initialization time computed from updated fields which is not the case between 281 !! two restarts => restartability issue. One can circumvent this, maybe, by assuming "interface separation", 282 !! e.g. a shift of the feedback interface inside child domain. 283 !! (ii) requires that all restart outputs of updated variables by agrif (e.g. passive tracers/tke/barotropic arrays) are done at the same 284 !! place. 285 !! 286 !!jc2: dynnxt must be the latest call. e3t(:,:,:,Nbb) are indeed updated in that routine 287 CALL tra_atf ( kstp, Nbb, Nnn, Naa, ts ) ! time filtering of "now" tracer arrays 288 CALL dyn_atf ( kstp, Nbb, Nnn, Naa, uu, vv, e3t, e3u, e3v ) ! time filtering of "now" velocities and scale factors 289 CALL ssh_atf ( kstp, Nbb, Nnn, Naa, ssh ) ! time filtering of "now" sea surface height 290 ! 242 243 !!an TO BE ADDED : dyn_nxt 244 !! CALL dyn_atf ( kstp, Nbb, Nnn, Naa, uu, vv, e3t, e3u, e3v ) ! time filtering of "now" velocities and scale factors 245 !!an TO BE ADDED : a simplifier 246 ! CALL ssh_atf ( kstp, Nbb, Nnn, Naa, ssh ) ! time filtering of "now" sea surface height 247 248 IF ( .NOT.( l_1st_euler ) ) THEN ! Only do time filtering for leapfrog timesteps 249 ! ! filtering "now" field 250 ssh(:,:,Nnn) = ssh(:,:,Nnn) + rn_atfp * ( ssh(:,:,Nbb) - 2 * ssh(:,:,Nnn) + ssh(:,:,Naa) ) 251 ENDIF 252 253 !!an 254 255 291 256 ! Swap time levels 292 257 Nrhs = Nbb … … 295 260 Naa = Nrhs 296 261 ! 297 IF(.NOT.ln_linssh) CALL dom_vvl_sf_update( kstp, Nbb, Nnn, Naa ) ! recompute vertical scale factors 298 ! 299 IF( ln_diahsb ) CALL dia_hsb ( kstp, Nbb, Nnn ) ! - ML - global conservation diagnostics 300 301 !!gm : This does not only concern the dynamics ==>>> add a new title 302 !!gm2: why ouput restart before AGRIF update? 303 !! 304 !!jc: That would be better, but see comment above 305 !! 262 CALL dom_vvl_sf_update( kstp, Nbb, Nnn, Naa ) ! recompute vertical scale factors 263 264 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 265 ! diagnostics and outputs 266 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 267 IF( ln_floats ) CALL flo_stp ( kstp, Nbb, Nnn ) ! drifting Floats 268 IF( ln_diacfl ) CALL dia_cfl ( kstp, Nnn ) ! Courant number diagnostics 269 270 CALL dia_wri ( kstp, Nnn ) ! ocean model: outputs 271 272 ! 306 273 IF( lrst_oce ) CALL rst_write ( kstp, Nbb, Nnn ) ! write output ocean restart file 307 IF( ln_sto_eos ) CALL sto_rst_write( kstp ) ! write restart file for stochastic parameters274 308 275 309 276 #if defined key_agrif … … 324 291 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 325 292 CALL stp_ctl ( kstp, Nbb, Nnn, indic ) 293 326 294 327 295 IF( kstp == nit000 ) THEN ! 1st time step only … … 331 299 ENDIF 332 300 333 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>334 ! Coupled mode335 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<336 !!gm why lk_oasis and not lk_cpl ????337 IF( lk_oasis ) CALL sbc_cpl_snd( kstp, Nbb, Nnn ) ! coupled mode : field exchanges338 301 ! 339 302 #if defined key_iomput … … 341 304 CALL iom_context_finalize( cxios_context ) ! needed for XIOS+AGRIF 342 305 IF(lrxios) CALL iom_context_finalize( crxios_context ) 343 IF( ln_crs ) CALL iom_context_finalize( trim(cxios_context)//"_crs" ) !344 306 ENDIF 345 307 #endif -
NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/stpctl.F90
r12529 r12614 67 67 INTEGER, DIMENSION(3) :: iu, is1, is2 ! min/max loc indices 68 68 REAL(wp) :: zzz ! local real 69 REAL(wp), DIMENSION( 9) :: zmax69 REAL(wp), DIMENSION(3) :: zmax 70 70 LOGICAL :: ll_wrtstp, ll_colruns, ll_wrtruns 71 71 CHARACTER(len=20) :: clname … … 91 91 istatus = NF90_DEF_VAR( idrun, 'abs_ssh_max', NF90_DOUBLE, (/ idtime /), idssh ) 92 92 istatus = NF90_DEF_VAR( idrun, 'abs_u_max', NF90_DOUBLE, (/ idtime /), idu ) 93 istatus = NF90_DEF_VAR( idrun, 's_min', NF90_DOUBLE, (/ idtime /), ids1 )94 istatus = NF90_DEF_VAR( idrun, 's_max', NF90_DOUBLE, (/ idtime /), ids2 )95 istatus = NF90_DEF_VAR( idrun, 't_min', NF90_DOUBLE, (/ idtime /), idt1 )96 istatus = NF90_DEF_VAR( idrun, 't_max', NF90_DOUBLE, (/ idtime /), idt2 )97 IF( ln_zad_Aimp ) THEN98 istatus = NF90_DEF_VAR( idrun, 'abs_wi_max', NF90_DOUBLE, (/ idtime /), idw1 )99 istatus = NF90_DEF_VAR( idrun, 'Cf_max', NF90_DOUBLE, (/ idtime /), idc1 )100 ENDIF101 93 istatus = NF90_ENDDEF(idrun) 102 zmax(8:9) = 0._wp ! initialise to zero in case ln_zad_Aimp option is not in use103 94 ENDIF 104 95 ENDIF … … 114 105 zmax(1) = MAXVAL( ABS( ssh(:,:,Kmm) + ssh_ref*tmask(:,:,1) ) ) ! ssh max 115 106 ELSE 116 zmax(1) = M AXVAL( ABS( ssh(:,:,Kmm) ) ) ! ssh max107 zmax(1) = MINVAL( e3t(:,:,1,Kmm) ) ! ssh min 117 108 ENDIF 118 109 zmax(2) = MAXVAL( ABS( uu(:,:,:,Kmm) ) ) ! velocity max (zonal only) 119 zmax(3) = MAXVAL( -ts(:,:,:,jp_sal,Kmm) , mask = tmask(:,:,:) == 1._wp ) ! minus salinity max 120 zmax(4) = MAXVAL( ts(:,:,:,jp_sal,Kmm) , mask = tmask(:,:,:) == 1._wp ) ! salinity max 121 zmax(5) = MAXVAL( -ts(:,:,:,jp_tem,Kmm) , mask = tmask(:,:,:) == 1._wp ) ! minus temperature max 122 zmax(6) = MAXVAL( ts(:,:,:,jp_tem,Kmm) , mask = tmask(:,:,:) == 1._wp ) ! temperature max 123 zmax(7) = REAL( nstop , wp ) ! stop indicator 124 IF( ln_zad_Aimp ) THEN 125 zmax(8) = MAXVAL( ABS( wi(:,:,:) ) , mask = wmask(:,:,:) == 1._wp ) ! implicit vertical vel. max 126 zmax(9) = MAXVAL( Cu_adv(:,:,:) , mask = tmask(:,:,:) == 1._wp ) ! partitioning coeff. max 127 ENDIF 110 zmax(3) = REAL( nstop , wp ) ! stop indicator 128 111 ! 129 112 IF( ll_colruns ) THEN 130 113 CALL mpp_max( "stpctl", zmax ) ! max over the global domain 131 nstop = NINT( zmax( 7) ) ! nstop indicator sheared among all local domains114 nstop = NINT( zmax(3) ) ! nstop indicator sheared among all local domains 132 115 ENDIF 133 116 ! !== run statistics ==! ("run.stat" files) 134 117 IF( ll_wrtruns ) THEN 135 WRITE(numrun,9500) kt, zmax(1), zmax(2) , -zmax(3), zmax(4)118 WRITE(numrun,9500) kt, zmax(1), zmax(2) 136 119 istatus = NF90_PUT_VAR( idrun, idssh, (/ zmax(1)/), (/kt/), (/1/) ) 137 120 istatus = NF90_PUT_VAR( idrun, idu, (/ zmax(2)/), (/kt/), (/1/) ) 138 istatus = NF90_PUT_VAR( idrun, ids1, (/-zmax(3)/), (/kt/), (/1/) )139 istatus = NF90_PUT_VAR( idrun, ids2, (/ zmax(4)/), (/kt/), (/1/) )140 istatus = NF90_PUT_VAR( idrun, idt1, (/-zmax(5)/), (/kt/), (/1/) )141 istatus = NF90_PUT_VAR( idrun, idt2, (/ zmax(6)/), (/kt/), (/1/) )142 IF( ln_zad_Aimp ) THEN143 istatus = NF90_PUT_VAR( idrun, idw1, (/ zmax(8)/), (/kt/), (/1/) )144 istatus = NF90_PUT_VAR( idrun, idc1, (/ zmax(9)/), (/kt/), (/1/) )145 ENDIF146 121 IF( MOD( kt , 100 ) == 0 ) istatus = NF90_SYNC(idrun) 147 122 IF( kt == nitend ) istatus = NF90_CLOSE(idrun) … … 149 124 ! !== error handling ==! 150 125 IF( ( sn_cfctl%l_glochk .OR. lsomeoce ) .AND. ( & ! domain contains some ocean points, check for sensible ranges 151 & zmax(1) > 20._wp .OR. & ! too large sea surface height ( > 20 m )126 & zmax(1) < 0._wp .OR. & ! negative sea surface height 152 127 & zmax(2) > 10._wp .OR. & ! too large velocity ( > 10 m/s) 153 & zmax(3) >= 0._wp .OR. & ! negative or zero sea surface salinity 154 & zmax(4) >= 100._wp .OR. & ! too large sea surface salinity ( > 100 ) 155 & zmax(4) < 0._wp .OR. & ! too large sea surface salinity (keep this line for sea-ice) 156 & ISNAN( zmax(1) + zmax(2) + zmax(3) ) ) ) THEN ! NaN encounter in the tests 128 & ISNAN( zmax(1) + zmax(2) ) ) ) THEN ! NaN encounter in the tests 157 129 IF( lk_mpp .AND. sn_cfctl%l_glochk ) THEN 158 130 ! have use mpp_max (because sn_cfctl%l_glochk=.T. and distributed) 159 131 CALL mpp_maxloc( 'stpctl', ABS(ssh(:,:,Kmm)) , ssmask(:,:) , zzz, ih ) 160 132 CALL mpp_maxloc( 'stpctl', ABS(uu(:,:,:,Kmm)) , umask (:,:,:), zzz, iu ) 161 CALL mpp_minloc( 'stpctl', ts(:,:,:,jp_sal,Kmm), tmask (:,:,:), zzz, is1 )162 CALL mpp_maxloc( 'stpctl', ts(:,:,:,jp_sal,Kmm), tmask (:,:,:), zzz, is2 )163 133 ELSE 164 134 ! find local min and max locations 165 135 ih(:) = MAXLOC( ABS( ssh(:,:,Kmm) ) ) + (/ nimpp - 1, njmpp - 1 /) 166 136 iu(:) = MAXLOC( ABS( uu (:,:,:,Kmm) ) ) + (/ nimpp - 1, njmpp - 1, 0 /) 167 is1(:) = MINLOC( ts(:,:,:,jp_sal,Kmm), mask = tmask(:,:,:) == 1._wp ) + (/ nimpp - 1, njmpp - 1, 0 /)168 is2(:) = MAXLOC( ts(:,:,:,jp_sal,Kmm), mask = tmask(:,:,:) == 1._wp ) + (/ nimpp - 1, njmpp - 1, 0 /)169 137 ENDIF 170 138 171 WRITE(ctmp1,*) ' stp_ctl: |ssh| > 20 m or |U| > 10 m/s or S <= 0 or S >= 100or NaN encounter in the tests'139 WRITE(ctmp1,*) ' stp_ctl: (e3t0) ssh < 0 m or |U| > 10 m/s or NaN encounter in the tests' 172 140 WRITE(ctmp2,9100) kt, zmax(1), ih(1) , ih(2) 173 141 WRITE(ctmp3,9200) kt, zmax(2), iu(1) , iu(2) , iu(3) 174 WRITE(ctmp4,9300) kt, - zmax(3), is1(1), is1(2), is1(3) 175 WRITE(ctmp5,9400) kt, zmax(4), is2(1), is2(2), is2(3) 176 WRITE(ctmp6,*) ' ===> output of last computed fields in output.abort.nc file' 142 WRITE(ctmp4,*) ' ===> output of last computed fields in output.abort.nc file' 177 143 178 144 CALL dia_wri_state( Kmm, 'output.abort' ) ! create an output.abort file … … 180 146 IF( .NOT. sn_cfctl%l_glochk ) THEN 181 147 WRITE(ctmp8,*) 'E R R O R message from sub-domain: ', narea 182 CALL ctl_stop( 'STOP', ctmp1, ' ', ctmp8, ' ', ctmp2, ctmp3, ctmp4 , ctmp5, ctmp6)148 CALL ctl_stop( 'STOP', ctmp1, ' ', ctmp8, ' ', ctmp2, ctmp3, ctmp4 ) 183 149 ELSE 184 CALL ctl_stop( ctmp1, ' ', ctmp2, ctmp3, ctmp4 , ctmp5, ' ', ctmp6, ' ')150 CALL ctl_stop( ctmp1, ' ', ctmp2, ctmp3, ctmp4 ) 185 151 ENDIF 186 152 … … 189 155 ENDIF 190 156 ! 191 9100 FORMAT (' kt=',i8,' |ssh| m ax: ',1pg11.4,', at i j : ',2i5)157 9100 FORMAT (' kt=',i8,' |ssh| min: ',1pg11.4,', at i j : ',2i5) 192 158 9200 FORMAT (' kt=',i8,' |U| max: ',1pg11.4,', at i j k: ',3i5) 193 9300 FORMAT (' kt=',i8,' S min: ',1pg11.4,', at i j k: ',3i5) 194 9400 FORMAT (' kt=',i8,' S max: ',1pg11.4,', at i j k: ',3i5) 195 9500 FORMAT(' it :', i8, ' |ssh|_max: ', D23.16, ' |U|_max: ', D23.16,' S_min: ', D23.16,' S_max: ', D23.16) 159 9500 FORMAT(' it :', i8, ' |ssh|_max: ', D23.16, ' |U|_max: ', D23.16) 196 160 ! 197 161 END SUBROUTINE stp_ctl
Note: See TracChangeset
for help on using the changeset viewer.
