- Timestamp:
- 2017-12-01T18:44:09+01:00 (6 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf_iso.F90
r6140 r8882 28 28 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 29 29 USE prtctl ! Print control 30 USE wrk_nemo ! Memory Allocation31 30 USE timing ! Timing 32 31 … … 37 36 PUBLIC dyn_ldf_iso_alloc ! called by nemogcm.F90 38 37 38 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: akzu, akzv !: vertical component of rotated lateral viscosity 39 39 40 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: zfuw, zdiu, zdju, zdj1u ! 2D workspace (dyn_ldf_iso) 40 41 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: zfvw, zdiv, zdjv, zdj1v ! - - … … 43 44 # include "vectopt_loop_substitute.h90" 44 45 !!---------------------------------------------------------------------- 45 !! NEMO/OPA 3.3 , NEMO Consortium (2011)46 !! NEMO/OPA 4.0 , NEMO Consortium (2017) 46 47 !! $Id$ 47 48 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) … … 53 54 !! *** ROUTINE dyn_ldf_iso_alloc *** 54 55 !!---------------------------------------------------------------------- 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 )56 ALLOCATE( akzu(jpi,jpj,jpk) , zfuw(jpi,jpk) , zdiu(jpi,jpk) , zdju(jpi,jpk) , zdj1u(jpi,jpk) , & 57 & akzv(jpi,jpj,jpk) , zfvw(jpi,jpk) , zdiv(jpi,jpk) , zdjv(jpi,jpk) , zdj1v(jpi,jpk) , STAT=dyn_ldf_iso_alloc ) 57 58 ! 58 59 IF( dyn_ldf_iso_alloc /= 0 ) CALL ctl_warn('dyn_ldf_iso_alloc: array allocate failed.') … … 99 100 !! 100 101 !! ** 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 102 !! -(ua,va) updated with the before geopotential harmonic mixing trend 103 !! -(akzu,akzv) to accompt for the diagonal vertical component 104 !! of the rotated operator in dynzdf module 105 105 !!---------------------------------------------------------------------- 106 106 INTEGER, INTENT( in ) :: kt ! ocean time-step index 107 107 ! 108 108 INTEGER :: ji, jj, jk ! dummy loop indices 109 REAL(wp) :: zabe1, zabe2, zcof1, zcof2 ! local scalars 110 REAL(wp) :: zmskt, zmskf ! - - 111 REAL(wp) :: zcoef0, zcoef3, zcoef4, zmkt, zmkf ! - - 112 REAL(wp) :: zuav, zvav, zuwslpi, zuwslpj, zvwslpi, zvwslpj ! - - 113 ! 114 REAL(wp), POINTER, DIMENSION(:,:) :: ziut, zjuf, zjvt, zivf, zdku, zdk1u, zdkv, zdk1v 109 REAL(wp) :: zabe1, zmskt, zmkt, zuav, zuwslpi, zuwslpj ! local scalars 110 REAL(wp) :: zabe2, zmskf, zmkf, zvav, zvwslpi, zvwslpj ! - - 111 REAL(wp) :: zcof0, zcof1, zcof2, zcof3, zcof4 ! - - 112 REAL(wp), DIMENSION(jpi,jpj) :: ziut, zivf, zdku, zdk1u ! 2D workspace 113 REAL(wp), DIMENSION(jpi,jpj) :: zjuf, zjvt, zdkv, zdk1v ! - - 115 114 !!---------------------------------------------------------------------- 116 115 ! 117 IF( nn_timing == 1 ) CALL timing_start('dyn_ldf_iso') 118 ! 119 CALL wrk_alloc( jpi, jpj, ziut, zjuf, zjvt, zivf, zdku, zdk1u, zdkv, zdk1v ) 116 IF( ln_timing ) CALL timing_start('dyn_ldf_iso') 120 117 ! 121 118 IF( kt == nit000 ) THEN … … 144 141 CALL lbc_lnk( uslp , 'U', -1. ) ; CALL lbc_lnk( vslp , 'V', -1. ) 145 142 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 143 ! 144 ENDIF 154 145 155 146 ! ! =============== … … 348 339 DO jk = 2, jpkm1 349 340 DO ji = 2, jpim1 350 zco ef0= 0.5* rn_aht_0 * umask(ji,jj,jk)341 zcof0 = 0.5_wp * rn_aht_0 * umask(ji,jj,jk) 351 342 ! 352 zuwslpi = zco ef0 * ( wslpi(ji+1,jj,jk) + wslpi(ji,jj,jk) )353 zuwslpj = zco ef0 * ( wslpj(ji+1,jj,jk) + wslpj(ji,jj,jk) )343 zuwslpi = zcof0 * ( wslpi(ji+1,jj,jk) + wslpi(ji,jj,jk) ) 344 zuwslpj = zcof0 * ( wslpj(ji+1,jj,jk) + wslpj(ji,jj,jk) ) 354 345 ! 355 zmkt = 1./MAX( tmask(ji,jj,jk-1)+tmask(ji+1,jj,jk-1) &356 + tmask(ji,jj,jk )+tmask(ji+1,jj,jk ) , 1. )357 zmkf = 1./MAX( fmask(ji,jj-1,jk-1) + fmask(ji,jj,jk-1) &358 + fmask(ji,jj-1,jk ) + fmask(ji,jj,jk ) , 1. )359 360 zco ef3 = - e2u(ji,jj) * zmkt * zuwslpi361 zco ef4 = - e1u(ji,jj) * zmkf * zuwslpj346 zmkt = 1./MAX( tmask(ji,jj,jk-1)+tmask(ji+1,jj,jk-1) & 347 + tmask(ji,jj,jk )+tmask(ji+1,jj,jk ) , 1. ) 348 zmkf = 1./MAX( fmask(ji,jj-1,jk-1) + fmask(ji,jj,jk-1) & 349 + fmask(ji,jj-1,jk ) + fmask(ji,jj,jk ) , 1. ) 350 351 zcof3 = - e2u(ji,jj) * zmkt * zuwslpi 352 zcof4 = - e1u(ji,jj) * zmkf * zuwslpj 362 353 ! vertical flux on u field 363 zfuw(ji,jk) = zco ef3 * ( zdiu (ji,jk-1) + zdiu (ji+1,jk-1)&364 +zdiu (ji,jk ) + zdiu (ji+1,jk )) &365 + zcoef4 * ( zdj1u(ji,jk-1) + zdju (ji ,jk-1)&366 +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_0354 zfuw(ji,jk) = zcof3 * ( zdiu (ji,jk-1) + zdiu (ji+1,jk-1) & 355 & + zdiu (ji,jk ) + zdiu (ji+1,jk ) ) & 356 & + zcof4 * ( zdj1u(ji,jk-1) + zdju (ji ,jk-1) & 357 & + zdj1u(ji,jk ) + zdju (ji ,jk ) ) 358 ! vertical mixing coefficient (akzu) 359 ! Note: zcof0 include rn_aht_0, so divided by rn_aht_0 to obtain slp^2 * rn_aht_0 360 akzu(ji,jj,jk) = ( zuwslpi * zuwslpi + zuwslpj * zuwslpj ) / rn_aht_0 370 361 END DO 371 362 END DO … … 374 365 DO jk = 2, jpkm1 375 366 DO ji = 2, jpim1 376 zco ef0 = 0.5* rn_aht_0 * vmask(ji,jj,jk)377 378 zvwslpi = zco ef0 * ( wslpi(ji,jj+1,jk) + wslpi(ji,jj,jk) )379 zvwslpj = zco ef0 * ( wslpj(ji,jj+1,jk) + wslpj(ji,jj,jk) )380 381 zmkf = 1./MAX( fmask(ji-1,jj,jk-1)+fmask(ji,jj,jk-1) &382 + fmask(ji-1,jj,jk )+fmask(ji,jj,jk ), 1. )383 zmkt = 1./MAX( tmask(ji,jj,jk-1)+tmask(ji,jj+1,jk-1) &384 + tmask(ji,jj,jk )+tmask(ji,jj+1,jk ), 1. )385 386 zco ef3 = - e2v(ji,jj) * zmkf * zvwslpi387 zco ef4 = - e1v(ji,jj) * zmkt * zvwslpj367 zcof0 = 0.5_wp * rn_aht_0 * vmask(ji,jj,jk) 368 ! 369 zvwslpi = zcof0 * ( wslpi(ji,jj+1,jk) + wslpi(ji,jj,jk) ) 370 zvwslpj = zcof0 * ( wslpj(ji,jj+1,jk) + wslpj(ji,jj,jk) ) 371 ! 372 zmkf = 1./MAX( fmask(ji-1,jj,jk-1)+fmask(ji,jj,jk-1) & 373 & + fmask(ji-1,jj,jk )+fmask(ji,jj,jk ) , 1. ) 374 zmkt = 1./MAX( tmask(ji,jj,jk-1)+tmask(ji,jj+1,jk-1) & 375 & + tmask(ji,jj,jk )+tmask(ji,jj+1,jk ) , 1. ) 376 377 zcof3 = - e2v(ji,jj) * zmkf * zvwslpi 378 zcof4 = - e1v(ji,jj) * zmkt * zvwslpj 388 379 ! vertical flux on v field 389 zfvw(ji,jk) = zco ef3 * ( zdiv (ji,jk-1) + zdiv (ji-1,jk-1)&390 & +zdiv (ji,jk ) + zdiv (ji-1,jk )) &391 & + zco ef4 * ( zdjv (ji,jk-1) + zdj1v(ji ,jk-1)&392 & +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_0380 zfvw(ji,jk) = zcof3 * ( zdiv (ji,jk-1) + zdiv (ji-1,jk-1) & 381 & + zdiv (ji,jk ) + zdiv (ji-1,jk ) ) & 382 & + zcof4 * ( zdjv (ji,jk-1) + zdj1v(ji ,jk-1) & 383 & + zdjv (ji,jk ) + zdj1v(ji ,jk ) ) 384 ! vertical mixing coefficient (akzv) 385 ! Note: zcof0 include rn_aht_0, so divided by rn_aht_0 to obtain slp^2 * rn_aht_0 386 akzv(ji,jj,jk) = ( zvwslpi * zvwslpi + zvwslpj * zvwslpj ) / rn_aht_0 396 387 END DO 397 388 END DO … … 409 400 END DO ! End of slab 410 401 ! ! =============== 411 CALL wrk_dealloc( jpi, jpj, ziut, zjuf, zjvt, zivf, zdku, zdk1u, zdkv, zdk1v )412 402 ! 413 IF( nn_timing == 1 )CALL timing_stop('dyn_ldf_iso')403 IF( ln_timing ) CALL timing_stop('dyn_ldf_iso') 414 404 ! 415 405 END SUBROUTINE dyn_ldf_iso
Note: See TracChangeset
for help on using the changeset viewer.
