Changeset 2926 for branches/2011/dev_r2802_NOCL_Smagorinsky/NEMOGCM/NEMO
- Timestamp:
- 2011-10-14T17:18:45+02:00 (13 years ago)
- Location:
- branches/2011/dev_r2802_NOCL_Smagorinsky/NEMOGCM/NEMO/OPA_SRC
- Files:
-
- 2 deleted
- 12 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2011/dev_r2802_NOCL_Smagorinsky/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf_bilap.F90
r2887 r2926 4 4 !! Ocean dynamics: lateral viscosity trend 5 5 !!====================================================================== 6 !! History : OPA ! 1990-09 (G. Madec) Original code 7 !! 4.0 ! 1993-03 (M. Guyon) symetrical conditions (M. Guyon) 8 !! 6.0 ! 1996-01 (G. Madec) statement function for e3 9 !! 8.0 ! 1997-07 (G. Madec) lbc calls 10 !! NEMO 1.0 ! 2002-08 (G. Madec) F90: Free form and module 11 !! 2.0 ! 2004-08 (C. Talandier) New trends organization 12 !!---------------------------------------------------------------------- 6 13 7 14 !!---------------------------------------------------------------------- … … 9 16 !! using an iso-level bilaplacian operator 10 17 !!---------------------------------------------------------------------- 11 !! * Modules used12 18 USE oce ! ocean dynamics and tracers 13 19 USE dom_oce ! ocean space and time domain … … 21 27 PRIVATE 22 28 23 !! * Routine accessibility 24 PUBLIC dyn_ldf_bilap ! called by step.F90 29 PUBLIC dyn_ldf_bilap ! called by step.F90 25 30 26 31 !! * Substitutions … … 31 36 !! NEMO/OPA 3.3 , NEMO Consortium (2010) 32 37 !! $Id$ 33 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 34 !!---------------------------------------------------------------------- 35 38 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 39 !!---------------------------------------------------------------------- 36 40 CONTAINS 37 41 … … 69 73 !! ** Action : - Update (ua,va) with the before iso-level biharmonic 70 74 !! mixing trend. 71 !!72 !! History :73 !! ! 90-09 (G. Madec) Original code74 !! ! 91-11 (G. Madec)75 !! ! 93-03 (M. Guyon) symetrical conditions (M. Guyon)76 !! ! 96-01 (G. Madec) statement function for e377 !! ! 97-07 (G. Madec) lbc calls78 !! 8.5 ! 02-08 (G. Madec) F90: Free form and module79 !! 9.0 ! 04-08 (C. Talandier) New trends organization80 75 !!---------------------------------------------------------------------- 81 !! * Arguments 82 INTEGER, INTENT( in ) :: kt ! ocean time-step index 83 84 !! * Local declarations 85 INTEGER :: ji, jj, jk ! dummy loop indices 86 REAL(wp) :: zua, zva, zbt, ze2u, ze2v ! temporary scalar 87 REAL(wp), DIMENSION(jpi,jpj) :: & 88 zcu, zcv ! temporary workspace 89 REAL(wp), DIMENSION(jpi,jpj,jpk) :: & 90 zuf, zut, zlu, zlv ! temporary workspace 76 USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released 77 USE wrk_nemo, ONLY: zcu => wrk_2d_1 , zcv => wrk_2d_2 ! 3D workspace 78 USE wrk_nemo, ONLY: zuf => wrk_3d_3 , zut => wrk_3d_4 ! 3D workspace 79 USE wrk_nemo, ONLY: zlu => wrk_3d_5 , zlv => wrk_3d_6 80 ! 81 INTEGER, INTENT(in) :: kt ! ocean time-step index 82 ! 83 INTEGER :: ji, jj, jk ! dummy loop indices 84 REAL(wp) :: zua, zva, zbt, ze2u, ze2v ! temporary scalar 91 85 !!---------------------------------------------------------------------- 92 !! OPA 8.5, LODYC-IPSL (2002) 93 !!---------------------------------------------------------------------- 94 95 IF( kt == nit000 ) THEN 96 IF(lwp) WRITE(numout,*) 97 IF(lwp) WRITE(numout,*) 'dyn_ldf_bilap : iso-level bilaplacian operator' 98 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~' 86 87 IF( wrk_in_use(2, 1,2) .OR. wrk_in_use(3, 3,4,5,6) ) THEN 88 CALL ctl_stop('dyn_ldf_bilap : requested workspace arrays unavailable') ; RETURN 89 ENDIF 90 91 IF( kt == nit000 .AND. lwp ) THEN 92 WRITE(numout,*) 93 WRITE(numout,*) 'dyn_ldf_bilap : iso-level bilaplacian operator' 94 WRITE(numout,*) '~~~~~~~~~~~~~' 99 95 ENDIF 100 96 … … 104 100 !!$ zlu(:,:,jpk) = 0.e0 105 101 !!$ zlv(:,:,jpk) = 0.e0 106 zuf(:,:,:) = 0. e0107 zut(:,:,:) = 0. e0108 zlu(:,:,:) = 0. e0109 zlv(:,:,:) = 0. e0102 zuf(:,:,:) = 0._wp 103 zut(:,:,:) = 0._wp 104 zlu(:,:,:) = 0._wp 105 zlv(:,:,:) = 0._wp 110 106 111 107 ! ! =============== … … 137 133 END DO 138 134 ENDIF 139 ENDDO 140 141 ! Boundary conditions on the laplacian (zlu,zlv) 142 CALL lbc_lnk( zlu, 'U', -1. ) 143 CALL lbc_lnk( zlv, 'V', -1. ) 144 135 END DO 136 CALL lbc_lnk( zlu, 'U', -1. ) ; CALL lbc_lnk( zlv, 'V', -1. ) ! Boundary conditions 137 145 138 146 139 DO jk = 1, jpkm1 … … 150 143 151 144 ! Multiply by the eddy viscosity coef. (at u- and v-points) 152 #if ! defined key_dynldf_smag153 145 zlu(:,:,jk) = zlu(:,:,jk) * fsahmu(:,:,jk) 154 146 zlv(:,:,jk) = zlv(:,:,jk) * fsahmv(:,:,jk) 155 #endif147 156 148 ! Contravariant "laplacian" 157 149 zcu(:,:) = e1u(:,:) * zlu(:,:,jk) … … 195 187 ! Bilaplacian 196 188 ! ----------- 197 #if !defined key_dynldf_smag 189 198 190 DO jj = 2, jpjm1 199 191 DO ji = fs_2, fs_jpim1 ! vector opt. … … 211 203 END DO 212 204 END DO 213 #else 214 215 DO jj = 2, jpjm1 216 DO ji = fs_2, fs_jpim1 ! vector opt. 217 ze2u = e2u(ji,jj) * fse3u(ji,jj,jk) 218 ze2v = e1v(ji,jj) * fse3v(ji,jj,jk) 219 ! horizontal biharmonic diffusive trends 220 zua = - ( zuf(ji ,jj,jk) - zuf(ji,jj-1,jk) ) / ze2u & 221 & + ( zut(ji+1,jj,jk) - zut(ji,jj ,jk) ) / e1u(ji,jj) 222 223 zva = + ( zuf(ji,jj ,jk) - zuf(ji-1,jj,jk) ) / ze2v & 224 & + ( zut(ji,jj+1,jk) - zut(ji ,jj,jk) ) / e2v(ji,jj) 225 ! add it to the general momentum trends 226 ua(ji,jj,jk) = ua(ji,jj,jk) + zua*fsahmu(ji,jj,jk) 227 va(ji,jj,jk) = va(ji,jj,jk) + zva*fsahmv(ji,jj,jk) 228 END DO 229 END DO 230 231 #endif 205 232 206 ! ! =============== 233 207 END DO ! End of slab 234 208 ! ! =============== 235 209 IF( wrk_not_released(2, 1,2) .OR. & 210 wrk_not_released(3, 3,4,5,6) ) CALL ctl_stop('dyn_ldf_bilap: failed to release workspace arrays') 211 ! 236 212 END SUBROUTINE dyn_ldf_bilap 237 213 -
branches/2011/dev_r2802_NOCL_Smagorinsky/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf_bilapg.F90
r2887 r2926 4 4 !! Ocean dynamics: lateral viscosity trend 5 5 !!====================================================================== 6 !! History : OPA ! 1997-07 (G. Madec) Original code 7 !! NEMO 1.0 ! 2002-08 (G. Madec) F90: Free form and module 8 !! 2.0 ! 2004-08 (C. Talandier) New trends organization 9 !!---------------------------------------------------------------------- 6 10 #if defined key_ldfslp || defined key_esopa 7 11 !!---------------------------------------------------------------------- … … 12 16 !! ldfguv : 13 17 !!---------------------------------------------------------------------- 14 !! * Modules used15 18 USE oce ! ocean dynamics and tracers 16 19 USE dom_oce ! ocean space and time domain 17 20 USE ldfdyn_oce ! ocean dynamics lateral physics 18 21 USE zdf_oce ! ocean vertical physics 19 USE in_out_manager ! I/O manager20 22 USE trdmod ! ocean dynamics trends 21 23 USE trdmod_oce ! ocean variables trends 22 24 USE ldfslp ! iso-neutral slopes available 25 USE in_out_manager ! I/O manager 26 USE lib_mpp ! MPP library 23 27 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 24 28 USE prtctl ! Print control … … 27 31 PRIVATE 28 32 29 !! * Routine accessibility 30 PUBLIC dyn_ldf_bilapg ! called by step.F90 33 PUBLIC dyn_ldf_bilapg ! called by step.F90 34 35 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: zfuw, zfvw , zdiu, zdiv ! 2D workspace (ldfguv) 36 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: zdju, zdj1u, zdjv, zdj1v ! 2D workspace (ldfguv) 31 37 32 38 !! * Substitutions … … 36 42 !! NEMO/OPA 3.3 , NEMO Consortium (2010) 37 43 !! $Id$ 38 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 39 !!---------------------------------------------------------------------- 40 44 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 45 !!---------------------------------------------------------------------- 41 46 CONTAINS 47 48 INTEGER FUNCTION dyn_ldf_bilapg_alloc() 49 !!---------------------------------------------------------------------- 50 !! *** ROUTINE dyn_ldf_bilapg_alloc *** 51 !!---------------------------------------------------------------------- 52 ALLOCATE( zfuw(jpi,jpk) , zfvw (jpi,jpk) , zdiu(jpi,jpk) , zdiv (jpi,jpk) , & 53 & zdju(jpi,jpk) , zdj1u(jpi,jpk) , zdjv(jpi,jpk) , zdj1v(jpi,jpk) , STAT=dyn_ldf_bilapg_alloc ) 54 ! 55 IF( dyn_ldf_bilapg_alloc /= 0 ) CALL ctl_warn('dyn_ldf_bilapg_alloc: failed to allocate arrays') 56 END FUNCTION dyn_ldf_bilapg_alloc 57 42 58 43 59 SUBROUTINE dyn_ldf_bilapg( kt ) … … 67 83 !! biharmonic mixing trend. 68 84 !! - save the trend in (zwk3,zwk4) ('key_trddyn') 69 !! 70 !! History : 71 !! 8.0 ! 97-07 (G. Madec) Original code 72 !! 8.5 ! 02-08 (G. Madec) F90: Free form and module 73 !! 9.0 ! 04-08 (C. Talandier) New trends organization 74 !!---------------------------------------------------------------------- 75 !! * Modules used 76 USE oce, ONLY : zwk3 => ta, & ! use ta as 3D workspace 77 zwk4 => sa ! use sa as 3D workspace 78 79 !! * Arguments 85 !!---------------------------------------------------------------------- 86 USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released 87 USE wrk_nemo, ONLY: zwk1 => wrk_3d_3 , zwk2 => wrk_3d_4 ! 3D workspace 88 USE oce , ONLY: zwk3 => ta , zwk4 => sa ! ta, sa used as 3D workspace 89 ! 80 90 INTEGER, INTENT( in ) :: kt ! ocean time-step index 81 82 !! * Local declarations 91 ! 83 92 INTEGER :: ji, jj, jk ! dummy loop indices 84 REAL(wp), DIMENSION(jpi,jpj,jpk) :: & 85 zwk1, zwk2 ! work array used for rotated biharmonic 86 ! ! operator on tracers and/or momentum 87 !!---------------------------------------------------------------------- 93 !!---------------------------------------------------------------------- 94 95 IF( wrk_in_use(3, 3,4) ) THEN 96 CALL ctl_stop('dyn_ldf_bilapg: requested workspace arrays unavailable') ; RETURN 97 ENDIF 88 98 89 99 IF( kt == nit000 ) THEN … … 93 103 zwk1(:,:,:) = 0.e0 ; zwk3(:,:,:) = 0.e0 94 104 zwk2(:,:,:) = 0.e0 ; zwk4(:,:,:) = 0.e0 105 ! ! allocate dyn_ldf_bilapg arrays 106 IF( dyn_ldf_bilapg_alloc() /= 0 ) CALL ctl_stop('STOP', 'dyn_ldf_bilapg: failed to allocate arrays') 95 107 ENDIF 96 108 97 109 ! Laplacian of (ub,vb) multiplied by ahm 98 110 ! -------------------------------------- 99 ! rotated harmonic operator applied to (ub,vb) 100 ! and multiply by ahmu, ahmv (output in (zwk1,zwk2) ) 101 102 CALL ldfguv ( ub, vb, zwk1, zwk2, 1 ) 103 104 105 ! Lateral boundary conditions on (zwk1,zwk2) 106 CALL lbc_lnk( zwk1, 'U', -1. ) 107 CALL lbc_lnk( zwk2, 'V', -1. ) 108 111 CALL ldfguv( ub, vb, zwk1, zwk2, 1 ) ! rotated harmonic operator applied to (ub,vb) 112 ! ! and multiply by ahmu, ahmv (output in (zwk1,zwk2) ) 113 CALL lbc_lnk( zwk1, 'U', -1. ) ; CALL lbc_lnk( zwk2, 'V', -1. ) ! Lateral boundary conditions 109 114 110 115 ! Bilaplacian of (ub,vb) 111 116 ! ---------------------- 112 ! rotated harmonic operator applied to (zwk1,zwk2) (output in (zwk3,zwk4) ) 113 114 CALL ldfguv ( zwk1, zwk2, zwk3, zwk4, 2 ) 115 116 117 ! Update the momentum trends (j-slab : 2, jpj-1) 117 CALL ldfguv( zwk1, zwk2, zwk3, zwk4, 2 ) ! rotated harmonic operator applied to (zwk1,zwk2) 118 ! ! (output in (zwk3,zwk4) ) 119 120 ! Update the momentum trends 118 121 ! -------------------------- 119 ! ! =============== 120 DO jj = 2, jpjm1 ! Vertical slab 121 ! ! =============== 122 DO jj = 2, jpjm1 ! add the diffusive trend to the general momentum trends 122 123 DO jk = 1, jpkm1 123 124 DO ji = 2, jpim1 124 ! add the diffusive trend to the general momentum trends125 125 ua(ji,jj,jk) = ua(ji,jj,jk) + zwk3(ji,jj,jk) 126 126 va(ji,jj,jk) = va(ji,jj,jk) + zwk4(ji,jj,jk) 127 127 END DO 128 128 END DO 129 ! ! ===============130 END DO ! End of slab131 ! ! ===============132 129 END DO 130 ! 131 IF( wrk_not_released(3, 3,4) ) CALL ctl_stop('dyn_ldf_bilapg: failed to release workspace arrays') 132 ! 133 133 END SUBROUTINE dyn_ldf_bilapg 134 134 … … 174 174 !! second order vertical derivative term) 175 175 !! 'key_trddyn' defined: the trend is saved for diagnostics. 176 !! 177 !! History : 178 !! 8.0 ! 97-07 (G. Madec) Original code 179 !! 8.5 ! 02-08 (G. Madec) F90: Free form and module 180 !!---------------------------------------------------------------------- 181 !! * Arguments 182 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( in ) :: & 183 pu, pv ! momentum fields (before u and v for the 1st call, and 184 ! ! laplacian of these fields multiplied by ahm for the 2nd 185 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( out ) :: & 186 plu, plv ! partial harmonic operator applied to 187 ! ! pu and pv (all the components except 188 ! ! second order vertical derivative term) 189 INTEGER, INTENT( in ) :: & 190 kahm ! =1 the laplacian is multiplied by the eddy diffusivity coef. 191 ! ! =2 no multiplication 192 193 !! * Local declarations 194 INTEGER :: ji, jj, jk ! dummy loop indices 195 REAL(wp) :: & 196 zabe1, zabe2, zcof1, zcof2, & ! temporary scalars 197 zcoef0, zcoef3, zcoef4 198 REAL(wp) :: & 199 zbur, zbvr, zmkt, zmkf, zuav, zvav, & 200 zuwslpi, zuwslpj, zvwslpi, zvwslpj 201 REAL(wp), DIMENSION(jpi,jpj) :: & 202 ziut, zjuf , zjvt, zivf, & ! workspace 203 zdku, zdk1u, zdkv, zdk1v 204 REAL(wp), DIMENSION(jpi,jpk) :: & 205 zfuw, zfvw, zdiu, zdiv, & ! workspace 206 zdju, zdj1u, zdjv, zdj1v 207 !!---------------------------------------------------------------------- 208 176 !!---------------------------------------------------------------------- 177 USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released 178 USE wrk_nemo, ONLY: ziut => wrk_2d_1 , zjuf => wrk_2d_2 , zjvt => wrk_2d_3 179 USE wrk_nemo, ONLY: zivf => wrk_2d_4 , zdku => wrk_2d_5 , zdk1u => wrk_2d_6 180 USE wrk_nemo, ONLY: zdkv => wrk_2d_7 , zdk1v => wrk_2d_8 181 !! 182 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in ) :: pu , pv ! 1st call: before horizontal velocity 183 ! ! 2nd call: ahm x these fields 184 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( out) :: plu, plv ! partial harmonic operator applied to 185 ! ! pu and pv (all the components except 186 ! ! second order vertical derivative term) 187 INTEGER , INTENT(in ) :: kahm ! =1 1st call ; =2 2nd call 188 ! 189 INTEGER :: ji, jj, jk ! dummy loop indices 190 REAL(wp) :: zabe1 , zabe2 , zcof1 , zcof2 ! local scalar 191 REAL(wp) :: zcoef0, zcoef3, zcoef4 ! - - 192 REAL(wp) :: zbur, zbvr, zmkt, zmkf, zuav, zvav ! - - 193 REAL(wp) :: zuwslpi, zuwslpj, zvwslpi, zvwslpj ! - - 194 !!---------------------------------------------------------------------- 195 196 IF( wrk_in_use(2, 1,2,3,4,5,6,7,8) ) THEN 197 CALL ctl_stop('dyn:ldfguv: requested workspace arrays unavailable') ; RETURN 198 END IF 209 199 ! ! ********** ! ! =============== 210 200 DO jk = 1, jpkm1 ! First step ! ! Horizontal slab … … 422 412 ! II.3 Divergence of vertical fluxes added to the horizontal divergence 423 413 ! --------------------------------------------------------------------- 424 #if defined key_dynldf_smag 414 425 415 IF( kahm == 1 ) THEN 426 416 ! multiply the laplacian by the eddy viscosity coefficient … … 458 448 STOP 'ldfguv' 459 449 ENDIF 460 #else461 IF( kahm == 2 ) THEN462 ! multiply the laplacian by the eddy viscosity coefficient463 DO jk = 1, jpkm1464 DO ji = 2, jpim1465 ! eddy coef. divided by the volume element466 zbur = fsahmu(ji,jj,jk) / ( e1u(ji,jj)*e2u(ji,jj)*fse3u(ji,jj,jk) )467 zbvr = fsahmv(ji,jj,jk) / ( e1v(ji,jj)*e2v(ji,jj)*fse3v(ji,jj,jk) )468 ! vertical divergence469 zuav = zfuw(ji,jk) - zfuw(ji,jk+1)470 zvav = zfvw(ji,jk) - zfvw(ji,jk+1)471 ! harmonic operator applied to (pu,pv) and multiply by ahm472 plu(ji,jj,jk) = ( plu(ji,jj,jk) + zuav ) * zbur473 plv(ji,jj,jk) = ( plv(ji,jj,jk) + zvav ) * zbvr474 END DO475 END DO476 ELSEIF( kahm == 1 ) THEN477 ! second call, no multiplication478 DO jk = 1, jpkm1479 DO ji = 2, jpim1480 ! inverse of the volume element481 zbur = 1. / ( e1u(ji,jj)*e2u(ji,jj)*fse3u(ji,jj,jk) )482 zbvr = 1. / ( e1v(ji,jj)*e2v(ji,jj)*fse3v(ji,jj,jk) )483 ! vertical divergence484 zuav = zfuw(ji,jk) - zfuw(ji,jk+1)485 zvav = zfvw(ji,jk) - zfvw(ji,jk+1)486 ! harmonic operator applied to (pu,pv)487 plu(ji,jj,jk) = ( plu(ji,jj,jk) + zuav ) * zbur488 plv(ji,jj,jk) = ( plv(ji,jj,jk) + zvav ) * zbvr489 END DO490 END DO491 ELSE492 IF(lwp)WRITE(numout,*) ' ldfguv: kahm= 1 or 2, here =', kahm493 IF(lwp)WRITE(numout,*) ' We stop'494 STOP 'ldfguv'495 ENDIF496 497 #endif498 450 ! ! =============== 499 451 END DO ! End of slab 500 452 ! ! =============== 453 454 IF( wrk_not_released(2, 1,2,3,4,5,6,7,8) ) CALL ctl_stop('dyn:ldfguv: failed to release workspace arrays') 455 ! 501 456 END SUBROUTINE ldfguv 502 457 … … 507 462 CONTAINS 508 463 SUBROUTINE dyn_ldf_bilapg( kt ) ! Dummy routine 464 INTEGER, INTENT(in) :: kt 509 465 WRITE(*,*) 'dyn_ldf_bilapg: You should not have seen this print! error?', kt 510 466 END SUBROUTINE dyn_ldf_bilapg -
branches/2011/dev_r2802_NOCL_Smagorinsky/NEMOGCM/NEMO/OPA_SRC/LDF/ldfdyn.F90
r2887 r2926 38 38 !! NEMO/OPA 3.3 , NEMO Consortium (2010) 39 39 !! $Id$ 40 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 41 !!---------------------------------------------------------------------- 42 40 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 41 !!---------------------------------------------------------------------- 43 42 CONTAINS 44 43 … … 63 62 !!---------------------------------------------------------------------- 64 63 INTEGER :: ioptio ! ??? 65 LOGICAL :: ll_print = .FALSE. ! Logical flag for printing viscosity coef.64 LOGICAL :: ll_print = .FALSE. ! Logical flag for printing viscosity coef. 66 65 !! 67 66 NAMELIST/namdyn_ldf/ ln_dynldf_lap , ln_dynldf_bilap, & 68 67 & ln_dynldf_level, ln_dynldf_hor , ln_dynldf_iso, & 69 & rn_ahm_0_lap , rn_ahmb_0 , rn_ahm_0_blp ,cmsmag1,cmsmag268 & rn_ahm_0_lap , rn_ahmb_0 , rn_ahm_0_blp 70 69 !!---------------------------------------------------------------------- 71 70 … … 86 85 WRITE(numout,*) ' background viscosity rn_ahmb_0 = ', rn_ahmb_0 87 86 WRITE(numout,*) ' horizontal bilaplacian eddy viscosity rn_ahm_0_blp = ', rn_ahm_0_blp 88 WRITE(numout,*) ' smagorinsky coefficient for laplacian cmsmag1 = ' , cmsmag189 WRITE(numout,*) ' smagorinsky coefficient for bilaplacian cmsmag2 = ' , cmsmag290 87 ENDIF 91 88 … … 131 128 ! Lateral eddy viscosity 132 129 ! ====================== 133 #if defined key_dynldf_smag && ! defined key_dynldf_c3d134 IF(lwp) WRITE(numout,*) ' key_dynldf_smag can not be used without key_dynldf_c3d'135 #endif136 137 130 #if defined key_dynldf_c3d 138 131 CALL ldf_dyn_c3d( ll_print ) ! ahm = 3D coef. = F( longitude, latitude, depth ) … … 213 206 REAL(wp), INTENT(in ) :: pwam ! width of inflection 214 207 REAL(wp), INTENT(in ) :: pbot ! bottom value (0<pbot<= 1) 215 REAL(wp), INTENT(in ), DIMENSION (jpk) :: pdep ! depth of the gridpoint (T, U, V, F)216 REAL(wp), INTENT(inout), DIMENSION (jpi,jpj,jpk) :: pah ! adimensional vertical profile208 REAL(wp), INTENT(in ), DIMENSION (:) :: pdep ! depth of the gridpoint (T, U, V, F) 209 REAL(wp), INTENT(inout), DIMENSION (:,:,:) :: pah ! adimensional vertical profile 217 210 !! 218 211 INTEGER :: jk ! dummy loop indices … … 255 248 REAL(wp), INTENT(in ) :: pwam ! width of inflection 256 249 REAL(wp), INTENT(in ) :: pbot ! bottom value (0<pbot<= 1) 257 REAL(wp), INTENT(in ), DIMENSION (jpi,jpj,jpk) :: pdep ! dep of the gridpoint (T, U, V, F)258 REAL(wp), INTENT(inout), DIMENSION (jpi,jpj,jpk) :: pah ! adimensional vertical profile250 REAL(wp), INTENT(in ), DIMENSION (:,:,:) :: pdep ! dep of the gridpoint (T, U, V, F) 251 REAL(wp), INTENT(inout), DIMENSION (:,:,:) :: pah ! adimensional vertical profile 259 252 !! 260 253 INTEGER :: jk ! dummy loop indices -
branches/2011/dev_r2802_NOCL_Smagorinsky/NEMOGCM/NEMO/OPA_SRC/LDF/ldfdyn_oce.F90
r2887 r2926 23 23 REAL(wp), PUBLIC :: rn_ahm_0_blp = 0._wp !: lateral bilaplacian eddy viscosity (m4/s) 24 24 REAL(wp), PUBLIC :: ahm0, ahmb0, ahm0_blp !: OLD namelist names 25 REAL(wp), PUBLIC :: cmsmag1 = 3._wp !: constant in laplacian Smagorinsky viscosity 26 REAL(wp), PUBLIC :: cmsmag2 = 3._wp !: constant in bilaplacian Smagorinsky viscosity 25 27 26 ! !!! eddy coeff. at U-,V-,W-pts [m2/s] 28 27 #if defined key_dynldf_c3d -
branches/2011/dev_r2802_NOCL_Smagorinsky/NEMOGCM/NEMO/OPA_SRC/LDF/ldfslp.F90
r2887 r2926 31 31 USE in_out_manager ! I/O manager 32 32 USE prtctl ! Print control 33 USE lib_mpp34 33 35 34 IMPLICIT NONE … … 692 691 ! !- bound the slopes: abs(zw.)<= 1/100 and zb..<0 693 692 ! kxz max= ah slope max =< e1 e3 /(pi**2 2 dt) 694 zbu = MIN( zbu , -100._wp* ABS( zau ) , -7.e+3_wp/fse3u(ji,jj,ik )* ABS( zau ) )695 zbv = MIN( zbv , -100._wp* ABS( zav ) , -7.e+3_wp/fse3v(ji,jj,ik )* ABS( zav ) )693 zbu = MIN( zbu , -100._wp* ABS( zau ) , -7.e+3_wp/fse3u(ji,jj,iku)* ABS( zau ) ) 694 zbv = MIN( zbv , -100._wp* ABS( zav ) , -7.e+3_wp/fse3v(ji,jj,ikv)* ABS( zav ) ) 696 695 ! !- Slope at u- & v-points (uslpml, vslpml) 697 uslpml(ji,jj) = zau / ( zbu - zeps ) * umask(ji,jj,ik )698 vslpml(ji,jj) = zav / ( zbv - zeps ) * vmask(ji,jj,ik )696 uslpml(ji,jj) = zau / ( zbu - zeps ) * umask(ji,jj,iku) 697 vslpml(ji,jj) = zav / ( zbv - zeps ) * vmask(ji,jj,ikv) 699 698 ! 700 699 ! !== i- & j-slopes at w-points just below the Mixed Layer ==! -
branches/2011/dev_r2802_NOCL_Smagorinsky/NEMOGCM/NEMO/OPA_SRC/LDF/ldftra.F90
r2887 r2926 68 68 & ln_traldf_grif , ln_traldf_gdia, & 69 69 & rn_aht_0 , rn_ahtb_0 , rn_aeiv_0, & 70 & rn_slpmax , chsmag70 & rn_slpmax 71 71 !!---------------------------------------------------------------------- 72 72 … … 157 157 ENDIF 158 158 #endif 159 #if defined key_traldf_smag && ! defined key_traldf_c3d160 CALL ctl_stop( 'key_traldf_smag can only be used with key_traldf_c3d' )161 #endif162 163 159 ! 164 160 END SUBROUTINE ldf_tra_init -
branches/2011/dev_r2802_NOCL_Smagorinsky/NEMOGCM/NEMO/OPA_SRC/LDF/ldftra_c3d.h90
r2887 r2926 89 89 ENDIF 90 90 91 # if defined key_traldf_eiv 91 # if defined key_traldf_eiv 92 92 aeiu(:,:,:) = aeiv0 ! set aeiu = aeiv at u- and v-points, 93 93 aeiv(:,:,:) = aeiv0 ! and aeiw at w-point … … 108 108 CALL lbc_lnk( aeiv, 'V', 1. ) 109 109 CALL lbc_lnk( aeiw, 'W', 1. ) 110 # endif 110 111 111 112 IF(lwp .AND. ld_print ) THEN … … 120 121 CALL prihre(aeiw(:,:,1),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 121 122 ENDIF 122 # endif123 123 ! 124 124 END SUBROUTINE ldf_tra_c3d -
branches/2011/dev_r2802_NOCL_Smagorinsky/NEMOGCM/NEMO/OPA_SRC/LDF/ldftra_oce.F90
r2887 r2926 30 30 REAL(wp), PUBLIC :: rn_aeiv_0 = 2000._wp !: eddy induced velocity coefficient (m2/s) 31 31 REAL(wp), PUBLIC :: rn_slpmax = 0.01_wp !: slope limit 32 REAL(wp), PUBLIC :: chsmag = 1._wp !: constant in Smagorinsky diffusivity 32 33 33 REAL(wp), PUBLIC :: aht0, ahtb0, aeiv0 !!: OLD namelist names 34 34 LOGICAL , PUBLIC :: l_triad_iso = .FALSE. !: calculate triads twice -
branches/2011/dev_r2802_NOCL_Smagorinsky/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_bilap.F90
r2887 r2926 28 28 USE diaptr ! poleward transport diagnostics 29 29 USE trc_oce ! share passive tracers/Ocean variables 30 USE lib_mpp ! MPP library 30 31 31 32 IMPLICIT NONE … … 73 74 !! biharmonic mixing trend. 74 75 !!---------------------------------------------------------------------- 75 USE oce , ztu => ua ! use ua as workspace 76 USE oce , ztv => va ! use va as workspace 76 USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released 77 USE oce , ONLY: ztu => ua , ztv => va ! (ua,va) used as workspace 78 USE wrk_nemo, ONLY: zeeu => wrk_2d_1 , zeev => wrk_2d_2 , zlt => wrk_2d_3 ! 2D workspace 77 79 !! 78 80 INTEGER , INTENT(in ) :: kt ! ocean time-step index … … 85 87 INTEGER :: ji, jj, jk, jn ! dummy loop indices 86 88 REAL(wp) :: zbtr, ztra ! local scalars 87 REAL(wp), DIMENSION(jpi,jpj) :: zeeu, zeev, zlt ! 2D workspace88 89 !!---------------------------------------------------------------------- 90 91 IF( wrk_in_use(2, 1,2,3) ) THEN 92 CALL ctl_stop('tra_ldf_bilap: requested workspace arrays unavailable') ; RETURN 93 ENDIF 89 94 90 95 IF( kt == nit000 ) THEN … … 123 128 END DO 124 129 ENDIF 125 #if ! defined key_traldf_smag126 130 DO jj = 2, jpjm1 ! Second derivative (divergence) time the eddy diffusivity coefficient 127 131 DO ji = fs_2, fs_jpim1 ! vector opt. … … 131 135 END DO 132 136 END DO 133 #else134 DO jj = 2, jpjm1 ! Second derivative (divergence) time the eddy diffusivity coefficient135 DO ji = fs_2, fs_jpim1 ! vector opt.136 zbtr = 1.0 / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )137 zlt(ji,jj) = zbtr * ( ztu(ji,jj,jk) - ztu(ji-1,jj,jk) &138 & + ztv(ji,jj,jk) - ztv(ji,jj-1,jk) )139 END DO140 END DO141 #endif142 137 CALL lbc_lnk( zlt, 'T', 1. ) ! Lateral boundary conditions (unchanged sgn) 143 138 … … 150 145 END DO 151 146 END DO 152 #if ! defined key_traldf_smag153 147 DO jj = 2, jpjm1 ! fourth derivative (divergence) and add to the general tracer trend 154 148 DO ji = fs_2, fs_jpim1 ! vector opt. … … 160 154 END DO 161 155 END DO 162 #else163 DO jj = 2, jpjm1 ! fourth derivative (divergence) and add to the general tracer trend164 DO ji = fs_2, fs_jpim1 ! vector opt.165 ! horizontal diffusive trends166 zbtr = 1.0 / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )167 ztra = zbtr * fsahtt(ji,jj,jk)* ( ztu(ji,jj,jk) - ztu(ji-1,jj,jk) + ztv(ji,jj,jk) - ztv(ji,jj-1,jk) )168 ! add it to the general tracer trends169 pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra170 END DO171 END DO172 173 #endif174 156 ! 175 157 END DO ! Horizontal slab … … 183 165 END DO ! tracer loop 184 166 ! ! =========== 167 IF( wrk_not_released(2, 1,2,3) ) CALL ctl_stop('tra_ldf_bilap: failed to release workspace arrays') 168 ! 185 169 END SUBROUTINE tra_ldf_bilap 186 170 -
branches/2011/dev_r2802_NOCL_Smagorinsky/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_bilapg.F90
r2887 r2926 24 24 USE diaptr ! poleward transport diagnostics 25 25 USE trc_oce ! share passive tracers/Ocean variables 26 USE lib_mpp ! MPP library 26 27 27 28 IMPLICIT NONE … … 65 66 !! biharmonic mixing trend. 66 67 !!---------------------------------------------------------------------- 68 USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released 69 USE wrk_nemo, ONLY: wk1 => wrk_4d_1 , wk2 => wrk_4d_2 ! 4D workspace 70 ! 67 71 INTEGER , INTENT(in ) :: kt ! ocean time-step index 68 72 CHARACTER(len=3), INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator) … … 70 74 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: ptb ! before and now tracer fields 71 75 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: pta ! tracer trend 72 !! 73 INTEGER :: ji, jj, jk, jn ! dummy loop indices 74 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt) :: wk1, wk2 ! 4D workspace 75 !!---------------------------------------------------------------------- 76 ! 77 INTEGER :: ji, jj, jk, jn ! dummy loop indices 78 !!---------------------------------------------------------------------- 79 80 IF( wrk_in_use(4, 1,2) ) THEN 81 CALL ctl_stop('tra_ldf_bilapg: requested workspace arrays unavailable') ; RETURN 82 ENDIF 76 83 77 84 IF( kt == nit000 ) THEN … … 107 114 END DO 108 115 END DO 116 ! 117 IF( wrk_not_released(4, 1,2) ) CALL ctl_stop('tra_ldf_bilapg : failed to release workspace arrays.') 109 118 ! 110 119 END SUBROUTINE tra_ldf_bilapg … … 149 158 !! 150 159 !!---------------------------------------------------------------------- 151 USE oce , zftv => ua ! use ua as workspace 152 !! 160 USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released, wrk_in_use_xz, wrk_not_released_xz 161 USE oce , ONLY: zftv => ua ! ua used as workspace 162 USE wrk_nemo, ONLY: zftu => wrk_2d_1 , zdkt => wrk_2d_2 , zdk1t => wrk_2d_3 163 USE wrk_nemo, ONLY: zftw => wrk_xz_1 , zdit => wrk_xz_2 164 USE wrk_nemo, ONLY: zdjt => wrk_xz_3 , zdj1t => wrk_xz_4 165 ! 153 166 INTEGER , INTENT(in ) :: kt ! ocean time-step index 154 167 CHARACTER(len=3), INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator) … … 166 179 REAL(wp) :: zbtr, ztah, ztav 167 180 REAL(wp) :: zcof0, zcof1, zcof2, zcof3, zcof4 168 REAL(wp), DIMENSION(jpi,jpj) :: zftu, zdkt, zdk1t ! workspace 169 REAL(wp), DIMENSION(jpi,jpk) :: zftw, zdit, zdjt, zdj1t ! 170 !!---------------------------------------------------------------------- 171 181 !!---------------------------------------------------------------------- 182 183 IF( wrk_in_use(2, 1,2,3) .OR. wrk_in_use_xz(1,2,3,4) )THEN 184 CALL ctl_stop('ldfght : requested workspace arrays unavailable') ; RETURN 185 ENDIF 172 186 ! 173 187 DO jn = 1, kjpt … … 285 299 286 300 ! II.3 Divergence of vertical fluxes added to the horizontal divergence 287 ! -------------------------------------------------------------------- 288 #if ! defined key_traldf_smag301 ! --------------------------------------------------------------------- 302 289 303 IF( kaht == 1 ) THEN 290 304 ! multiply the laplacian by the eddy diffusivity coefficient … … 320 334 ! ! =============== 321 335 END DO 322 #else 323 IF( kaht == 2 ) THEN 324 ! multiply the laplacian by the eddy diffusivity coefficient 325 DO jk = 1, jpkm1 326 DO ji = 2, jpim1 327 ! eddy coef. divided by the volume element 328 zbtr = 1.0 / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 329 ! vertical divergence 330 ztav = fsahtt(ji,jj,jk) * ( zftw(ji,jk) - zftw(ji,jk+1) ) 331 ! harmonic operator applied to (pt,ps) and multiply by aht 332 plt(ji,jj,jk,jn) = ( plt(ji,jj,jk,jn) + ztav ) * zbtr 333 END DO 334 END DO 335 ELSEIF( kaht == 1 ) THEN 336 ! second call, no multiplication 337 DO jk = 1, jpkm1 338 DO ji = 2, jpim1 339 ! inverse of the volume element 340 zbtr = 1.0 / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 341 ! vertical divergence 342 ztav = zftw(ji,jk) - zftw(ji,jk+1) 343 ! harmonic operator applied to (pt,ps) 344 plt(ji,jj,jk,jn) = ( plt(ji,jj,jk,jn) + ztav ) * zbtr 345 END DO 346 END DO 347 ELSE 348 IF(lwp) WRITE(numout,*) ' ldfght: kaht= 1 or 2, here =', kaht 349 IF(lwp) WRITE(numout,*) ' We stop' 350 STOP 'ldfght' 351 ENDIF 352 ! ! =============== 353 END DO ! End of slab 354 ! ! =============== 355 END DO 356 357 358 #endif 336 ! 337 IF( wrk_not_released(2, 1,2,3) .OR. & 338 wrk_not_released_xz(1,2,3,4) ) CALL ctl_stop('ldfght : failed to release workspace arrays.') 339 ! 359 340 END SUBROUTINE ldfght 360 341 -
branches/2011/dev_r2802_NOCL_Smagorinsky/NEMOGCM/NEMO/OPA_SRC/step.F90
r2887 r2926 156 156 IF( lk_traldf_eiv ) CALL ldf_eiv( kstp ) ! eddy induced velocity coefficient 157 157 #endif 158 #if defined key_traldf_c3d && key_traldf_smag159 CALL ldf_tra_smag( kstp ) ! Smagorinsky diffusivity160 # endif161 #if defined key_dynldf_c3d && key_dynldf_smag162 CALL ldf_dyn_smag( kstp ) ! Smagorinsky viscosity163 # endif164 165 158 166 159 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> -
branches/2011/dev_r2802_NOCL_Smagorinsky/NEMOGCM/NEMO/OPA_SRC/step_oce.F90
r2887 r2926 60 60 USE ldfslp ! iso-neutral slopes (ldf_slp routine) 61 61 USE ldfeiv ! eddy induced velocity coef. (ldf_eiv routine) 62 USE ldftra_smag ! smagorinsky diffusivity . (ldf_tra_smag routine)63 USE ldfdyn_smag ! smagorinsky viscosity (ldf_dyn_smag routine)64 62 65 63 USE zdftmx ! tide-induced vertical mixing (zdf_tmx routine)
Note: See TracChangeset
for help on using the changeset viewer.