- Timestamp:
- 2020-07-01T16:09:00+02:00 (4 years ago)
- Location:
- NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/tests/CANAL/MY_SRC
- Files:
-
- 1 deleted
- 8 edited
-
diawri.F90 (modified) (36 diffs)
-
domvvl.F90 (modified) (19 diffs)
-
sbcmod.F90 (deleted)
-
stpctl.F90 (modified) (5 diffs)
-
trazdf.F90 (modified) (6 diffs)
-
usrdef_hgr.F90 (modified) (2 diffs)
-
usrdef_istate.F90 (modified) (6 diffs)
-
usrdef_sbc.F90 (modified) (3 diffs)
-
usrdef_zgr.F90 (modified) (1 diff)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/tests/CANAL/MY_SRC/diawri.F90
r12489 r13197 26 26 !!---------------------------------------------------------------------- 27 27 USE oce ! ocean dynamics and tracers 28 USE isf_oce 29 USE isfcpl 30 USE abl ! abl variables in case ln_abl = .true. 28 31 USE dom_oce ! ocean space and time domain 29 32 USE phycst ! physical constants … … 65 68 PUBLIC dia_wri_state 66 69 PUBLIC dia_wri_alloc ! Called by nemogcm module 67 70 #if ! defined key_iomput 71 PUBLIC dia_wri_alloc_abl ! Called by sbcabl module (if ln_abl = .true.) 72 #endif 68 73 INTEGER :: nid_T, nz_T, nh_T, ndim_T, ndim_hT ! grid_T file 69 74 INTEGER :: nb_T , ndim_bT ! grid_T file … … 71 76 INTEGER :: nid_V, nz_V, nh_V, ndim_V, ndim_hV ! grid_V file 72 77 INTEGER :: nid_W, nz_W, nh_W ! grid_W file 78 INTEGER :: nid_A, nz_A, nh_A, ndim_A, ndim_hA ! grid_ABL file 73 79 INTEGER :: ndex(1) ! ??? 74 80 INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_hT, ndex_hU, ndex_hV 81 INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_hA, ndex_A ! ABL 75 82 INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_T, ndex_U, ndex_V 76 83 INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_bT 77 84 85 !! * Substitutions 86 # include "do_loop_substitute.h90" 78 87 !!---------------------------------------------------------------------- 79 88 !! NEMO/OCE 4.0 , NEMO Consortium (2018) … … 147 156 CALL iom_put( "sst", ts(:,:,1,jp_tem,Kmm) ) ! surface temperature 148 157 IF ( iom_use("sbt") ) THEN 149 DO jj = 1, jpj 150 DO ji = 1, jpi 151 ikbot = mbkt(ji,jj) 152 z2d(ji,jj) = ts(ji,jj,ikbot,jp_tem,Kmm) 153 END DO 154 END DO 158 DO_2D_11_11 159 ikbot = mbkt(ji,jj) 160 z2d(ji,jj) = ts(ji,jj,ikbot,jp_tem,Kmm) 161 END_2D 155 162 CALL iom_put( "sbt", z2d ) ! bottom temperature 156 163 ENDIF … … 159 166 CALL iom_put( "sss", ts(:,:,1,jp_sal,Kmm) ) ! surface salinity 160 167 IF ( iom_use("sbs") ) THEN 161 DO jj = 1, jpj 162 DO ji = 1, jpi 163 ikbot = mbkt(ji,jj) 164 z2d(ji,jj) = ts(ji,jj,ikbot,jp_sal,Kmm) 165 END DO 166 END DO 168 DO_2D_11_11 169 ikbot = mbkt(ji,jj) 170 z2d(ji,jj) = ts(ji,jj,ikbot,jp_sal,Kmm) 171 END_2D 167 172 CALL iom_put( "sbs", z2d ) ! bottom salinity 168 173 ENDIF … … 171 176 zztmp = rho0 * 0.25 172 177 z2d(:,:) = 0._wp 173 DO jj = 2, jpjm1 174 DO ji = fs_2, fs_jpim1 ! vector opt. 175 zztmp2 = ( ( rCdU_bot(ji+1,jj)+rCdU_bot(ji ,jj) ) * uu(ji ,jj,mbku(ji ,jj),Kmm) )**2 & 176 & + ( ( rCdU_bot(ji ,jj)+rCdU_bot(ji-1,jj) ) * uu(ji-1,jj,mbku(ji-1,jj),Kmm) )**2 & 177 & + ( ( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj ) ) * vv(ji,jj ,mbkv(ji,jj ),Kmm) )**2 & 178 & + ( ( rCdU_bot(ji,jj )+rCdU_bot(ji,jj-1) ) * vv(ji,jj-1,mbkv(ji,jj-1),Kmm) )**2 179 z2d(ji,jj) = zztmp * SQRT( zztmp2 ) * tmask(ji,jj,1) 180 ! 181 END DO 182 END DO 178 DO_2D_00_00 179 zztmp2 = ( ( rCdU_bot(ji+1,jj)+rCdU_bot(ji ,jj) ) * uu(ji ,jj,mbku(ji ,jj),Kmm) )**2 & 180 & + ( ( rCdU_bot(ji ,jj)+rCdU_bot(ji-1,jj) ) * uu(ji-1,jj,mbku(ji-1,jj),Kmm) )**2 & 181 & + ( ( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj ) ) * vv(ji,jj ,mbkv(ji,jj ),Kmm) )**2 & 182 & + ( ( rCdU_bot(ji,jj )+rCdU_bot(ji,jj-1) ) * vv(ji,jj-1,mbkv(ji,jj-1),Kmm) )**2 183 z2d(ji,jj) = zztmp * SQRT( zztmp2 ) * tmask(ji,jj,1) 184 ! 185 END_2D 183 186 CALL lbc_lnk( 'diawri', z2d, 'T', 1. ) 184 187 CALL iom_put( "taubot", z2d ) … … 188 191 CALL iom_put( "ssu", uu(:,:,1,Kmm) ) ! surface i-current 189 192 IF ( iom_use("sbu") ) THEN 190 DO jj = 1, jpj 191 DO ji = 1, jpi 192 ikbot = mbku(ji,jj) 193 z2d(ji,jj) = uu(ji,jj,ikbot,Kmm) 194 END DO 195 END DO 193 DO_2D_11_11 194 ikbot = mbku(ji,jj) 195 z2d(ji,jj) = uu(ji,jj,ikbot,Kmm) 196 END_2D 196 197 CALL iom_put( "sbu", z2d ) ! bottom i-current 197 198 ENDIF … … 200 201 CALL iom_put( "ssv", vv(:,:,1,Kmm) ) ! surface j-current 201 202 IF ( iom_use("sbv") ) THEN 202 DO jj = 1, jpj 203 DO ji = 1, jpi 204 ikbot = mbkv(ji,jj) 205 z2d(ji,jj) = vv(ji,jj,ikbot,Kmm) 206 END DO 207 END DO 203 DO_2D_11_11 204 ikbot = mbkv(ji,jj) 205 z2d(ji,jj) = vv(ji,jj,ikbot,Kmm) 206 END_2D 208 207 CALL iom_put( "sbv", z2d ) ! bottom j-current 209 208 ENDIF 210 209 210 IF( ln_zad_Aimp ) ww = ww + wi ! Recombine explicit and implicit parts of vertical velocity for diagnostic output 211 ! 211 212 CALL iom_put( "woce", ww ) ! vertical velocity 212 213 IF( iom_use('w_masstr') .OR. iom_use('w_masstr2') ) THEN ! vertical mass transport & its square value … … 219 220 IF( iom_use('w_masstr2') ) CALL iom_put( "w_masstr2", z3d(:,:,:) * z3d(:,:,:) ) 220 221 ENDIF 222 ! 223 IF( ln_zad_Aimp ) ww = ww - wi ! Remove implicit part of vertical velocity that was added for diagnostic output 221 224 222 225 CALL iom_put( "avt" , avt ) ! T vert. eddy diff. coef. … … 227 230 IF( iom_use('logavs') ) CALL iom_put( "logavs", LOG( MAX( 1.e-20_wp, avs(:,:,:) ) ) ) 228 231 229 IF ( iom_use("salgrad") .OR. iom_use("salgrad2") ) THEN230 z3d(:,:,jpk) = 0.231 DO jk = 1, jpkm1232 DO jj = 2, jpjm1 ! sal gradient233 DO ji = fs_2, fs_jpim1 ! vector opt.234 zztmp = ts(ji,jj,jk,jp_sal,Kmm)235 zztmpx = ( ts(ji+1,jj,jk,jp_sal,Kmm) - zztmp ) * r1_e1u(ji,jj) + ( zztmp - ts(ji-1,jj ,jk,jp_sal,Kmm) ) * r1_e1u(ji-1,jj)236 zztmpy = ( ts(ji,jj+1,jk,jp_sal,Kmm) - zztmp ) * r1_e2v(ji,jj) + ( zztmp - ts(ji ,jj-1,jk,jp_sal,Kmm) ) * r1_e2v(ji,jj-1)237 z3d(ji,jj,jk) = 0.25 * ( zztmpx * zztmpx + zztmpy * zztmpy ) &238 & * umask(ji,jj,jk) * umask(ji-1,jj,jk) * vmask(ji,jj,jk) * umask(ji,jj-1,jk)239 END DO240 END DO241 END DO242 CALL lbc_lnk( 'diawri', z3d, 'T', 1. )243 CALL iom_put( "salgrad2", z3d ) ! square of module of sal gradient244 z3d(:,:,:) = SQRT( z3d(:,:,:) )245 CALL iom_put( "salgrad" , z3d ) ! module of sal gradient246 ENDIF247 248 232 IF ( iom_use("sstgrad") .OR. iom_use("sstgrad2") ) THEN 249 DO jj = 2, jpjm1 ! sst gradient 250 DO ji = fs_2, fs_jpim1 ! vector opt. 251 zztmp = ts(ji,jj,1,jp_tem,Kmm) 252 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) 253 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) 254 z2d(ji,jj) = 0.25 * ( zztmpx * zztmpx + zztmpy * zztmpy ) & 255 & * umask(ji,jj,1) * umask(ji-1,jj,1) * vmask(ji,jj,1) * umask(ji,jj-1,1) 256 END DO 257 END DO 233 DO_2D_00_00 234 zztmp = ts(ji,jj,1,jp_tem,Kmm) 235 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) 236 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) 237 z2d(ji,jj) = 0.25 * ( zztmpx * zztmpx + zztmpy * zztmpy ) & 238 & * umask(ji,jj,1) * umask(ji-1,jj,1) * vmask(ji,jj,1) * umask(ji,jj-1,1) 239 END_2D 258 240 CALL lbc_lnk( 'diawri', z2d, 'T', 1. ) 259 241 CALL iom_put( "sstgrad2", z2d ) ! square of module of sst gradient … … 265 247 IF( iom_use("heatc") ) THEN 266 248 z2d(:,:) = 0._wp 267 DO jk = 1, jpkm1 268 DO jj = 1, jpj 269 DO ji = 1, jpi 270 z2d(ji,jj) = z2d(ji,jj) + e3t(ji,jj,jk,Kmm) * ts(ji,jj,jk,jp_tem,Kmm) * tmask(ji,jj,jk) 271 END DO 272 END DO 273 END DO 249 DO_3D_11_11( 1, jpkm1 ) 250 z2d(ji,jj) = z2d(ji,jj) + e3t(ji,jj,jk,Kmm) * ts(ji,jj,jk,jp_tem,Kmm) * tmask(ji,jj,jk) 251 END_3D 274 252 CALL iom_put( "heatc", rho0_rcp * z2d ) ! vertically integrated heat content (J/m2) 275 253 ENDIF … … 277 255 IF( iom_use("saltc") ) THEN 278 256 z2d(:,:) = 0._wp 279 DO jk = 1, jpkm1 280 DO jj = 1, jpj 281 DO ji = 1, jpi 282 z2d(ji,jj) = z2d(ji,jj) + e3t(ji,jj,jk,Kmm) * ts(ji,jj,jk,jp_sal,Kmm) * tmask(ji,jj,jk) 283 END DO 284 END DO 285 END DO 257 DO_3D_11_11( 1, jpkm1 ) 258 z2d(ji,jj) = z2d(ji,jj) + e3t(ji,jj,jk,Kmm) * ts(ji,jj,jk,jp_sal,Kmm) * tmask(ji,jj,jk) 259 END_3D 286 260 CALL iom_put( "saltc", rho0 * z2d ) ! vertically integrated salt content (PSU*kg/m2) 287 261 ENDIF … … 289 263 IF( iom_use("salt2c") ) THEN 290 264 z2d(:,:) = 0._wp 291 DO jk = 1, jpkm1 292 DO jj = 1, jpj 293 DO ji = 1, jpi 294 z2d(ji,jj) = z2d(ji,jj) + e3t(ji,jj,jk,Kmm) * ts(ji,jj,jk,jp_sal,Kmm) * ts(ji,jj,jk,jp_sal,Kmm) * tmask(ji,jj,jk) 295 END DO 296 END DO 297 END DO 265 DO_3D_11_11( 1, jpkm1 ) 266 z2d(ji,jj) = z2d(ji,jj) + e3t(ji,jj,jk,Kmm) * ts(ji,jj,jk,jp_sal,Kmm) * ts(ji,jj,jk,jp_sal,Kmm) * tmask(ji,jj,jk) 267 END_3D 298 268 CALL iom_put( "salt2c", rho0 * z2d ) ! vertically integrated salt content (PSU*kg/m2) 299 269 ENDIF … … 301 271 IF ( iom_use("eken") ) THEN 302 272 z3d(:,:,jpk) = 0._wp 303 DO jk = 1, jpkm1 304 DO jj = 2, jpjm1 305 DO ji = fs_2, fs_jpim1 ! vector opt. 306 zztmp = 0.25_wp * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm) 307 z3d(ji,jj,jk) = zztmp * ( uu(ji-1,jj,jk,Kmm)**2 * e2u(ji-1,jj) * e3u(ji-1,jj,jk,Kmm) & 308 & + uu(ji ,jj,jk,Kmm)**2 * e2u(ji ,jj) * e3u(ji ,jj,jk,Kmm) & 309 & + vv(ji,jj-1,jk,Kmm)**2 * e1v(ji,jj-1) * e3v(ji,jj-1,jk,Kmm) & 310 & + vv(ji,jj ,jk,Kmm)**2 * e1v(ji,jj ) * e3v(ji,jj ,jk,Kmm) ) 311 END DO 312 END DO 313 END DO 273 DO_3D_00_00( 1, jpkm1 ) 274 zztmp = 0.25_wp * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm) 275 z3d(ji,jj,jk) = zztmp * ( uu(ji-1,jj,jk,Kmm)**2 * e2u(ji-1,jj) * e3u(ji-1,jj,jk,Kmm) & 276 & + uu(ji ,jj,jk,Kmm)**2 * e2u(ji ,jj) * e3u(ji ,jj,jk,Kmm) & 277 & + vv(ji,jj-1,jk,Kmm)**2 * e1v(ji,jj-1) * e3v(ji,jj-1,jk,Kmm) & 278 & + vv(ji,jj ,jk,Kmm)**2 * e1v(ji,jj ) * e3v(ji,jj ,jk,Kmm) ) 279 END_3D 314 280 CALL lbc_lnk( 'diawri', z3d, 'T', 1. ) 315 281 CALL iom_put( "eken", z3d ) ! kinetic energy … … 321 287 z3d(1,:, : ) = 0._wp 322 288 z3d(:,1, : ) = 0._wp 323 DO jk = 1, jpkm1 324 DO jj = 2, jpj 325 DO ji = 2, jpi 326 z3d(ji,jj,jk) = 0.25_wp * ( uu(ji ,jj,jk,Kmm) * uu(ji ,jj,jk,Kmm) * e1e2u(ji ,jj) * e3u(ji ,jj,jk,Kmm) & 327 & + uu(ji-1,jj,jk,Kmm) * uu(ji-1,jj,jk,Kmm) * e1e2u(ji-1,jj) * e3u(ji-1,jj,jk,Kmm) & 328 & + vv(ji,jj ,jk,Kmm) * vv(ji,jj ,jk,Kmm) * e1e2v(ji,jj ) * e3v(ji,jj ,jk,Kmm) & 329 & + vv(ji,jj-1,jk,Kmm) * vv(ji,jj-1,jk,Kmm) * e1e2v(ji,jj-1) * e3v(ji,jj-1,jk,Kmm) ) & 330 & * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm) * tmask(ji,jj,jk) 331 END DO 332 END DO 333 END DO 334 289 DO_3D_00_00( 1, jpkm1 ) 290 z3d(ji,jj,jk) = 0.25_wp * ( uu(ji ,jj,jk,Kmm) * uu(ji ,jj,jk,Kmm) * e1e2u(ji ,jj) * e3u(ji ,jj,jk,Kmm) & 291 & + uu(ji-1,jj,jk,Kmm) * uu(ji-1,jj,jk,Kmm) * e1e2u(ji-1,jj) * e3u(ji-1,jj,jk,Kmm) & 292 & + vv(ji,jj ,jk,Kmm) * vv(ji,jj ,jk,Kmm) * e1e2v(ji,jj ) * e3v(ji,jj ,jk,Kmm) & 293 & + vv(ji,jj-1,jk,Kmm) * vv(ji,jj-1,jk,Kmm) * e1e2v(ji,jj-1) * e3v(ji,jj-1,jk,Kmm) ) & 294 & * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm) * tmask(ji,jj,jk) 295 END_3D 335 296 CALL lbc_lnk( 'diawri', z3d, 'T', 1. ) 336 297 CALL iom_put( "ke", z3d ) ! kinetic energy 337 298 338 299 z2d(:,:) = 0._wp 339 DO jk = 1, jpkm1 340 DO jj = 1, jpj 341 DO ji = 1, jpi 342 z2d(ji,jj) = z2d(ji,jj) + e3t(ji,jj,jk,Kmm) * z3d(ji,jj,jk) * tmask(ji,jj,jk) 343 END DO 344 END DO 345 END DO 300 DO_3D_11_11( 1, jpkm1 ) 301 z2d(ji,jj) = z2d(ji,jj) + e3t(ji,jj,jk,Kmm) * z3d(ji,jj,jk) * tmask(ji,jj,jk) 302 END_3D 346 303 CALL iom_put( "ke_zint", z2d ) ! vertically integrated kinetic energy 347 304 … … 353 310 354 311 z3d(:,:,jpk) = 0._wp 355 DO jk = 1, jpkm1 356 DO jj = 1, jpjm1 357 DO ji = 1, fs_jpim1 ! vector opt. 358 z3d(ji,jj,jk) = ( e2v(ji+1,jj ) * vv(ji+1,jj ,jk,Kmm) - e2v(ji,jj) * vv(ji,jj,jk,Kmm) & 359 & - e1u(ji ,jj+1) * uu(ji ,jj+1,jk,Kmm) + e1u(ji,jj) * uu(ji,jj,jk,Kmm) ) * r1_e1e2f(ji,jj) 360 END DO 361 END DO 362 END DO 312 DO_3D_00_00( 1, jpkm1 ) 313 z3d(ji,jj,jk) = ( e2v(ji+1,jj ) * vv(ji+1,jj ,jk,Kmm) - e2v(ji,jj) * vv(ji,jj,jk,Kmm) & 314 & - e1u(ji ,jj+1) * uu(ji ,jj+1,jk,Kmm) + e1u(ji,jj) * uu(ji,jj,jk,Kmm) ) * r1_e1e2f(ji,jj) 315 END_3D 363 316 CALL lbc_lnk( 'diawri', z3d, 'F', 1. ) 364 317 CALL iom_put( "relvor", z3d ) ! relative vorticity 365 318 366 DO jk = 1, jpkm1 367 DO jj = 1, jpj 368 DO ji = 1, jpi 369 z3d(ji,jj,jk) = ff_f(ji,jj) + z3d(ji,jj,jk) 370 END DO 371 END DO 372 END DO 319 DO_3D_11_11( 1, jpkm1 ) 320 z3d(ji,jj,jk) = ff_f(ji,jj) + z3d(ji,jj,jk) 321 END_3D 373 322 CALL iom_put( "absvor", z3d ) ! absolute vorticity 374 323 375 DO jk = 1, jpkm1 376 DO jj = 1, jpjm1 377 DO ji = 1, fs_jpim1 ! vector opt. 378 ze3 = ( e3t(ji,jj+1,jk,Kmm)*tmask(ji,jj+1,jk) + e3t(ji+1,jj+1,jk,Kmm)*tmask(ji+1,jj+1,jk) & 379 & + e3t(ji,jj ,jk,Kmm)*tmask(ji,jj ,jk) + e3t(ji+1,jj ,jk,Kmm)*tmask(ji+1,jj ,jk) ) 380 IF( ze3 /= 0._wp ) THEN ; ze3 = 4._wp / ze3 381 ELSE ; ze3 = 0._wp 382 ENDIF 383 z3d(ji,jj,jk) = ze3 * z3d(ji,jj,jk) 384 END DO 385 END DO 386 END DO 324 DO_3D_00_00( 1, jpkm1 ) 325 ze3 = ( e3t(ji,jj+1,jk,Kmm)*tmask(ji,jj+1,jk) + e3t(ji+1,jj+1,jk,Kmm)*tmask(ji+1,jj+1,jk) & 326 & + e3t(ji,jj ,jk,Kmm)*tmask(ji,jj ,jk) + e3t(ji+1,jj ,jk,Kmm)*tmask(ji+1,jj ,jk) ) 327 IF( ze3 /= 0._wp ) THEN ; ze3 = 4._wp / ze3 328 ELSE ; ze3 = 0._wp 329 ENDIF 330 z3d(ji,jj,jk) = ze3 * z3d(ji,jj,jk) 331 END_3D 387 332 CALL lbc_lnk( 'diawri', z3d, 'F', 1. ) 388 333 CALL iom_put( "potvor", z3d ) ! potential vorticity 389 334 390 335 ENDIF 391 392 336 ! 393 337 IF( iom_use("u_masstr") .OR. iom_use("u_masstr_vint") .OR. iom_use("u_heattr") .OR. iom_use("u_salttr") ) THEN … … 404 348 IF( iom_use("u_heattr") ) THEN 405 349 z2d(:,:) = 0._wp 406 DO jk = 1, jpkm1 407 DO jj = 2, jpjm1 408 DO ji = fs_2, fs_jpim1 ! vector opt. 409 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) ) 410 END DO 411 END DO 412 END DO 350 DO_3D_00_00( 1, jpkm1 ) 351 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) ) 352 END_3D 413 353 CALL lbc_lnk( 'diawri', z2d, 'U', -1. ) 414 354 CALL iom_put( "u_heattr", 0.5*rcp * z2d ) ! heat transport in i-direction … … 417 357 IF( iom_use("u_salttr") ) THEN 418 358 z2d(:,:) = 0.e0 419 DO jk = 1, jpkm1 420 DO jj = 2, jpjm1 421 DO ji = fs_2, fs_jpim1 ! vector opt. 422 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) ) 423 END DO 424 END DO 425 END DO 359 DO_3D_00_00( 1, jpkm1 ) 360 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) ) 361 END_3D 426 362 CALL lbc_lnk( 'diawri', z2d, 'U', -1. ) 427 363 CALL iom_put( "u_salttr", 0.5 * z2d ) ! heat transport in i-direction … … 439 375 IF( iom_use("v_heattr") ) THEN 440 376 z2d(:,:) = 0.e0 441 DO jk = 1, jpkm1 442 DO jj = 2, jpjm1 443 DO ji = fs_2, fs_jpim1 ! vector opt. 444 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) ) 445 END DO 446 END DO 447 END DO 377 DO_3D_00_00( 1, jpkm1 ) 378 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) ) 379 END_3D 448 380 CALL lbc_lnk( 'diawri', z2d, 'V', -1. ) 449 381 CALL iom_put( "v_heattr", 0.5*rcp * z2d ) ! heat transport in j-direction … … 452 384 IF( iom_use("v_salttr") ) THEN 453 385 z2d(:,:) = 0._wp 454 DO jk = 1, jpkm1 455 DO jj = 2, jpjm1 456 DO ji = fs_2, fs_jpim1 ! vector opt. 457 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) ) 458 END DO 459 END DO 460 END DO 386 DO_3D_00_00( 1, jpkm1 ) 387 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) ) 388 END_3D 461 389 CALL lbc_lnk( 'diawri', z2d, 'V', -1. ) 462 390 CALL iom_put( "v_salttr", 0.5 * z2d ) ! heat transport in j-direction … … 465 393 IF( iom_use("tosmint") ) THEN 466 394 z2d(:,:) = 0._wp 467 DO jk = 1, jpkm1 468 DO jj = 2, jpjm1 469 DO ji = fs_2, fs_jpim1 ! vector opt. 470 z2d(ji,jj) = z2d(ji,jj) + e3t(ji,jj,jk,Kmm) * ts(ji,jj,jk,jp_tem,Kmm) 471 END DO 472 END DO 473 END DO 395 DO_3D_00_00( 1, jpkm1 ) 396 z2d(ji,jj) = z2d(ji,jj) + e3t(ji,jj,jk,Kmm) * ts(ji,jj,jk,jp_tem,Kmm) 397 END_3D 474 398 CALL lbc_lnk( 'diawri', z2d, 'T', -1. ) 475 399 CALL iom_put( "tosmint", rho0 * z2d ) ! Vertical integral of temperature … … 477 401 IF( iom_use("somint") ) THEN 478 402 z2d(:,:)=0._wp 479 DO jk = 1, jpkm1 480 DO jj = 2, jpjm1 481 DO ji = fs_2, fs_jpim1 ! vector opt. 482 z2d(ji,jj) = z2d(ji,jj) + e3t(ji,jj,jk,Kmm) * ts(ji,jj,jk,jp_sal,Kmm) 483 END DO 484 END DO 485 END DO 403 DO_3D_00_00( 1, jpkm1 ) 404 z2d(ji,jj) = z2d(ji,jj) + e3t(ji,jj,jk,Kmm) * ts(ji,jj,jk,jp_sal,Kmm) 405 END_3D 486 406 CALL lbc_lnk( 'diawri', z2d, 'T', -1. ) 487 407 CALL iom_put( "somint", rho0 * z2d ) ! Vertical integral of salinity … … 490 410 CALL iom_put( "bn2", rn2 ) ! Brunt-Vaisala buoyancy frequency (N^2) 491 411 ! 492 412 493 413 IF (ln_dia25h) CALL dia_25h( kt, Kmm ) ! 25h averaging 494 414 … … 506 426 INTEGER, DIMENSION(2) :: ierr 507 427 !!---------------------------------------------------------------------- 508 ierr = 0 509 ALLOCATE( ndex_hT(jpi*jpj) , ndex_T(jpi*jpj*jpk) , & 510 & ndex_hU(jpi*jpj) , ndex_U(jpi*jpj*jpk) , & 511 & ndex_hV(jpi*jpj) , ndex_V(jpi*jpj*jpk) , STAT=ierr(1) ) 428 IF( nn_write == -1 ) THEN 429 dia_wri_alloc = 0 430 ELSE 431 ierr = 0 432 ALLOCATE( ndex_hT(jpi*jpj) , ndex_T(jpi*jpj*jpk) , & 433 & ndex_hU(jpi*jpj) , ndex_U(jpi*jpj*jpk) , & 434 & ndex_hV(jpi*jpj) , ndex_V(jpi*jpj*jpk) , STAT=ierr(1) ) 512 435 ! 513 dia_wri_alloc = MAXVAL(ierr) 514 CALL mpp_sum( 'diawri', dia_wri_alloc ) 436 dia_wri_alloc = MAXVAL(ierr) 437 CALL mpp_sum( 'diawri', dia_wri_alloc ) 438 ! 439 ENDIF 515 440 ! 516 441 END FUNCTION dia_wri_alloc 442 443 INTEGER FUNCTION dia_wri_alloc_abl() 444 !!---------------------------------------------------------------------- 445 ALLOCATE( ndex_hA(jpi*jpj), ndex_A (jpi*jpj*jpkam1), STAT=dia_wri_alloc_abl) 446 CALL mpp_sum( 'diawri', dia_wri_alloc_abl ) 447 ! 448 END FUNCTION dia_wri_alloc_abl 517 449 518 450 … … 538 470 INTEGER :: ierr ! error code return from allocation 539 471 INTEGER :: iimi, iima, ipk, it, itmod, ijmi, ijma ! local integers 472 INTEGER :: ipka ! ABL 540 473 INTEGER :: jn, ierror ! local integers 541 474 REAL(wp) :: zsto, zout, zmax, zjulian ! local scalars … … 543 476 REAL(wp), DIMENSION(jpi,jpj) :: zw2d ! 2D workspace 544 477 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zw3d ! 3D workspace 478 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: zw3d_abl ! ABL 3D workspace 545 479 !!---------------------------------------------------------------------- 546 480 ! … … 576 510 ijmi = 1 ; ijma = jpj 577 511 ipk = jpk 512 IF(ln_abl) ipka = jpkam1 578 513 579 514 ! define time axis … … 678 613 & "m", ipk, gdepw_1d, nz_W, "down" ) 679 614 615 IF( ln_abl ) THEN 616 ! Define the ABL grid FILE ( nid_A ) 617 CALL dia_nam( clhstnam, nn_write, 'grid_ABL' ) 618 IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam ! filename 619 CALL histbeg( clhstnam, jpi, glamt, jpj, gphit, & ! Horizontal grid: glamt and gphit 620 & iimi, iima-iimi+1, ijmi, ijma-ijmi+1, & 621 & nit000-1, zjulian, rn_Dt, nh_A, nid_A, domain_id=nidom, snc4chunks=snc4set ) 622 CALL histvert( nid_A, "ght_abl", "Vertical T levels", & ! Vertical grid: gdept 623 & "m", ipka, ght_abl(2:jpka), nz_A, "up" ) 624 ! ! Index of ocean points 625 ALLOCATE( zw3d_abl(jpi,jpj,ipka) ) 626 zw3d_abl(:,:,:) = 1._wp 627 CALL wheneq( jpi*jpj*ipka, zw3d_abl, 1, 1., ndex_A , ndim_A ) ! volume 628 CALL wheneq( jpi*jpj , zw3d_abl, 1, 1., ndex_hA, ndim_hA ) ! surface 629 DEALLOCATE(zw3d_abl) 630 ENDIF 680 631 681 632 ! Declare all the output fields as NETCDF variables … … 727 678 CALL histdef( nid_T, "sowindsp", "wind speed at 10m" , "m/s" , & ! wndm 728 679 & jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout ) 729 ! 680 ! 681 IF( ln_abl ) THEN 682 CALL histdef( nid_A, "t_abl", "Potential Temperature" , "K" , & ! t_abl 683 & jpi, jpj, nh_A, ipka, 1, ipka, nz_A, 32, clop, zsto, zout ) 684 CALL histdef( nid_A, "q_abl", "Humidity" , "kg/kg" , & ! q_abl 685 & jpi, jpj, nh_A, ipka, 1, ipka, nz_A, 32, clop, zsto, zout ) 686 CALL histdef( nid_A, "u_abl", "Atmospheric U-wind " , "m/s" , & ! u_abl 687 & jpi, jpj, nh_A, ipka, 1, ipka, nz_A, 32, clop, zsto, zout ) 688 CALL histdef( nid_A, "v_abl", "Atmospheric V-wind " , "m/s" , & ! v_abl 689 & jpi, jpj, nh_A, ipka, 1, ipka, nz_A, 32, clop, zsto, zout ) 690 CALL histdef( nid_A, "tke_abl", "Atmospheric TKE " , "m2/s2" , & ! tke_abl 691 & jpi, jpj, nh_A, ipka, 1, ipka, nz_A, 32, clop, zsto, zout ) 692 CALL histdef( nid_A, "avm_abl", "Atmospheric turbulent viscosity", "m2/s" , & ! avm_abl 693 & jpi, jpj, nh_A, ipka, 1, ipka, nz_A, 32, clop, zsto, zout ) 694 CALL histdef( nid_A, "avt_abl", "Atmospheric turbulent diffusivity", "m2/s2", & ! avt_abl 695 & jpi, jpj, nh_A, ipka, 1, ipka, nz_A, 32, clop, zsto, zout ) 696 CALL histdef( nid_A, "pblh", "Atmospheric boundary layer height " , "m", & ! pblh 697 & jpi, jpj, nh_A, 1 , 1, 1 , -99 , 32, clop, zsto, zout ) 698 #if defined key_si3 699 CALL histdef( nid_A, "oce_frac", "Fraction of open ocean" , " ", & ! ato_i 700 & jpi, jpj, nh_A, 1 , 1, 1 , -99 , 32, clop, zsto, zout ) 701 #endif 702 CALL histend( nid_A, snc4chunks=snc4set ) 703 ENDIF 704 ! 730 705 IF( ln_icebergs ) THEN 731 706 CALL histdef( nid_T, "calving" , "calving mass input" , "kg/s" , & … … 885 860 CALL histwrite( nid_T, "soicecov", it, fr_i , ndim_hT, ndex_hT ) ! ice fraction 886 861 CALL histwrite( nid_T, "sowindsp", it, wndm , ndim_hT, ndex_hT ) ! wind speed 887 ! 862 ! 863 IF( ln_abl ) THEN 864 ALLOCATE( zw3d_abl(jpi,jpj,jpka) ) 865 IF( ln_mskland ) THEN 866 DO jk=1,jpka 867 zw3d_abl(:,:,jk) = tmask(:,:,1) 868 END DO 869 ELSE 870 zw3d_abl(:,:,:) = 1._wp 871 ENDIF 872 CALL histwrite( nid_A, "pblh" , it, pblh(:,:) *zw3d_abl(:,:,1 ), ndim_hA, ndex_hA ) ! pblh 873 CALL histwrite( nid_A, "u_abl" , it, u_abl (:,:,2:jpka,nt_n )*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A ) ! u_abl 874 CALL histwrite( nid_A, "v_abl" , it, v_abl (:,:,2:jpka,nt_n )*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A ) ! v_abl 875 CALL histwrite( nid_A, "t_abl" , it, tq_abl (:,:,2:jpka,nt_n,1)*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A ) ! t_abl 876 CALL histwrite( nid_A, "q_abl" , it, tq_abl (:,:,2:jpka,nt_n,2)*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A ) ! q_abl 877 CALL histwrite( nid_A, "tke_abl", it, tke_abl (:,:,2:jpka,nt_n )*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A ) ! tke_abl 878 CALL histwrite( nid_A, "avm_abl", it, avm_abl (:,:,2:jpka )*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A ) ! avm_abl 879 CALL histwrite( nid_A, "avt_abl", it, avt_abl (:,:,2:jpka )*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A ) ! avt_abl 880 #if defined key_si3 881 CALL histwrite( nid_A, "oce_frac" , it, ato_i(:,:) , ndim_hA, ndex_hA ) ! ato_i 882 #endif 883 DEALLOCATE(zw3d_abl) 884 ENDIF 885 ! 888 886 IF( ln_icebergs ) THEN 889 887 ! … … 931 929 CALL histwrite( nid_V, "sometauy", it, vtau , ndim_hV, ndex_hV ) ! j-wind stress 932 930 933 CALL histwrite( nid_W, "vovecrtz", it, ww , ndim_T, ndex_T ) ! vert. current 931 IF( ln_zad_Aimp ) THEN 932 CALL histwrite( nid_W, "vovecrtz", it, ww + wi , ndim_T, ndex_T ) ! vert. current 933 ELSE 934 CALL histwrite( nid_W, "vovecrtz", it, ww , ndim_T, ndex_T ) ! vert. current 935 ENDIF 934 936 CALL histwrite( nid_W, "votkeavt", it, avt , ndim_T, ndex_T ) ! T vert. eddy diff. coef. 935 937 CALL histwrite( nid_W, "votkeavm", it, avm , ndim_T, ndex_T ) ! T vert. eddy visc. coef. … … 951 953 CALL histclo( nid_V ) 952 954 CALL histclo( nid_W ) 955 IF(ln_abl) CALL histclo( nid_A ) 953 956 ENDIF 954 957 ! … … 974 977 CHARACTER (len=* ), INTENT( in ) :: cdfile_name ! name of the file created 975 978 !! 976 INTEGER :: inum 979 INTEGER :: inum, jk 977 980 !!---------------------------------------------------------------------- 978 981 ! … … 981 984 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~ and forcing fields file created ' 982 985 IF(lwp) WRITE(numout,*) ' and named :', cdfile_name, '...nc' 983 984 #if defined key_si3 985 CALL iom_open( TRIM(cdfile_name), inum, ldwrt = .TRUE., kdlev = jpl ) 986 #else 987 CALL iom_open( TRIM(cdfile_name), inum, ldwrt = .TRUE. ) 988 #endif 989 986 ! 987 CALL iom_open( TRIM(cdfile_name), inum, ldwrt = .TRUE. ) 988 ! 990 989 CALL iom_rstput( 0, 0, inum, 'votemper', ts(:,:,:,jp_tem,Kmm) ) ! now temperature 991 990 CALL iom_rstput( 0, 0, inum, 'vosaline', ts(:,:,:,jp_sal,Kmm) ) ! now salinity … … 993 992 CALL iom_rstput( 0, 0, inum, 'vozocrtx', uu(:,:,:,Kmm) ) ! now i-velocity 994 993 CALL iom_rstput( 0, 0, inum, 'vomecrty', vv(:,:,:,Kmm) ) ! now j-velocity 995 CALL iom_rstput( 0, 0, inum, 'vovecrtz', ww ) ! now k-velocity 994 IF( ln_zad_Aimp ) THEN 995 CALL iom_rstput( 0, 0, inum, 'vovecrtz', ww + wi ) ! now k-velocity 996 ELSE 997 CALL iom_rstput( 0, 0, inum, 'vovecrtz', ww ) ! now k-velocity 998 ENDIF 999 CALL iom_rstput( 0, 0, inum, 'risfdep', risfdep ) ! now k-velocity 1000 CALL iom_rstput( 0, 0, inum, 'ht' , ht ) ! now water column height 1001 ! 1002 IF ( ln_isf ) THEN 1003 IF (ln_isfcav_mlt) THEN 1004 CALL iom_rstput( 0, 0, inum, 'fwfisf_cav', fwfisf_cav ) ! now k-velocity 1005 CALL iom_rstput( 0, 0, inum, 'rhisf_cav_tbl', rhisf_tbl_cav ) ! now k-velocity 1006 CALL iom_rstput( 0, 0, inum, 'rfrac_cav_tbl', rfrac_tbl_cav ) ! now k-velocity 1007 CALL iom_rstput( 0, 0, inum, 'misfkb_cav', REAL(misfkb_cav,wp) ) ! now k-velocity 1008 CALL iom_rstput( 0, 0, inum, 'misfkt_cav', REAL(misfkt_cav,wp) ) ! now k-velocity 1009 CALL iom_rstput( 0, 0, inum, 'mskisf_cav', REAL(mskisf_cav,wp), ktype = jp_i1 ) 1010 END IF 1011 IF (ln_isfpar_mlt) THEN 1012 CALL iom_rstput( 0, 0, inum, 'isfmsk_par', REAL(mskisf_par,wp) ) ! now k-velocity 1013 CALL iom_rstput( 0, 0, inum, 'fwfisf_par', fwfisf_par ) ! now k-velocity 1014 CALL iom_rstput( 0, 0, inum, 'rhisf_par_tbl', rhisf_tbl_par ) ! now k-velocity 1015 CALL iom_rstput( 0, 0, inum, 'rfrac_par_tbl', rfrac_tbl_par ) ! now k-velocity 1016 CALL iom_rstput( 0, 0, inum, 'misfkb_par', REAL(misfkb_par,wp) ) ! now k-velocity 1017 CALL iom_rstput( 0, 0, inum, 'misfkt_par', REAL(misfkt_par,wp) ) ! now k-velocity 1018 CALL iom_rstput( 0, 0, inum, 'mskisf_par', REAL(mskisf_par,wp), ktype = jp_i1 ) 1019 END IF 1020 END IF 1021 ! 996 1022 IF( ALLOCATED(ahtu) ) THEN 997 1023 CALL iom_rstput( 0, 0, inum, 'ahtu', ahtu ) ! aht at u-point … … 1017 1043 CALL iom_rstput( 0, 0, inum, 'sdvecrtz', wsd ) ! now StokesDrift k-velocity 1018 1044 ENDIF 1019 1045 IF ( ln_abl ) THEN 1046 CALL iom_rstput ( 0, 0, inum, "uz1_abl", u_abl(:,:,2,nt_a ) ) ! now first level i-wind 1047 CALL iom_rstput ( 0, 0, inum, "vz1_abl", v_abl(:,:,2,nt_a ) ) ! now first level j-wind 1048 CALL iom_rstput ( 0, 0, inum, "tz1_abl", tq_abl(:,:,2,nt_a,1) ) ! now first level temperature 1049 CALL iom_rstput ( 0, 0, inum, "qz1_abl", tq_abl(:,:,2,nt_a,2) ) ! now first level humidity 1050 ENDIF 1051 ! 1052 CALL iom_close( inum ) 1053 ! 1020 1054 #if defined key_si3 1021 1055 IF( nn_ice == 2 ) THEN ! condition needed in case agrif + ice-model but no-ice in child grid 1056 CALL iom_open( TRIM(cdfile_name)//'_ice', inum, ldwrt = .TRUE., kdlev = jpl, cdcomp = 'ICE' ) 1022 1057 CALL ice_wri_state( inum ) 1058 CALL iom_close( inum ) 1023 1059 ENDIF 1024 1060 #endif 1025 ! 1026 CALL iom_close( inum ) 1027 ! 1061 1028 1062 END SUBROUTINE dia_wri_state 1029 1063 -
NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/tests/CANAL/MY_SRC/domvvl.F90
r12489 r13197 37 37 38 38 PUBLIC dom_vvl_init ! called by domain.F90 39 PUBLIC dom_vvl_zgr ! called by isfcpl.F90 39 40 PUBLIC dom_vvl_sf_nxt ! called by step.F90 40 41 PUBLIC dom_vvl_sf_update ! called by step.F90 … … 62 63 REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:) :: frq_rst_hdv ! retoring period for low freq. divergence 63 64 65 !! * Substitutions 66 # include "do_loop_substitute.h90" 64 67 !!---------------------------------------------------------------------- 65 68 !! NEMO/OCE 4.0 , NEMO Consortium (2018) … … 116 119 INTEGER, INTENT(in) :: Kbb, Kmm, Kaa 117 120 ! 121 IF(lwp) WRITE(numout,*) 122 IF(lwp) WRITE(numout,*) 'dom_vvl_init : Variable volume activated' 123 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~' 124 ! 125 CALL dom_vvl_ctl ! choose vertical coordinate (z_star, z_tilde or layer) 126 ! 127 ! ! Allocate module arrays 128 IF( dom_vvl_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'dom_vvl_init : unable to allocate arrays' ) 129 ! 130 ! ! Read or initialize e3t_(b/n), tilde_e3t_(b/n) and hdiv_lf 131 CALL dom_vvl_rst( nit000, Kbb, Kmm, 'READ' ) 132 e3t(:,:,jpk,Kaa) = e3t_0(:,:,jpk) ! last level always inside the sea floor set one for all 133 ! 134 CALL dom_vvl_zgr(Kbb, Kmm, Kaa) ! interpolation scale factor, depth and water column 135 ! 136 END SUBROUTINE dom_vvl_init 137 ! 138 SUBROUTINE dom_vvl_zgr(Kbb, Kmm, Kaa) 139 !!---------------------------------------------------------------------- 140 !! *** ROUTINE dom_vvl_init *** 141 !! 142 !! ** Purpose : Interpolation of all scale factors, 143 !! depths and water column heights 144 !! 145 !! ** Method : - interpolate scale factors 146 !! 147 !! ** Action : - e3t_(n/b) and tilde_e3t_(n/b) 148 !! - Regrid: e3(u/v)_n 149 !! e3(u/v)_b 150 !! e3w_n 151 !! e3(u/v)w_b 152 !! e3(u/v)w_n 153 !! gdept_n, gdepw_n and gde3w_n 154 !! - h(t/u/v)_0 155 !! - frq_rst_e3t and frq_rst_hdv 156 !! 157 !! Reference : Leclair, M., and G. Madec, 2011, Ocean Modelling. 158 !!---------------------------------------------------------------------- 159 INTEGER, INTENT(in) :: Kbb, Kmm, Kaa 160 !!---------------------------------------------------------------------- 118 161 INTEGER :: ji, jj, jk 119 162 INTEGER :: ii0, ii1, ij0, ij1 120 163 REAL(wp):: zcoef 121 164 !!---------------------------------------------------------------------- 122 !123 IF(lwp) WRITE(numout,*)124 IF(lwp) WRITE(numout,*) 'dom_vvl_init : Variable volume activated'125 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'126 !127 CALL dom_vvl_ctl ! choose vertical coordinate (z_star, z_tilde or layer)128 !129 ! ! Allocate module arrays130 IF( dom_vvl_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'dom_vvl_init : unable to allocate arrays' )131 !132 ! ! Read or initialize e3t_(b/n), tilde_e3t_(b/n) and hdiv_lf133 CALL dom_vvl_rst( nit000, Kbb, Kmm, 'READ' )134 e3t(:,:,jpk,Kaa) = e3t_0(:,:,jpk) ! last level always inside the sea floor set one for all135 165 ! 136 166 ! !== Set of all other vertical scale factors ==! (now and before) … … 160 190 gdept(:,:,1,Kbb) = 0.5_wp * e3w(:,:,1,Kbb) 161 191 gdepw(:,:,1,Kbb) = 0.0_wp 162 DO jk = 2, jpk ! vertical sum 163 DO jj = 1,jpj 164 DO ji = 1,jpi 165 ! zcoef = tmask - wmask ! 0 everywhere tmask = wmask, ie everywhere expect at jk = mikt 166 ! ! 1 everywhere from mbkt to mikt + 1 or 1 (if no isf) 167 ! ! 0.5 where jk = mikt 192 DO_3D_11_11( 2, jpk ) 193 ! zcoef = tmask - wmask ! 0 everywhere tmask = wmask, ie everywhere expect at jk = mikt 194 ! ! 1 everywhere from mbkt to mikt + 1 or 1 (if no isf) 195 ! ! 0.5 where jk = mikt 168 196 !!gm ??????? BUG ? gdept(:,:,:,Kmm) as well as gde3w does not include the thickness of ISF ?? 169 zcoef = ( tmask(ji,jj,jk) - wmask(ji,jj,jk) ) 170 gdepw(ji,jj,jk,Kmm) = gdepw(ji,jj,jk-1,Kmm) + e3t(ji,jj,jk-1,Kmm) 171 gdept(ji,jj,jk,Kmm) = zcoef * ( gdepw(ji,jj,jk ,Kmm) + 0.5 * e3w(ji,jj,jk,Kmm)) & 172 & + (1-zcoef) * ( gdept(ji,jj,jk-1,Kmm) + e3w(ji,jj,jk,Kmm)) 173 gde3w(ji,jj,jk) = gdept(ji,jj,jk,Kmm) - ssh(ji,jj,Kmm) 174 gdepw(ji,jj,jk,Kbb) = gdepw(ji,jj,jk-1,Kbb) + e3t(ji,jj,jk-1,Kbb) 175 gdept(ji,jj,jk,Kbb) = zcoef * ( gdepw(ji,jj,jk ,Kbb) + 0.5 * e3w(ji,jj,jk,Kbb)) & 176 & + (1-zcoef) * ( gdept(ji,jj,jk-1,Kbb) + e3w(ji,jj,jk,Kbb)) 177 END DO 178 END DO 179 END DO 197 zcoef = ( tmask(ji,jj,jk) - wmask(ji,jj,jk) ) 198 gdepw(ji,jj,jk,Kmm) = gdepw(ji,jj,jk-1,Kmm) + e3t(ji,jj,jk-1,Kmm) 199 gdept(ji,jj,jk,Kmm) = zcoef * ( gdepw(ji,jj,jk ,Kmm) + 0.5 * e3w(ji,jj,jk,Kmm)) & 200 & + (1-zcoef) * ( gdept(ji,jj,jk-1,Kmm) + e3w(ji,jj,jk,Kmm)) 201 gde3w(ji,jj,jk) = gdept(ji,jj,jk,Kmm) - ssh(ji,jj,Kmm) 202 gdepw(ji,jj,jk,Kbb) = gdepw(ji,jj,jk-1,Kbb) + e3t(ji,jj,jk-1,Kbb) 203 gdept(ji,jj,jk,Kbb) = zcoef * ( gdepw(ji,jj,jk ,Kbb) + 0.5 * e3w(ji,jj,jk,Kbb)) & 204 & + (1-zcoef) * ( gdept(ji,jj,jk-1,Kbb) + e3w(ji,jj,jk,Kbb)) 205 END_3D 180 206 ! 181 207 ! !== thickness of the water column !! (ocean portion only) … … 212 238 ENDIF 213 239 IF ( ln_vvl_zstar_at_eqtor ) THEN ! use z-star in vicinity of the Equator 214 DO jj = 1, jpj 215 DO ji = 1, jpi 240 DO_2D_11_11 216 241 !!gm case |gphi| >= 6 degrees is useless initialized just above by default 217 IF( ABS(gphit(ji,jj)) >= 6.) THEN 218 ! values outside the equatorial band and transition zone (ztilde) 219 frq_rst_e3t(ji,jj) = 2.0_wp * rpi / ( MAX( rn_rst_e3t , rsmall ) * 86400.e0_wp ) 220 frq_rst_hdv(ji,jj) = 2.0_wp * rpi / ( MAX( rn_lf_cutoff, rsmall ) * 86400.e0_wp ) 221 ELSEIF( ABS(gphit(ji,jj)) <= 2.5) THEN ! Equator strip ==> z-star 222 ! values inside the equatorial band (ztilde as zstar) 223 frq_rst_e3t(ji,jj) = 0.0_wp 224 frq_rst_hdv(ji,jj) = 1.0_wp / rn_Dt 225 ELSE ! transition band (2.5 to 6 degrees N/S) 226 ! ! (linearly transition from z-tilde to z-star) 227 frq_rst_e3t(ji,jj) = 0.0_wp + (frq_rst_e3t(ji,jj)-0.0_wp)*0.5_wp & 228 & * ( 1.0_wp - COS( rad*(ABS(gphit(ji,jj))-2.5_wp) & 229 & * 180._wp / 3.5_wp ) ) 230 frq_rst_hdv(ji,jj) = (1.0_wp / rn_Dt) & 231 & + ( frq_rst_hdv(ji,jj)-(1.e0_wp / rn_Dt) )*0.5_wp & 232 & * ( 1._wp - COS( rad*(ABS(gphit(ji,jj))-2.5_wp) & 233 & * 180._wp / 3.5_wp ) ) 234 ENDIF 235 END DO 236 END DO 242 IF( ABS(gphit(ji,jj)) >= 6.) THEN 243 ! values outside the equatorial band and transition zone (ztilde) 244 frq_rst_e3t(ji,jj) = 2.0_wp * rpi / ( MAX( rn_rst_e3t , rsmall ) * 86400.e0_wp ) 245 frq_rst_hdv(ji,jj) = 2.0_wp * rpi / ( MAX( rn_lf_cutoff, rsmall ) * 86400.e0_wp ) 246 ELSEIF( ABS(gphit(ji,jj)) <= 2.5) THEN ! Equator strip ==> z-star 247 ! values inside the equatorial band (ztilde as zstar) 248 frq_rst_e3t(ji,jj) = 0.0_wp 249 frq_rst_hdv(ji,jj) = 1.0_wp / rn_Dt 250 ELSE ! transition band (2.5 to 6 degrees N/S) 251 ! ! (linearly transition from z-tilde to z-star) 252 frq_rst_e3t(ji,jj) = 0.0_wp + (frq_rst_e3t(ji,jj)-0.0_wp)*0.5_wp & 253 & * ( 1.0_wp - COS( rad*(ABS(gphit(ji,jj))-2.5_wp) & 254 & * 180._wp / 3.5_wp ) ) 255 frq_rst_hdv(ji,jj) = (1.0_wp / rn_Dt) & 256 & + ( frq_rst_hdv(ji,jj)-(1.e0_wp / rn_Dt) )*0.5_wp & 257 & * ( 1._wp - COS( rad*(ABS(gphit(ji,jj))-2.5_wp) & 258 & * 180._wp / 3.5_wp ) ) 259 ENDIF 260 END_2D 237 261 IF( cn_cfg == "orca" .OR. cn_cfg == "ORCA" ) THEN 238 262 IF( nn_cfg == 3 ) THEN ! ORCA2: Suppress ztilde in the Foxe Basin for ORCA2 … … 264 288 ENDIF 265 289 ! 266 END SUBROUTINE dom_vvl_ init290 END SUBROUTINE dom_vvl_zgr 267 291 268 292 … … 329 353 END DO 330 354 ! 331 IF( ln_vvl_ztilde .OR. ln_vvl_layer.AND. ll_do_bclinic ) THEN ! z_tilde or layer coordinate !332 ! ! ------baroclinic part------ !355 IF( (ln_vvl_ztilde .OR. ln_vvl_layer) .AND. ll_do_bclinic ) THEN ! z_tilde or layer coordinate ! 356 ! ! ------baroclinic part------ ! 333 357 ! I - initialization 334 358 ! ================== … … 383 407 zwu(:,:) = 0._wp 384 408 zwv(:,:) = 0._wp 385 DO jk = 1, jpkm1 ! a - first derivative: diffusive fluxes 386 DO jj = 1, jpjm1 387 DO ji = 1, fs_jpim1 ! vector opt. 388 un_td(ji,jj,jk) = rn_ahe3 * umask(ji,jj,jk) * e2_e1u(ji,jj) & 389 & * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji+1,jj ,jk) ) 390 vn_td(ji,jj,jk) = rn_ahe3 * vmask(ji,jj,jk) * e1_e2v(ji,jj) & 391 & * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji ,jj+1,jk) ) 392 zwu(ji,jj) = zwu(ji,jj) + un_td(ji,jj,jk) 393 zwv(ji,jj) = zwv(ji,jj) + vn_td(ji,jj,jk) 394 END DO 395 END DO 396 END DO 397 DO jj = 1, jpj ! b - correction for last oceanic u-v points 398 DO ji = 1, jpi 399 un_td(ji,jj,mbku(ji,jj)) = un_td(ji,jj,mbku(ji,jj)) - zwu(ji,jj) 400 vn_td(ji,jj,mbkv(ji,jj)) = vn_td(ji,jj,mbkv(ji,jj)) - zwv(ji,jj) 401 END DO 402 END DO 403 DO jk = 1, jpkm1 ! c - second derivative: divergence of diffusive fluxes 404 DO jj = 2, jpjm1 405 DO ji = fs_2, fs_jpim1 ! vector opt. 406 tilde_e3t_a(ji,jj,jk) = tilde_e3t_a(ji,jj,jk) + ( un_td(ji-1,jj ,jk) - un_td(ji,jj,jk) & 407 & + vn_td(ji ,jj-1,jk) - vn_td(ji,jj,jk) & 408 & ) * r1_e1e2t(ji,jj) 409 END DO 410 END DO 411 END DO 409 DO_3D_10_10( 1, jpkm1 ) 410 un_td(ji,jj,jk) = rn_ahe3 * umask(ji,jj,jk) * e2_e1u(ji,jj) & 411 & * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji+1,jj ,jk) ) 412 vn_td(ji,jj,jk) = rn_ahe3 * vmask(ji,jj,jk) * e1_e2v(ji,jj) & 413 & * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji ,jj+1,jk) ) 414 zwu(ji,jj) = zwu(ji,jj) + un_td(ji,jj,jk) 415 zwv(ji,jj) = zwv(ji,jj) + vn_td(ji,jj,jk) 416 END_3D 417 DO_2D_11_11 418 un_td(ji,jj,mbku(ji,jj)) = un_td(ji,jj,mbku(ji,jj)) - zwu(ji,jj) 419 vn_td(ji,jj,mbkv(ji,jj)) = vn_td(ji,jj,mbkv(ji,jj)) - zwv(ji,jj) 420 END_2D 421 DO_3D_00_00( 1, jpkm1 ) 422 tilde_e3t_a(ji,jj,jk) = tilde_e3t_a(ji,jj,jk) + ( un_td(ji-1,jj ,jk) - un_td(ji,jj,jk) & 423 & + vn_td(ji ,jj-1,jk) - vn_td(ji,jj,jk) & 424 & ) * r1_e1e2t(ji,jj) 425 END_3D 412 426 ! ! d - thickness diffusion transport: boundary conditions 413 427 ! (stored for tracer advction and continuity equation) … … 416 430 ! 4 - Time stepping of baroclinic scale factors 417 431 ! --------------------------------------------- 418 ! Leapfrog time stepping419 ! ~~~~~~~~~~~~~~~~~~~~~~420 432 CALL lbc_lnk( 'domvvl', tilde_e3t_a(:,:,:), 'T', 1._wp ) 421 433 tilde_e3t_a(:,:,:) = tilde_e3t_b(:,:,:) + rDt * tmask(:,:,:) * tilde_e3t_a(:,:,:) … … 613 625 tilde_e3t_n(:,:,:) = tilde_e3t_a(:,:,:) 614 626 ENDIF 615 gdept(:,:,:,Kbb) = gdept(:,:,:,Kmm)616 gdepw(:,:,:,Kbb) = gdepw(:,:,:,Kmm)617 618 e3t(:,:,:,Kmm) = e3t(:,:,:,Kaa)619 e3u(:,:,:,Kmm) = e3u(:,:,:,Kaa)620 e3v(:,:,:,Kmm) = e3v(:,:,:,Kaa)621 627 622 628 ! Compute all missing vertical scale factor and depths … … 641 647 gdepw(:,:,1,Kmm) = 0.0_wp 642 648 gde3w(:,:,1) = gdept(:,:,1,Kmm) - ssh(:,:,Kmm) 643 DO jk = 2, jpk 644 DO jj = 1,jpj 645 DO ji = 1,jpi 646 ! zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk)) ! 0 everywhere tmask = wmask, ie everywhere expect at jk = mikt 647 ! 1 for jk = mikt 648 zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk)) 649 gdepw(ji,jj,jk,Kmm) = gdepw(ji,jj,jk-1,Kmm) + e3t(ji,jj,jk-1,Kmm) 650 gdept(ji,jj,jk,Kmm) = zcoef * ( gdepw(ji,jj,jk ,Kmm) + 0.5 * e3w(ji,jj,jk,Kmm) ) & 651 & + (1-zcoef) * ( gdept(ji,jj,jk-1,Kmm) + e3w(ji,jj,jk,Kmm) ) 652 gde3w(ji,jj,jk) = gdept(ji,jj,jk,Kmm) - ssh(ji,jj,Kmm) 653 END DO 654 END DO 655 END DO 649 DO_3D_11_11( 2, jpk ) 650 ! zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk)) ! 0 everywhere tmask = wmask, ie everywhere expect at jk = mikt 651 ! 1 for jk = mikt 652 zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk)) 653 gdepw(ji,jj,jk,Kmm) = gdepw(ji,jj,jk-1,Kmm) + e3t(ji,jj,jk-1,Kmm) 654 gdept(ji,jj,jk,Kmm) = zcoef * ( gdepw(ji,jj,jk ,Kmm) + 0.5 * e3w(ji,jj,jk,Kmm) ) & 655 & + (1-zcoef) * ( gdept(ji,jj,jk-1,Kmm) + e3w(ji,jj,jk,Kmm) ) 656 gde3w(ji,jj,jk) = gdept(ji,jj,jk,Kmm) - ssh(ji,jj,Kmm) 657 END_3D 656 658 657 659 ! Local depth and Inverse of the local depth of the water … … 700 702 ! 701 703 CASE( 'U' ) !* from T- to U-point : hor. surface weighted mean 702 DO jk = 1, jpk 703 DO jj = 1, jpjm1 704 DO ji = 1, fs_jpim1 ! vector opt. 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 DO 709 END DO 710 END DO 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 711 709 CALL lbc_lnk( 'domvvl', pe3_out(:,:,:), 'U', 1._wp ) 712 710 pe3_out(:,:,:) = pe3_out(:,:,:) + e3u_0(:,:,:) 713 711 ! 714 712 CASE( 'V' ) !* from T- to V-point : hor. surface weighted mean 715 DO jk = 1, jpk 716 DO jj = 1, jpjm1 717 DO ji = 1, fs_jpim1 ! vector opt. 718 pe3_out(ji,jj,jk) = 0.5_wp * ( vmask(ji,jj,jk) * (1.0_wp - zlnwd) + zlnwd ) * r1_e1e2v(ji,jj) & 719 & * ( e1e2t(ji,jj ) * ( pe3_in(ji,jj ,jk) - e3t_0(ji,jj ,jk) ) & 720 & + e1e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) ) 721 END DO 722 END DO 723 END DO 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 724 718 CALL lbc_lnk( 'domvvl', pe3_out(:,:,:), 'V', 1._wp ) 725 719 pe3_out(:,:,:) = pe3_out(:,:,:) + e3v_0(:,:,:) 726 720 ! 727 721 CASE( 'F' ) !* from U-point to F-point : hor. surface weighted mean 728 DO jk = 1, jpk 729 DO jj = 1, jpjm1 730 DO ji = 1, fs_jpim1 ! vector opt. 731 pe3_out(ji,jj,jk) = 0.5_wp * ( umask(ji,jj,jk) * umask(ji,jj+1,jk) * (1.0_wp - zlnwd) + zlnwd ) & 732 & * r1_e1e2f(ji,jj) & 733 & * ( e1e2u(ji,jj ) * ( pe3_in(ji,jj ,jk) - e3u_0(ji,jj ,jk) ) & 734 & + e1e2u(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3u_0(ji,jj+1,jk) ) ) 735 END DO 736 END DO 737 END DO 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 738 728 CALL lbc_lnk( 'domvvl', pe3_out(:,:,:), 'F', 1._wp ) 739 729 pe3_out(:,:,:) = pe3_out(:,:,:) + e3f_0(:,:,:) … … 810 800 id4 = iom_varid( numror, 'tilde_e3t_n', ldstop = .FALSE. ) 811 801 id5 = iom_varid( numror, 'hdiv_lf', ldstop = .FALSE. ) 802 ! 812 803 ! ! --------- ! 813 804 ! ! all cases ! 814 805 ! ! --------- ! 806 ! 815 807 IF( MIN( id1, id2 ) > 0 ) THEN ! all required arrays exist 816 808 CALL iom_get( numror, jpdom_autoglo, 'e3t_b', e3t(:,:,:,Kbb), ldxios = lrxios ) … … 828 820 IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t(:,:,:,Kmm) not found in restart files' 829 821 IF(lwp) write(numout,*) 'e3t_n set equal to e3t_b.' 830 IF(lwp) write(numout,*) 'l_1st_euler is forced to .true.'822 IF(lwp) write(numout,*) 'l_1st_euler is forced to true' 831 823 CALL iom_get( numror, jpdom_autoglo, 'e3t_b', e3t(:,:,:,Kbb), ldxios = lrxios ) 832 824 e3t(:,:,:,Kmm) = e3t(:,:,:,Kbb) … … 835 827 IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t(:,:,:,Kbb) not found in restart files' 836 828 IF(lwp) write(numout,*) 'e3t_b set equal to e3t_n.' 837 IF(lwp) write(numout,*) 'l_1st_euler is forced to .true.'829 IF(lwp) write(numout,*) 'l_1st_euler is forced to true' 838 830 CALL iom_get( numror, jpdom_autoglo, 'e3t_n', e3t(:,:,:,Kmm), ldxios = lrxios ) 839 831 e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm) … … 842 834 IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t(:,:,:,Kmm) not found in restart file' 843 835 IF(lwp) write(numout,*) 'Compute scale factor from sshn' 844 IF(lwp) write(numout,*) 'l_1st_euler is forced to .true.'836 IF(lwp) write(numout,*) 'l_1st_euler is forced to true' 845 837 DO jk = 1, jpk 846 838 e3t(:,:,jk,Kmm) = e3t_0(:,:,jk) * ( ht_0(:,:) + ssh(:,:,Kmm) ) & … … 895 887 ssh(:,:,Kbb) = -ssh_ref 896 888 897 DO jj = 1, jpj 898 DO ji = 1, jpi 899 IF( ht_0(ji,jj)-ssh_ref < rn_wdmin1 ) THEN ! if total depth is less than min depth 900 ssh(ji,jj,Kbb) = rn_wdmin1 - (ht_0(ji,jj) ) 901 ssh(ji,jj,Kmm) = rn_wdmin1 - (ht_0(ji,jj) ) 902 ENDIF 903 ENDDO 904 ENDDO 889 DO_2D_11_11 890 IF( ht_0(ji,jj)-ssh_ref < rn_wdmin1 ) THEN ! if total depth is less than min depth 891 ssh(ji,jj,Kbb) = rn_wdmin1 - (ht_0(ji,jj) ) 892 ssh(ji,jj,Kmm) = rn_wdmin1 - (ht_0(ji,jj) ) 893 ENDIF 894 END_2D 905 895 ENDIF !If test case else 906 896 … … 913 903 e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm) 914 904 915 DO ji = 1, jpi 916 DO jj = 1, jpj 917 IF ( ht_0(ji,jj) .LE. 0.0 .AND. NINT( ssmask(ji,jj) ) .EQ. 1) THEN 918 CALL ctl_stop( 'dom_vvl_rst: ht_0 must be positive at potentially wet points' ) 919 ENDIF 920 END DO 921 END DO 905 DO_2D_11_11 906 IF ( ht_0(ji,jj) .LE. 0.0 .AND. NINT( ssmask(ji,jj) ) .EQ. 1) THEN 907 CALL ctl_stop( 'dom_vvl_rst: ht_0 must be positive at potentially wet points' ) 908 ENDIF 909 END_2D 922 910 ! 923 911 ELSE 924 912 ! 925 ! usr_def_istate called here only to get sshb, that is needed to initialize e3t(Kbb) and e3t(Kmm) 926 CALL usr_def_istate( gdept_0, tmask, ts(:,:,:,:,Kbb), uu(:,:,:,Kbb), vv(:,:,:,Kbb), ssh(:,:,Kbb) ) 927 ! usr_def_istate will be called again in istate_init to initialize ts(bn), ssh(bn), u(bn) and v(bn) 913 ! usr_def_istate called here only to get ssh(Kbb) needed to initialize e3t(Kbb) and e3t(Kmm) 914 ! 915 CALL usr_def_istate( gdept_0, tmask, ts(:,:,:,:,Kbb), uu(:,:,:,Kbb), vv(:,:,:,Kbb), ssh(:,:,Kbb) ) 916 ! 917 ! usr_def_istate will be called again in istate_init to initialize ts, ssh, u and v 928 918 ! 929 919 DO jk=1,jpk 930 e3t(:,:,jk,K mm) = e3t_0(:,:,jk) * ( ht_0(:,:) + ssh(:,:,Kbb) ) &931 & 932 & + e3t_0(:,:,jk) * ( 1._wp - tmask(:,:,jk) ) ! make sure e3t_b!= 0 on land points920 e3t(:,:,jk,Kbb) = e3t_0(:,:,jk) * ( ht_0(:,:) + ssh(:,:,Kbb) ) & 921 & / ( ht_0(:,:) + 1._wp - ssmask(:,:) ) * tmask(:,:,jk) & 922 & + e3t_0(:,:,jk) * ( 1._wp - tmask(:,:,jk) ) ! make sure e3t(:,:,:,Kbb) != 0 on land points 933 923 END DO 934 924 e3t(:,:,:,Kmm) = e3t(:,:,:,Kbb) 935 ssh(:,: ,Kmm) = ssh(:,: ,Kbb)! needed later for gde3w925 ssh(:,:,Kmm) = ssh(:,:,Kbb) ! needed later for gde3w 936 926 ! 937 927 END IF ! end of ll_wd edits … … 1025 1015 ! 1026 1016 IF( ioptio /= 1 ) CALL ctl_stop( 'Choose ONE vertical coordinate in namelist nam_vvl' ) 1027 IF( .NOT. ln_vvl_zstar .AND. ln_isf ) CALL ctl_stop( 'Only vvl_zstar has been tested with ice shelf cavity' )1028 1017 ! 1029 1018 IF(lwp) THEN ! Print the choice -
NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/tests/CANAL/MY_SRC/stpctl.F90
r13159 r13197 19 19 USE dom_oce ! ocean space and time domain variables 20 20 USE c1d ! 1D vertical configuration 21 USE zdf_oce , ONLY : ln_zad_Aimp ! ocean vertical physics variables22 USE wet_dry, ONLY : ll_wd, ssh_ref ! reference depth for negative bathy23 !24 21 USE diawri ! Standard run outputs (dia_wri_state routine) 22 ! 25 23 USE in_out_manager ! I/O manager 26 24 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 27 25 USE lib_mpp ! distributed memory computing 28 ! 26 USE zdf_oce , ONLY : ln_zad_Aimp ! ocean vertical physics variables 27 USE wet_dry, ONLY : ll_wd, ssh_ref ! reference depth for negative bathy 28 29 29 USE netcdf ! NetCDF library 30 30 IMPLICIT NONE … … 33 33 PUBLIC stp_ctl ! routine called by step.F90 34 34 35 INTEGER :: nrunid ! netcdf file id36 INTEGER, DIMENSION(8) :: nvarid ! netcdf variable id35 INTEGER :: idrun, idtime, idssh, idu, ids1, ids2, idt1, idt2, idc1, idw1, istatus 36 LOGICAL :: lsomeoce 37 37 !!---------------------------------------------------------------------- 38 38 !! NEMO/OCE 4.0 , NEMO Consortium (2018) … … 42 42 CONTAINS 43 43 44 SUBROUTINE stp_ctl( kt, K mm)44 SUBROUTINE stp_ctl( kt, Kbb, Kmm, kindic ) 45 45 !!---------------------------------------------------------------------- 46 46 !! *** ROUTINE stp_ctl *** … … 50 50 !! ** Method : - Save the time step in numstp 51 51 !! - Print it each 50 time steps 52 !! - Stop the run IF problem encountered by setting nstop > 052 !! - Stop the run IF problem encountered by setting indic=-3 53 53 !! Problems checked: |ssh| maximum larger than 10 m 54 54 !! |U| maximum larger than 10 m/s … … 57 57 !! ** Actions : "time.step" file = last ocean time-step 58 58 !! "run.stat" file = run statistics 59 !! nstop indicator sheared among all local domain59 !! nstop indicator sheared among all local domain (lk_mpp=T) 60 60 !!---------------------------------------------------------------------- 61 61 INTEGER, INTENT(in ) :: kt ! ocean time-step index 62 INTEGER, INTENT(in ) :: Kmm ! ocean time level index 62 INTEGER, INTENT(in ) :: Kbb, Kmm ! ocean time level index 63 INTEGER, INTENT(inout) :: kindic ! error indicator 63 64 !! 64 INTEGER :: ji ! dummy loop indices 65 INTEGER :: idtime, istatus 66 INTEGER , DIMENSION(9) :: iareasum, iareamin, iareamax 67 INTEGER , DIMENSION(3,4) :: iloc ! min/max loc indices 68 REAL(wp) :: zzz ! local real 69 REAL(wp), DIMENSION(9) :: zmax, zmaxlocal 70 LOGICAL :: ll_wrtstp, ll_colruns, ll_wrtruns 71 LOGICAL, DIMENSION(jpi,jpj,jpk) :: llmsk 72 CHARACTER(len=20) :: clname 65 INTEGER :: ji, jj, jk ! dummy loop indices 66 INTEGER, DIMENSION(2) :: ih ! min/max loc indices 67 INTEGER, DIMENSION(3) :: iu, is1, is2 ! min/max loc indices 68 REAL(wp) :: zzz ! local real 69 REAL(wp), DIMENSION(9) :: zmax 70 LOGICAL :: ll_wrtstp, ll_colruns, ll_wrtruns 71 CHARACTER(len=20) :: clname 73 72 !!---------------------------------------------------------------------- 74 IF( nstop > 0 .AND. ngrdstop > -1 ) RETURN ! stpctl was already called by a child grid 75 ! 76 ll_wrtstp = ( MOD( kt-nit000, sn_cfctl%ptimincr ) == 0 ) .OR. ( kt == nitend ) 77 ll_colruns = ll_wrtstp .AND. sn_cfctl%l_runstat .AND. jpnij > 1 78 ll_wrtruns = ( ll_colruns .OR. jpnij == 1 ) .AND. lwm 79 ! 80 IF( kt == nit000 ) THEN 81 ! 82 IF( lwp ) THEN 83 WRITE(numout,*) 84 WRITE(numout,*) 'stp_ctl : time-stepping control' 85 WRITE(numout,*) '~~~~~~~' 86 ENDIF 87 ! ! open time.step ascii file, done only by 1st subdomain 88 IF( lwm ) CALL ctl_opn( numstp, 'time.step', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp, narea ) 89 ! 90 IF( ll_wrtruns ) THEN 91 ! ! open run.stat ascii file, done only by 1st subdomain 73 ! 74 ll_wrtstp = ( MOD( kt, sn_cfctl%ptimincr ) == 0 ) .OR. ( kt == nitend ) 75 ll_colruns = ll_wrtstp .AND. ( sn_cfctl%l_runstat ) 76 ll_wrtruns = ll_colruns .AND. lwm 77 IF( kt == nit000 .AND. lwp ) THEN 78 WRITE(numout,*) 79 WRITE(numout,*) 'stp_ctl : time-stepping control' 80 WRITE(numout,*) '~~~~~~~' 81 ! ! open time.step file 82 IF( lwm ) CALL ctl_opn( numstp, 'time.step', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp, narea ) 83 ! ! open run.stat file(s) at start whatever 84 ! ! the value of sn_cfctl%ptimincr 85 IF( lwm .AND. ( sn_cfctl%l_runstat ) ) THEN 92 86 CALL ctl_opn( numrun, 'run.stat', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp, narea ) 93 ! ! open run.stat.nc netcdf file, done only by 1st subdomain94 87 clname = 'run.stat.nc' 95 88 IF( .NOT. Agrif_Root() ) clname = TRIM(Agrif_CFixed())//"_"//TRIM(clname) 96 istatus = NF90_CREATE( TRIM(clname), NF90_CLOBBER, nrunid)97 istatus = NF90_DEF_DIM( nrunid, 'time', NF90_UNLIMITED, idtime )98 istatus = NF90_DEF_VAR( nrunid, 'abs_ssh_max', NF90_DOUBLE, (/ idtime /), nvarid(1))99 istatus = NF90_DEF_VAR( nrunid, 'abs_u_max', NF90_DOUBLE, (/ idtime /), nvarid(2))100 istatus = NF90_DEF_VAR( nrunid, 's_min', NF90_DOUBLE, (/ idtime /), nvarid(3))101 istatus = NF90_DEF_VAR( nrunid, 's_max', NF90_DOUBLE, (/ idtime /), nvarid(4))102 istatus = NF90_DEF_VAR( nrunid, 't_min', NF90_DOUBLE, (/ idtime /), nvarid(5))103 istatus = NF90_DEF_VAR( nrunid, 't_max', NF90_DOUBLE, (/ idtime /), nvarid(6))89 istatus = NF90_CREATE( TRIM(clname), NF90_CLOBBER, idrun ) 90 istatus = NF90_DEF_DIM( idrun, 'time', NF90_UNLIMITED, idtime ) 91 istatus = NF90_DEF_VAR( idrun, 'abs_ssh_max', NF90_DOUBLE, (/ idtime /), idssh ) 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 ) 104 97 IF( ln_zad_Aimp ) THEN 105 istatus = NF90_DEF_VAR( nrunid, 'Cf_max', NF90_DOUBLE, (/ idtime /), nvarid(7))106 istatus = NF90_DEF_VAR( nrunid,'abs_wi_max',NF90_DOUBLE, (/ idtime /), nvarid(8))98 istatus = NF90_DEF_VAR( idrun, 'abs_wi_max', NF90_DOUBLE, (/ idtime /), idw1 ) 99 istatus = NF90_DEF_VAR( idrun, 'Cf_max', NF90_DOUBLE, (/ idtime /), idc1 ) 107 100 ENDIF 108 istatus = NF90_ENDDEF(nrunid) 109 ENDIF 110 ! 111 ENDIF 112 ! 113 ! !== write current time step ==! 114 ! !== done only by 1st subdomain at writting timestep ==! 115 IF( lwm .AND. ll_wrtstp ) THEN 101 istatus = NF90_ENDDEF(idrun) 102 zmax(8:9) = 0._wp ! initialise to zero in case ln_zad_Aimp option is not in use 103 ENDIF 104 ENDIF 105 IF( kt == nit000 ) lsomeoce = COUNT( ssmask(:,:) == 1._wp ) > 0 106 ! 107 IF(lwm .AND. ll_wrtstp) THEN !== current time step ==! ("time.step" file) 116 108 WRITE ( numstp, '(1x, i8)' ) kt 117 109 REWIND( numstp ) 118 110 ENDIF 119 ! !== test of local extrema ==! 120 ! !== done by all processes at every time step ==! 121 ! 122 ! define zmax default value. needed for land processors 123 IF( ll_colruns ) THEN ! default value: must not be kept when calling mpp_max -> must be as small as possible 124 zmax(:) = -HUGE(1._wp) 125 ELSE ! default value: must not give true for any of the tests bellow (-> avoid manipulating HUGE...) 126 zmax(:) = 0._wp 127 zmax(3) = -1._wp ! avoid salinity minimum at 0. 128 ENDIF 129 ! 130 llmsk(:,:,1) = ssmask(:,:) == 1._wp 131 IF( COUNT( llmsk(:,:,1) ) > 0 ) THEN ! avoid huge values sent back for land processors... 132 IF( ll_wd ) THEN 133 zmax(1) = MAXVAL( ABS( ssh(:,:,Kmm) + ssh_ref ), mask = llmsk(:,:,1) ) ! ssh max 134 ELSE 135 zmax(1) = MAXVAL( ABS( ssh(:,:,Kmm) ), mask = llmsk(:,:,1) ) ! ssh max 136 ENDIF 137 ENDIF 138 zmax(2) = MAXVAL( ABS( uu(:,:,:,Kmm) ) ) ! velocity max (zonal only) 139 llmsk(:,:,:) = tmask(:,:,:) == 1._wp 140 IF( COUNT( llmsk(:,:,:) ) > 0 ) THEN ! avoid huge values sent back for land processors... 141 zmax(3) = MAXVAL( -ts(:,:,:,jp_sal,Kmm), mask = llmsk ) ! minus salinity max 142 zmax(4) = MAXVAL( ts(:,:,:,jp_sal,Kmm), mask = llmsk ) ! salinity max 143 IF( ll_colruns .OR. jpnij == 1 ) THEN ! following variables are used only in the netcdf file 144 zmax(5) = MAXVAL( -ts(:,:,:,jp_tem,Kmm), mask = llmsk ) ! minus temperature max 145 zmax(6) = MAXVAL( ts(:,:,:,jp_tem,Kmm), mask = llmsk ) ! temperature max 146 IF( ln_zad_Aimp ) THEN 147 zmax(7) = MAXVAL( Cu_adv(:,:,:) , mask = llmsk ) ! partitioning coeff. max 148 llmsk(:,:,:) = wmask(:,:,:) == 1._wp 149 IF( COUNT( llmsk(:,:,:) ) > 0 ) THEN ! avoid huge values sent back for land processors... 150 zmax(8) = MAXVAL(ABS( wi(:,:,:) ), mask = llmsk ) ! implicit vertical vel. max 151 ENDIF 152 ENDIF 153 ENDIF 154 ENDIF 155 zmax(9) = REAL( nstop, wp ) ! stop indicator 156 ! !== get global extrema ==! 157 ! !== done by all processes if writting run.stat ==! 111 ! 112 ! !== test of extrema ==! 113 IF( ll_wd ) THEN 114 zmax(1) = MAXVAL( ABS( ssh(:,:,Kmm) + ssh_ref*tmask(:,:,1) ) ) ! ssh max 115 ELSE 116 zmax(1) = MAXVAL( ABS( ssh(:,:,Kmm) ) ) ! ssh max 117 ENDIF 118 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 128 ! 158 129 IF( ll_colruns ) THEN 159 zmaxlocal(:) = zmax(:)160 130 CALL mpp_max( "stpctl", zmax ) ! max over the global domain 161 nstop = NINT( zmax(9) ) ! update nstop indicator (now sheared among all local domains) 162 ENDIF 163 ! !== write "run.stat" files ==! 164 ! !== done only by 1st subdomain at writting timestep ==! 131 nstop = NINT( zmax(7) ) ! nstop indicator sheared among all local domains 132 ENDIF 133 ! !== run statistics ==! ("run.stat" files) 165 134 IF( ll_wrtruns ) THEN 166 135 WRITE(numrun,9500) kt, zmax(1), zmax(2), -zmax(3), zmax(4) 167 istatus = NF90_PUT_VAR( nrunid, nvarid(1), (/ zmax(1)/), (/kt/), (/1/) )168 istatus = NF90_PUT_VAR( nrunid, nvarid(2), (/ zmax(2)/), (/kt/), (/1/) )169 istatus = NF90_PUT_VAR( nrunid, nvarid(3), (/-zmax(3)/), (/kt/), (/1/) )170 istatus = NF90_PUT_VAR( nrunid, nvarid(4), (/ zmax(4)/), (/kt/), (/1/) )171 istatus = NF90_PUT_VAR( nrunid, nvarid(5), (/-zmax(5)/), (/kt/), (/1/) )172 istatus = NF90_PUT_VAR( nrunid, nvarid(6), (/ zmax(6)/), (/kt/), (/1/) )136 istatus = NF90_PUT_VAR( idrun, idssh, (/ zmax(1)/), (/kt/), (/1/) ) 137 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/) ) 173 142 IF( ln_zad_Aimp ) THEN 174 istatus = NF90_PUT_VAR( nrunid, nvarid(7), (/ zmax(7)/), (/kt/), (/1/) ) 175 istatus = NF90_PUT_VAR( nrunid, nvarid(8), (/ zmax(8)/), (/kt/), (/1/) ) 176 ENDIF 177 IF( kt == nitend ) istatus = NF90_CLOSE(nrunid) 143 istatus = NF90_PUT_VAR( idrun, idw1, (/ zmax(8)/), (/kt/), (/1/) ) 144 istatus = NF90_PUT_VAR( idrun, idc1, (/ zmax(9)/), (/kt/), (/1/) ) 145 ENDIF 146 IF( MOD( kt , 100 ) == 0 ) istatus = NF90_SYNC(idrun) 147 IF( kt == nitend ) istatus = NF90_CLOSE(idrun) 178 148 END IF 179 ! !== error handling ==! 180 ! !== done by all processes at every time step ==! 181 ! 182 IF( zmax(1) > 20._wp .OR. & ! too large sea surface height ( > 20 m ) 183 & zmax(2) > 10._wp .OR. & ! too large velocity ( > 10 m/s) 184 !!$ & zmax(3) >= 0._wp .OR. & ! negative or zero sea surface salinity 185 !!$ & zmax(4) >= 100._wp .OR. & ! too large sea surface salinity ( > 100 ) 186 !!$ & zmax(4) < 0._wp .OR. & ! too large sea surface salinity (keep this line for sea-ice) 187 & ISNAN( zmax(1) + zmax(2) + zmax(3) ) .OR. & ! NaN encounter in the tests 188 & ABS( zmax(1) + zmax(2) + zmax(3) ) > HUGE(1._wp) ) THEN ! Infinity encounter in the tests 149 ! !== error handling ==! 150 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 ) 152 & 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 157 IF( lk_mpp .AND. sn_cfctl%l_glochk ) THEN 158 ! have use mpp_max (because sn_cfctl%l_glochk=.T. and distributed) 159 CALL mpp_maxloc( 'stpctl', ABS(ssh(:,:,Kmm)) , ssmask(:,:) , zzz, ih ) 160 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 ELSE 164 ! find local min and max locations 165 ih(:) = MAXLOC( ABS( ssh(:,:,Kmm) ) ) + (/ nimpp - 1, njmpp - 1 /) 166 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 ENDIF 170 171 WRITE(ctmp1,*) ' stp_ctl: |ssh| > 20 m or |U| > 10 m/s or S <= 0 or S >= 100 or NaN encounter in the tests' 172 WRITE(ctmp2,9100) kt, zmax(1), ih(1) , ih(2) 173 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' 177 178 CALL dia_wri_state( Kmm, 'output.abort' ) ! create an output.abort file 179 180 IF( .NOT. sn_cfctl%l_glochk ) THEN 181 WRITE(ctmp8,*) 'E R R O R message from sub-domain: ', narea 182 CALL ctl_stop( 'STOP', ctmp1, ' ', ctmp8, ' ', ctmp2, ctmp3, ctmp4, ctmp5, ctmp6 ) 183 ELSE 184 CALL ctl_stop( ctmp1, ' ', ctmp2, ctmp3, ctmp4, ctmp5, ' ', ctmp6, ' ' ) 185 ENDIF 186 187 kindic = -3 189 188 ! 190 iloc(:,:) = 0 191 IF( ll_colruns ) THEN ! zmax is global, so it is the same on all subdomains -> no dead lock with mpp_maxloc 192 ! first: close the netcdf file, so we can read it 193 IF( lwm .AND. kt /= nitend ) istatus = NF90_CLOSE(nrunid) 194 ! get global loc on the min/max 195 CALL mpp_maxloc( 'stpctl', ABS(ssh(:,:, Kmm)), ssmask(:,: ), zzz, iloc(1:2,1) ) ! mpp_maxloc ok if mask = F 196 CALL mpp_maxloc( 'stpctl', ABS( uu(:,:,:, Kmm)), umask(:,:,:), zzz, iloc(1:3,2) ) 197 CALL mpp_minloc( 'stpctl', ts(:,:,:,jp_sal,Kmm) , tmask(:,:,:), zzz, iloc(1:3,3) ) 198 CALL mpp_maxloc( 'stpctl', ts(:,:,:,jp_sal,Kmm) , tmask(:,:,:), zzz, iloc(1:3,4) ) 199 ! find which subdomain has the max. 200 iareamin(:) = jpnij+1 ; iareamax(:) = 0 ; iareasum(:) = 0 201 DO ji = 1, 9 202 IF( zmaxlocal(ji) == zmax(ji) ) THEN 203 iareamin(ji) = narea ; iareamax(ji) = narea ; iareasum(ji) = 1 204 ENDIF 205 END DO 206 CALL mpp_min( "stpctl", iareamin ) ! min over the global domain 207 CALL mpp_max( "stpctl", iareamax ) ! max over the global domain 208 CALL mpp_sum( "stpctl", iareasum ) ! sum over the global domain 209 ELSE ! find local min and max locations: 210 ! if we are here, this means that the subdomain contains some oce points -> no need to test the mask used in maxloc 211 iloc(1:2,1) = MAXLOC( ABS( ssh(:,:, Kmm)), mask = ssmask(:,: ) == 1._wp ) + (/ nimpp - 1, njmpp - 1 /) 212 iloc(1:3,2) = MAXLOC( ABS( uu(:,:,:, Kmm)), mask = umask(:,:,:) == 1._wp ) + (/ nimpp - 1, njmpp - 1, 0 /) 213 iloc(1:3,3) = MINLOC( ts(:,:,:,jp_sal,Kmm) , mask = tmask(:,:,:) == 1._wp ) + (/ nimpp - 1, njmpp - 1, 0 /) 214 iloc(1:3,4) = MAXLOC( ts(:,:,:,jp_sal,Kmm) , mask = tmask(:,:,:) == 1._wp ) + (/ nimpp - 1, njmpp - 1, 0 /) 215 iareamin(:) = narea ; iareamax(:) = narea ; iareasum(:) = 1 ! this is local information 216 ENDIF 217 ! 218 WRITE(ctmp1,*) ' stp_ctl: |ssh| > 20 m or |U| > 10 m/s or S <= 0 or S >= 100 or NaN encounter in the tests' 219 CALL wrt_line( ctmp2, kt, '|ssh| max', zmax(1), iloc(:,1), iareasum(1), iareamin(1), iareamax(1) ) 220 CALL wrt_line( ctmp3, kt, '|U| max', zmax(2), iloc(:,2), iareasum(2), iareamin(2), iareamax(2) ) 221 CALL wrt_line( ctmp4, kt, 'Sal min', -zmax(3), iloc(:,3), iareasum(3), iareamin(3), iareamax(3) ) 222 CALL wrt_line( ctmp5, kt, 'Sal max', zmax(4), iloc(:,4), iareasum(4), iareamin(4), iareamax(4) ) 223 IF( Agrif_Root() ) THEN 224 WRITE(ctmp6,*) ' ===> output of last computed fields in output.abort* files' 225 ELSE 226 WRITE(ctmp6,*) ' ===> output of last computed fields in '//TRIM(Agrif_CFixed())//'_output.abort* files' 227 ENDIF 228 ! 229 CALL dia_wri_state( Kmm, 'output.abort' ) ! create an output.abort file 230 ! 231 IF( ll_colruns .or. jpnij == 1 ) THEN ! all processes synchronized -> use lwp to print in opened ocean.output files 232 IF(lwp) THEN ; CALL ctl_stop( ctmp1, ' ', ctmp2, ctmp3, ctmp4, ctmp5, ' ', ctmp6 ) 233 ELSE ; nstop = MAX(1, nstop) ! make sure nstop > 0 (automatically done when calling ctl_stop) 234 ENDIF 235 ELSE ! only mpi subdomains with errors are here -> STOP now 236 CALL ctl_stop( 'STOP', ctmp1, ' ', ctmp2, ctmp3, ctmp4, ctmp5, ' ', ctmp6 ) 237 ENDIF 238 ! 239 ENDIF 240 ! 241 IF( nstop > 0 ) THEN ! an error was detected and we did not abort yet... 242 ngrdstop = Agrif_Fixed() ! store which grid got this error 243 IF( .NOT. ll_colruns .AND. jpnij > 1 ) CALL ctl_stop( 'STOP' ) ! we must abort here to avoid MPI deadlock 244 ENDIF 245 ! 189 ENDIF 190 ! 191 9100 FORMAT (' kt=',i8,' |ssh| max: ',1pg11.4,', at i j : ',2i5) 192 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) 246 195 9500 FORMAT(' it :', i8, ' |ssh|_max: ', D23.16, ' |U|_max: ', D23.16,' S_min: ', D23.16,' S_max: ', D23.16) 247 196 ! 248 197 END SUBROUTINE stp_ctl 249 250 251 SUBROUTINE wrt_line( cdline, kt, cdprefix, pval, kloc, ksum, kmin, kmax )252 !!----------------------------------------------------------------------253 !! *** ROUTINE wrt_line ***254 !!255 !! ** Purpose : write information line256 !!257 !!----------------------------------------------------------------------258 CHARACTER(len=*), INTENT( out) :: cdline259 CHARACTER(len=*), INTENT(in ) :: cdprefix260 REAL(wp), INTENT(in ) :: pval261 INTEGER, DIMENSION(3), INTENT(in ) :: kloc262 INTEGER, INTENT(in ) :: kt, ksum, kmin, kmax263 !264 CHARACTER(len=80) :: clsuff265 CHARACTER(len=9 ) :: clkt, clsum, clmin, clmax266 CHARACTER(len=9 ) :: cli, clj, clk267 CHARACTER(len=1 ) :: clfmt268 CHARACTER(len=4 ) :: cl4 ! needed to be able to compile with Agrif, I don't know why269 INTEGER :: ifmtk270 !!----------------------------------------------------------------------271 WRITE(clkt , '(i9)') kt272 273 WRITE(clfmt, '(i1)') INT(LOG10(REAL(jpnij ,wp))) + 1 ! how many digits to we need to write ? (we decide max = 9)274 !!! WRITE(clsum, '(i'//clfmt//')') ksum ! this is creating a compilation error with AGRIF275 cl4 = '(i'//clfmt//')' ; WRITE(clsum, cl4) ksum276 WRITE(clfmt, '(i1)') INT(LOG10(REAL(MAX(1,jpnij-1),wp))) + 1 ! how many digits to we need to write ? (we decide max = 9)277 cl4 = '(i'//clfmt//')' ; WRITE(clmin, cl4) kmin-1278 WRITE(clmax, cl4) kmax-1279 !280 WRITE(clfmt, '(i1)') INT(LOG10(REAL(jpiglo,wp))) + 1 ! how many digits to we need to write jpiglo? (we decide max = 9)281 cl4 = '(i'//clfmt//')' ; WRITE(cli, cl4) kloc(1) ! this is ok with AGRIF282 WRITE(clfmt, '(i1)') INT(LOG10(REAL(jpjglo,wp))) + 1 ! how many digits to we need to write jpjglo? (we decide max = 9)283 cl4 = '(i'//clfmt//')' ; WRITE(clj, cl4) kloc(2) ! this is ok with AGRIF284 !285 IF( ksum == 1 ) THEN ; WRITE(clsuff,9100) TRIM(clmin)286 ELSE ; WRITE(clsuff,9200) TRIM(clsum), TRIM(clmin), TRIM(clmax)287 ENDIF288 IF(kloc(3) == 0) THEN289 ifmtk = INT(LOG10(REAL(jpk,wp))) + 1 ! how many digits to we need to write jpk? (we decide max = 9)290 clk = REPEAT(' ', ifmtk) ! create the equivalent in blank string291 WRITE(cdline,9300) TRIM(ADJUSTL(clkt)), TRIM(ADJUSTL(cdprefix)), pval, TRIM(cli), TRIM(clj), clk(1:ifmtk), TRIM(clsuff)292 ELSE293 WRITE(clfmt, '(i1)') INT(LOG10(REAL(jpk,wp))) + 1 ! how many digits to we need to write jpk? (we decide max = 9)294 !!! WRITE(clk, '(i'//clfmt//')') kloc(3) ! this is creating a compilation error with AGRIF295 cl4 = '(i'//clfmt//')' ; WRITE(clk, cl4) kloc(3) ! this is ok with AGRIF296 WRITE(cdline,9400) TRIM(ADJUSTL(clkt)), TRIM(ADJUSTL(cdprefix)), pval, TRIM(cli), TRIM(clj), TRIM(clk), TRIM(clsuff)297 ENDIF298 !299 9100 FORMAT('MPI rank ', a)300 9200 FORMAT('found in ', a, ' MPI tasks, spread out among ranks ', a, ' to ', a)301 9300 FORMAT('kt ', a, ' ', a, ' ', 1pg11.4, ' at i j ', a, ' ', a, ' ', a, ' ', a)302 9400 FORMAT('kt ', a, ' ', a, ' ', 1pg11.4, ' at i j k ', a, ' ', a, ' ', a, ' ', a)303 !304 END SUBROUTINE wrt_line305 306 198 307 199 !!====================================================================== -
NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/tests/CANAL/MY_SRC/trazdf.F90
r12489 r13197 35 35 PUBLIC tra_zdf_imp ! called by trczdf.F90 36 36 37 !! * Substitutions 38 # include "do_loop_substitute.h90" 37 39 !!---------------------------------------------------------------------- 38 40 !! NEMO/OCE 4.0 , NEMO Consortium (2018) … … 77 79 ! JMM avoid negative salinities near river outlet ! Ugly fix 78 80 ! JMM : restore negative salinities to small salinities: 79 !!$ WHERE( pts(:,:,:,jp_sal,Kaa) < 0._wp ) pts(:,:,:,jp_sal,Kaa) = 0.1_wp81 !!$ WHERE( pts(:,:,:,jp_sal,Kaa) < 0._wp ) pts(:,:,:,jp_sal,Kaa) = 0.1_wp 80 82 !!gm 81 83 … … 95 97 ENDIF 96 98 ! ! print mean trends (used for debugging) 97 IF( ln_ctl) CALL prt_ctl( tab3d_1=pts(:,:,:,jp_tem,Kaa), clinfo1=' zdf - Ta: ', mask1=tmask, &98 & tab3d_2=pts(:,:,:,jp_sal,Kaa), clinfo2= ' Sa: ', mask2=tmask, clinfo3='tra' )99 IF(sn_cfctl%l_prtctl) CALL prt_ctl( tab3d_1=pts(:,:,:,jp_tem,Kaa), clinfo1=' zdf - Ta: ', mask1=tmask, & 100 & tab3d_2=pts(:,:,:,jp_sal,Kaa), clinfo2= ' Sa: ', mask2=tmask, clinfo3='tra' ) 99 101 ! 100 102 IF( ln_timing ) CALL timing_stop('tra_zdf') … … 154 156 IF( l_ldfslp ) THEN ! isoneutral diffusion: add the contribution 155 157 IF( ln_traldf_msc ) THEN ! MSC iso-neutral operator 156 DO jk = 2, jpkm1 157 DO jj = 2, jpjm1 158 DO ji = fs_2, fs_jpim1 ! vector opt. 159 zwt(ji,jj,jk) = zwt(ji,jj,jk) + akz(ji,jj,jk) 160 END DO 161 END DO 162 END DO 158 DO_3D_00_00( 2, jpkm1 ) 159 zwt(ji,jj,jk) = zwt(ji,jj,jk) + akz(ji,jj,jk) 160 END_3D 163 161 ELSE ! standard or triad iso-neutral operator 164 DO jk = 2, jpkm1 165 DO jj = 2, jpjm1 166 DO ji = fs_2, fs_jpim1 ! vector opt. 167 zwt(ji,jj,jk) = zwt(ji,jj,jk) + ah_wslp2(ji,jj,jk) 168 END DO 169 END DO 170 END DO 162 DO_3D_00_00( 2, jpkm1 ) 163 zwt(ji,jj,jk) = zwt(ji,jj,jk) + ah_wslp2(ji,jj,jk) 164 END_3D 171 165 ENDIF 172 166 ENDIF … … 174 168 ! Diagonal, lower (i), upper (s) (including the bottom boundary condition since avt is masked) 175 169 IF( ln_zad_Aimp ) THEN ! Adaptive implicit vertical advection 176 DO jk = 1, jpkm1 177 DO jj = 2, jpjm1 178 DO ji = fs_2, fs_jpim1 ! vector opt. (ensure same order of calculation as below if wi=0.) 179 zzwi = - p2dt * zwt(ji,jj,jk ) / e3w(ji,jj,jk ,Kmm) 180 zzws = - p2dt * zwt(ji,jj,jk+1) / e3w(ji,jj,jk+1,Kmm) 181 zwd(ji,jj,jk) = e3t(ji,jj,jk,Kaa) - zzwi - zzws & 182 & + p2dt * ( MAX( wi(ji,jj,jk ) , 0._wp ) - MIN( wi(ji,jj,jk+1) , 0._wp ) ) 183 zwi(ji,jj,jk) = zzwi + p2dt * MIN( wi(ji,jj,jk ) , 0._wp ) 184 zws(ji,jj,jk) = zzws - p2dt * MAX( wi(ji,jj,jk+1) , 0._wp ) 185 END DO 186 END DO 187 END DO 170 DO_3D_00_00( 1, jpkm1 ) 171 zzwi = - p2dt * zwt(ji,jj,jk ) / e3w(ji,jj,jk ,Kmm) 172 zzws = - p2dt * zwt(ji,jj,jk+1) / e3w(ji,jj,jk+1,Kmm) 173 zwd(ji,jj,jk) = e3t(ji,jj,jk,Kaa) - zzwi - zzws & 174 & + p2dt * ( MAX( wi(ji,jj,jk ) , 0._wp ) - MIN( wi(ji,jj,jk+1) , 0._wp ) ) 175 zwi(ji,jj,jk) = zzwi + p2dt * MIN( wi(ji,jj,jk ) , 0._wp ) 176 zws(ji,jj,jk) = zzws - p2dt * MAX( wi(ji,jj,jk+1) , 0._wp ) 177 END_3D 188 178 ELSE 189 DO jk = 1, jpkm1 190 DO jj = 2, jpjm1 191 DO ji = fs_2, fs_jpim1 ! vector opt. 192 zwi(ji,jj,jk) = - p2dt * zwt(ji,jj,jk ) / e3w(ji,jj,jk,Kmm) 193 zws(ji,jj,jk) = - p2dt * zwt(ji,jj,jk+1) / e3w(ji,jj,jk+1,Kmm) 194 zwd(ji,jj,jk) = e3t(ji,jj,jk,Kaa) - zwi(ji,jj,jk) - zws(ji,jj,jk) 195 END DO 196 END DO 197 END DO 179 DO_3D_00_00( 1, jpkm1 ) 180 zwi(ji,jj,jk) = - p2dt * zwt(ji,jj,jk ) / e3w(ji,jj,jk,Kmm) 181 zws(ji,jj,jk) = - p2dt * zwt(ji,jj,jk+1) / e3w(ji,jj,jk+1,Kmm) 182 zwd(ji,jj,jk) = e3t(ji,jj,jk,Kaa) - zwi(ji,jj,jk) - zws(ji,jj,jk) 183 END_3D 198 184 ENDIF 199 185 ! … … 217 203 ! used as a work space array: its value is modified. 218 204 ! 219 DO jj = 2, jpjm1 !* 1st recurrence: Tk = Dk - Ik Sk-1 / Tk-1 (increasing k) 220 DO ji = fs_2, fs_jpim1 ! done one for all passive tracers (so included in the IF instruction) 221 zwt(ji,jj,1) = zwd(ji,jj,1) 222 END DO 223 END DO 224 DO jk = 2, jpkm1 225 DO jj = 2, jpjm1 226 DO ji = fs_2, fs_jpim1 227 zwt(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwt(ji,jj,jk-1) 228 END DO 229 END DO 230 END DO 205 DO_2D_00_00 206 zwt(ji,jj,1) = zwd(ji,jj,1) 207 END_2D 208 DO_3D_00_00( 2, jpkm1 ) 209 zwt(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwt(ji,jj,jk-1) 210 END_3D 231 211 ! 232 212 ENDIF 233 213 ! 234 DO jj = 2, jpjm1 !* 2nd recurrence: Zk = Yk - Ik / Tk-1 Zk-1 235 DO ji = fs_2, fs_jpim1 236 pt(ji,jj,1,jn,Kaa) = e3t(ji,jj,1,Kbb) * pt(ji,jj,1,jn,Kbb) + p2dt * e3t(ji,jj,1,Kmm) * pt(ji,jj,1,jn,Krhs) 237 END DO 238 END DO 239 DO jk = 2, jpkm1 240 DO jj = 2, jpjm1 241 DO ji = fs_2, fs_jpim1 242 zrhs = e3t(ji,jj,jk,Kbb) * pt(ji,jj,jk,jn,Kbb) + p2dt * e3t(ji,jj,jk,Kmm) * pt(ji,jj,jk,jn,Krhs) ! zrhs=right hand side 243 pt(ji,jj,jk,jn,Kaa) = zrhs - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) * pt(ji,jj,jk-1,jn,Kaa) 244 END DO 245 END DO 246 END DO 214 DO_2D_00_00 215 pt(ji,jj,1,jn,Kaa) = e3t(ji,jj,1,Kbb) * pt(ji,jj,1,jn,Kbb) + p2dt * e3t(ji,jj,1,Kmm) * pt(ji,jj,1,jn,Krhs) 216 END_2D 217 DO_3D_00_00( 2, jpkm1 ) 218 zrhs = e3t(ji,jj,jk,Kbb) * pt(ji,jj,jk,jn,Kbb) + p2dt * e3t(ji,jj,jk,Kmm) * pt(ji,jj,jk,jn,Krhs) ! zrhs=right hand side 219 pt(ji,jj,jk,jn,Kaa) = zrhs - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) * pt(ji,jj,jk-1,jn,Kaa) 220 END_3D 247 221 ! 248 DO jj = 2, jpjm1 !* 3d recurrence: Xk = (Zk - Sk Xk+1 ) / Tk (result is the after tracer) 249 DO ji = fs_2, fs_jpim1 250 pt(ji,jj,jpkm1,jn,Kaa) = pt(ji,jj,jpkm1,jn,Kaa) / zwt(ji,jj,jpkm1) * tmask(ji,jj,jpkm1) 251 END DO 252 END DO 253 DO jk = jpk-2, 1, -1 254 DO jj = 2, jpjm1 255 DO ji = fs_2, fs_jpim1 256 pt(ji,jj,jk,jn,Kaa) = ( pt(ji,jj,jk,jn,Kaa) - zws(ji,jj,jk) * pt(ji,jj,jk+1,jn,Kaa) ) & 257 & / zwt(ji,jj,jk) * tmask(ji,jj,jk) 258 END DO 259 END DO 260 END DO 222 DO_2D_00_00 223 pt(ji,jj,jpkm1,jn,Kaa) = pt(ji,jj,jpkm1,jn,Kaa) / zwt(ji,jj,jpkm1) * tmask(ji,jj,jpkm1) 224 END_2D 225 DO_3DS_00_00( jpk-2, 1, -1 ) 226 pt(ji,jj,jk,jn,Kaa) = ( pt(ji,jj,jk,jn,Kaa) - zws(ji,jj,jk) * pt(ji,jj,jk+1,jn,Kaa) ) & 227 & / zwt(ji,jj,jk) * tmask(ji,jj,jk) 228 END_3D 261 229 ! ! ================= ! 262 230 END DO ! end tracer loop ! -
NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/tests/CANAL/MY_SRC/usrdef_hgr.F90
r10074 r13197 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) … … 88 90 #endif 89 91 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 * zti 96 plamu(ji,jj) = zlam0 + rn_dx * zui 97 plamv(ji,jj) = plamt(ji,jj) 98 plamf(ji,jj) = plamu(ji,jj) 99 100 pphit(ji,jj) = zphi0 + rn_dy * ztj 101 pphiv(ji,jj) = zphi0 + rn_dy * zvj 102 pphiu(ji,jj) = pphit(ji,jj) 103 pphif(ji,jj) = pphiv(ji,jj) 104 END DO 105 END DO 92 DO_2D_11_11 93 zti = FLOAT( ji - 1 + nimpp - 1 ) ; ztj = FLOAT( jj - 1 + njmpp - 1 ) 94 zui = FLOAT( ji - 1 + nimpp - 1 ) + 0.5_wp ; zvj = FLOAT( jj - 1 + njmpp - 1 ) + 0.5_wp 95 96 plamt(ji,jj) = zlam0 + rn_dx * zti 97 plamu(ji,jj) = zlam0 + rn_dx * zui 98 plamv(ji,jj) = plamt(ji,jj) 99 plamf(ji,jj) = plamu(ji,jj) 100 101 pphit(ji,jj) = zphi0 + rn_dy * ztj 102 pphiv(ji,jj) = zphi0 + rn_dy * zvj 103 pphiu(ji,jj) = pphit(ji,jj) 104 pphif(ji,jj) = pphiv(ji,jj) 105 END_2D 106 106 ! 107 107 ! Horizontal scale factors (in meters) -
NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/tests/CANAL/MY_SRC/usrdef_istate.F90
r12489 r13197 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) … … 164 166 pssh(:,1) = - ff_t(:,1) / grav * pu(:,1,1) * e2t(:,1) 165 167 DO jl=1, jpnj 166 DO jj=nldj, nlej 167 DO ji=nldi, nlei 168 pssh(ji,jj) = pssh(ji,jj-1) - ff_t(ji,jj) / grav * pu(ji,jj,1) * e2t(ji,jj) 169 END DO 170 END DO 168 DO_2D_00_00 169 pssh(ji,jj) = pssh(ji,jj-1) - ff_t(ji,jj) / grav * pu(ji,jj,1) * e2t(ji,jj) 170 END_2D 171 171 CALL lbc_lnk( 'usrdef_istate', pssh, 'T', 1. ) 172 172 END DO … … 183 183 CASE(4) ! geostrophic zonal pulse 184 184 185 DO jj=1, jpj 186 DO ji=1, jpi 187 IF ( ABS(glamt(ji,jj)) <= zjetx ) THEN 188 zdu = rn_uzonal 189 ELSEIF ( ABS(glamt(ji,jj)) <= zjetx + 100. ) THEN 190 zdu = rn_uzonal * ( ( zjetx-ABS(glamt(ji,jj)) )/100. + 1. ) 191 ELSE 192 zdu = 0. 193 END IF 194 IF ( ABS(gphit(ji,jj)) <= zjety ) THEN 195 pssh(ji,jj) = - ff_t(ji,jj) * zdu * gphit(ji,jj) * 1.e3 / grav 196 pu(ji,jj,:) = zdu 197 pts(ji,jj,:,jp_sal) = zdu / rn_uzonal + 1. 198 ELSE 199 pssh(ji,jj) = - ff_t(ji,jj) * zdu * SIGN(zjety,gphit(ji,jj)) * 1.e3 / grav 200 pu(ji,jj,:) = 0. 201 pts(ji,jj,:,jp_sal) = 1. 202 END IF 203 END DO 204 END DO 185 DO_2D_11_11 186 IF ( ABS(glamt(ji,jj)) <= zjetx ) THEN 187 zdu = rn_uzonal 188 ELSEIF ( ABS(glamt(ji,jj)) <= zjetx + 100. ) THEN 189 zdu = rn_uzonal * ( ( zjetx-ABS(glamt(ji,jj)) )/100. + 1. ) 190 ELSE 191 zdu = 0. 192 END IF 193 IF ( ABS(gphit(ji,jj)) <= zjety ) THEN 194 pssh(ji,jj) = - ff_t(ji,jj) * zdu * gphit(ji,jj) * 1.e3 / grav 195 pu(ji,jj,:) = zdu 196 pts(ji,jj,:,jp_sal) = zdu / rn_uzonal + 1. 197 ELSE 198 pssh(ji,jj) = - ff_t(ji,jj) * zdu * SIGN(zjety,gphit(ji,jj)) * 1.e3 / grav 199 pu(ji,jj,:) = 0. 200 pts(ji,jj,:,jp_sal) = 1. 201 END IF 202 END_2D 205 203 206 204 ! temperature: 207 205 pts(:,:,:,jp_tem) = 10._wp * ptmask(:,:,:) 208 206 pv(:,:,:) = 0. 209 210 207 211 208 CASE(5) ! vortex … … 220 217 zP0 = rho0 * zf0 * zumax * zlambda * SQRT(EXP(1._wp)/2._wp) 221 218 ! 222 DO jj=1, jpj 223 DO ji=1, jpi 224 zx = glamt(ji,jj) * 1.e3 225 zy = gphit(ji,jj) * 1.e3 226 ! Surface pressure: P(x,y,z) = F(z) * Psurf(x,y) 227 zpsurf = zP0 * EXP(-(zx**2+zy**2)*zr_lambda2) - rho0 * ff_t(ji,jj) * rn_uzonal * zy 228 ! Sea level: 229 pssh(ji,jj) = 0. 230 DO jl=1,5 231 zdt = pssh(ji,jj) 232 zdzF = (1._wp - EXP(zdt-zH)) / (zH - 1._wp + EXP(-zH)) ! F'(z) 233 zrho1 = rho0 * (1._wp + zn2*zdt/grav) - zdzF * zpsurf / grav ! -1/g Dz(P) = -1/g * F'(z) * Psurf(x,y) 234 pssh(ji,jj) = zpsurf / (zrho1*grav) * ptmask(ji,jj,1) ! ssh = Psurf / (Rho*g) 235 END DO 236 ! temperature: 237 DO jk=1,jpk 238 zdt = pdept(ji,jj,jk) 239 zrho1 = rho0 * (1._wp + zn2*zdt/grav) 240 IF (zdt < zH) THEN 241 zdzF = (1._wp-EXP(zdt-zH)) / (zH-1._wp + EXP(-zH)) ! F'(z) 242 zrho1 = zrho1 - zdzF * zpsurf / grav ! -1/g Dz(P) = -1/g * F'(z) * Psurf(x,y) 243 ENDIF 244 ! pts(ji,jj,jk,jp_tem) = (20._wp + (rho0-zrho1) / 0.28_wp) * ptmask(ji,jj,jk) 245 pts(ji,jj,jk,jp_tem) = (10._wp + (rho0-zrho1) / 0.28_wp) * ptmask(ji,jj,jk) 246 END DO 247 END DO 248 END DO 219 DO_2D_11_11 220 zx = glamt(ji,jj) * 1.e3 221 zy = gphit(ji,jj) * 1.e3 222 ! Surface pressure: P(x,y,z) = F(z) * Psurf(x,y) 223 zpsurf = zP0 * EXP(-(zx**2+zy**2)*zr_lambda2) - rho0 * ff_t(ji,jj) * rn_uzonal * zy 224 ! Sea level: 225 pssh(ji,jj) = 0. 226 DO jl=1,5 227 zdt = pssh(ji,jj) 228 zdzF = (1._wp - EXP(zdt-zH)) / (zH - 1._wp + EXP(-zH)) ! F'(z) 229 zrho1 = rho0 * (1._wp + zn2*zdt/grav) - zdzF * zpsurf / grav ! -1/g Dz(P) = -1/g * F'(z) * Psurf(x,y) 230 pssh(ji,jj) = zpsurf / (zrho1*grav) * ptmask(ji,jj,1) ! ssh = Psurf / (Rho*g) 231 END DO 232 ! temperature: 233 DO jk=1,jpk 234 zdt = pdept(ji,jj,jk) 235 zrho1 = rho0 * (1._wp + zn2*zdt/grav) 236 IF (zdt < zH) THEN 237 zdzF = (1._wp-EXP(zdt-zH)) / (zH-1._wp + EXP(-zH)) ! F'(z) 238 zrho1 = zrho1 - zdzF * zpsurf / grav ! -1/g Dz(P) = -1/g * F'(z) * Psurf(x,y) 239 ENDIF 240 ! pts(ji,jj,jk,jp_tem) = (20._wp + (rho0-zrho1) / 0.28_wp) * ptmask(ji,jj,jk) 241 pts(ji,jj,jk,jp_tem) = (10._wp + (rho0-zrho1) / 0.28_wp) * ptmask(ji,jj,jk) 242 END DO 243 END_2D 249 244 ! 250 245 ! salinity: … … 253 248 ! velocities: 254 249 za = 2._wp * zP0 / zlambda**2 255 DO jj=1, jpj 256 DO ji=1, jpim1 257 zx = glamu(ji,jj) * 1.e3 258 zy = gphiu(ji,jj) * 1.e3 259 DO jk=1, jpk 260 zdu = 0.5_wp * (pdept(ji,jj,jk) + pdept(ji+1,jj,jk)) 261 IF (zdu < zH) THEN 262 zf = (zH-1._wp-zdu+EXP(zdu-zH)) / (zH-1._wp+EXP(-zH)) 263 zdyPs = - za * zy * EXP(-(zx**2+zy**2)*zr_lambda2) - rho0 * ff_t(ji,jj) * rn_uzonal 264 pu(ji,jj,jk) = - zf / ( rho0 * ff_t(ji,jj) ) * zdyPs * ptmask(ji,jj,jk) * ptmask(ji+1,jj,jk) 265 ELSE 266 pu(ji,jj,jk) = 0._wp 267 ENDIF 268 END DO 269 END DO 270 END DO 271 ! 272 DO jj=1, jpjm1 273 DO ji=1, jpi 274 zx = glamv(ji,jj) * 1.e3 275 zy = gphiv(ji,jj) * 1.e3 276 DO jk=1, jpk 277 zdv = 0.5_wp * (pdept(ji,jj,jk) + pdept(ji,jj+1,jk)) 278 IF (zdv < zH) THEN 279 zf = (zH-1._wp-zdv+EXP(zdv-zH)) / (zH-1._wp+EXP(-zH)) 280 zdxPs = - za * zx * EXP(-(zx**2+zy**2)*zr_lambda2) 281 pv(ji,jj,jk) = zf / ( rho0 * ff_f(ji,jj) ) * zdxPs * ptmask(ji,jj,jk) * ptmask(ji,jj+1,jk) 282 ELSE 283 pv(ji,jj,jk) = 0._wp 284 ENDIF 285 END DO 286 END DO 287 END DO 250 DO_2D_00_00 251 zx = glamu(ji,jj) * 1.e3 252 zy = gphiu(ji,jj) * 1.e3 253 DO jk=1, jpk 254 zdu = 0.5_wp * (pdept(ji,jj,jk) + pdept(ji+1,jj,jk)) 255 IF (zdu < zH) THEN 256 zf = (zH-1._wp-zdu+EXP(zdu-zH)) / (zH-1._wp+EXP(-zH)) 257 zdyPs = - za * zy * EXP(-(zx**2+zy**2)*zr_lambda2) - rho0 * ff_t(ji,jj) * rn_uzonal 258 pu(ji,jj,jk) = - zf / ( rho0 * ff_t(ji,jj) ) * zdyPs * ptmask(ji,jj,jk) * ptmask(ji+1,jj,jk) 259 ELSE 260 pu(ji,jj,jk) = 0._wp 261 ENDIF 262 END DO 263 END_2D 264 ! 265 DO_2D_00_00 266 zx = glamv(ji,jj) * 1.e3 267 zy = gphiv(ji,jj) * 1.e3 268 DO jk=1, jpk 269 zdv = 0.5_wp * (pdept(ji,jj,jk) + pdept(ji,jj+1,jk)) 270 IF (zdv < zH) THEN 271 zf = (zH-1._wp-zdv+EXP(zdv-zH)) / (zH-1._wp+EXP(-zH)) 272 zdxPs = - za * zx * EXP(-(zx**2+zy**2)*zr_lambda2) 273 pv(ji,jj,jk) = zf / ( rho0 * ff_f(ji,jj) ) * zdxPs * ptmask(ji,jj,jk) * ptmask(ji,jj+1,jk) 274 ELSE 275 pv(ji,jj,jk) = 0._wp 276 ENDIF 277 END DO 278 END_2D 288 279 ! 289 280 END SELECT 290 281 291 282 IF (ln_sshnoise) THEN 292 283 CALL RANDOM_NUMBER(zrandom) … … 294 285 END IF 295 286 CALL lbc_lnk( 'usrdef_istate', pssh, 'T', 1. ) 296 CALL lbc_lnk( 'usrdef_istate', pts, 'T', 1. ) 297 CALL lbc_lnk( 'usrdef_istate', pu, 'U', -1. ) 298 CALL lbc_lnk( 'usrdef_istate', pv, 'V', -1. ) 287 CALL lbc_lnk( 'usrdef_istate', pts , 'T', 1. ) 288 CALL lbc_lnk_multi( 'usrdef_istate', pu, 'U', -1., pv, 'V', -1. ) 299 289 300 290 END SUBROUTINE usr_def_istate -
NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/tests/CANAL/MY_SRC/usrdef_sbc.F90
r12377 r13197 38 38 CONTAINS 39 39 40 SUBROUTINE usrdef_sbc_oce( kt, K mm, Kbb )40 SUBROUTINE usrdef_sbc_oce( kt, Kbb ) 41 41 !!--------------------------------------------------------------------- 42 42 !! *** ROUTINE usr_def_sbc *** … … 53 53 !!---------------------------------------------------------------------- 54 54 INTEGER, INTENT(in) :: kt ! ocean time step 55 INTEGER, INTENT(in) :: Kbb , Kmm! ocean time index55 INTEGER, INTENT(in) :: Kbb ! ocean time index 56 56 INTEGER :: ji, jj ! dummy loop indices 57 57 REAL(wp) :: zrhoair = 1.22 ! approximate air density [Kg/m3] … … 86 86 87 87 WHERE( ABS(gphit) <= rn_windszy/2. ) 88 zwndrel(:,:) = rn_u10 - rn_uofac * uu(:,:,1,K mm)88 zwndrel(:,:) = rn_u10 - rn_uofac * uu(:,:,1,Kbb) 89 89 ELSEWHERE 90 zwndrel(:,:) = - rn_uofac * uu(:,:,1,K mm)90 zwndrel(:,:) = - rn_uofac * uu(:,:,1,Kbb) 91 91 END WHERE 92 92 utau(:,:) = zrhocd * zwndrel(:,:) * zwndrel(:,:) 93 93 94 zwndrel(:,:) = - rn_uofac * vv(:,:,1,K mm)94 zwndrel(:,:) = - rn_uofac * vv(:,:,1,Kbb) 95 95 vtau(:,:) = zrhocd * zwndrel(:,:) * zwndrel(:,:) 96 96 -
NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/tests/CANAL/MY_SRC/usrdef_zgr.F90
r12377 r13197 204 204 CALL lbc_lnk( 'usrdef_zgr', z2d, 'T', 1. ) ! set surrounding land to zero (here jperio=0 ==>> closed) 205 205 ! 206 k_bot(:,:) = INT( z2d(:,:) )! =jpkm1 over the ocean point, =0 elsewhere206 k_bot(:,:) = NINT( z2d(:,:) ) ! =jpkm1 over the ocean point, =0 elsewhere 207 207 ! 208 208 k_top(:,:) = MIN( 1 , k_bot(:,:) ) ! = 1 over the ocean point, =0 elsewhere
Note: See TracChangeset
for help on using the changeset viewer.
