- Timestamp:
- 2016-07-19T10:38:35+02:00 (8 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/NERC/dev_r5549_BDY_ZEROGRAD/NEMOGCM/NEMO/OPA_SRC/LDF/ldfdyn.F90
r4624 r6808 6 6 !! History : OPA ! 1997-07 (G. Madec) multi dimensional coefficients 7 7 !! NEMO 1.0 ! 2002-09 (G. Madec) F90: Free form and module 8 !! ----------------------------------------------------------------------9 10 !!---------------------------------------------------------------------- 11 !! ldf_dyn_init : initialization, namelist read, and parameters control 12 !! ldf_dyn_c3d : 3D eddy viscosity coefficient initialization13 !! ldf_dyn_ c2d : 2D eddy viscosity coefficient initialization14 !! ldf_dyn _c1d : 1D eddy viscosity coefficient initialization8 !! 3.7 ! 2014-01 (F. Lemarie, G. Madec) restructuration/simplification of ahm specification, 9 !! ! add velocity dependent coefficient and optional read in file 10 !!---------------------------------------------------------------------- 11 12 !!---------------------------------------------------------------------- 13 !! ldf_dyn_init : initialization, namelist read, and parameters control 14 !! ldf_dyn : update lateral eddy viscosity coefficients at each time step 15 15 !!---------------------------------------------------------------------- 16 16 USE oce ! ocean dynamics and tracers 17 17 USE dom_oce ! ocean space and time domain 18 USE ldfdyn_oce ! ocean dynamics lateral physics19 18 USE phycst ! physical constants 20 USE ldf slp ! ???21 USE ioipsl19 USE ldfc1d_c2d ! lateral diffusion: 1D and 2D cases 20 ! 22 21 USE in_out_manager ! I/O manager 22 USE iom ! I/O module for ehanced bottom friction file 23 USE timing ! Timing 23 24 USE lib_mpp ! distribued memory computing library 24 25 USE lbclnk ! ocean lateral boundary conditions (or mpp link) … … 28 29 PRIVATE 29 30 30 PUBLIC ldf_dyn_init ! called by opa.F90 31 32 INTERFACE ldf_zpf 33 MODULE PROCEDURE ldf_zpf_1d, ldf_zpf_1d_3d, ldf_zpf_3d 34 END INTERFACE 31 PUBLIC ldf_dyn_init ! called by nemogcm.F90 32 PUBLIC ldf_dyn ! called by step.F90 33 34 ! !!* Namelist namdyn_ldf : lateral mixing on momentum * 35 LOGICAL , PUBLIC :: ln_dynldf_lap !: laplacian operator 36 LOGICAL , PUBLIC :: ln_dynldf_blp !: bilaplacian operator 37 LOGICAL , PUBLIC :: ln_dynldf_lev !: iso-level direction 38 LOGICAL , PUBLIC :: ln_dynldf_hor !: horizontal (geopotential) direction 39 LOGICAL , PUBLIC :: ln_dynldf_iso !: iso-neutral direction 40 INTEGER , PUBLIC :: nn_ahm_ijk_t !: choice of time & space variations of the lateral eddy viscosity coef. 41 REAL(wp), PUBLIC :: rn_ahm_0 !: lateral laplacian eddy viscosity [m2/s] 42 REAL(wp), PUBLIC :: rn_ahm_b !: lateral laplacian background eddy viscosity [m2/s] 43 REAL(wp), PUBLIC :: rn_bhm_0 !: lateral bilaplacian eddy viscosity [m4/s] 44 45 LOGICAL , PUBLIC :: l_ldfdyn_time !: flag for time variation of the lateral eddy viscosity coef. 46 47 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: ahmt, ahmf !: eddy diffusivity coef. at U- and V-points [m2/s or m4/s] 48 49 REAL(wp) :: r1_12 = 1._wp / 12._wp ! =1/12 50 REAL(wp) :: r1_4 = 0.25_wp ! =1/4 51 REAL(wp) :: r1_288 = 1._wp / 288._wp ! =1/( 12^2 * 2 ) 35 52 36 53 !! * Substitutions 37 # include " domzgr_substitute.h90"38 !!---------------------------------------------------------------------- 39 !! NEMO/OPA 3. 3 , NEMO Consortium (2010)54 # include "vectopt_loop_substitute.h90" 55 !!---------------------------------------------------------------------- 56 !! NEMO/OPA 3.7 , NEMO Consortium (2014) 40 57 !! $Id$ 41 58 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) … … 49 66 !! ** Purpose : set the horizontal ocean dynamics physics 50 67 !! 51 !! ** Method : 52 !! - default option : ahm = constant coef. = rn_ahm_0 (namelist) 53 !! - 'key_dynldf_c1d': ahm = F(depth) see ldf_dyn_c1d.h90 54 !! - 'key_dynldf_c2d': ahm = F(latitude,longitude) see ldf_dyn_c2d.h90 55 !! - 'key_dynldf_c3d': ahm = F(latitude,longitude,depth) see ldf_dyn_c3d.h90 56 !! 57 !! N.B. User defined include files. By default, 3d and 2d coef. 58 !! are set to a constant value given in the namelist and the 1d 59 !! coefficients are initialized to a hyperbolic tangent vertical 60 !! profile. 61 !! 62 !! Reference : Madec, G. and M. Imbard, 1996: Climate Dynamics, 12, 381-388. 63 !!---------------------------------------------------------------------- 64 INTEGER :: ioptio ! ??? 65 INTEGER :: ios ! Local : output status for namelist read 66 LOGICAL :: ll_print = .FALSE. ! Logical flag for printing viscosity coef. 67 !! 68 NAMELIST/namdyn_ldf/ ln_dynldf_lap , ln_dynldf_bilap, & 69 & ln_dynldf_level, ln_dynldf_hor , ln_dynldf_iso, & 70 & rn_ahm_0_lap , rn_ahmb_0 , rn_ahm_0_blp , & 71 & rn_cmsmag_1 , rn_cmsmag_2 , rn_cmsh, & 72 & rn_ahm_m_lap , rn_ahm_m_blp 73 74 !!---------------------------------------------------------------------- 75 68 !! ** Method : the eddy viscosity coef. specification depends on: 69 !! - the operator: 70 !! ln_dynldf_lap = T laplacian operator 71 !! ln_dynldf_blp = T bilaplacian operator 72 !! - the parameter nn_ahm_ijk_t: 73 !! nn_ahm_ijk_t = 0 => = constant 74 !! = 10 => = F(z) : = constant with a reduction of 1/4 with depth 75 !! =-20 => = F(i,j) = shape read in 'eddy_viscosity.nc' file 76 !! = 20 = F(i,j) = F(e1,e2) or F(e1^3,e2^3) (lap or bilap case) 77 !! =-30 => = F(i,j,k) = shape read in 'eddy_viscosity.nc' file 78 !! = 30 = F(i,j,k) = 2D (case 20) + decrease with depth (case 10) 79 !! = 31 = F(i,j,k,t) = F(local velocity) ( |u|e /12 laplacian operator 80 !! or |u|e^3/12 bilaplacian operator ) 81 !!---------------------------------------------------------------------- 82 INTEGER :: jk ! dummy loop indices 83 INTEGER :: ierr, inum, ios ! local integer 84 REAL(wp) :: zah0 ! local scalar 85 ! 86 NAMELIST/namdyn_ldf/ ln_dynldf_lap, ln_dynldf_blp, & 87 & ln_dynldf_lev, ln_dynldf_hor, ln_dynldf_iso, & 88 & nn_ahm_ijk_t , rn_ahm_0, rn_ahm_b, rn_bhm_0 89 !!---------------------------------------------------------------------- 90 ! 76 91 REWIND( numnam_ref ) ! Namelist namdyn_ldf in reference namelist : Lateral physics 77 92 READ ( numnam_ref, namdyn_ldf, IOSTAT = ios, ERR = 901) … … 88 103 WRITE(numout,*) '~~~~~~~' 89 104 WRITE(numout,*) ' Namelist namdyn_ldf : set lateral mixing parameters' 90 WRITE(numout,*) ' laplacian operator ln_dynldf_lap = ', ln_dynldf_lap 91 WRITE(numout,*) ' bilaplacian operator ln_dynldf_bilap = ', ln_dynldf_bilap 92 WRITE(numout,*) ' iso-level ln_dynldf_level = ', ln_dynldf_level 93 WRITE(numout,*) ' horizontal (geopotential) ln_dynldf_hor = ', ln_dynldf_hor 94 WRITE(numout,*) ' iso-neutral ln_dynldf_iso = ', ln_dynldf_iso 95 WRITE(numout,*) ' horizontal laplacian eddy viscosity rn_ahm_0_lap = ', rn_ahm_0_lap 96 WRITE(numout,*) ' background viscosity rn_ahmb_0 = ', rn_ahmb_0 97 WRITE(numout,*) ' horizontal bilaplacian eddy viscosity rn_ahm_0_blp = ', rn_ahm_0_blp 98 WRITE(numout,*) ' upper limit for laplacian eddy visc rn_ahm_m_lap = ', rn_ahm_m_lap 99 WRITE(numout,*) ' upper limit for bilap eddy viscosity rn_ahm_m_blp = ', rn_ahm_m_blp 100 101 ENDIF 102 103 ahm0 = rn_ahm_0_lap ! OLD namelist variables defined from DOCTOR namelist variables 104 ahmb0 = rn_ahmb_0 105 ahm0_blp = rn_ahm_0_blp 106 107 ! ... check of lateral diffusive operator on tracers 108 ! ==> will be done in trazdf module 109 110 ! ... Space variation of eddy coefficients 111 ioptio = 0 112 #if defined key_dynldf_c3d 113 IF(lwp) WRITE(numout,*) ' momentum mixing coef. = F( latitude, longitude, depth)' 114 ioptio = ioptio+1 115 #endif 116 #if defined key_dynldf_c2d 117 IF(lwp) WRITE(numout,*) ' momentum mixing coef. = F( latitude, longitude)' 118 ioptio = ioptio+1 119 #endif 120 #if defined key_dynldf_c1d 121 IF(lwp) WRITE(numout,*) ' momentum mixing coef. = F( depth )' 122 ioptio = ioptio+1 123 IF( ln_sco ) CALL ctl_stop( 'key_dynldf_c1d cannot be used in s-coordinate (ln_sco)' ) 124 #endif 125 IF( ioptio == 0 ) THEN 126 IF(lwp) WRITE(numout,*) ' momentum mixing coef. = constant (default option)' 127 ELSEIF( ioptio > 1 ) THEN 128 CALL ctl_stop( 'use only one of the following keys: key_dynldf_c3d, key_dynldf_c2d, key_dynldf_c1d' ) 129 ENDIF 130 131 132 IF( ln_dynldf_bilap ) THEN 133 IF(lwp) WRITE(numout,*) ' biharmonic momentum diffusion' 134 IF( .NOT. ln_dynldf_lap ) ahm0 = ahm0_blp ! Allow spatially varying coefs, which use ahm0 as input 135 IF( ahm0_blp > 0 .AND. .NOT. lk_esopa ) CALL ctl_stop( 'The horizontal viscosity coef. ahm0 must be negative' ) 136 ELSE 137 IF(lwp) WRITE(numout,*) ' harmonic momentum diff. (default)' 138 IF( ahm0 < 0 .AND. .NOT. lk_esopa ) CALL ctl_stop( 'The horizontal viscosity coef. ahm0 must be positive' ) 139 ENDIF 140 141 142 ! Lateral eddy viscosity 143 ! ====================== 144 #if defined key_dynldf_c3d 145 CALL ldf_dyn_c3d( ll_print ) ! ahm = 3D coef. = F( longitude, latitude, depth ) 146 #elif defined key_dynldf_c2d 147 CALL ldf_dyn_c2d( ll_print ) ! ahm = 1D coef. = F( longitude, latitude ) 148 #elif defined key_dynldf_c1d 149 CALL ldf_dyn_c1d( ll_print ) ! ahm = 1D coef. = F( depth ) 150 #else 151 ! Constant coefficients 152 IF(lwp) WRITE(numout,*) 153 IF(lwp) WRITE(numout,*) 'inildf: constant eddy viscosity coef. ' 154 IF(lwp) WRITE(numout,*) '~~~~~~' 155 IF(lwp) WRITE(numout,*) ' ahm1 = ahm2 = ahm0 = ',ahm0 156 #endif 157 nkahm_smag = 0 158 #if defined key_dynldf_smag 159 nkahm_smag = 1 160 #endif 161 105 ! 106 WRITE(numout,*) ' type :' 107 WRITE(numout,*) ' laplacian operator ln_dynldf_lap = ', ln_dynldf_lap 108 WRITE(numout,*) ' bilaplacian operator ln_dynldf_blp = ', ln_dynldf_blp 109 ! 110 WRITE(numout,*) ' direction of action :' 111 WRITE(numout,*) ' iso-level ln_dynldf_lev = ', ln_dynldf_lev 112 WRITE(numout,*) ' horizontal (geopotential) ln_dynldf_hor = ', ln_dynldf_hor 113 WRITE(numout,*) ' iso-neutral ln_dynldf_iso = ', ln_dynldf_iso 114 ! 115 WRITE(numout,*) ' coefficients :' 116 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' 118 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' 120 ENDIF 121 122 ! ! Parameter control 123 IF( .NOT.ln_dynldf_lap .AND. .NOT.ln_dynldf_blp ) THEN 124 IF(lwp) WRITE(numout,*) ' No viscous operator selected. ahmt and ahmf are not allocated' 125 l_ldfdyn_time = .FALSE. 126 RETURN 127 ENDIF 128 ! 129 IF( ln_dynldf_blp .AND. ln_dynldf_iso ) THEN ! iso-neutral bilaplacian not implemented 130 CALL ctl_stop( 'dyn_ldf_init: iso-neutral bilaplacian not coded yet' ) 131 ENDIF 132 133 ! ... Space/Time variation of eddy coefficients 134 ! ! allocate the ahm arrays 135 ALLOCATE( ahmt(jpi,jpj,jpk) , ahmf(jpi,jpj,jpk) , STAT=ierr ) 136 IF( ierr /= 0 ) CALL ctl_stop( 'STOP', 'ldf_dyn_init: failed to allocate arrays') 137 ! 138 ahmt(:,:,jpk) = 0._wp ! last level always 0 139 ahmf(:,:,jpk) = 0._wp 140 ! 141 ! ! value of eddy mixing coef. 142 IF ( ln_dynldf_lap ) THEN ; zah0 = rn_ahm_0 ! laplacian operator 143 ELSEIF( ln_dynldf_blp ) THEN ; zah0 = ABS( rn_bhm_0 ) ! bilaplacian operator 144 ELSE ! NO viscous operator 145 CALL ctl_warn( 'ldf_dyn_init: No lateral viscous operator used ' ) 146 ENDIF 147 ! 148 l_ldfdyn_time = .FALSE. ! no time variation except in case defined below 149 ! 150 IF( ln_dynldf_lap .OR. ln_dynldf_blp ) THEN ! only if a lateral diffusion operator is used 151 ! 152 SELECT CASE( nn_ahm_ijk_t ) ! Specification of space time variations of ahmt, ahmf 153 ! 154 CASE( 0 ) !== constant ==! 155 IF(lwp) WRITE(numout,*) ' momentum mixing coef. = constant ' 156 ahmt(:,:,:) = zah0 * tmask(:,:,:) 157 ahmf(:,:,:) = zah0 * fmask(:,:,:) 158 ! 159 CASE( 10 ) !== fixed profile ==! 160 IF(lwp) WRITE(numout,*) ' momentum mixing coef. = F( depth )' 161 ahmt(:,:,1) = zah0 * tmask(:,:,1) ! constant surface value 162 ahmf(:,:,1) = zah0 * fmask(:,:,1) 163 CALL ldf_c1d( 'DYN', r1_4, ahmt(:,:,1), ahmf(:,:,1), ahmt, ahmf ) 164 ! 165 CASE ( -20 ) !== fixed horizontal shape read in file ==! 166 IF(lwp) WRITE(numout,*) ' momentum mixing coef. = F(i,j) read in eddy_viscosity.nc file' 167 CALL iom_open( 'eddy_viscosity_2D.nc', inum ) 168 CALL iom_get ( inum, jpdom_data, 'ahmt_2d', ahmt(:,:,1) ) 169 CALL iom_get ( inum, jpdom_data, 'ahmf_2d', ahmf(:,:,1) ) 170 CALL iom_close( inum ) 171 !!gm Question : info for LAP or BLP case to take into account the SQRT in the bilaplacian case ??? 172 !! do we introduce a scaling by the max value of the array, and then multiply by zah0 ???? 173 !! better: check that the max is <=1 i.e. it is a shape from 0 to 1, not a coef that has physical dimension 174 DO jk = 2, jpkm1 175 ahmt(:,:,jk) = ahmt(:,:,1) * tmask(:,:,jk) 176 ahmf(:,:,jk) = ahmf(:,:,1) * fmask(:,:,jk) 177 END DO 178 ! 179 CASE( 20 ) !== fixed horizontal shape ==! 180 IF(lwp) WRITE(numout,*) ' momentum mixing coef. = F( e1, e2 ) or F( e1^3, e2^3 ) (lap. or blp. case)' 181 IF( ln_dynldf_lap ) CALL ldf_c2d( 'DYN', 'LAP', zah0, ahmt, ahmf ) ! surface value proportional to scale factor 182 IF( ln_dynldf_blp ) CALL ldf_c2d( 'DYN', 'BLP', zah0, ahmt, ahmf ) ! surface value proportional to scale factor^3 183 ! 184 CASE( -30 ) !== fixed 3D shape read in file ==! 185 IF(lwp) WRITE(numout,*) ' momentum mixing coef. = F(i,j,k) read in eddy_diffusivity_3D.nc file' 186 CALL iom_open( 'eddy_viscosity_3D.nc', inum ) 187 CALL iom_get ( inum, jpdom_data, 'ahmt_3d', ahmt ) 188 CALL iom_get ( inum, jpdom_data, 'ahmf_3d', ahmf ) 189 CALL iom_close( inum ) 190 !!gm Question : info for LAP or BLP case to take into account the SQRT in the bilaplacian case ???? 191 !! do we introduce a scaling by the max value of the array, and then multiply by zah0 ???? 192 DO jk = 1, jpkm1 193 ahmt(:,:,jk) = ahmt(:,:,jk) * tmask(:,:,jk) 194 ahmf(:,:,jk) = ahmf(:,:,jk) * fmask(:,:,jk) 195 END DO 196 ! 197 CASE( 30 ) !== fixed 3D shape ==! 198 IF(lwp) WRITE(numout,*) ' momentum mixing coef. = F( latitude, longitude, depth )' 199 IF( ln_dynldf_lap ) CALL ldf_c2d( 'DYN', 'LAP', zah0, ahmt, ahmf ) ! surface value proportional to scale factor 200 IF( ln_dynldf_blp ) CALL ldf_c2d( 'DYN', 'BLP', zah0, ahmt, ahmf ) ! surface value proportional to scale factor 201 ! ! reduction with depth 202 CALL ldf_c1d( 'DYN', r1_4, ahmt(:,:,1), ahmf(:,:,1), ahmt, ahmf ) 203 ! 204 CASE( 31 ) !== time varying 3D field ==! 205 IF(lwp) WRITE(numout,*) ' momentum mixing coef. = F( latitude, longitude, depth , time )' 206 IF(lwp) WRITE(numout,*) ' proportional to the velocity : |u|e/12 or |u|e^3/12' 207 ! 208 l_ldfdyn_time = .TRUE. ! will be calculated by call to ldf_dyn routine in step.F90 209 ! 210 CASE DEFAULT 211 CALL ctl_stop('ldf_dyn_init: wrong choice for nn_ahm_ijk_t, the type of space-time variation of ahm') 212 END SELECT 213 ! 214 IF( ln_dynldf_blp .AND. .NOT. l_ldfdyn_time ) THEN ! bilapcian and no time variation: 215 ahmt(:,:,:) = SQRT( ahmt(:,:,:) ) ! take the square root of the coefficient 216 ahmf(:,:,:) = SQRT( ahmf(:,:,:) ) 217 ENDIF 218 ! 219 ENDIF 162 220 ! 163 221 END SUBROUTINE ldf_dyn_init 164 222 165 #if defined key_dynldf_c3d 166 # include "ldfdyn_c3d.h90" 167 #elif defined key_dynldf_c2d 168 # include "ldfdyn_c2d.h90" 169 #elif defined key_dynldf_c1d 170 # include "ldfdyn_c1d.h90" 171 #endif 172 173 174 SUBROUTINE ldf_zpf_1d( ld_print, pdam, pwam, pbot, pdep, pah ) 175 !!---------------------------------------------------------------------- 176 !! *** ROUTINE ldf_zpf *** 177 !! 178 !! ** Purpose : vertical adimensional profile for eddy coefficient 179 !! 180 !! ** Method : 1D eddy viscosity coefficients ( depth ) 181 !!---------------------------------------------------------------------- 182 LOGICAL , INTENT(in ) :: ld_print ! If true, output arrays on numout 183 REAL(wp), INTENT(in ) :: pdam ! depth of the inflection point 184 REAL(wp), INTENT(in ) :: pwam ! width of inflection 185 REAL(wp), INTENT(in ) :: pbot ! bottom value (0<pbot<= 1) 186 REAL(wp), INTENT(in ), DIMENSION(jpk) :: pdep ! depth of the gridpoint (T, U, V, F) 187 REAL(wp), INTENT(inout), DIMENSION(jpk) :: pah ! adimensional vertical profile 188 !! 189 INTEGER :: jk ! dummy loop indices 190 REAL(wp) :: zm00, zm01, zmhb, zmhs ! temporary scalars 191 !!---------------------------------------------------------------------- 192 193 zm00 = TANH( ( pdam - gdept_1d(1 ) ) / pwam ) 194 zm01 = TANH( ( pdam - gdept_1d(jpkm1) ) / pwam ) 195 zmhs = zm00 / zm01 196 zmhb = ( 1.e0 - pbot ) / ( 1.e0 - zmhs ) / zm01 197 198 DO jk = 1, jpk 199 pah(jk) = 1.e0 + zmhb * ( zm00 - TANH( ( pdam - pdep(jk) ) / pwam ) ) 200 END DO 201 202 IF(lwp .AND. ld_print ) THEN ! Control print 203 WRITE(numout,*) 204 WRITE(numout,*) ' ahm profile : ' 205 WRITE(numout,*) 206 WRITE(numout,'(" jk ahm "," depth t-level " )') 207 DO jk = 1, jpk 208 WRITE(numout,'(i6,2f12.4,3x,2f12.4)') jk, pah(jk), pdep(jk) 209 END DO 210 ENDIF 211 ! 212 END SUBROUTINE ldf_zpf_1d 213 214 215 SUBROUTINE ldf_zpf_1d_3d( ld_print, pdam, pwam, pbot, pdep, pah ) 216 !!---------------------------------------------------------------------- 217 !! *** ROUTINE ldf_zpf *** 218 !! 219 !! ** Purpose : vertical adimensional profile for eddy coefficient 220 !! 221 !! ** Method : 1D eddy viscosity coefficients ( depth ) 222 !!---------------------------------------------------------------------- 223 LOGICAL , INTENT(in ) :: ld_print ! If true, output arrays on numout 224 REAL(wp), INTENT(in ) :: pdam ! depth of the inflection point 225 REAL(wp), INTENT(in ) :: pwam ! width of inflection 226 REAL(wp), INTENT(in ) :: pbot ! bottom value (0<pbot<= 1) 227 REAL(wp), INTENT(in ), DIMENSION (:) :: pdep ! depth of the gridpoint (T, U, V, F) 228 REAL(wp), INTENT(inout), DIMENSION (:,:,:) :: pah ! adimensional vertical profile 229 !! 230 INTEGER :: jk ! dummy loop indices 231 REAL(wp) :: zm00, zm01, zmhb, zmhs, zcf ! temporary scalars 232 !!---------------------------------------------------------------------- 233 234 zm00 = TANH( ( pdam - gdept_1d(1 ) ) / pwam ) 235 zm01 = TANH( ( pdam - gdept_1d(jpkm1) ) / pwam ) 236 zmhs = zm00 / zm01 237 zmhb = ( 1.e0 - pbot ) / ( 1.e0 - zmhs ) / zm01 238 239 DO jk = 1, jpk 240 zcf = 1.e0 + zmhb * ( zm00 - TANH( ( pdam - pdep(jk) ) / pwam ) ) 241 pah(:,:,jk) = zcf 242 END DO 243 244 IF(lwp .AND. ld_print ) THEN ! Control print 245 WRITE(numout,*) 246 WRITE(numout,*) ' ahm profile : ' 247 WRITE(numout,*) 248 WRITE(numout,'(" jk ahm "," depth t-level " )') 249 DO jk = 1, jpk 250 WRITE(numout,'(i6,2f12.4,3x,2f12.4)') jk, pah(1,1,jk), pdep(jk) 251 END DO 252 ENDIF 253 ! 254 END SUBROUTINE ldf_zpf_1d_3d 255 256 257 SUBROUTINE ldf_zpf_3d( ld_print, pdam, pwam, pbot, pdep, pah ) 258 !!---------------------------------------------------------------------- 259 !! *** ROUTINE ldf_zpf *** 260 !! 261 !! ** Purpose : vertical adimensional profile for eddy coefficient 262 !! 263 !! ** Method : 3D for partial step or s-coordinate 264 !!---------------------------------------------------------------------- 265 LOGICAL , INTENT(in ) :: ld_print ! If true, output arrays on numout 266 REAL(wp), INTENT(in ) :: pdam ! depth of the inflection point 267 REAL(wp), INTENT(in ) :: pwam ! width of inflection 268 REAL(wp), INTENT(in ) :: pbot ! bottom value (0<pbot<= 1) 269 REAL(wp), INTENT(in ), DIMENSION (:,:,:) :: pdep ! dep of the gridpoint (T, U, V, F) 270 REAL(wp), INTENT(inout), DIMENSION (:,:,:) :: pah ! adimensional vertical profile 271 !! 272 INTEGER :: jk ! dummy loop indices 273 REAL(wp) :: zm00, zm01, zmhb, zmhs ! temporary scalars 274 !!---------------------------------------------------------------------- 275 276 zm00 = TANH( ( pdam - gdept_1d(1 ) ) / pwam ) 277 zm01 = TANH( ( pdam - gdept_1d(jpkm1) ) / pwam ) 278 zmhs = zm00 / zm01 279 zmhb = ( 1.e0 - pbot ) / ( 1.e0 - zmhs ) / zm01 280 281 DO jk = 1, jpk 282 pah(:,:,jk) = 1.e0 + zmhb * ( zm00 - TANH( ( pdam - pdep(:,:,jk) ) / pwam ) ) 283 END DO 284 285 IF(lwp .AND. ld_print ) THEN ! Control print 286 WRITE(numout,*) 287 WRITE(numout,*) ' ahm profile : ' 288 WRITE(numout,*) 289 WRITE(numout,'(" jk ahm "," depth t-level " )') 290 DO jk = 1, jpk 291 WRITE(numout,'(i6,2f12.4,3x,2f12.4)') jk, pah(1,1,jk), pdep(1,1,jk) 292 END DO 293 ENDIF 294 ! 295 END SUBROUTINE ldf_zpf_3d 223 224 SUBROUTINE ldf_dyn( kt ) 225 !!---------------------------------------------------------------------- 226 !! *** ROUTINE ldf_dyn *** 227 !! 228 !! ** Purpose : update at kt the momentum lateral mixing coeff. (ahmt and ahmf) 229 !! 230 !! ** Method : time varying eddy viscosity coefficients: 231 !! 232 !! nn_ahm_ijk_t = 31 ahmt, ahmf = F(i,j,k,t) = F(local velocity) 233 !! ( |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 laplacian 235 !! 236 !! ** action : ahmt, ahmf update at each time step 237 !!---------------------------------------------------------------------- 238 INTEGER, INTENT(in) :: kt ! time step index 239 ! 240 INTEGER :: ji, jj, jk ! dummy loop indices 241 REAL(wp) :: zu2pv2_ij_p1, zu2pv2_ij, zu2pv2_ij_m1, zetmax, zefmax ! local scalar 242 !!---------------------------------------------------------------------- 243 ! 244 IF( nn_timing == 1 ) CALL timing_start('ldf_dyn') 245 ! 246 SELECT CASE( nn_ahm_ijk_t ) !== Eddy vicosity coefficients ==! 247 ! 248 CASE( 31 ) !== time varying 3D field ==! = F( local velocity ) 249 ! 250 IF( ln_dynldf_lap ) THEN ! laplacian operator : |u| e /12 = |u/144| e 251 DO jk = 1, jpkm1 252 DO jj = 2, jpjm1 253 DO ji = fs_2, fs_jpim1 254 zu2pv2_ij_p1 = ub(ji ,jj+1,jk) * ub(ji ,jj+1,jk) + vb(ji+1,jj ,jk) * vb(ji+1,jj ,jk) 255 zu2pv2_ij = ub(ji ,jj ,jk) * ub(ji ,jj ,jk) + vb(ji ,jj ,jk) * vb(ji ,jj ,jk) 256 zu2pv2_ij_m1 = ub(ji-1,jj ,jk) * ub(ji-1,jj ,jk) + vb(ji ,jj-1,jk) * vb(ji ,jj-1,jk) 257 zetmax = MAX( e1t(ji,jj) , e2t(ji,jj) ) 258 zefmax = MAX( e1f(ji,jj) , e2f(ji,jj) ) 259 ahmt(ji,jj,jk) = SQRT( (zu2pv2_ij + zu2pv2_ij_m1) * r1_288 ) * zetmax * tmask(ji,jj,jk) ! 288= 12*12 * 2 260 ahmf(ji,jj,jk) = SQRT( (zu2pv2_ij + zu2pv2_ij_p1) * r1_288 ) * zefmax * fmask(ji,jj,jk) 261 END DO 262 END DO 263 END DO 264 ELSEIF( ln_dynldf_blp ) THEN ! bilaplacian operator : sqrt( |u| e^3 /12 ) = sqrt( |u/144| e ) * e 265 DO jk = 1, jpkm1 266 DO jj = 2, jpjm1 267 DO ji = fs_2, fs_jpim1 268 zu2pv2_ij_p1 = ub(ji ,jj+1,jk) * ub(ji ,jj+1,jk) + vb(ji+1,jj ,jk) * vb(ji+1,jj ,jk) 269 zu2pv2_ij = ub(ji ,jj ,jk) * ub(ji ,jj ,jk) + vb(ji ,jj ,jk) * vb(ji ,jj ,jk) 270 zu2pv2_ij_m1 = ub(ji-1,jj ,jk) * ub(ji-1,jj ,jk) + vb(ji ,jj-1,jk) * vb(ji ,jj-1,jk) 271 zetmax = MAX( e1t(ji,jj) , e2t(ji,jj) ) 272 zefmax = MAX( e1f(ji,jj) , e2f(ji,jj) ) 273 ahmt(ji,jj,jk) = SQRT( SQRT( (zu2pv2_ij + zu2pv2_ij_m1) * r1_288 ) * zetmax ) * zetmax * tmask(ji,jj,jk) 274 ahmf(ji,jj,jk) = SQRT( SQRT( (zu2pv2_ij + zu2pv2_ij_p1) * r1_288 ) * zefmax ) * zefmax * fmask(ji,jj,jk) 275 END DO 276 END DO 277 END DO 278 ENDIF 279 ! 280 CALL lbc_lnk( ahmt, 'T', 1. ) ; CALL lbc_lnk( ahmf, 'F', 1. ) 281 ! 282 END SELECT 283 ! 284 CALL iom_put( "ahmt_2d", ahmt(:,:,1) ) ! surface u-eddy diffusivity coeff. 285 CALL iom_put( "ahmf_2d", ahmf(:,:,1) ) ! surface v-eddy diffusivity coeff. 286 CALL iom_put( "ahmt_3d", ahmt(:,:,:) ) ! 3D u-eddy diffusivity coeff. 287 CALL iom_put( "ahmf_3d", ahmf(:,:,:) ) ! 3D v-eddy diffusivity coeff. 288 ! 289 IF( nn_timing == 1 ) CALL timing_stop('ldf_dyn') 290 ! 291 END SUBROUTINE ldf_dyn 296 292 297 293 !!======================================================================
Note: See TracChangeset
for help on using the changeset viewer.
