Changeset 7646 for trunk/NEMOGCM/NEMO/OPA_SRC/LDF
- Timestamp:
- 2017-02-06T10:25:03+01:00 (7 years ago)
- Location:
- trunk/NEMOGCM/NEMO/OPA_SRC/LDF
- Files:
-
- 3 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMOGCM/NEMO/OPA_SRC/LDF/ldfc1d_c2d.F90
r6140 r7646 59 59 REAL(wp) :: zw , zdep2 ! - - 60 60 !!---------------------------------------------------------------------- 61 62 IF(lwp) THEN 63 WRITE(numout,*) 64 WRITE(numout,*) ' ldf_c1d : set a given profile to eddy diffusivity/viscosity coefficients' 65 WRITE(numout,*) ' ~~~~~~~' 66 ENDIF 61 67 62 68 ! initialization of the profile … … 130 136 ! 131 137 IF(lwp) THEN 132 WRITE(numout,*) 'ldf_c2d : aht = rn_aht0 * max(e1,e2)/e_equator ( laplacian) ' 133 WRITE(numout,*) '~~~~~~~ or = rn_bht0 * [max(e1,e2)/e_equator]**3 (bilaplacian)' 138 WRITE(numout,*) 139 WRITE(numout,*) ' ldf_c2d : aht = rn_aht0 * max(e1,e2)/e_equator ( laplacian) ' 140 WRITE(numout,*) ' ~~~~~~~ or = rn_bht0 * [max(e1,e2)/e_equator]**3 (bilaplacian)' 134 141 WRITE(numout,*) 135 142 ENDIF -
trunk/NEMOGCM/NEMO/OPA_SRC/LDF/ldfdyn.F90
r6140 r7646 42 42 REAL(wp), PUBLIC :: rn_ahm_b !: lateral laplacian background eddy viscosity [m2/s] 43 43 REAL(wp), PUBLIC :: rn_bhm_0 !: lateral bilaplacian eddy viscosity [m4/s] 44 !! If nn_ahm_ijk_t = 32 a time and space varying Smagorinsky viscosity 45 !! will be computed. 46 REAL(wp), PUBLIC :: rn_csmc !: Smagorinsky constant of proportionality 47 REAL(wp), PUBLIC :: rn_minfac !: Multiplicative factor of theorectical minimum Smagorinsky viscosity 48 REAL(wp), PUBLIC :: rn_maxfac !: Multiplicative factor of theorectical maximum Smagorinsky viscosity 44 49 45 50 LOGICAL , PUBLIC :: l_ldfdyn_time !: flag for time variation of the lateral eddy viscosity coef. 46 51 47 52 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: ahmt, ahmf !: eddy diffusivity coef. at U- and V-points [m2/s or m4/s] 53 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: dtensq !: horizontal tension squared (Smagorinsky only) 54 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: dshesq !: horizontal shearing strain squared (Smagorinsky only) 55 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: esqt, esqf !: Square of the local gridscale (e1e2/(e1+e2))**2 48 56 49 57 REAL(wp) :: r1_12 = 1._wp / 12._wp ! =1/12 50 58 REAL(wp) :: r1_4 = 0.25_wp ! =1/4 59 REAL(wp) :: r1_8 = 0.125_wp ! =1/8 51 60 REAL(wp) :: r1_288 = 1._wp / 288._wp ! =1/( 12^2 * 2 ) 52 61 … … 79 88 !! = 31 = F(i,j,k,t) = F(local velocity) ( |u|e /12 laplacian operator 80 89 !! or |u|e^3/12 bilaplacian operator ) 81 !!---------------------------------------------------------------------- 82 INTEGER :: jk ! dummy loop indices 90 !! = 32 = F(i,j,k,t) = F(local deformation rate and gridscale) (D and L) (Smagorinsky) 91 !! ( L^2|D| laplacian operator 92 !! or L^4|D|/8 bilaplacian operator ) 93 !!---------------------------------------------------------------------- 94 INTEGER :: ji, jj, jk ! dummy loop indices 83 95 INTEGER :: ierr, inum, ios ! local integer 84 96 REAL(wp) :: zah0 ! local scalar … … 86 98 NAMELIST/namdyn_ldf/ ln_dynldf_lap, ln_dynldf_blp, & 87 99 & ln_dynldf_lev, ln_dynldf_hor, ln_dynldf_iso, & 88 & nn_ahm_ijk_t , rn_ahm_0, rn_ahm_b, rn_bhm_0 100 & nn_ahm_ijk_t , rn_ahm_0, rn_ahm_b, rn_bhm_0, & 101 & rn_csmc , rn_minfac, rn_maxfac 89 102 !!---------------------------------------------------------------------- 90 103 ! … … 115 128 WRITE(numout,*) ' coefficients :' 116 129 WRITE(numout,*) ' type of time-space variation nn_ahm_ijk_t = ', nn_ahm_ijk_t 117 WRITE(numout,*) ' lateral laplacian eddy viscosity rn_ahm_0 _lap= ', rn_ahm_0, ' m2/s'130 WRITE(numout,*) ' lateral laplacian eddy viscosity rn_ahm_0 = ', rn_ahm_0, ' m2/s' 118 131 WRITE(numout,*) ' background viscosity (iso case) rn_ahm_b = ', rn_ahm_b, ' m2/s' 119 WRITE(numout,*) ' lateral bilaplacian eddy viscosity rn_ahm_0_blp = ', rn_bhm_0, ' m4/s' 132 WRITE(numout,*) ' lateral bilaplacian eddy viscosity rn_bhm_0 = ', rn_bhm_0, ' m4/s' 133 WRITE(numout,*) ' smagorinsky settings (nn_ahm_ijk_t = 32) :' 134 WRITE(numout,*) ' Smagorinsky coefficient rn_csmc = ', rn_csmc 135 WRITE(numout,*) ' factor multiplier for theorectical lower limit for ' 136 WRITE(numout,*) ' Smagorinsky eddy visc (def. 1.0) rn_minfac = ', rn_minfac 137 WRITE(numout,*) ' factor multiplier for theorectical lower upper for ' 138 WRITE(numout,*) ' Smagorinsky eddy visc (def. 1.0) rn_maxfac = ', rn_maxfac 120 139 ENDIF 121 140 … … 208 227 l_ldfdyn_time = .TRUE. ! will be calculated by call to ldf_dyn routine in step.F90 209 228 ! 229 CASE( 32 ) !== time varying 3D field ==! 230 IF(lwp) WRITE(numout,*) ' momentum mixing coef. = F( latitude, longitude, depth , time )' 231 IF(lwp) WRITE(numout,*) ' proportional to the local deformation rate and gridscale (Smagorinsky)' 232 IF(lwp) WRITE(numout,*) ' : L^2|D| or L^4|D|/8' 233 ! 234 l_ldfdyn_time = .TRUE. ! will be calculated by call to ldf_dyn routine in step.F90 235 ! 236 ! allocate arrays used in ldf_dyn. 237 ALLOCATE( dtensq(jpi,jpj) , dshesq(jpi,jpj) , esqt(jpi,jpj) , esqf(jpi,jpj) , STAT=ierr ) 238 IF( ierr /= 0 ) CALL ctl_stop( 'STOP', 'ldf_dyn_init: failed to allocate Smagorinsky arrays') 239 ! 240 ! Set local gridscale values 241 DO jj = 2, jpjm1 242 DO ji = fs_2, fs_jpim1 243 esqt(ji,jj) = ( e1e2t(ji,jj) /( e1t(ji,jj) + e2t(ji,jj) ) )**2 244 esqf(ji,jj) = ( e1e2f(ji,jj) /( e1f(ji,jj) + e2f(ji,jj) ) )**2 245 END DO 246 END DO 247 ! 210 248 CASE DEFAULT 211 249 CALL ctl_stop('ldf_dyn_init: wrong choice for nn_ahm_ijk_t, the type of space-time variation of ahm') … … 232 270 !! nn_ahm_ijk_t = 31 ahmt, ahmf = F(i,j,k,t) = F(local velocity) 233 271 !! ( |u|e /12 or |u|e^3/12 for laplacian or bilaplacian operator ) 234 !! BLP case : sqrt of the eddy coef, since bilaplacian is en re-entrant laplacian235 272 !! 236 !! ** action : ahmt, ahmf update at each time step 273 !! nn_ahm_ijk_t = 32 ahmt, ahmf = F(i,j,k,t) = F(local deformation rate and gridscale) (D and L) (Smagorinsky) 274 !! ( L^2|D| or L^4|D|/8 for laplacian or bilaplacian operator ) 275 !! 276 !! ** note : in BLP cases the sqrt of the eddy coef is returned, since bilaplacian is en re-entrant laplacian 277 !! ** action : ahmt, ahmf updated at each time step 237 278 !!---------------------------------------------------------------------- 238 279 INTEGER, INTENT(in) :: kt ! time step index … … 240 281 INTEGER :: ji, jj, jk ! dummy loop indices 241 282 REAL(wp) :: zu2pv2_ij_p1, zu2pv2_ij, zu2pv2_ij_m1, zetmax, zefmax ! local scalar 283 REAL(wp) :: zcmsmag, zstabf_lo, zstabf_up, zdelta, zdb ! local scalar 242 284 !!---------------------------------------------------------------------- 243 285 ! … … 248 290 CASE( 31 ) !== time varying 3D field ==! = F( local velocity ) 249 291 ! 250 IF( ln_dynldf_lap ) THEN !laplacian operator : |u| e /12 = |u/144| e292 IF( ln_dynldf_lap ) THEN ! laplacian operator : |u| e /12 = |u/144| e 251 293 DO jk = 1, jpkm1 252 294 DO jj = 2, jpjm1 … … 280 322 CALL lbc_lnk( ahmt, 'T', 1. ) ; CALL lbc_lnk( ahmf, 'F', 1. ) 281 323 ! 324 ! 325 CASE( 32 ) !== time varying 3D field ==! = F( local deformation rate and gridscale ) (Smagorinsky) 326 ! 327 IF( ln_dynldf_lap .OR. ln_dynldf_blp ) THEN ! laplacian operator : (C_smag/pi)^2 L^2 |D| 328 ! 329 zcmsmag = (rn_csmc/rpi)**2 ! (C_smag/pi)^2 330 zstabf_lo = rn_minfac * rn_minfac / ( 2._wp * 4._wp * zcmsmag ) ! lower limit stability factor scaling 331 zstabf_up = rn_maxfac / ( 4._wp * zcmsmag * 2._wp * rdt ) ! upper limit stability factor scaling 332 IF( ln_dynldf_blp ) zstabf_lo = ( 16._wp / 9._wp ) * zstabf_lo ! provide |U|L^3/12 lower limit instead 333 ! ! of |U|L^3/16 in blp case 334 DO jk = 1, jpkm1 335 ! 336 DO jj = 2, jpj 337 DO ji = 2, jpi 338 zdb = ( ( ub(ji,jj,jk) * r1_e2u(ji,jj) - ub(ji-1,jj,jk) * r1_e2u(ji-1,jj) ) & 339 & * r1_e1t(ji,jj) * e2t(ji,jj) & 340 & - ( vb(ji,jj,jk) * r1_e1v(ji,jj) - vb(ji,jj-1,jk) * r1_e1v(ji,jj-1) ) & 341 & * r1_e2t(ji,jj) * e1t(ji,jj) ) * tmask(ji,jj,jk) 342 dtensq(ji,jj) = zdb*zdb 343 END DO 344 END DO 345 ! 346 DO jj = 1, jpjm1 347 DO ji = 1, jpim1 348 zdb = ( ( ub(ji,jj+1,jk) * r1_e1u(ji,jj+1) - ub(ji,jj,jk) * r1_e1u(ji,jj) ) & 349 & * r1_e2f(ji,jj) * e1f(ji,jj) & 350 & + ( vb(ji+1,jj,jk) * r1_e2v(ji+1,jj) - vb(ji,jj,jk) * r1_e2v(ji,jj) ) & 351 & * r1_e1f(ji,jj) * e2f(ji,jj) ) * fmask(ji,jj,jk) 352 dshesq(ji,jj) = zdb*zdb 353 END DO 354 END DO 355 ! 356 DO jj = 2, jpjm1 357 DO ji = fs_2, fs_jpim1 358 ! 359 zu2pv2_ij_p1 = ub(ji ,jj+1,jk) * ub(ji ,jj+1,jk) + vb(ji+1,jj ,jk) * vb(ji+1,jj ,jk) 360 zu2pv2_ij = ub(ji ,jj ,jk) * ub(ji ,jj ,jk) + vb(ji ,jj ,jk) * vb(ji ,jj ,jk) 361 zu2pv2_ij_m1 = ub(ji-1,jj ,jk) * ub(ji-1,jj ,jk) + vb(ji ,jj-1,jk) * vb(ji ,jj-1,jk) 362 ! T-point value 363 zdelta = zcmsmag * esqt(ji,jj) ! L^2 * (C_smag/pi)^2 364 ahmt(ji,jj,jk) = zdelta * sqrt( dtensq(ji,jj) + & 365 & r1_4 * ( dshesq(ji,jj) + dshesq(ji,jj-1) + & 366 & dshesq(ji-1,jj) + dshesq(ji-1,jj-1) ) ) 367 ahmt(ji,jj,jk) = MAX( ahmt(ji,jj,jk), & 368 & SQRT( (zu2pv2_ij + zu2pv2_ij_m1) * zdelta * zstabf_lo ) ) ! Impose lower limit == minfac * |U|L/2 369 ahmt(ji,jj,jk) = MIN( ahmt(ji,jj,jk), zdelta * zstabf_up ) ! Impose upper limit == maxfac * L^2/(4*2dt) 370 ! F-point value 371 zdelta = zcmsmag * esqf(ji,jj) ! L^2 * (C_smag/pi)^2 372 ahmf(ji,jj,jk) = zdelta * sqrt( dshesq(ji,jj) + & 373 & r1_4 * ( dtensq(ji,jj) + dtensq(ji,jj+1) + & 374 & dtensq(ji+1,jj) + dtensq(ji+1,jj+1) ) ) 375 ahmf(ji,jj,jk) = MAX( ahmf(ji,jj,jk), & 376 & SQRT( (zu2pv2_ij + zu2pv2_ij_p1) * zdelta * zstabf_lo ) ) ! Impose lower limit == minfac * |U|L/2 377 ahmf(ji,jj,jk) = MIN( ahmf(ji,jj,jk), zdelta * zstabf_up ) ! Impose upper limit == maxfac * L^2/(4*2dt) 378 ! 379 END DO 380 END DO 381 END DO 382 ! 383 ENDIF 384 ! 385 IF( ln_dynldf_blp ) THEN ! bilaplacian operator : sqrt( (C_smag/pi)^2 L^4 |D|/8) 386 ! = sqrt( A_lap_smag L^2/8 ) 387 ! stability limits already applied to laplacian values 388 ! effective default limits are |U|L^3/12 < B_hm < L^4/(32*2dt) 389 ! 390 DO jk = 1, jpkm1 391 ! 392 DO jj = 2, jpjm1 393 DO ji = fs_2, fs_jpim1 394 ! 395 ahmt(ji,jj,jk) = sqrt( r1_8 * esqt(ji,jj) * ahmt(ji,jj,jk) ) 396 ahmf(ji,jj,jk) = sqrt( r1_8 * esqf(ji,jj) * ahmf(ji,jj,jk) ) 397 ! 398 END DO 399 END DO 400 END DO 401 ! 402 ENDIF 403 ! 404 CALL lbc_lnk( ahmt, 'T', 1. ) ; CALL lbc_lnk( ahmf, 'F', 1. ) 405 ! 282 406 END SELECT 283 407 ! -
trunk/NEMOGCM/NEMO/OPA_SRC/LDF/ldftra.F90
r6140 r7646 24 24 USE ldfslp ! lateral diffusion: slope of iso-neutral surfaces 25 25 USE ldfc1d_c2d ! lateral diffusion: 1D & 2D cases 26 USE dia ar5, ONLY: lk_diaar526 USE diaptr 27 27 ! 28 USE trc_oce, ONLY: lk_offline ! offline flag29 28 USE in_out_manager ! I/O manager 30 29 USE iom ! I/O module for ehanced bottom friction file … … 298 297 ! 299 298 INTEGER :: ji, jj, jk ! dummy loop indices 300 REAL(wp) :: zaht, zah t_min, z1_f20 ! local scalar301 !!---------------------------------------------------------------------- 302 ! 303 IF( nn_aei_ijk_t == 21 ) THEN ! eddy induced velocity coefficients299 REAL(wp) :: zaht, zahf, zaht_min, z1_f20 ! local scalar 300 !!---------------------------------------------------------------------- 301 ! 302 IF( ln_ldfeiv .AND. nn_aei_ijk_t == 21 ) THEN ! eddy induced velocity coefficients 304 303 ! ! =F(growth rate of baroclinic instability) 305 304 ! ! max value rn_aeiv_0 ; decreased to 0 within 20N-20S 306 305 CALL ldf_eiv( kt, rn_aeiv_0, aeiu, aeiv ) 307 IF(lwp .AND. kt<=nit000+20 ) WRITE(numout,*) ' kt , ldf_eiv appel', kt308 306 ENDIF 309 307 ! … … 314 312 ! ! max value rn_aht_0 (rn_aeiv_0 if nn_aei_ijk_t=21) 315 313 ! ! increase to rn_aht_0 within 20N-20S 316 IF( nn_aei_ijk_t /= 21 ) THEN 317 CALL ldf_eiv( kt, rn_aht_0, ahtu, ahtv ) 318 IF(lwp .AND. kt<=nit000+20 ) WRITE(numout,*) ' kt , ldf_eiv appel 2', kt 319 ELSE 314 IF( ln_ldfeiv .AND. nn_aei_ijk_t == 21 ) THEN ! use the already computed aei. 320 315 ahtu(:,:,1) = aeiu(:,:,1) 321 316 ahtv(:,:,1) = aeiv(:,:,1) 322 IF(lwp .AND. kt<=nit000+20 ) WRITE(numout,*) ' kt , ahtu=aeiu', kt 317 ELSE ! compute aht. 318 CALL ldf_eiv( kt, rn_aht_0, ahtu, ahtv ) 323 319 ENDIF 324 320 ! … … 327 323 DO jj = 1, jpj 328 324 DO ji = 1, jpi 329 zaht = ( 1._wp - MIN( 1._wp , ABS( ff(ji,jj) * z1_f20 ) ) ) * ( rn_aht_0 - zaht_min ) 325 !!gm CAUTION : here we assume lat/lon grid in 20deg N/S band (like all ORCA cfg) 326 !! ==>>> The Coriolis value is identical for t- & u_points, and for v- and f-points 327 zaht = ( 1._wp - MIN( 1._wp , ABS( ff_t(ji,jj) * z1_f20 ) ) ) * ( rn_aht_0 - zaht_min ) 328 zahf = ( 1._wp - MIN( 1._wp , ABS( ff_f(ji,jj) * z1_f20 ) ) ) * ( rn_aht_0 - zaht_min ) 330 329 ahtu(ji,jj,1) = ( MAX( zaht_min, ahtu(ji,jj,1) ) + zaht ) * umask(ji,jj,1) ! min value zaht_min 331 ahtv(ji,jj,1) = ( MAX( zaht_min, ahtv(ji,jj,1) ) + zah t) * vmask(ji,jj,1) ! increase within 20S-20N330 ahtv(ji,jj,1) = ( MAX( zaht_min, ahtv(ji,jj,1) ) + zahf ) * vmask(ji,jj,1) ! increase within 20S-20N 332 331 END DO 333 332 END DO … … 352 351 END SELECT 353 352 ! 354 IF( .NOT.lk_offline ) THEN 355 CALL iom_put( "ahtu_2d", ahtu(:,:,1) ) ! surface u-eddy diffusivity coeff. 356 CALL iom_put( "ahtv_2d", ahtv(:,:,1) ) ! surface v-eddy diffusivity coeff. 357 CALL iom_put( "ahtu_3d", ahtu(:,:,:) ) ! 3D u-eddy diffusivity coeff. 358 CALL iom_put( "ahtv_3d", ahtv(:,:,:) ) ! 3D v-eddy diffusivity coeff. 359 ! 353 CALL iom_put( "ahtu_2d", ahtu(:,:,1) ) ! surface u-eddy diffusivity coeff. 354 CALL iom_put( "ahtv_2d", ahtv(:,:,1) ) ! surface v-eddy diffusivity coeff. 355 CALL iom_put( "ahtu_3d", ahtu(:,:,:) ) ! 3D u-eddy diffusivity coeff. 356 CALL iom_put( "ahtv_3d", ahtv(:,:,:) ) ! 3D v-eddy diffusivity coeff. 357 ! 360 358 !!gm : THE IF below is to be checked (comes from Seb) 361 IF( ln_ldfeiv ) THEN 362 CALL iom_put( "aeiu_2d", aeiu(:,:,1) ) ! surface u-EIV coeff. 363 CALL iom_put( "aeiv_2d", aeiv(:,:,1) ) ! surface v-EIV coeff. 364 CALL iom_put( "aeiu_3d", aeiu(:,:,:) ) ! 3D u-EIV coeff. 365 CALL iom_put( "aeiv_3d", aeiv(:,:,:) ) ! 3D v-EIV coeff. 366 ENDIF 359 IF( ln_ldfeiv ) THEN 360 CALL iom_put( "aeiu_2d", aeiu(:,:,1) ) ! surface u-EIV coeff. 361 CALL iom_put( "aeiv_2d", aeiv(:,:,1) ) ! surface v-EIV coeff. 362 CALL iom_put( "aeiu_3d", aeiu(:,:,:) ) ! 3D u-EIV coeff. 363 CALL iom_put( "aeiv_3d", aeiv(:,:,:) ) ! 3D v-EIV coeff. 367 364 ENDIF 368 365 ! … … 555 552 END DO 556 553 557 !!gm IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN ! ORCA R2558 !!gm DO jj = 2, jpjm1559 !!gm DO ji = fs_2, fs_jpim1 ! vector opt.560 !!gm ! Take the minimum between aeiw and 1000 m2/s over shelves (depth shallower than 650 m)561 !!gm IF( mbkt(ji,jj) <= 20 ) zaeiw(ji,jj) = MIN( zaeiw(ji,jj), 1000. )562 !!gm END DO563 !!gm END DO564 !!gm ENDIF565 566 554 ! !== Bound on eiv coeff. ==! 567 555 z1_f20 = 1._wp / ( 2._wp * omega * sin( rad * 20._wp ) ) 568 556 DO jj = 2, jpjm1 569 557 DO ji = fs_2, fs_jpim1 ! vector opt. 570 zzaei = MIN( 1._wp, ABS( ff (ji,jj) * z1_f20 ) ) * zaeiw(ji,jj) ! tropical decrease558 zzaei = MIN( 1._wp, ABS( ff_t(ji,jj) * z1_f20 ) ) * zaeiw(ji,jj) ! tropical decrease 571 559 zaeiw(ji,jj) = MIN( zzaei , paei0 ) ! Max value = paei0 572 560 END DO … … 730 718 CALL iom_put( "woce_eiv", zw3d ) 731 719 ! 720 ! 721 ! 722 CALL wrk_alloc( jpi,jpj, zw2d ) 723 ! 724 zztmp = 0.5_wp * rau0 * rcp 725 IF( iom_use('ueiv_heattr') .OR. iom_use('ueiv_heattr3d') ) THEN 726 zw2d(:,:) = 0._wp 727 zw3d(:,:,:) = 0._wp 728 DO jk = 1, jpkm1 729 DO jj = 2, jpjm1 730 DO ji = fs_2, fs_jpim1 ! vector opt. 731 zw3d(ji,jj,jk) = zw3d(ji,jj,jk) + ( psi_uw(ji,jj,jk+1) - psi_uw(ji,jj,jk) ) & 732 & * ( tsn (ji,jj,jk,jp_tem) + tsn (ji+1,jj,jk,jp_tem) ) 733 zw2d(ji,jj) = zw2d(ji,jj) + zw3d(ji,jj,jk) 734 END DO 735 END DO 736 END DO 737 CALL lbc_lnk( zw2d, 'U', -1. ) 738 CALL lbc_lnk( zw3d, 'U', -1. ) 739 CALL iom_put( "ueiv_heattr" , zztmp * zw2d ) ! heat transport in i-direction 740 CALL iom_put( "ueiv_heattr3d", zztmp * zw3d ) ! heat transport in i-direction 741 ENDIF 742 zw2d(:,:) = 0._wp 743 zw3d(:,:,:) = 0._wp 744 DO jk = 1, jpkm1 745 DO jj = 2, jpjm1 746 DO ji = fs_2, fs_jpim1 ! vector opt. 747 zw3d(ji,jj,jk) = zw3d(ji,jj,jk) + ( psi_vw(ji,jj,jk+1) - psi_vw(ji,jj,jk) ) & 748 & * ( tsn (ji,jj,jk,jp_tem) + tsn (ji,jj+1,jk,jp_tem) ) 749 zw2d(ji,jj) = zw2d(ji,jj) + zw3d(ji,jj,jk) 750 END DO 751 END DO 752 END DO 753 CALL lbc_lnk( zw2d, 'V', -1. ) 754 CALL iom_put( "veiv_heattr", zztmp * zw2d ) ! heat transport in j-direction 755 CALL iom_put( "veiv_heattr", zztmp * zw3d ) ! heat transport in j-direction 756 ! 757 IF( ln_diaptr ) CALL dia_ptr_hst( jp_tem, 'eiv', 0.5 * zw3d ) 758 ! 759 zztmp = 0.5_wp * 0.5 760 IF( iom_use('ueiv_salttr') .OR. iom_use('ueiv_salttr3d')) THEN 761 zw2d(:,:) = 0._wp 762 zw3d(:,:,:) = 0._wp 763 DO jk = 1, jpkm1 764 DO jj = 2, jpjm1 765 DO ji = fs_2, fs_jpim1 ! vector opt. 766 zw3d(ji,jj,jk) = zw3d(ji,jj,jk) * ( psi_uw(ji,jj,jk+1) - psi_uw(ji,jj,jk) ) & 767 & * ( tsn (ji,jj,jk,jp_sal) + tsn (ji+1,jj,jk,jp_sal) ) 768 zw2d(ji,jj) = zw2d(ji,jj) + zw3d(ji,jj,jk) 769 END DO 770 END DO 771 END DO 772 CALL lbc_lnk( zw2d, 'U', -1. ) 773 CALL lbc_lnk( zw3d, 'U', -1. ) 774 CALL iom_put( "ueiv_salttr", zztmp * zw2d ) ! salt transport in i-direction 775 CALL iom_put( "ueiv_salttr3d", zztmp * zw3d ) ! salt transport in i-direction 776 ENDIF 777 zw2d(:,:) = 0._wp 778 zw3d(:,:,:) = 0._wp 779 DO jk = 1, jpkm1 780 DO jj = 2, jpjm1 781 DO ji = fs_2, fs_jpim1 ! vector opt. 782 zw3d(ji,jj,jk) = zw3d(ji,jj,jk) + ( psi_vw(ji,jj,jk+1) - psi_vw(ji,jj,jk) ) & 783 & * ( tsn (ji,jj,jk,jp_sal) + tsn (ji,jj+1,jk,jp_sal) ) 784 zw2d(ji,jj) = zw2d(ji,jj) + zw3d(ji,jj,jk) 785 END DO 786 END DO 787 END DO 788 CALL lbc_lnk( zw2d, 'V', -1. ) 789 CALL iom_put( "veiv_salttr", zztmp * zw2d ) ! salt transport in j-direction 790 CALL iom_put( "veiv_salttr", zztmp * zw3d ) ! salt transport in j-direction 791 ! 792 IF( ln_diaptr ) CALL dia_ptr_hst( jp_sal, 'eiv', 0.5 * zw3d ) 793 ! 794 CALL wrk_dealloc( jpi,jpj, zw2d ) 732 795 CALL wrk_dealloc( jpi,jpj,jpk, zw3d ) 733 !734 !735 IF( lk_diaar5 ) THEN !== eiv heat transport: calculate and output ==!736 CALL wrk_alloc( jpi,jpj, zw2d )737 !738 zztmp = 0.5_wp * rau0 * rcp739 zw2d(:,:) = 0._wp740 DO jk = 1, jpkm1741 DO jj = 2, jpjm1742 DO ji = fs_2, fs_jpim1 ! vector opt.743 zw2d(ji,jj) = zw2d(ji,jj) + zztmp * ( psi_uw(ji,jj,jk+1) - psi_uw(ji,jj,jk) ) &744 & * ( tsn (ji,jj,jk,jp_tem) + tsn (ji+1,jj,jk,jp_tem) )745 END DO746 END DO747 END DO748 CALL lbc_lnk( zw2d, 'U', -1. )749 CALL iom_put( "ueiv_heattr", zw2d ) ! heat transport in i-direction750 zw2d(:,:) = 0._wp751 DO jk = 1, jpkm1752 DO jj = 2, jpjm1753 DO ji = fs_2, fs_jpim1 ! vector opt.754 zw2d(ji,jj) = zw2d(ji,jj) + zztmp * ( psi_vw(ji,jj,jk+1) - psi_vw(ji,jj,jk) ) &755 & * ( tsn (ji,jj,jk,jp_tem) + tsn (ji,jj+1,jk,jp_tem) )756 END DO757 END DO758 END DO759 CALL lbc_lnk( zw2d, 'V', -1. )760 CALL iom_put( "veiv_heattr", zw2d ) ! heat transport in i-direction761 !762 CALL wrk_dealloc( jpi,jpj, zw2d )763 ENDIF764 796 ! 765 797 IF( nn_timing == 1 ) CALL timing_stop( 'ldf_eiv_dia')
Note: See TracChangeset
for help on using the changeset viewer.