Changeset 5948 for branches/2014/dev_r4650_UKMO12_CFL_diags_take2/NEMOGCM/NEMO/OPA_SRC/LDF/ldfdyn.F90
- Timestamp:
- 2015-11-30T11:47:24+01:00 (8 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2014/dev_r4650_UKMO12_CFL_diags_take2/NEMOGCM/NEMO/OPA_SRC/LDF/ldfdyn.F90
r5947 r5948 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 54 # include "domzgr_substitute.h90" 38 !!---------------------------------------------------------------------- 39 !! NEMO/OPA 3.3 , NEMO Consortium (2010) 55 # include "vectopt_loop_substitute.h90" 56 !!---------------------------------------------------------------------- 57 !! NEMO/OPA 3.7 , NEMO Consortium (2014) 40 58 !! $Id$ 41 59 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) … … 49 67 !! ** Purpose : set the horizontal ocean dynamics physics 50 68 !! 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 69 !! ** Method : the eddy viscosity coef. specification depends on: 70 !! - the operator: 71 !! ln_dynldf_lap = T laplacian operator 72 !! ln_dynldf_blp = T bilaplacian operator 73 !! - the parameter nn_ahm_ijk_t: 74 !! nn_ahm_ijk_t = 0 => = constant 75 !! = 10 => = F(z) : = constant with a reduction of 1/4 with depth 76 !! =-20 => = F(i,j) = shape read in 'eddy_viscosity.nc' file 77 !! = 20 = F(i,j) = F(e1,e2) or F(e1^3,e2^3) (lap or bilap case) 78 !! =-30 => = F(i,j,k) = shape read in 'eddy_viscosity.nc' file 79 !! = 30 = F(i,j,k) = 2D (case 20) + decrease with depth (case 10) 80 !! = 31 = F(i,j,k,t) = F(local velocity) ( |u|e /12 laplacian operator 81 !! or |u|e^3/12 bilaplacian operator ) 82 !!---------------------------------------------------------------------- 83 INTEGER :: jk ! dummy loop indices 84 INTEGER :: ierr, inum, ios ! local integer 85 REAL(wp) :: zah0 ! local scalar 86 ! 87 NAMELIST/namdyn_ldf/ ln_dynldf_lap, ln_dynldf_blp, & 88 & ln_dynldf_lev, ln_dynldf_hor, ln_dynldf_iso, & 89 & nn_ahm_ijk_t , rn_ahm_0, rn_ahm_b, rn_bhm_0 90 !!---------------------------------------------------------------------- 91 ! 76 92 REWIND( numnam_ref ) ! Namelist namdyn_ldf in reference namelist : Lateral physics 77 93 READ ( numnam_ref, namdyn_ldf, IOSTAT = ios, ERR = 901) … … 88 104 WRITE(numout,*) '~~~~~~~' 89 105 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 106 ! 107 WRITE(numout,*) ' type :' 108 WRITE(numout,*) ' laplacian operator ln_dynldf_lap = ', ln_dynldf_lap 109 WRITE(numout,*) ' bilaplacian operator ln_dynldf_blp = ', ln_dynldf_blp 110 ! 111 WRITE(numout,*) ' direction of action :' 112 WRITE(numout,*) ' iso-level ln_dynldf_lev = ', ln_dynldf_lev 113 WRITE(numout,*) ' horizontal (geopotential) ln_dynldf_hor = ', ln_dynldf_hor 114 WRITE(numout,*) ' iso-neutral ln_dynldf_iso = ', ln_dynldf_iso 115 ! 116 WRITE(numout,*) ' coefficients :' 117 WRITE(numout,*) ' type of time-space variation nn_ahm_ijk_t = ', nn_ahm_ijk_t 118 WRITE(numout,*) ' lateral laplacian eddy viscosity rn_ahm_0_lap = ', rn_ahm_0, ' m2/s' 119 WRITE(numout,*) ' background viscosity (iso case) rn_ahm_b = ', rn_ahm_b, ' m2/s' 120 WRITE(numout,*) ' lateral bilaplacian eddy viscosity rn_ahm_0_blp = ', rn_bhm_0, ' m4/s' 121 ENDIF 122 123 ! ! Parameter control 124 IF( .NOT.ln_dynldf_lap .AND. .NOT.ln_dynldf_blp ) THEN 125 IF(lwp) WRITE(numout,*) ' No viscous operator selected. ahmt and ahmf are not allocated' 126 l_ldfdyn_time = .FALSE. 127 RETURN 128 ENDIF 129 ! 130 IF( ln_dynldf_blp .AND. ln_dynldf_iso ) THEN ! iso-neutral bilaplacian not implemented 131 CALL ctl_stop( 'dyn_ldf_init: iso-neutral bilaplacian not coded yet' ) 132 ENDIF 133 134 ! ... Space/Time variation of eddy coefficients 135 ! ! allocate the ahm arrays 136 ALLOCATE( ahmt(jpi,jpj,jpk) , ahmf(jpi,jpj,jpk) , STAT=ierr ) 137 IF( ierr /= 0 ) CALL ctl_stop( 'STOP', 'ldf_dyn_init: failed to allocate arrays') 138 ! 139 ahmt(:,:,jpk) = 0._wp ! last level always 0 140 ahmf(:,:,jpk) = 0._wp 141 ! 142 ! ! value of eddy mixing coef. 143 IF ( ln_dynldf_lap ) THEN ; zah0 = rn_ahm_0 ! laplacian operator 144 ELSEIF( ln_dynldf_blp ) THEN ; zah0 = ABS( rn_bhm_0 ) ! bilaplacian operator 145 ELSE ! NO viscous operator 146 CALL ctl_warn( 'ldf_dyn_init: No lateral viscous operator used ' ) 147 ENDIF 148 ! 149 l_ldfdyn_time = .FALSE. ! no time variation except in case defined below 150 ! 151 IF( ln_dynldf_lap .OR. ln_dynldf_blp ) THEN ! only if a lateral diffusion operator is used 152 ! 153 SELECT CASE( nn_ahm_ijk_t ) ! Specification of space time variations of ahmt, ahmf 154 ! 155 CASE( 0 ) !== constant ==! 156 IF(lwp) WRITE(numout,*) ' momentum mixing coef. = constant ' 157 ahmt(:,:,:) = zah0 * tmask(:,:,:) 158 ahmf(:,:,:) = zah0 * fmask(:,:,:) 159 ! 160 CASE( 10 ) !== fixed profile ==! 161 IF(lwp) WRITE(numout,*) ' momentum mixing coef. = F( depth )' 162 ahmt(:,:,1) = zah0 * tmask(:,:,1) ! constant surface value 163 ahmf(:,:,1) = zah0 * fmask(:,:,1) 164 CALL ldf_c1d( 'DYN', r1_4, ahmt(:,:,1), ahmf(:,:,1), ahmt, ahmf ) 165 ! 166 CASE ( -20 ) !== fixed horizontal shape read in file ==! 167 IF(lwp) WRITE(numout,*) ' momentum mixing coef. = F(i,j) read in eddy_viscosity.nc file' 168 CALL iom_open( 'eddy_viscosity_2D.nc', inum ) 169 CALL iom_get ( inum, jpdom_data, 'ahmt_2d', ahmt(:,:,1) ) 170 CALL iom_get ( inum, jpdom_data, 'ahmf_2d', ahmf(:,:,1) ) 171 CALL iom_close( inum ) 172 !!gm Question : info for LAP or BLP case to take into account the SQRT in the bilaplacian case ??? 173 !! do we introduce a scaling by the max value of the array, and then multiply by zah0 ???? 174 !! better: check that the max is <=1 i.e. it is a shape from 0 to 1, not a coef that has physical dimension 175 DO jk = 2, jpkm1 176 ahmt(:,:,jk) = ahmt(:,:,1) * tmask(:,:,jk) 177 ahmf(:,:,jk) = ahmf(:,:,1) * fmask(:,:,jk) 178 END DO 179 ! 180 CASE( 20 ) !== fixed horizontal shape ==! 181 IF(lwp) WRITE(numout,*) ' momentum mixing coef. = F( e1, e2 ) or F( e1^3, e2^3 ) (lap. or blp. case)' 182 IF( ln_dynldf_lap ) CALL ldf_c2d( 'DYN', 'LAP', zah0, ahmt, ahmf ) ! surface value proportional to scale factor 183 IF( ln_dynldf_blp ) CALL ldf_c2d( 'DYN', 'BLP', zah0, ahmt, ahmf ) ! surface value proportional to scale factor^3 184 ! 185 CASE( -30 ) !== fixed 3D shape read in file ==! 186 IF(lwp) WRITE(numout,*) ' momentum mixing coef. = F(i,j,k) read in eddy_diffusivity_3D.nc file' 187 CALL iom_open( 'eddy_viscosity_3D.nc', inum ) 188 CALL iom_get ( inum, jpdom_data, 'ahmt_3d', ahmt ) 189 CALL iom_get ( inum, jpdom_data, 'ahmf_3d', ahmf ) 190 CALL iom_close( inum ) 191 !!gm Question : info for LAP or BLP case to take into account the SQRT in the bilaplacian case ???? 192 !! do we introduce a scaling by the max value of the array, and then multiply by zah0 ???? 193 DO jk = 1, jpkm1 194 ahmt(:,:,jk) = ahmt(:,:,jk) * tmask(:,:,jk) 195 ahmf(:,:,jk) = ahmf(:,:,jk) * fmask(:,:,jk) 196 END DO 197 ! 198 CASE( 30 ) !== fixed 3D shape ==! 199 IF(lwp) WRITE(numout,*) ' momentum mixing coef. = F( latitude, longitude, depth )' 200 IF( ln_dynldf_lap ) CALL ldf_c2d( 'DYN', 'LAP', zah0, ahmt, ahmf ) ! surface value proportional to scale factor 201 IF( ln_dynldf_blp ) CALL ldf_c2d( 'DYN', 'BLP', zah0, ahmt, ahmf ) ! surface value proportional to scale factor 202 ! ! reduction with depth 203 CALL ldf_c1d( 'DYN', r1_4, ahmt(:,:,1), ahmf(:,:,1), ahmt, ahmf ) 204 ! 205 CASE( 31 ) !== time varying 3D field ==! 206 IF(lwp) WRITE(numout,*) ' momentum mixing coef. = F( latitude, longitude, depth , time )' 207 IF(lwp) WRITE(numout,*) ' proportional to the velocity : |u|e/12 or |u|e^3/12' 208 ! 209 l_ldfdyn_time = .TRUE. ! will be calculated by call to ldf_dyn routine in step.F90 210 ! 211 CASE DEFAULT 212 CALL ctl_stop('ldf_dyn_init: wrong choice for nn_ahm_ijk_t, the type of space-time variation of ahm') 213 END SELECT 214 ! 215 IF( ln_dynldf_blp .AND. .NOT. l_ldfdyn_time ) THEN ! bilapcian and no time variation: 216 ahmt(:,:,:) = SQRT( ahmt(:,:,:) ) ! take the square root of the coefficient 217 ahmf(:,:,:) = SQRT( ahmf(:,:,:) ) 218 ENDIF 219 ! 220 ENDIF 162 221 ! 163 222 END SUBROUTINE ldf_dyn_init 164 223 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 224 225 SUBROUTINE ldf_dyn( kt ) 226 !!---------------------------------------------------------------------- 227 !! *** ROUTINE ldf_dyn *** 228 !! 229 !! ** Purpose : update at kt the momentum lateral mixing coeff. (ahmt and ahmf) 230 !! 231 !! ** Method : time varying eddy viscosity coefficients: 232 !! 233 !! nn_ahm_ijk_t = 31 ahmt, ahmf = F(i,j,k,t) = F(local velocity) 234 !! ( |u|e /12 or |u|e^3/12 for laplacian or bilaplacian operator ) 235 !! BLP case : sqrt of the eddy coef, since bilaplacian is en re-entrant laplacian 236 !! 237 !! ** action : ahmt, ahmf update at each time step 238 !!---------------------------------------------------------------------- 239 INTEGER, INTENT(in) :: kt ! time step index 240 ! 241 INTEGER :: ji, jj, jk ! dummy loop indices 242 REAL(wp) :: zu2pv2_ij_p1, zu2pv2_ij, zu2pv2_ij_m1, zetmax, zefmax ! local scalar 243 !!---------------------------------------------------------------------- 244 ! 245 IF( nn_timing == 1 ) CALL timing_start('ldf_dyn') 246 ! 247 SELECT CASE( nn_ahm_ijk_t ) !== Eddy vicosity coefficients ==! 248 ! 249 CASE( 31 ) !== time varying 3D field ==! = F( local velocity ) 250 ! 251 IF( ln_dynldf_lap ) THEN ! laplacian operator : |u| e /12 = |u/144| e 252 DO jk = 1, jpkm1 253 DO jj = 2, jpjm1 254 DO ji = fs_2, fs_jpim1 255 zu2pv2_ij_p1 = ub(ji ,jj+1,jk) * ub(ji ,jj+1,jk) + vb(ji+1,jj ,jk) * vb(ji+1,jj ,jk) 256 zu2pv2_ij = ub(ji ,jj ,jk) * ub(ji ,jj ,jk) + vb(ji ,jj ,jk) * vb(ji ,jj ,jk) 257 zu2pv2_ij_m1 = ub(ji-1,jj ,jk) * ub(ji-1,jj ,jk) + vb(ji ,jj-1,jk) * vb(ji ,jj-1,jk) 258 zetmax = MAX( e1t(ji,jj) , e2t(ji,jj) ) 259 zefmax = MAX( e1f(ji,jj) , e2f(ji,jj) ) 260 ahmt(ji,jj,jk) = SQRT( (zu2pv2_ij + zu2pv2_ij_m1) * r1_288 ) * zetmax * tmask(ji,jj,jk) ! 288= 12*12 * 2 261 ahmf(ji,jj,jk) = SQRT( (zu2pv2_ij + zu2pv2_ij_p1) * r1_288 ) * zefmax * fmask(ji,jj,jk) 262 END DO 263 END DO 264 END DO 265 ELSEIF( ln_dynldf_blp ) THEN ! bilaplacian operator : sqrt( |u| e^3 /12 ) = sqrt( |u/144| e ) * e 266 DO jk = 1, jpkm1 267 DO jj = 2, jpjm1 268 DO ji = fs_2, fs_jpim1 269 zu2pv2_ij_p1 = ub(ji ,jj+1,jk) * ub(ji ,jj+1,jk) + vb(ji+1,jj ,jk) * vb(ji+1,jj ,jk) 270 zu2pv2_ij = ub(ji ,jj ,jk) * ub(ji ,jj ,jk) + vb(ji ,jj ,jk) * vb(ji ,jj ,jk) 271 zu2pv2_ij_m1 = ub(ji-1,jj ,jk) * ub(ji-1,jj ,jk) + vb(ji ,jj-1,jk) * vb(ji ,jj-1,jk) 272 zetmax = MAX( e1t(ji,jj) , e2t(ji,jj) ) 273 zefmax = MAX( e1f(ji,jj) , e2f(ji,jj) ) 274 ahmt(ji,jj,jk) = SQRT( SQRT( (zu2pv2_ij + zu2pv2_ij_m1) * r1_288 ) * zetmax ) * zetmax * tmask(ji,jj,jk) 275 ahmf(ji,jj,jk) = SQRT( SQRT( (zu2pv2_ij + zu2pv2_ij_p1) * r1_288 ) * zefmax ) * zefmax * fmask(ji,jj,jk) 276 END DO 277 END DO 278 END DO 279 ENDIF 280 ! 281 CALL lbc_lnk( ahmt, 'T', 1. ) ; CALL lbc_lnk( ahmf, 'F', 1. ) 282 ! 283 END SELECT 284 ! 285 CALL iom_put( "ahmt_2d", ahmt(:,:,1) ) ! surface u-eddy diffusivity coeff. 286 CALL iom_put( "ahmf_2d", ahmf(:,:,1) ) ! surface v-eddy diffusivity coeff. 287 CALL iom_put( "ahmt_3d", ahmt(:,:,:) ) ! 3D u-eddy diffusivity coeff. 288 CALL iom_put( "ahmf_3d", ahmf(:,:,:) ) ! 3D v-eddy diffusivity coeff. 289 ! 290 IF( nn_timing == 1 ) CALL timing_stop('ldf_dyn') 291 ! 292 END SUBROUTINE ldf_dyn 296 293 297 294 !!======================================================================
Note: See TracChangeset
for help on using the changeset viewer.
