- Timestamp:
- 2017-06-15T08:44:30+02:00 (7 years ago)
- Location:
- branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC
- Files:
-
- 6 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90
r8143 r8178 25 25 !! dia_wri_state : create an output NetCDF file for a single instantaeous ocean state and forcing fields 26 26 !!---------------------------------------------------------------------- 27 USE oce ! ocean dynamics and tracers 28 USE dom_oce ! ocean space and time domain 29 USE dynadv, ONLY: ln_dynadv_vec 30 USE zdf_oce ! ocean vertical physics 31 USE zdfdrg ! ocean vertical physics: top/bottom friction 32 USE ldftra ! lateral physics: eddy diffusivity coef. 33 USE ldfdyn ! lateral physics: eddy viscosity coef. 34 USE sbc_oce ! Surface boundary condition: ocean fields 35 USE sbc_ice ! Surface boundary condition: ice fields 36 USE icb_oce ! Icebergs 37 USE icbdia ! Iceberg budgets 38 USE sbcssr ! restoring term toward SST/SSS climatology 39 USE phycst ! physical constants 40 USE zdfmxl ! mixed layer 41 USE dianam ! build name of file (routine) 42 ! USE zdfddm ! vertical physics: double diffusion 43 USE diahth ! thermocline diagnostics 44 USE wet_dry ! wetting and drying 45 USE sbcwave ! wave parameters 27 USE oce ! ocean dynamics and tracers 28 USE dom_oce ! ocean space and time domain 29 USE phycst ! physical constants 30 USE dianam ! build name of file (routine) 31 USE diahth ! thermocline diagnostics 32 USE dynadv , ONLY: ln_dynadv_vec 33 USE icb_oce ! Icebergs 34 USE icbdia ! Iceberg budgets 35 USE ldftra ! lateral physics: eddy diffusivity coef. 36 USE ldfdyn ! lateral physics: eddy viscosity coef. 37 USE sbc_oce ! Surface boundary condition: ocean fields 38 USE sbc_ice ! Surface boundary condition: ice fields 39 USE sbcssr ! restoring term toward SST/SSS climatology 40 USE sbcwave ! wave parameters 41 USE wet_dry ! wetting and drying 42 USE zdf_oce ! ocean vertical physics 43 USE zdfdrg ! ocean vertical physics: top/bottom friction 44 USE zdfmxl ! mixed layer 46 45 ! 47 USE lbclnk 48 USE in_out_manager 49 USE diatmb 50 USE dia25h 51 USE iom 52 USE ioipsl 46 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 47 USE in_out_manager ! I/O manager 48 USE diatmb ! Top,middle,bottom output 49 USE dia25h ! 25h Mean output 50 USE iom ! 51 USE ioipsl ! 53 52 54 53 #if defined key_lim2 … … 194 193 ENDIF 195 194 196 CALL iom_put( "uoce", un(:,:,:) )! 3D i-current197 CALL iom_put( "ssu", un(:,:,1) )! surface i-current195 CALL iom_put( "uoce", un(:,:,:) ) ! 3D i-current 196 CALL iom_put( "ssu", un(:,:,1) ) ! surface i-current 198 197 IF ( iom_use("sbu") ) THEN 199 198 DO jj = 1, jpj … … 206 205 ENDIF 207 206 208 CALL iom_put( "voce", vn(:,:,:) )! 3D j-current209 CALL iom_put( "ssv", vn(:,:,1) )! surface j-current207 CALL iom_put( "voce", vn(:,:,:) ) ! 3D j-current 208 CALL iom_put( "ssv", vn(:,:,1) ) ! surface j-current 210 209 IF ( iom_use("sbv") ) THEN 211 210 DO jj = 1, jpj … … 229 228 ENDIF 230 229 231 CALL iom_put( "avt" , avt )! T vert. eddy diff. coef.232 CALL iom_put( "avs" , avs )! S vert. eddy diff. coef.233 CALL iom_put( "avm" , avm u )! T vert. eddy visc. coef.230 CALL iom_put( "avt" , avt ) ! T vert. eddy diff. coef. 231 CALL iom_put( "avs" , avs ) ! S vert. eddy diff. coef. 232 CALL iom_put( "avm" , avm ) ! T vert. eddy visc. coef. 234 233 235 234 IF( iom_use('logavt') ) CALL iom_put( "logavt", LOG( MAX( 1.e-20_wp, avt(:,:,:) ) ) ) … … 247 246 END DO 248 247 CALL lbc_lnk( z2d, 'T', 1. ) 249 CALL iom_put( "sstgrad2", z2d )! square of module of sst gradient248 CALL iom_put( "sstgrad2", z2d ) ! square of module of sst gradient 250 249 z2d(:,:) = SQRT( z2d(:,:) ) 251 CALL iom_put( "sstgrad" , z2d )! module of sst gradient250 CALL iom_put( "sstgrad" , z2d ) ! module of sst gradient 252 251 ENDIF 253 252 … … 262 261 END DO 263 262 END DO 264 CALL iom_put( "heatc", (rau0 * rcp) * z2d )! vertically integrated heat content (J/m2)263 CALL iom_put( "heatc", rau0_rcp * z2d ) ! vertically integrated heat content (J/m2) 265 264 ENDIF 266 265 … … 274 273 END DO 275 274 END DO 276 CALL iom_put( "saltc", rau0 * z2d ) ! vertically integrated salt content (PSU*kg/m2)275 CALL iom_put( "saltc", rau0 * z2d ) ! vertically integrated salt content (PSU*kg/m2) 277 276 ENDIF 278 277 ! 279 278 IF ( iom_use("eken") ) THEN 280 z3d(:,:,jk) = 0._wp ! kinetic energy279 z3d(:,:,jk) = 0._wp 281 280 DO jk = 1, jpkm1 282 281 DO jj = 2, jpjm1 283 282 DO ji = fs_2, fs_jpim1 ! vector opt. 284 zztmp = r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) 285 z3d(ji,jj,jk) = 0.25_wp * zztmp * ( & 286 & un(ji-1,jj,jk)**2 * e2u(ji-1,jj) * e3u_n(ji-1,jj,jk) & 287 & + un(ji ,jj,jk)**2 * e2u(ji ,jj) * e3u_n(ji ,jj,jk) & 288 & + vn(ji,jj-1,jk)**2 * e1v(ji,jj-1) * e3v_n(ji,jj-1,jk) & 289 & + vn(ji,jj ,jk)**2 * e1v(ji,jj ) * e3v_n(ji,jj ,jk) ) 283 zztmp = 0.25_wp * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) 284 z3d(ji,jj,jk) = zztmp * ( un(ji-1,jj,jk)**2 * e2u(ji-1,jj) * e3u_n(ji-1,jj,jk) & 285 & + un(ji ,jj,jk)**2 * e2u(ji ,jj) * e3u_n(ji ,jj,jk) & 286 & + vn(ji,jj-1,jk)**2 * e1v(ji,jj-1) * e3v_n(ji,jj-1,jk) & 287 & + vn(ji,jj ,jk)**2 * e1v(ji,jj ) * e3v_n(ji,jj ,jk) ) 290 288 END DO 291 289 END DO 292 290 END DO 293 291 CALL lbc_lnk( z3d, 'T', 1. ) 294 CALL iom_put( "eken", z3d ) 292 CALL iom_put( "eken", z3d ) ! kinetic energy 295 293 ENDIF 296 294 ! … … 304 302 z2d(:,:) = z2d(:,:) + z3d(:,:,jk) 305 303 END DO 306 CALL iom_put( "u_masstr" , z3d )! mass transport in i-direction307 CALL iom_put( "u_masstr_vint", z2d ) 304 CALL iom_put( "u_masstr" , z3d ) ! mass transport in i-direction 305 CALL iom_put( "u_masstr_vint", z2d ) ! mass transport in i-direction vertical sum 308 306 ENDIF 309 307 310 308 IF( iom_use("u_heattr") ) THEN 311 z2d(:,:) = 0. e0309 z2d(:,:) = 0._wp 312 310 DO jk = 1, jpkm1 313 311 DO jj = 2, jpjm1 … … 318 316 END DO 319 317 CALL lbc_lnk( z2d, 'U', -1. ) 320 CALL iom_put( "u_heattr", (0.5 * rcp)* z2d ) ! heat transport in i-direction318 CALL iom_put( "u_heattr", 0.5*rcp * z2d ) ! heat transport in i-direction 321 319 ENDIF 322 320 … … 331 329 END DO 332 330 CALL lbc_lnk( z2d, 'U', -1. ) 333 CALL iom_put( "u_salttr", 0.5 * z2d ) 331 CALL iom_put( "u_salttr", 0.5 * z2d ) ! heat transport in i-direction 334 332 ENDIF 335 333 … … 340 338 z3d(:,:,jk) = rau0 * vn(:,:,jk) * e1v(:,:) * e3v_n(:,:,jk) * vmask(:,:,jk) 341 339 END DO 342 CALL iom_put( "v_masstr", z3d ) 340 CALL iom_put( "v_masstr", z3d ) ! mass transport in j-direction 343 341 ENDIF 344 342 … … 353 351 END DO 354 352 CALL lbc_lnk( z2d, 'V', -1. ) 355 CALL iom_put( "v_heattr", (0.5 * rcp)* z2d ) ! heat transport in j-direction353 CALL iom_put( "v_heattr", 0.5*rcp * z2d ) ! heat transport in j-direction 356 354 ENDIF 357 355 358 356 IF( iom_use("v_salttr") ) THEN 359 z2d(:,:) = 0. e0357 z2d(:,:) = 0._wp 360 358 DO jk = 1, jpkm1 361 359 DO jj = 2, jpjm1 … … 366 364 END DO 367 365 CALL lbc_lnk( z2d, 'V', -1. ) 368 CALL iom_put( "v_salttr", 0.5 * z2d ) ! heat transport in j-direction 369 ENDIF 370 371 ! Vertical integral of temperature 366 CALL iom_put( "v_salttr", 0.5 * z2d ) ! heat transport in j-direction 367 ENDIF 368 372 369 IF( iom_use("tosmint") ) THEN 373 z2d(:,:) =0._wp370 z2d(:,:) = 0._wp 374 371 DO jk = 1, jpkm1 375 372 DO jj = 2, jpjm1 376 373 DO ji = fs_2, fs_jpim1 ! vector opt. 377 z2d(ji,jj) = z2d(ji,jj) + rau0 *e3t_n(ji,jj,jk) * tsn(ji,jj,jk,jp_tem)374 z2d(ji,jj) = z2d(ji,jj) + e3t_n(ji,jj,jk) * tsn(ji,jj,jk,jp_tem) 378 375 END DO 379 376 END DO 380 377 END DO 381 378 CALL lbc_lnk( z2d, 'T', -1. ) 382 CALL iom_put( "tosmint", z2d ) 383 ENDIF 384 385 ! Vertical integral of salinity 379 CALL iom_put( "tosmint", rau0 * z2d ) ! Vertical integral of temperature 380 ENDIF 386 381 IF( iom_use("somint") ) THEN 387 382 z2d(:,:)=0._wp … … 389 384 DO jj = 2, jpjm1 390 385 DO ji = fs_2, fs_jpim1 ! vector opt. 391 z2d(ji,jj) = z2d(ji,jj) + rau0 *e3t_n(ji,jj,jk) * tsn(ji,jj,jk,jp_sal)386 z2d(ji,jj) = z2d(ji,jj) + e3t_n(ji,jj,jk) * tsn(ji,jj,jk,jp_sal) 392 387 END DO 393 388 END DO 394 389 END DO 395 390 CALL lbc_lnk( z2d, 'T', -1. ) 396 CALL iom_put( "somint", z2d ) 397 ENDIF 398 399 CALL iom_put( "bn2", rn2 ) !Brunt-Vaisala buoyancy frequency (N^2) 400 ! 401 402 IF (ln_diatmb) THEN ! If we want tmb values 403 CALL dia_tmb 404 ENDIF 405 IF (ln_dia25h) THEN 406 CALL dia_25h( kt ) 407 ENDIF 391 CALL iom_put( "somint", rau0 * z2d ) ! Vertical integral of salinity 392 ENDIF 393 394 CALL iom_put( "bn2", rn2 ) ! Brunt-Vaisala buoyancy frequency (N^2) 395 ! 396 397 IF (ln_diatmb) CALL dia_tmb ! tmb values 398 399 IF (ln_dia25h) CALL dia_25h( kt ) ! 25h averaging 408 400 409 401 IF( nn_timing == 1 ) CALL timing_stop('dia_wri') … … 729 721 CALL histdef( nid_W, "votkeavt", "Vertical Eddy Diffusivity" , "m2/s" , & ! avt 730 722 & jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout ) 731 CALL histdef( nid_W, "votkeavm", "Vertical Eddy Viscosity" , "m2/s" , & ! avm u723 CALL histdef( nid_W, "votkeavm", "Vertical Eddy Viscosity" , "m2/s" , & ! avm 732 724 & jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout ) 733 725 … … 856 848 CALL histwrite( nid_W, "vovecrtz", it, wn , ndim_T, ndex_T ) ! vert. current 857 849 CALL histwrite( nid_W, "votkeavt", it, avt , ndim_T, ndex_T ) ! T vert. eddy diff. coef. 858 CALL histwrite( nid_W, "votkeavm", it, avm u, ndim_T, ndex_T ) ! T vert. eddy visc. coef.850 CALL histwrite( nid_W, "votkeavm", it, avm , ndim_T, ndex_T ) ! T vert. eddy visc. coef. 859 851 IF( ln_zdfddm ) THEN 860 852 CALL histwrite( nid_W, "voddmavs", it, avs , ndim_T, ndex_T ) ! S vert. eddy diff. coef. -
branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf.F90
r7753 r8178 37 37 38 38 ! ! Parameter to control the type of lateral viscous operator 39 INTEGER, PARAMETER, PUBLIC :: np_ERROR =-10 ! error in setting the operator40 INTEGER, PARAMETER, PUBLIC :: np_no_ldf = 00 ! without operator (i.e. no lateral viscous trend)39 INTEGER, PARAMETER, PUBLIC :: np_ERROR =-10 !: error in setting the operator 40 INTEGER, PARAMETER, PUBLIC :: np_no_ldf = 00 !: without operator (i.e. no lateral viscous trend) 41 41 ! !! laplacian ! bilaplacian ! 42 INTEGER, PARAMETER, PUBLIC :: np_lap = 10 , np_blp = 20 ! iso-level operator43 INTEGER, PARAMETER, PUBLIC :: np_lap_i = 11 ! iso-neutral or geopotential operator42 INTEGER, PARAMETER, PUBLIC :: np_lap = 10 , np_blp = 20 !: iso-level operator 43 INTEGER, PARAMETER, PUBLIC :: np_lap_i = 11 !: iso-neutral or geopotential operator 44 44 45 INTEGER :: nldf !type of lateral diffusion used defined from ln_dynldf_... (namlist logicals)45 INTEGER, PUBLIC :: nldf !: type of lateral diffusion used defined from ln_dynldf_... (namlist logicals) 46 46 47 47 !! * Substitutions -
branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf_iso.F90
r6140 r8178 37 37 PUBLIC dyn_ldf_iso_alloc ! called by nemogcm.F90 38 38 39 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: akzu, akzv !: vertical component of rotated lateral viscosity 40 39 41 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: zfuw, zdiu, zdju, zdj1u ! 2D workspace (dyn_ldf_iso) 40 42 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: zfvw, zdiv, zdjv, zdj1v ! - - … … 53 55 !! *** ROUTINE dyn_ldf_iso_alloc *** 54 56 !!---------------------------------------------------------------------- 55 ALLOCATE( zfuw(jpi,jpk) , zdiu(jpi,jpk) , zdju(jpi,jpk) , zdj1u(jpi,jpk) , &56 & zfvw(jpi,jpk) , zdiv(jpi,jpk) , zdjv(jpi,jpk) , zdj1v(jpi,jpk) , STAT=dyn_ldf_iso_alloc )57 ALLOCATE( akzu(jpi,jpj,jpk) , zfuw(jpi,jpk) , zdiu(jpi,jpk) , zdju(jpi,jpk) , zdj1u(jpi,jpk) , & 58 & akzv(jpi,jpj,jpk) , zfvw(jpi,jpk) , zdiv(jpi,jpk) , zdjv(jpi,jpk) , zdj1v(jpi,jpk) , STAT=dyn_ldf_iso_alloc ) 57 59 ! 58 60 IF( dyn_ldf_iso_alloc /= 0 ) CALL ctl_warn('dyn_ldf_iso_alloc: array allocate failed.') … … 99 101 !! 100 102 !! ** Action : 101 !! Update (ua,va) arrays with the before geopotential biharmonic 102 !! mixing trend. 103 !! Update (avmu,avmv) to accompt for the diagonal vertical component 104 !! of the rotated operator in dynzdf module 103 !! -(ua,va) updated with the before geopotential harmonic mixing trend 104 !! -(akzu,akzv) to accompt for the diagonal vertical component 105 !! of the rotated operator in dynzdf module 105 106 !!---------------------------------------------------------------------- 106 107 INTEGER, INTENT( in ) :: kt ! ocean time-step index … … 144 145 CALL lbc_lnk( uslp , 'U', -1. ) ; CALL lbc_lnk( vslp , 'V', -1. ) 145 146 CALL lbc_lnk( wslpi, 'W', -1. ) ; CALL lbc_lnk( wslpj, 'W', -1. ) 146 147 !!bug 148 IF( kt == nit000 ) then 149 IF(lwp) WRITE(numout,*) ' max slop: u', SQRT( MAXVAL(uslp*uslp)), ' v ', SQRT(MAXVAL(vslp)), & 150 & ' wi', sqrt(MAXVAL(wslpi)) , ' wj', sqrt(MAXVAL(wslpj)) 151 endif 152 !!end 153 ENDIF 147 ! 148 ENDIF 154 149 155 150 ! ! =============== … … 365 360 + zcoef4 * ( zdj1u(ji,jk-1) + zdju (ji ,jk-1) & 366 361 +zdj1u(ji,jk ) + zdju (ji ,jk ) ) 367 ! update avmu (add isopycnal vertical coefficient to avmu)368 ! Caution: zcoef0 include rn_aht_0, so divided by rn_aht_0 to obtain slp^2 * rn_aht_0369 a vmu(ji,jj,jk) = avmu(ji,jj,jk) +( zuwslpi * zuwslpi + zuwslpj * zuwslpj ) / rn_aht_0362 ! vertical mixing coefficient (akzu) 363 ! Note: zcoef0 include rn_aht_0, so divided by rn_aht_0 to obtain slp^2 * rn_aht_0 364 akzu(ji,jj,jk) = ( zuwslpi * zuwslpi + zuwslpj * zuwslpj ) / rn_aht_0 370 365 END DO 371 366 END DO … … 391 386 & + zcoef4 * ( zdjv (ji,jk-1) + zdj1v(ji ,jk-1) & 392 387 & +zdjv (ji,jk ) + zdj1v(ji ,jk ) ) 393 ! update avmv (add isopycnal vertical coefficient to avmv)394 ! Caution: zcoef0 include rn_aht_0, so divided by rn_aht_0 to obtain slp^2 * rn_aht_0395 a vmv(ji,jj,jk) = avmv(ji,jj,jk) +( zvwslpi * zvwslpi + zvwslpj * zvwslpj ) / rn_aht_0388 ! vertical mixing coefficient (akzv) 389 ! Note: zcoef0 include rn_aht_0, so divided by rn_aht_0 to obtain slp^2 * rn_aht_0 390 akzv(ji,jj,jk) = ( zvwslpi * zvwslpi + zvwslpj * zvwslpj ) / rn_aht_0 396 391 END DO 397 392 END DO -
branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf.F90
r8160 r8178 6 6 !! History : 1.0 ! 2005-11 (G. Madec) Original code 7 7 !! 3.3 ! 2010-10 (C. Ethe, G. Madec) reorganisation of initialisation phase 8 !! 4.0 ! 2017-06 (G. Madec) remove the explicit time-stepping option 9 !!---------------------------------------------------------------------- 10 11 !!---------------------------------------------------------------------- 12 !! dyn_zdf : Update the momentum trend with the vertical diffusion8 !! 4.0 ! 2017-06 (G. Madec) remove the explicit time-stepping option + avm at t-point 9 !!---------------------------------------------------------------------- 10 11 !!---------------------------------------------------------------------- 12 !! dyn_zdf : compute the after velocity through implicit calculation of vertical mixing 13 13 !!---------------------------------------------------------------------- 14 14 USE oce ! ocean dynamics and tracers variables … … 18 18 USE zdf_oce ! ocean vertical physics variables 19 19 USE zdfdrg ! vertical physics: top/bottom drag coef. 20 USE dynadv , ONLY: ln_dynadv_vec ! dynamics: advection form 20 USE dynadv ,ONLY: ln_dynadv_vec ! dynamics: advection form 21 USE dynldf ,ONLY: nldf, np_lap_i ! dynamics: type of lateral mixing 22 USE dynldf_iso,ONLY: akzu, akzv ! dynamics: vertical component of rotated lateral mixing 21 23 USE ldfdyn ! lateral diffusion: eddy viscosity coef. 22 24 USE trd_oce ! trends: ocean variables … … 26 28 USE lib_mpp ! MPP library 27 29 USE prtctl ! Print control 28 USE wrk_nemo ! Memory Allocation29 30 USE timing ! Timing 30 31 … … 49 50 !! *** ROUTINE dyn_zdf *** 50 51 !! 51 !! ** Purpose : compute the vertical ocean dynamics physics. 52 !!--------------------------------------------------------------------- 53 INTEGER, INTENT( in ) :: kt ! ocean time-step index 54 ! 55 REAL(wp), POINTER, DIMENSION(:,:,:) :: ztrdu, ztrdv 56 !!--------------------------------------------------------------------- 57 ! 58 IF( nn_timing == 1 ) CALL timing_start('dyn_zdf') 59 ! 60 ! !* set time step 61 IF( neuler == 0 .AND. kt == nit000 ) THEN ; r2dt = rdt ! = rdt (restart with Euler time stepping) 62 ELSEIF( kt <= nit000 + 1 ) THEN ; r2dt = 2. * rdt ! = 2 rdt (leapfrog) 63 ENDIF 64 65 IF( l_trddyn ) THEN !* temporary save of ta and sa trends 66 CALL wrk_alloc( jpi, jpj, jpk, ztrdu, ztrdv ) 67 ztrdu(:,:,:) = ua(:,:,:) 68 ztrdv(:,:,:) = va(:,:,:) 69 ENDIF 70 ! !* compute lateral mixing trend and add it to the general trend 71 ! 72 CALL dyn_zdf_imp( kt, r2dt ) 73 ! 74 IF( l_trddyn ) THEN ! save the vertical diffusive trends for further diagnostics 75 ztrdu(:,:,:) = ( ua(:,:,:) - ub(:,:,:) ) / r2dt - ztrdu(:,:,:) 76 ztrdv(:,:,:) = ( va(:,:,:) - vb(:,:,:) ) / r2dt - ztrdv(:,:,:) 77 CALL trd_dyn( ztrdu, ztrdv, jpdyn_zdf, kt ) 78 CALL wrk_dealloc( jpi, jpj, jpk, ztrdu, ztrdv ) 79 ENDIF 80 ! ! print mean trends (used for debugging) 81 IF(ln_ctl) CALL prt_ctl( tab3d_1=ua, clinfo1=' zdf - Ua: ', mask1=umask, & 82 & tab3d_2=va, clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' ) 83 ! 84 IF( nn_timing == 1 ) CALL timing_stop('dyn_zdf') 85 ! 86 END SUBROUTINE dyn_zdf 87 88 89 SUBROUTINE dyn_zdf_imp( kt, p2dt ) 90 !!---------------------------------------------------------------------- 91 !! *** ROUTINE dyn_zdf_imp *** 92 !! 93 !! ** Purpose : Compute the trend due to the vert. momentum diffusion 52 !! ** Purpose : compute the trend due to the vert. momentum diffusion 94 53 !! together with the Leap-Frog time stepping using an 95 54 !! implicit scheme. … … 100 59 !! - update the after velocity with the implicit vertical mixing. 101 60 !! This requires to solver the following system: 102 !! ua = ua + 1/e3u_a dk+1[ avmu/ e3uw_a dk[ua] ]61 !! ua = ua + 1/e3u_a dk+1[ mi(avm) / e3uw_a dk[ua] ] 103 62 !! with the following surface/top/bottom boundary condition: 104 63 !! surface: wind stress input (averaged over kt-1/2 & kt+1/2) 105 64 !! top & bottom : top stress (iceshelf-ocean) & bottom stress (cf zdfdrg.F90) 106 65 !! 107 !! ** Action : (ua,va) after velocity66 !! ** Action : (ua,va) after velocity 108 67 !!--------------------------------------------------------------------- 109 68 INTEGER , INTENT(in) :: kt ! ocean time-step index 110 REAL(wp), INTENT(in) :: p2dt ! vertical profile of tracer time-step111 !112 INTEGER :: ji, jj, jk ! dummy loop indices113 INTEGER :: ikbu, ikbv ! local integers114 REAL(wp) :: zzw i, ze3ua ! local scalars115 REAL(wp) :: zzws, ze3va ! - -116 REAL(wp), DIMENSION( jpi,jpj,jpk) :: zwi, zwd, zws117 !!--------------------------------------------------------------------- -118 ! 119 IF( nn_timing == 1 ) CALL timing_start('dyn_zdf_imp')120 ! 121 IF( kt == nit000 ) THEN 69 ! 70 INTEGER :: ji, jj, jk ! dummy loop indices 71 INTEGER :: iku, ikv ! local integers 72 REAL(wp) :: zzwi, ze3ua, zdt ! local scalars 73 REAL(wp) :: zzws, ze3va ! - - 74 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zwi, zwd, zws ! 3D workspace 75 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ztrdu, ztrdv ! - - 76 !!--------------------------------------------------------------------- 77 ! 78 IF( nn_timing == 1 ) CALL timing_start('dyn_zdf') 79 ! 80 IF( kt == nit000 ) THEN !* initialization 122 81 IF(lwp) WRITE(numout,*) 123 82 IF(lwp) WRITE(numout,*) 'dyn_zdf_imp : vertical momentum diffusion implicit operator' … … 128 87 ENDIF 129 88 ENDIF 130 ! 131 ! !== Time step dynamics ==! 132 ! 133 IF( ln_dynadv_vec .OR. ln_linssh ) THEN ! applied on velocity 134 DO jk = 1, jpkm1 135 ua(:,:,jk) = ( ub(:,:,jk) + p2dt * ua(:,:,jk) ) * umask(:,:,jk) 136 va(:,:,jk) = ( vb(:,:,jk) + p2dt * va(:,:,jk) ) * vmask(:,:,jk) 137 END DO 138 ELSE ! applied on thickness weighted velocity 89 ! !* set time step 90 IF( neuler == 0 .AND. kt == nit000 ) THEN ; r2dt = rdt ! = rdt (restart with Euler time stepping) 91 ELSEIF( kt <= nit000 + 1 ) THEN ; r2dt = 2. * rdt ! = 2 rdt (leapfrog) 92 ENDIF 93 94 IF( l_trddyn ) THEN !* temporary save of ta and sa trends 95 ALLOCATE( ztrdu(jpi,jpj,jpk), ztrdv(jpi,jpj,jpk) ) 96 ztrdu(:,:,:) = ua(:,:,:) 97 ztrdv(:,:,:) = va(:,:,:) 98 ENDIF 99 ! 100 ! !== RHS: Leap-Frog time stepping on all trends but the vertical mixing ==! (put in ua,va) 101 ! 102 ! ! time stepping except vertical diffusion 103 IF( ln_dynadv_vec .OR. ln_linssh ) THEN ! applied on velocity 104 DO jk = 1, jpkm1 105 ua(:,:,jk) = ( ub(:,:,jk) + r2dt * ua(:,:,jk) ) * umask(:,:,jk) 106 va(:,:,jk) = ( vb(:,:,jk) + r2dt * va(:,:,jk) ) * vmask(:,:,jk) 107 END DO 108 ELSE ! applied on thickness weighted velocity 139 109 DO jk = 1, jpkm1 140 110 ua(:,:,jk) = ( e3u_b(:,:,jk) * ub(:,:,jk) & 141 & + p2dt * e3u_n(:,:,jk) * ua(:,:,jk) ) / e3u_a(:,:,jk) * umask(:,:,jk)111 & + r2dt * e3u_n(:,:,jk) * ua(:,:,jk) ) / e3u_a(:,:,jk) * umask(:,:,jk) 142 112 va(:,:,jk) = ( e3v_b(:,:,jk) * vb(:,:,jk) & 143 & + p2dt * e3v_n(:,:,jk) * va(:,:,jk) ) / e3v_a(:,:,jk) * vmask(:,:,jk) 144 END DO 145 ENDIF 146 ! 147 ! !== Apply semi-implicit bottom friction ==! 148 ! 149 ! Only needed for semi-implicit bottom friction setup. The explicit 150 ! bottom friction has been included in "u(v)a" which act as the R.H.S 151 ! column vector of the tri-diagonal matrix equation 152 ! 153 IF( ln_drgimp ) THEN 154 DO jj = 2, jpjm1 155 DO ji = 2, jpim1 156 ikbu = mbku(ji,jj) ! ocean bottom level at u- and v-points 157 ikbv = mbkv(ji,jj) ! (deepest ocean u- and v-points) 158 avmu(ji,jj,ikbu+1) = -0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) * e3uw_n(ji,jj,ikbu+1) 159 avmv(ji,jj,ikbv+1) = -0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) * e3vw_n(ji,jj,ikbv+1) 160 END DO 161 END DO 162 IF ( ln_isfcav ) THEN 163 DO jj = 2, jpjm1 164 DO ji = 2, jpim1 165 ikbu = miku(ji,jj) ! ocean top level at u- and v-points 166 ikbv = mikv(ji,jj) ! (first wet ocean u- and v-points) 167 ! top Cd is masked (=0 outside cavities) no need of test on mik>=2 168 avmu(ji,jj,ikbu) = -0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * e3uw_n(ji,jj,ikbu) 169 avmv(ji,jj,ikbv) = -0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * e3vw_n(ji,jj,ikbv) 170 END DO 171 END DO 172 END IF 173 ENDIF 174 ! 175 ! With split-explicit free surface, barotropic stress is treated explicitly 176 ! Update velocities at the bottom. 177 ! J. Chanut: The bottom stress is computed considering after barotropic velocities, which does 178 ! not lead to the effective stress seen over the whole barotropic loop. 179 ! G. Madec : in linear free surface, e3u_a = e3u_n = e3u_0, so systematic use of e3u_a 113 & + r2dt * e3v_n(:,:,jk) * va(:,:,jk) ) / e3v_a(:,:,jk) * vmask(:,:,jk) 114 END DO 115 ENDIF 116 ! ! add top/bottom friction 117 ! With split-explicit free surface, barotropic stress is treated explicitly Update velocities at the bottom. 118 ! J. Chanut: The bottom stress is computed considering after barotropic velocities, which does 119 ! not lead to the effective stress seen over the whole barotropic loop. 120 ! G. Madec : in linear free surface, e3u_a = e3u_n = e3u_0, so systematic use of e3u_a 180 121 IF( ln_drgimp .AND. ln_dynspg_ts ) THEN 181 122 DO jk = 1, jpkm1 ! remove barotropic velocities … … 185 126 DO jj = 2, jpjm1 ! Add bottom/top stress due to barotropic component only 186 127 DO ji = fs_2, fs_jpim1 ! vector opt. 187 ik bu = mbku(ji,jj) ! ocean bottom level at u- and v-points188 ik bv = mbkv(ji,jj) ! (deepest ocean u- and v-points)189 ze3ua = ( 1._wp - r_vvl ) * e3u_n(ji,jj,ik bu) + r_vvl * e3u_a(ji,jj,ikbu)190 ze3va = ( 1._wp - r_vvl ) * e3v_n(ji,jj,ik bv) + r_vvl * e3v_a(ji,jj,ikbv)191 ua(ji,jj,ik bu) = ua(ji,jj,ikbu) + p2dt * 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) * ua_b(ji,jj) / ze3ua192 va(ji,jj,ik bv) = va(ji,jj,ikbv) + p2dt * 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) * va_b(ji,jj) / ze3va128 iku = mbku(ji,jj) ! ocean bottom level at u- and v-points 129 ikv = mbkv(ji,jj) ! (deepest ocean u- and v-points) 130 ze3ua = ( 1._wp - r_vvl ) * e3u_n(ji,jj,iku) + r_vvl * e3u_a(ji,jj,iku) 131 ze3va = ( 1._wp - r_vvl ) * e3v_n(ji,jj,ikv) + r_vvl * e3v_a(ji,jj,ikv) 132 ua(ji,jj,iku) = ua(ji,jj,iku) + r2dt * 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) * ua_b(ji,jj) / ze3ua 133 va(ji,jj,ikv) = va(ji,jj,ikv) + r2dt * 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) * va_b(ji,jj) / ze3va 193 134 END DO 194 135 END DO … … 196 137 DO jj = 2, jpjm1 197 138 DO ji = fs_2, fs_jpim1 ! vector opt. 198 ik bu = miku(ji,jj) ! top ocean level at u- and v-points199 ik bv = mikv(ji,jj) ! (first wet ocean u- and v-points)200 ze3ua = ( 1._wp - r_vvl ) * e3u_n(ji,jj,ik bu) + r_vvl * e3u_a(ji,jj,ikbu)201 ze3va = ( 1._wp - r_vvl ) * e3v_n(ji,jj,ik bv) + r_vvl * e3v_a(ji,jj,ikbv)202 ua(ji,jj,ik bu) = ua(ji,jj,ikbu) + p2dt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * ua_b(ji,jj) / ze3ua203 va(ji,jj,ik bv) = va(ji,jj,ikbv) + p2dt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * va_b(ji,jj) / ze3va139 iku = miku(ji,jj) ! top ocean level at u- and v-points 140 ikv = mikv(ji,jj) ! (first wet ocean u- and v-points) 141 ze3ua = ( 1._wp - r_vvl ) * e3u_n(ji,jj,iku) + r_vvl * e3u_a(ji,jj,iku) 142 ze3va = ( 1._wp - r_vvl ) * e3v_n(ji,jj,ikv) + r_vvl * e3v_a(ji,jj,ikv) 143 ua(ji,jj,iku) = ua(ji,jj,iku) + r2dt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * ua_b(ji,jj) / ze3ua 144 va(ji,jj,ikv) = va(ji,jj,ikv) + r2dt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * va_b(ji,jj) / ze3va 204 145 END DO 205 146 END DO … … 209 150 ! !== Vertical diffusion on u ==! 210 151 ! 211 ! Matrix and second member construction 212 ! bottom boundary condition: both zwi and zws must be masked as avmu can take 213 ! non zero value at the ocean bottom depending on the bottom friction used. 214 ! 215 DO jk = 1, jpkm1 ! Matrix 216 DO jj = 2, jpjm1 217 DO ji = fs_2, fs_jpim1 ! vector opt. 218 ze3ua = ( 1._wp - r_vvl ) * e3u_n(ji,jj,jk) + r_vvl * e3u_a(ji,jj,jk) ! after scale factor at T-point 219 zzwi = - p2dt * avmu(ji,jj,jk ) / ( ze3ua * e3uw_n(ji,jj,jk ) ) 220 zzws = - p2dt * avmu(ji,jj,jk+1) / ( ze3ua * e3uw_n(ji,jj,jk+1) ) 221 zwi(ji,jj,jk) = zzwi * wumask(ji,jj,jk ) 222 zws(ji,jj,jk) = zzws * wumask(ji,jj,jk+1) 223 zwd(ji,jj,jk) = 1._wp - zzwi - zzws 224 END DO 225 END DO 226 END DO 227 DO jj = 2, jpjm1 ! Surface boundary conditions 152 ! !* Matrix construction 153 zdt = r2dt * 0.5 154 IF( nldf == np_lap_i ) THEN ! rotated lateral mixing: add its vertical mixing (akzu) 155 DO jk = 1, jpkm1 156 DO jj = 2, jpjm1 157 DO ji = fs_2, fs_jpim1 ! vector opt. 158 ze3ua = ( 1._wp - r_vvl ) * e3u_n(ji,jj,jk) + r_vvl * e3u_a(ji,jj,jk) ! after scale factor at T-point 159 zzwi = - zdt * ( avm(ji+1,jj,jk ) + avm(ji,jj,jk ) + akzu(ji,jj,jk ) ) & 160 & / ( ze3ua * e3uw_n(ji,jj,jk ) ) * wumask(ji,jj,jk ) 161 zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) + akzu(ji,jj,jk+1) ) & 162 & / ( ze3ua * e3uw_n(ji,jj,jk+1) ) * wumask(ji,jj,jk+1) 163 zwi(ji,jj,jk) = zzwi 164 zws(ji,jj,jk) = zzws 165 zwd(ji,jj,jk) = 1._wp - zzwi - zzws 166 END DO 167 END DO 168 END DO 169 ELSE ! standard case 170 DO jk = 1, jpkm1 171 DO jj = 2, jpjm1 172 DO ji = fs_2, fs_jpim1 ! vector opt. 173 ze3ua = ( 1._wp - r_vvl ) * e3u_n(ji,jj,jk) + r_vvl * e3u_a(ji,jj,jk) ! after scale factor at T-point 174 zzwi = - zdt * ( avm(ji+1,jj,jk ) + avm(ji,jj,jk ) ) / ( ze3ua * e3uw_n(ji,jj,jk ) ) * wumask(ji,jj,jk ) 175 zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) ) / ( ze3ua * e3uw_n(ji,jj,jk+1) ) * wumask(ji,jj,jk+1) 176 zwi(ji,jj,jk) = zzwi 177 zws(ji,jj,jk) = zzws 178 zwd(ji,jj,jk) = 1._wp - zzwi - zzws 179 END DO 180 END DO 181 END DO 182 ENDIF 183 ! 184 DO jj = 2, jpjm1 !* Surface boundary conditions 228 185 DO ji = fs_2, fs_jpim1 ! vector opt. 229 186 zwi(ji,jj,1) = 0._wp … … 231 188 END DO 232 189 END DO 233 190 ! 191 ! !== Apply semi-implicit bottom friction ==! 192 ! 193 ! Only needed for semi-implicit bottom friction setup. The explicit 194 ! bottom friction has been included in "u(v)a" which act as the R.H.S 195 ! column vector of the tri-diagonal matrix equation 196 ! 197 IF ( ln_drgimp ) THEN ! implicit bottom friction 198 DO jj = 2, jpjm1 199 DO ji = 2, jpim1 200 iku = mbku(ji,jj) ! ocean bottom level at u- and v-points 201 ze3ua = ( 1._wp - r_vvl ) * e3u_n(ji,jj,iku) + r_vvl * e3u_a(ji,jj,iku) ! after scale factor at T-point 202 zwd(ji,jj,iku) = zwd(ji,jj,iku) - r2dt * 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) / ze3ua 203 END DO 204 END DO 205 IF ( ln_isfcav ) THEN ! top friction (always implicit) 206 DO jj = 2, jpjm1 207 DO ji = 2, jpim1 208 !!gm top Cd is masked (=0 outside cavities) no need of test on mik>=2 ==>> it has been suppressed 209 iku = miku(ji,jj) ! ocean top level at u- and v-points 210 ze3ua = ( 1._wp - r_vvl ) * e3u_n(ji,jj,iku) + r_vvl * e3u_a(ji,jj,iku) ! after scale factor at T-point 211 zwd(ji,jj,iku) = zwd(ji,jj,iku) - r2dt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) / ze3ua 212 END DO 213 END DO 214 END IF 215 ENDIF 216 ! 234 217 ! Matrix inversion starting from the first level 235 218 !----------------------------------------------------------------------- … … 258 241 DO ji = fs_2, fs_jpim1 ! vector opt. 259 242 ze3ua = ( 1._wp - r_vvl ) * e3u_n(ji,jj,1) + r_vvl * e3u_a(ji,jj,1) 260 ua(ji,jj,1) = ua(ji,jj,1) + p2dt * 0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) ) &243 ua(ji,jj,1) = ua(ji,jj,1) + r2dt * 0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) ) & 261 244 & / ( ze3ua * rau0 ) * umask(ji,jj,1) 262 245 END DO … … 285 268 ! !== Vertical diffusion on v ==! 286 269 ! 287 ! Matrix and second member construction 288 ! bottom boundary condition: both zwi and zws must be masked as avmv can take 289 ! non zero value at the ocean bottom depending on the bottom friction used 290 ! 291 DO jk = 1, jpkm1 ! Matrix 292 DO jj = 2, jpjm1 293 DO ji = fs_2, fs_jpim1 ! vector opt. 294 ze3va = ( 1._wp - r_vvl ) * e3v_n(ji,jj,jk) + r_vvl * e3v_a(ji,jj,jk) ! after scale factor at T-point 295 zzwi = - p2dt * avmv (ji,jj,jk ) / ( ze3va * e3vw_n(ji,jj,jk ) ) 296 zzws = - p2dt * avmv (ji,jj,jk+1) / ( ze3va * e3vw_n(ji,jj,jk+1) ) 297 zwi(ji,jj,jk) = zzwi * wvmask(ji,jj,jk ) 298 zws(ji,jj,jk) = zzws * wvmask(ji,jj,jk+1) 299 zwd(ji,jj,jk) = 1._wp - zzwi - zzws 300 END DO 301 END DO 302 END DO 303 DO jj = 2, jpjm1 ! Surface boundary conditions 270 ! !* Matrix construction 271 zdt = r2dt * 0.5 272 IF( nldf == np_lap_i ) THEN ! rotated lateral mixing: add its vertical mixing (akzu) 273 DO jk = 1, jpkm1 274 DO jj = 2, jpjm1 275 DO ji = fs_2, fs_jpim1 ! vector opt. 276 ze3va = ( 1._wp - r_vvl ) * e3v_n(ji,jj,jk) + r_vvl * e3v_a(ji,jj,jk) ! after scale factor at T-point 277 zzwi = - zdt * ( avm(ji,jj+1,jk )+ avm(ji,jj,jk ) + akzv(ji,jj,jk ) ) & 278 & / ( ze3va * e3vw_n(ji,jj,jk ) ) * wvmask(ji,jj,jk ) 279 zzws = - zdt * ( avm(ji,jj+1,jk+1)+ avm(ji,jj,jk+1) + akzv(ji,jj,jk+1) ) & 280 & / ( ze3va * e3vw_n(ji,jj,jk+1) ) * wvmask(ji,jj,jk+1) 281 zwi(ji,jj,jk) = zzwi * wvmask(ji,jj,jk ) 282 zws(ji,jj,jk) = zzws * wvmask(ji,jj,jk+1) 283 zwd(ji,jj,jk) = 1._wp - zzwi - zzws 284 END DO 285 END DO 286 END DO 287 ELSE ! standard case 288 DO jk = 1, jpkm1 289 DO jj = 2, jpjm1 290 DO ji = fs_2, fs_jpim1 ! vector opt. 291 ze3va = ( 1._wp - r_vvl ) * e3v_n(ji,jj,jk) + r_vvl * e3v_a(ji,jj,jk) ! after scale factor at T-point 292 zzwi = - zdt * ( avm(ji,jj+1,jk )+ avm(ji,jj,jk ) ) / ( ze3va * e3vw_n(ji,jj,jk ) ) * wvmask(ji,jj,jk ) 293 zzws = - zdt * ( avm(ji,jj+1,jk+1)+ avm(ji,jj,jk+1) ) / ( ze3va * e3vw_n(ji,jj,jk+1) ) * wvmask(ji,jj,jk+1) 294 zwi(ji,jj,jk) = zzwi * wvmask(ji,jj,jk ) 295 zws(ji,jj,jk) = zzws * wvmask(ji,jj,jk+1) 296 zwd(ji,jj,jk) = 1._wp - zzwi - zzws 297 END DO 298 END DO 299 END DO 300 ENDIF 301 ! 302 DO jj = 2, jpjm1 !* Surface boundary conditions 304 303 DO ji = fs_2, fs_jpim1 ! vector opt. 305 304 zwi(ji,jj,1) = 0._wp … … 307 306 END DO 308 307 END DO 308 ! !== Apply semi-implicit top/bottom friction ==! 309 ! 310 ! Only needed for semi-implicit bottom friction setup. The explicit 311 ! bottom friction has been included in "u(v)a" which act as the R.H.S 312 ! column vector of the tri-diagonal matrix equation 313 ! 314 IF( ln_drgimp ) THEN 315 DO jj = 2, jpjm1 316 DO ji = 2, jpim1 317 ikv = mbkv(ji,jj) ! (deepest ocean u- and v-points) 318 ze3va = ( 1._wp - r_vvl ) * e3v_n(ji,jj,ikv) + r_vvl * e3v_a(ji,jj,ikv) ! after scale factor at T-point 319 zwd(ji,jj,ikv) = zwd(ji,jj,ikv) - r2dt * 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) / ze3va 320 END DO 321 END DO 322 IF ( ln_isfcav ) THEN 323 DO jj = 2, jpjm1 324 DO ji = 2, jpim1 325 ikv = mikv(ji,jj) ! (first wet ocean u- and v-points) 326 ze3va = ( 1._wp - r_vvl ) * e3v_n(ji,jj,ikv) + r_vvl * e3v_a(ji,jj,ikv) ! after scale factor at T-point 327 zwd(ji,jj,iku) = zwd(ji,jj,iku) - r2dt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) / ze3va 328 END DO 329 END DO 330 ENDIF 331 ENDIF 309 332 310 333 ! Matrix inversion … … 334 357 DO ji = fs_2, fs_jpim1 ! vector opt. 335 358 ze3va = ( 1._wp - r_vvl ) * e3v_n(ji,jj,1) + r_vvl * e3v_a(ji,jj,1) 336 va(ji,jj,1) = va(ji,jj,1) + p2dt * 0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) ) &359 va(ji,jj,1) = va(ji,jj,1) + r2dt * 0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) ) & 337 360 & / ( ze3va * rau0 ) * vmask(ji,jj,1) 338 361 END DO … … 359 382 END DO 360 383 ! 361 IF( nn_timing == 1 ) CALL timing_stop('dyn_zdf_imp') 362 ! 363 END SUBROUTINE dyn_zdf_imp 384 IF( l_trddyn ) THEN ! save the vertical diffusive trends for further diagnostics 385 ztrdu(:,:,:) = ( ua(:,:,:) - ub(:,:,:) ) / r2dt - ztrdu(:,:,:) 386 ztrdv(:,:,:) = ( va(:,:,:) - vb(:,:,:) ) / r2dt - ztrdv(:,:,:) 387 CALL trd_dyn( ztrdu, ztrdv, jpdyn_zdf, kt ) 388 DEALLOCATE( ztrdu, ztrdv ) 389 ENDIF 390 ! ! print mean trends (used for debugging) 391 IF(ln_ctl) CALL prt_ctl( tab3d_1=ua, clinfo1=' zdf - Ua: ', mask1=umask, & 392 & tab3d_2=va, clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' ) 393 ! 394 IF( nn_timing == 1 ) CALL timing_stop('dyn_zdf') 395 ! 396 END SUBROUTINE dyn_zdf 364 397 365 398 !!============================================================================== -
branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdf_oce.F90
r8160 r8178 6 6 !! history : 1.0 ! 2002-06 (G. Madec) Original code 7 7 !! 3.2 ! 2009-07 (G. Madec) addition of avm 8 !! 4.0 ! 2017-05 (G. Madec) drag coef. defined at t-point (zdfdrg.F90)8 !! 4.0 ! 2017-05 (G. Madec) avm and drag coef. defined at t-point 9 9 !!---------------------------------------------------------------------- 10 10 USE par_oce ! ocean parameters … … 44 44 45 45 46 REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:) :: avm, avt, avs !: vertical mixing coefficient (w-point) [m2/s] 47 !!gm 48 REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:) :: avmu , avmv !: vertical viscosity coef at uw- & vw-pts [m2/s] 49 !!gm 50 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: avm_k, avt_k !: avm, avt computed by turbulent closure alone 51 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: en !: now turbulent kinetic energy [m2/s2] 52 REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:) :: avtb_2d !: horizontal shape of background Kz profile 53 REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:) :: avmb , avtb !: background profile of avm and avt 46 REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:) :: avm, avt, avs !: vertical mixing coefficients (w-point) [m2/s] 47 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: avm_k , avt_k !: Kz computed by turbulent closure alone [m2/s] 48 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: en !: now turbulent kinetic energy [m2/s2] 49 REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:) :: avmb , avtb !: background profile of avm and avt [m2/s] 50 REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:) :: avtb_2d !: horizontal shape of background Kz profile [-] 54 51 55 52 !!---------------------------------------------------------------------- … … 65 62 !!---------------------------------------------------------------------- 66 63 ! 67 ALLOCATE( avm (jpi,jpj,jpk) , avt (jpi,jpj,jpk) , avs(jpi,jpj,jpk) , & 68 & avm_k(jpi,jpj,jpk) , avt_k(jpi,jpj,jpk) , en (jpi,jpj,jpk) , & 69 & avmb(jpk) , avtb(jpk) , avtb_2d(jpi,jpj) , & 70 & avmu(jpi,jpj,jpk), avmv(jpi,jpj,jpk) , STAT = zdf_oce_alloc ) 64 ALLOCATE( avm (jpi,jpj,jpk) , avm_k(jpi,jpj,jpk) , avs(jpi,jpj,jpk) , & 65 & avt (jpi,jpj,jpk) , avt_k(jpi,jpj,jpk) , en (jpi,jpj,jpk) , & 66 & avmb(jpk) , avtb(jpk) , avtb_2d(jpi,jpj) , STAT = zdf_oce_alloc ) 71 67 ! 72 68 IF( zdf_oce_alloc /= 0 ) CALL ctl_warn('zdf_oce_alloc: failed to allocate arrays') -
branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/step.F90
r8160 r8178 29 29 !! 3.7 ! 2014-10 (G. Madec) LDF simplication 30 30 !! - ! 2014-12 (G. Madec) remove KPP scheme 31 !! - ! 2015-11 (J. Chanut) free surface simplification 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) 32 33 !!---------------------------------------------------------------------- 33 34 … … 36 37 !!---------------------------------------------------------------------- 37 38 USE step_oce ! time stepping definition modules 38 !!gm to be removed when removing avmu, avmv39 USE zdf_oce ! ocean vertical physics variables40 !!gm41 39 ! 42 40 USE iom ! xIOs server … … 48 46 49 47 !!---------------------------------------------------------------------- 50 !! NEMO/OPA 3.7 , NEMO Consortium (2015)48 !! NEMO/OPA 4.0 , NEMO Consortium (2017) 51 49 !! $Id$ 52 50 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) … … 131 129 CALL zdf_phy( kstp ) ! vertical physics update (top/bot drag, avt, avs, avm + MLD) 132 130 133 134 !!gm ===>>> TO BE REMOVED (require changes in zdf_oce, dynzdf(_imp,_exp), dynldf_iso, diawri)135 DO jk = 1, jpkm1 !* vertical eddy viscosity at wu- and wv-points136 DO jj = 2, jpjm1137 DO ji = 2, jpim1138 avmu(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk) + avm(ji+1,jj ,jk) ) * wumask(ji,jj,jk)139 avmv(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk) + avm(ji ,jj+1,jk) ) * wvmask(ji,jj,jk)140 END DO141 END DO142 END DO143 CALL lbc_lnk( avmu, 'U', 1. ) ; CALL lbc_lnk( avmv, 'V', 1. ) ! Lateral boundary conditions144 !!gm end145 146 147 131 ! LATERAL PHYSICS 148 132 !
Note: See TracChangeset
for help on using the changeset viewer.