- Timestamp:
- 2018-04-23T10:44:07+02:00 (6 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/LDF/ldftra.F90
r9169 r9490 50 50 LOGICAL , PUBLIC :: ln_traldf_hor !: horizontal (geopotential) direction 51 51 ! LOGICAL , PUBLIC :: ln_traldf_iso !: iso-neutral direction (see ldfslp) 52 ! != iso-neutral options =! 52 53 ! LOGICAL , PUBLIC :: ln_traldf_triad !: griffies triad scheme (see ldfslp) 53 54 LOGICAL , PUBLIC :: ln_traldf_msc !: Method of Stabilizing Correction … … 58 59 ! != Coefficients =! 59 60 INTEGER , PUBLIC :: nn_aht_ijk_t !: choice of time & space variations of the lateral eddy diffusivity coef. 60 REAL(wp), PUBLIC :: rn_aht_0 !: laplacian lateral eddy diffusivity [m2/s] 61 REAL(wp), PUBLIC :: rn_bht_0 !: bilaplacian lateral eddy diffusivity [m4/s] 62 63 ! !!* Namelist namtra_ldfeiv : eddy induced velocity param. * 61 ! ! time invariant coefficients: aht_0 = 1/2 Ud*Ld (lap case) 62 ! ! bht_0 = 1/12 Ud*Ld^3 (blp case) 63 REAL(wp), PUBLIC :: rn_Ud !: lateral diffusive velocity [m/s] 64 REAL(wp), PUBLIC :: rn_Ld !: lateral diffusive length [m] 65 66 ! !!* Namelist namtra_eiv : eddy induced velocity param. * 64 67 ! != Use/diagnose eiv =! 65 68 LOGICAL , PUBLIC :: ln_ldfeiv !: eddy induced velocity flag … … 67 70 ! != Coefficients =! 68 71 INTEGER , PUBLIC :: nn_aei_ijk_t !: choice of time/space variation of the eiv coeff. 69 REAL(wp), PUBLIC :: rn_aeiv_0 !: eddy induced velocity coefficient [m2/s] 72 REAL(wp), PUBLIC :: rn_Ue !: lateral diffusive velocity [m/s] 73 REAL(wp), PUBLIC :: rn_Le !: lateral diffusive length [m] 70 74 75 ! ! Flag to control the type of lateral diffusive operator 76 INTEGER, PARAMETER, PUBLIC :: np_ERROR =-10 ! error in specification of lateral diffusion 77 INTEGER, PARAMETER, PUBLIC :: np_no_ldf = 00 ! without operator (i.e. no lateral diffusive trend) 78 ! !! laplacian ! bilaplacian ! 79 INTEGER, PARAMETER, PUBLIC :: np_lap = 10 , np_blp = 20 ! iso-level operator 80 INTEGER, PARAMETER, PUBLIC :: np_lap_i = 11 , np_blp_i = 21 ! standard iso-neutral or geopotential operator 81 INTEGER, PARAMETER, PUBLIC :: np_lap_it = 12 , np_blp_it = 22 ! triad iso-neutral or geopotential operator 82 83 INTEGER , PUBLIC :: nldf_tra = 0 !: type of lateral diffusion used defined from ln_traldf_... (namlist logicals) 71 84 LOGICAL , PUBLIC :: l_ldftra_time = .FALSE. !: flag for time variation of the lateral eddy diffusivity coef. 72 LOGICAL , PUBLIC :: l_ldfeiv_time = .FALSE. ! flag for time variation of the eiv coef.85 LOGICAL , PUBLIC :: l_ldfeiv_time = .FALSE. !: flag for time variation of the eiv coef. 73 86 74 87 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: ahtu, ahtv !: eddy diffusivity coef. at U- and V-points [m2/s] 75 88 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: aeiu, aeiv !: eddy induced velocity coeff. [m2/s] 76 89 90 REAL(wp) :: aht0, aei0 ! constant eddy coefficients (deduced from namelist values) [m2/s] 91 REAL(wp) :: r1_2 = 0.5_wp ! =1/2 77 92 REAL(wp) :: r1_4 = 0.25_wp ! =1/4 78 93 REAL(wp) :: r1_12 = 1._wp / 12._wp ! =1/12 … … 108 123 !! =-30 => = F(i,j,k) = shape read in 'eddy_diffusivity.nc' file 109 124 !! = 30 = F(i,j,k) = 2D (case 20) + decrease with depth (case 10) 110 !! = 31 = F(i,j,k,t) = F(local velocity) ( |u|e /12laplacian operator111 !! or |u|e^3/12bilaplacian operator )125 !! = 31 = F(i,j,k,t) = F(local velocity) ( 1/2 |u|e laplacian operator 126 !! or 1/12 |u|e^3 bilaplacian operator ) 112 127 !! * initialisation of the eddy induced velocity coefficient by a call to ldf_eiv_init 113 128 !! 114 !! ** action : ahtu, ahtv initialized once for all or l_ldftra_time set to true 115 !! aeiu, aeiv initialized once for all or l_ldfeiv_time set to true 116 !!---------------------------------------------------------------------- 117 INTEGER :: jk ! dummy loop indices 118 INTEGER :: ierr, inum, ios ! local integer 119 REAL(wp) :: zah0 ! local scalar 129 !! ** action : ahtu, ahtv initialized one for all or l_ldftra_time set to true 130 !! aeiu, aeiv initialized one for all or l_ldfeiv_time set to true 131 !!---------------------------------------------------------------------- 132 INTEGER :: jk ! dummy loop indices 133 INTEGER :: ioptio, ierr, inum, ios, inn ! local integer 134 REAL(wp) :: zah_max, zUfac ! - - 135 CHARACTER(len=5) :: cl_Units ! units (m2/s or m4/s) 120 136 !! 121 137 NAMELIST/namtra_ldf/ ln_traldf_NONE, ln_traldf_lap , ln_traldf_blp , & ! type of operator … … 123 139 & ln_traldf_iso , ln_traldf_msc , rn_slpmax , & ! option for iso-neutral operator 124 140 & ln_triad_iso , ln_botmix_triad, rn_sw_triad , & ! option for triad operator 125 & rn_aht_0 , rn_bht_0 , nn_aht_ijk_t! lateral eddy coefficient141 & nn_aht_ijk_t , rn_Ud , rn_Ld ! lateral eddy coefficient 126 142 !!---------------------------------------------------------------------- 127 143 ! 128 144 IF(lwp) THEN ! control print 129 145 WRITE(numout,*) 130 WRITE(numout,*) 'ldf_tra_init : lateral tracer physics'146 WRITE(numout,*) 'ldf_tra_init : lateral tracer diffusion' 131 147 WRITE(numout,*) '~~~~~~~~~~~~ ' 132 148 ENDIF 149 133 150 ! 134 151 ! Choice of lateral tracer physics … … 154 171 WRITE(numout,*) ' iso-neutral Madec operator ln_traldf_iso = ', ln_traldf_iso 155 172 WRITE(numout,*) ' iso-neutral triad operator ln_traldf_triad = ', ln_traldf_triad 156 WRITE(numout,*) ' iso-neutral (Method of Stab. Corr.)ln_traldf_msc = ', ln_traldf_msc173 WRITE(numout,*) ' use the Method of Stab. Correction ln_traldf_msc = ', ln_traldf_msc 157 174 WRITE(numout,*) ' maximum isoppycnal slope rn_slpmax = ', rn_slpmax 158 175 WRITE(numout,*) ' pure lateral mixing in ML ln_triad_iso = ', ln_triad_iso … … 160 177 WRITE(numout,*) ' lateral mixing on bottom ln_botmix_triad = ', ln_botmix_triad 161 178 WRITE(numout,*) ' coefficients :' 162 WRITE(numout,*) ' lateral eddy diffusivity (lap case) rn_aht_0 = ', rn_aht_0163 WRITE(numout,*) ' lateral eddy diffusivity (bilap case) rn_bht_0 = ', rn_bht_0164 179 WRITE(numout,*) ' type of time-space variation nn_aht_ijk_t = ', nn_aht_ijk_t 165 ENDIF 166 ! 167 ! ! Parameter control 168 ! 169 IF( ln_traldf_NONE ) THEN 170 IF(lwp) WRITE(numout,*) ' ==>>> No diffusive operator selected. ahtu and ahtv are not allocated' 171 l_ldftra_time = .FALSE. 172 RETURN 173 ENDIF 180 WRITE(numout,*) ' lateral diffusive velocity (if cst) rn_Ud = ', rn_Ud, ' m/s' 181 WRITE(numout,*) ' lateral diffusive length (if cst) rn_Ld = ', rn_Ld, ' m' 182 ENDIF 183 ! 184 ! 185 ! Operator and its acting direction (set nldf_tra) 186 ! ================================= 187 ! 188 nldf_tra = np_ERROR 189 ioptio = 0 190 IF( ln_traldf_NONE ) THEN ; nldf_tra = np_no_ldf ; ioptio = ioptio + 1 ; ENDIF 191 IF( ln_traldf_lap ) THEN ; ioptio = ioptio + 1 ; ENDIF 192 IF( ln_traldf_blp ) THEN ; ioptio = ioptio + 1 ; ENDIF 193 IF( ioptio /= 1 ) CALL ctl_stop( 'tra_ldf_init: use ONE of the 3 operator options (NONE/lap/blp)' ) 194 ! 195 IF( .NOT.ln_traldf_NONE ) THEN !== direction ==>> type of operator ==! 196 ioptio = 0 197 IF( ln_traldf_lev ) ioptio = ioptio + 1 198 IF( ln_traldf_hor ) ioptio = ioptio + 1 199 IF( ln_traldf_iso ) ioptio = ioptio + 1 200 IF( ioptio /= 1 ) CALL ctl_stop( 'tra_ldf_init: use ONE direction (level/hor/iso)' ) 201 ! 202 ! ! defined the type of lateral diffusion from ln_traldf_... logicals 203 ierr = 0 204 IF ( ln_traldf_lap ) THEN ! laplacian operator 205 IF ( ln_zco ) THEN ! z-coordinate 206 IF ( ln_traldf_lev ) nldf_tra = np_lap ! iso-level = horizontal (no rotation) 207 IF ( ln_traldf_hor ) nldf_tra = np_lap ! iso-level = horizontal (no rotation) 208 IF ( ln_traldf_iso ) nldf_tra = np_lap_i ! iso-neutral: standard ( rotation) 209 IF ( ln_traldf_triad ) nldf_tra = np_lap_it ! iso-neutral: triad ( rotation) 210 ENDIF 211 IF ( ln_zps ) THEN ! z-coordinate with partial step 212 IF ( ln_traldf_lev ) ierr = 1 ! iso-level not allowed 213 IF ( ln_traldf_hor ) nldf_tra = np_lap ! horizontal (no rotation) 214 IF ( ln_traldf_iso ) nldf_tra = np_lap_i ! iso-neutral: standard (rotation) 215 IF ( ln_traldf_triad ) nldf_tra = np_lap_it ! iso-neutral: triad (rotation) 216 ENDIF 217 IF ( ln_sco ) THEN ! s-coordinate 218 IF ( ln_traldf_lev ) nldf_tra = np_lap ! iso-level (no rotation) 219 IF ( ln_traldf_hor ) nldf_tra = np_lap_i ! horizontal ( rotation) 220 IF ( ln_traldf_iso ) nldf_tra = np_lap_i ! iso-neutral: standard ( rotation) 221 IF ( ln_traldf_triad ) nldf_tra = np_lap_it ! iso-neutral: triad ( rotation) 222 ENDIF 223 ENDIF 224 ! 225 IF( ln_traldf_blp ) THEN ! bilaplacian operator 226 IF ( ln_zco ) THEN ! z-coordinate 227 IF ( ln_traldf_lev ) nldf_tra = np_blp ! iso-level = horizontal (no rotation) 228 IF ( ln_traldf_hor ) nldf_tra = np_blp ! iso-level = horizontal (no rotation) 229 IF ( ln_traldf_iso ) nldf_tra = np_blp_i ! iso-neutral: standard ( rotation) 230 IF ( ln_traldf_triad ) nldf_tra = np_blp_it ! iso-neutral: triad ( rotation) 231 ENDIF 232 IF ( ln_zps ) THEN ! z-coordinate with partial step 233 IF ( ln_traldf_lev ) ierr = 1 ! iso-level not allowed 234 IF ( ln_traldf_hor ) nldf_tra = np_blp ! horizontal (no rotation) 235 IF ( ln_traldf_iso ) nldf_tra = np_blp_i ! iso-neutral: standard ( rotation) 236 IF ( ln_traldf_triad ) nldf_tra = np_blp_it ! iso-neutral: triad ( rotation) 237 ENDIF 238 IF ( ln_sco ) THEN ! s-coordinate 239 IF ( ln_traldf_lev ) nldf_tra = np_blp ! iso-level (no rotation) 240 IF ( ln_traldf_hor ) nldf_tra = np_blp_it ! horizontal ( rotation) 241 IF ( ln_traldf_iso ) nldf_tra = np_blp_i ! iso-neutral: standard ( rotation) 242 IF ( ln_traldf_triad ) nldf_tra = np_blp_it ! iso-neutral: triad ( rotation) 243 ENDIF 244 ENDIF 245 IF ( ierr == 1 ) CALL ctl_stop( 'iso-level in z-partial step, not allowed' ) 246 ENDIF 247 ! 248 IF( ln_ldfeiv .AND. .NOT.( ln_traldf_iso .OR. ln_traldf_triad ) ) & 249 & CALL ctl_stop( 'ln_ldfeiv=T requires iso-neutral laplacian diffusion' ) 250 IF( ln_isfcav .AND. ln_traldf_triad ) & 251 & CALL ctl_stop( ' ice shelf cavity and traldf_triad not tested' ) 252 ! 253 IF( nldf_tra == np_lap_i .OR. nldf_tra == np_lap_it .OR. & 254 & nldf_tra == np_blp_i .OR. nldf_tra == np_blp_it ) l_ldfslp = .TRUE. ! slope of neutral surfaces required 174 255 ! 175 256 IF( ln_traldf_blp .AND. ( ln_traldf_iso .OR. ln_traldf_triad) ) THEN ! iso-neutral bilaplacian need MSC … … 177 258 ENDIF 178 259 ! 260 IF(lwp) THEN 261 WRITE(numout,*) 262 SELECT CASE( nldf_tra ) 263 CASE( np_no_ldf ) ; WRITE(numout,*) ' ==>>> NO lateral diffusion' 264 CASE( np_lap ) ; WRITE(numout,*) ' ==>>> laplacian iso-level operator' 265 CASE( np_lap_i ) ; WRITE(numout,*) ' ==>>> Rotated laplacian operator (standard)' 266 CASE( np_lap_it ) ; WRITE(numout,*) ' ==>>> Rotated laplacian operator (triad)' 267 CASE( np_blp ) ; WRITE(numout,*) ' ==>>> bilaplacian iso-level operator' 268 CASE( np_blp_i ) ; WRITE(numout,*) ' ==>>> Rotated bilaplacian operator (standard)' 269 CASE( np_blp_it ) ; WRITE(numout,*) ' ==>>> Rotated bilaplacian operator (triad)' 270 END SELECT 271 WRITE(numout,*) 272 ENDIF 273 274 ! 179 275 ! Space/time variation of eddy coefficients 180 276 ! =========================================== 181 ! ! allocate the aht arrays 182 ALLOCATE( ahtu(jpi,jpj,jpk) , ahtv(jpi,jpj,jpk) , STAT=ierr ) 183 IF( ierr /= 0 ) CALL ctl_stop( 'STOP', 'ldf_tra_init: failed to allocate arrays') 184 ! 185 ahtu(:,:,jpk) = 0._wp ! last level always 0 186 ahtv(:,:,jpk) = 0._wp 187 ! 188 ! ! value of eddy mixing coef. 189 IF ( ln_traldf_lap ) THEN ; zah0 = rn_aht_0 ! laplacian operator 190 ELSEIF( ln_traldf_blp ) THEN ; zah0 = ABS( rn_bht_0 ) ! bilaplacian operator 191 ENDIF 192 ! 193 l_ldftra_time = .FALSE. ! no time variation except in case defined below 194 ! 195 IF( ln_traldf_lap .OR. ln_traldf_blp ) THEN ! only if a lateral diffusion operator is used 196 ! 197 SELECT CASE( nn_aht_ijk_t ) ! Specification of space time variations of ehtu, ahtv 277 ! 278 l_ldftra_time = .FALSE. ! no time variation except in case defined below 279 ! 280 IF( ln_traldf_NONE ) THEN !== no explicit diffusive operator ==! 281 ! 282 IF(lwp) WRITE(numout,*) ' ==>>> No diffusive operator selected. ahtu and ahtv are not allocated' 283 RETURN 284 ! 285 ELSE !== a lateral diffusion operator is used ==! 286 ! 287 ! ! allocate the aht arrays 288 ALLOCATE( ahtu(jpi,jpj,jpk) , ahtv(jpi,jpj,jpk) , STAT=ierr ) 289 IF( ierr /= 0 ) CALL ctl_stop( 'STOP', 'ldf_tra_init: failed to allocate arrays') 290 ! 291 ahtu(:,:,jpk) = 0._wp ! last level always 0 292 ahtv(:,:,jpk) = 0._wp 293 !. 294 ! ! value of lap/blp eddy mixing coef. 295 IF( ln_traldf_lap ) THEN ; zUfac = r1_2 *rn_Ud ; inn = 1 ; cl_Units = ' m2/s' ! laplacian 296 ELSEIF( ln_traldf_blp ) THEN ; zUfac = r1_12*rn_Ud ; inn = 3 ; cl_Units = ' m4/s' ! bilaplacian 297 ENDIF 298 aht0 = zUfac * rn_Ld**inn ! mixing coefficient 299 zah_max = zUfac * (ra*rad)**inn ! maximum reachable coefficient (value at the Equator for e1=1 degree) 300 ! 301 ! 302 SELECT CASE( nn_aht_ijk_t ) !* Specification of space-time variations of ahtu, ahtv 198 303 ! 199 304 CASE( 0 ) !== constant ==! 200 IF(lwp) WRITE(numout,*) ' ==>>> tracer mixing coef. = constant = ', rn_aht_0201 ahtu(:,:, :) = zah0 * umask(:,:,:)202 ahtv(:,:, :) = zah0 * vmask(:,:,:)305 IF(lwp) WRITE(numout,*) ' ==>>> eddy diffusivity = constant = ', aht0, cl_Units 306 ahtu(:,:,1:jpkm1) = aht0 307 ahtv(:,:,1:jpkm1) = aht0 203 308 ! 204 309 CASE( 10 ) !== fixed profile ==! 205 IF(lwp) WRITE(numout,*) ' ==>>> tracer mixing coef. = F( depth )' 206 ahtu(:,:,1) = zah0 * umask(:,:,1) ! constant surface value 207 ahtv(:,:,1) = zah0 * vmask(:,:,1) 208 CALL ldf_c1d( 'TRA', r1_4, ahtu(:,:,1), ahtv(:,:,1), ahtu, ahtv ) 209 ! 210 CASE ( -20 ) !== fixed horizontal shape read in file ==! 211 IF(lwp) WRITE(numout,*) ' ==>>> tracer mixing coef. = F(i,j) read in eddy_diffusivity.nc file' 310 IF(lwp) WRITE(numout,*) ' ==>>> eddy diffusivity = F( depth )' 311 IF(lwp) WRITE(numout,*) ' surface eddy diffusivity = constant = ', aht0, cl_Units 312 ahtu(:,:,1) = aht0 ! constant surface value 313 ahtv(:,:,1) = aht0 314 CALL ldf_c1d( 'TRA', ahtu(:,:,1), ahtv(:,:,1), ahtu, ahtv ) 315 ! 316 CASE ( -20 ) !== fixed horizontal shape and magnitude read in file ==! 317 IF(lwp) WRITE(numout,*) ' ==>>> eddy diffusivity = F(i,j) read in eddy_diffusivity.nc file' 212 318 CALL iom_open( 'eddy_diffusivity_2D.nc', inum ) 213 319 CALL iom_get ( inum, jpdom_data, 'ahtu_2D', ahtu(:,:,1) ) … … 215 321 CALL iom_close( inum ) 216 322 DO jk = 2, jpkm1 217 ahtu(:,:,jk) = ahtu(:,:,1) * umask(:,:,jk)218 ahtv(:,:,jk) = ahtv(:,:,1) * vmask(:,:,jk)323 ahtu(:,:,jk) = ahtu(:,:,1) 324 ahtv(:,:,jk) = ahtv(:,:,1) 219 325 END DO 220 326 ! 221 327 CASE( 20 ) !== fixed horizontal shape ==! 222 IF(lwp) WRITE(numout,*) ' ==>>> tracer mixing coef. = F( e1, e2 ) or F( e1^3, e2^3 ) (lap or blp case)' 223 IF( ln_traldf_lap ) CALL ldf_c2d( 'TRA', 'LAP', zah0, ahtu, ahtv ) ! surface value proportional to scale factor 224 IF( ln_traldf_blp ) CALL ldf_c2d( 'TRA', 'BLP', zah0, ahtu, ahtv ) ! surface value proportional to scale factor 328 IF(lwp) WRITE(numout,*) ' ==>>> eddy diffusivity = F( e1, e2 ) or F( e1^3, e2^3 ) (lap or blp case)' 329 IF(lwp) WRITE(numout,*) ' using a fixed diffusive velocity = ', rn_Ud,' m/s and Ld = Max(e1,e2)' 330 IF(lwp) WRITE(numout,*) ' maximum reachable coefficient (at the Equator) = ', zah_max, cl_Units, ' for e1=1°)' 331 CALL ldf_c2d( 'TRA', zUfac , inn , ahtu, ahtv ) ! value proportional to scale factor^inn 225 332 ! 226 333 CASE( 21 ) !== time varying 2D field ==! 227 IF(lwp) WRITE(numout,*) ' ==>>> tracer mixing coef. = F( latitude, longitude, time )' 228 IF(lwp) WRITE(numout,*) ' = F( growth rate of baroclinic instability )' 229 IF(lwp) WRITE(numout,*) ' min value = 0.1 * rn_aht_0' 230 IF(lwp) WRITE(numout,*) ' max value = rn_aht_0 (rn_aeiv_0 if nn_aei_ijk_t=21)' 231 IF(lwp) WRITE(numout,*) ' increased to rn_aht_0 within 20N-20S' 334 IF(lwp) WRITE(numout,*) ' ==>>> eddy diffusivity = F( latitude, longitude, time )' 335 IF(lwp) WRITE(numout,*) ' = F( growth rate of baroclinic instability )' 336 IF(lwp) WRITE(numout,*) ' min value = 0.2 * aht0 (with aht0= 1/2 rn_Ud*rn_Ld)' 337 IF(lwp) WRITE(numout,*) ' max value = aei0 (with aei0=1/2 rn_Ue*Le increased to aht0 within 20N-20S' 232 338 ! 233 339 l_ldftra_time = .TRUE. ! will be calculated by call to ldf_tra routine in step.F90 234 340 ! 235 IF( ln_traldf_blp ) THEN 236 CALL ctl_stop( 'ldf_tra_init: aht=F(growth rate of baroc. insta.) incompatible with bilaplacian operator' ) 237 ENDIF 341 IF( ln_traldf_blp ) CALL ctl_stop( 'ldf_tra_init: aht=F( growth rate of baroc. insta .)', & 342 & ' incompatible with bilaplacian operator' ) 238 343 ! 239 344 CASE( -30 ) !== fixed 3D shape read in file ==! 240 IF(lwp) WRITE(numout,*) ' ==>>> tracer mixing coef.= F(i,j,k) read in eddy_diffusivity.nc file'345 IF(lwp) WRITE(numout,*) ' ==>>> eddy diffusivity = F(i,j,k) read in eddy_diffusivity.nc file' 241 346 CALL iom_open( 'eddy_diffusivity_3D.nc', inum ) 242 347 CALL iom_get ( inum, jpdom_data, 'ahtu_3D', ahtu ) 243 348 CALL iom_get ( inum, jpdom_data, 'ahtv_3D', ahtv ) 244 349 CALL iom_close( inum ) 245 DO jk = 1, jpkm1246 ahtu(:,:,jk) = ahtu(:,:,jk) * umask(:,:,jk)247 ahtv(:,:,jk) = ahtv(:,:,jk) * vmask(:,:,jk)248 END DO249 350 ! 250 351 CASE( 30 ) !== fixed 3D shape ==! 251 IF(lwp) WRITE(numout,*) ' ==>>> tracer mixing coef.= F( latitude, longitude, depth )'252 IF( ln_traldf_lap ) CALL ldf_c2d( 'TRA', 'LAP', zah0, ahtu, ahtv ) ! surface value proportional to scale factor253 IF( ln_traldf_blp ) CALL ldf_c2d( 'TRA', 'BLP', zah0, ahtu, ahtv ) ! surface value proportional to scale factor254 ! ! reduction with depth255 CALL ldf_c1d( 'TRA', r1_4, ahtu(:,:,1), ahtv(:,:,1), ahtu, ahtv )352 IF(lwp) WRITE(numout,*) ' ==>>> eddy diffusivity = F( latitude, longitude, depth )' 353 IF(lwp) WRITE(numout,*) ' using a fixed diffusive velocity = ', rn_Ud,' m/s and Ld = Max(e1,e2)' 354 IF(lwp) WRITE(numout,*) ' maximum reachable coefficient (at the Equator) = ', zah_max, cl_Units, ' for e1=1°)' 355 CALL ldf_c2d( 'TRA', zUfac , inn , ahtu, ahtv ) ! surface value proportional to scale factor^inn 356 CALL ldf_c1d( 'TRA', ahtu(:,:,1), ahtv(:,:,1), ahtu, ahtv ) ! reduction with depth 256 357 ! 257 358 CASE( 31 ) !== time varying 3D field ==! 258 IF(lwp) WRITE(numout,*) ' ==>>> tracer mixing coef.= F( latitude, longitude, depth , time )'259 IF(lwp) WRITE(numout,*) ' proportional to the velocity : |u|e/12 or |u|e^3/12'359 IF(lwp) WRITE(numout,*) ' ==>>> eddy diffusivity = F( latitude, longitude, depth , time )' 360 IF(lwp) WRITE(numout,*) ' proportional to the velocity : 1/2 |u|e or 1/12 |u|e^3' 260 361 ! 261 362 l_ldftra_time = .TRUE. ! will be calculated by call to ldf_tra routine in step.F90 … … 265 366 END SELECT 266 367 ! 267 IF( ln_traldf_blp .AND. .NOT. l_ldftra_time ) THEN 268 ahtu(:,:,:) = SQRT( ahtu(:,:,:) ) 269 ahtv(:,:,:) = SQRT( ahtv(:,:,:) ) 368 IF( .NOT.l_ldftra_time ) THEN !* No time variation 369 IF( ln_traldf_lap ) THEN ! laplacian operator (mask only) 370 ahtu(:,:,1:jpkm1) = ahtu(:,:,1:jpkm1) * umask(:,:,1:jpkm1) 371 ahtv(:,:,1:jpkm1) = ahtv(:,:,1:jpkm1) * vmask(:,:,1:jpkm1) 372 ELSEIF( ln_traldf_blp ) THEN ! bilaplacian operator (square root + mask) 373 ahtu(:,:,1:jpkm1) = SQRT( ahtu(:,:,1:jpkm1) ) * umask(:,:,1:jpkm1) 374 ahtv(:,:,1:jpkm1) = SQRT( ahtv(:,:,1:jpkm1) ) * vmask(:,:,1:jpkm1) 375 ENDIF 270 376 ENDIF 271 377 ! … … 281 387 !! ** Purpose : update at kt the tracer lateral mixing coeff. (aht and aeiv) 282 388 !! 283 !! ** Method : 389 !! ** Method : * time varying eddy diffusivity coefficients: 284 390 !! 285 391 !! nn_aei_ijk_t = 21 aeiu, aeiv = F(i,j, t) = F(growth rate of baroclinic instability) … … 290 396 !! or |u|e^3/12 bilaplacian operator ) 291 397 !! 398 !! * time varying EIV coefficients: call to ldf_eiv routine 399 !! 292 400 !! ** action : ahtu, ahtv update at each time step 293 401 !! aeiu, aeiv - - - - (if ln_ldfeiv=T) … … 296 404 ! 297 405 INTEGER :: ji, jj, jk ! dummy loop indices 298 REAL(wp) :: zaht, zahf, zaht_min, z 1_f20! local scalar406 REAL(wp) :: zaht, zahf, zaht_min, zDaht, z1_f20 ! local scalar 299 407 !!---------------------------------------------------------------------- 300 408 ! 301 409 IF( ln_ldfeiv .AND. nn_aei_ijk_t == 21 ) THEN ! eddy induced velocity coefficients 302 410 ! ! =F(growth rate of baroclinic instability) 303 ! ! max value rn_aeiv_0 ; decreased to 0 within 20N-20S304 CALL ldf_eiv( kt, rn_aeiv_0, aeiu, aeiv )411 ! ! max value aeiv_0 ; decreased to 0 within 20N-20S 412 CALL ldf_eiv( kt, aei0, aeiu, aeiv ) 305 413 ENDIF 306 414 ! … … 308 416 ! 309 417 CASE( 21 ) !== time varying 2D field ==! = F( growth rate of baroclinic instability ) 310 ! ! min value rn_aht_0 / 10311 ! ! max value rn_aht_0 (rn_aeiv_0 if nn_aei_ijk_t=21)312 ! ! increase to rn_aht_0 within 20N-20S418 ! ! min value 0.2*aht0 419 ! ! max value aht0 (aei0 if nn_aei_ijk_t=21) 420 ! ! increase to aht0 within 20N-20S 313 421 IF( ln_ldfeiv .AND. nn_aei_ijk_t == 21 ) THEN ! use the already computed aei. 314 422 ahtu(:,:,1) = aeiu(:,:,1) 315 423 ahtv(:,:,1) = aeiv(:,:,1) 316 424 ELSE ! compute aht. 317 CALL ldf_eiv( kt, rn_aht_0, ahtu, ahtv )425 CALL ldf_eiv( kt, aht0, ahtu, ahtv ) 318 426 ENDIF 319 427 ! 320 z1_f20 = 1._wp / ( 2._wp * omega * SIN( rad * 20._wp ) ) ! 1 / ff(20 degrees) 321 zaht_min = 0.2_wp * rn_aht_0 ! minimum value for aht 428 z1_f20 = 1._wp / ( 2._wp * omega * SIN( rad * 20._wp ) ) ! 1 / ff(20 degrees) 429 zaht_min = 0.2_wp * aht0 ! minimum value for aht 430 zDaht = aht0 - zaht_min 322 431 DO jj = 1, jpj 323 432 DO ji = 1, jpi 324 433 !!gm CAUTION : here we assume lat/lon grid in 20deg N/S band (like all ORCA cfg) 325 434 !! ==>>> The Coriolis value is identical for t- & u_points, and for v- and f-points 326 zaht = ( 1._wp - MIN( 1._wp , ABS( ff_t(ji,jj) * z1_f20 ) ) ) * ( rn_aht_0 - zaht_min )327 zahf = ( 1._wp - MIN( 1._wp , ABS( ff_f(ji,jj) * z1_f20 ) ) ) * ( rn_aht_0 - zaht_min )328 ahtu(ji,jj,1) = ( MAX( zaht_min, ahtu(ji,jj,1) ) + zaht ) * umask(ji,jj,1)! min value zaht_min329 ahtv(ji,jj,1) = ( MAX( zaht_min, ahtv(ji,jj,1) ) + zahf ) * vmask(ji,jj,1)! increase within 20S-20N330 END DO 331 END DO 332 DO jk = 2, jpkm1 ! deeper value = surface value435 zaht = ( 1._wp - MIN( 1._wp , ABS( ff_t(ji,jj) * z1_f20 ) ) ) * zDaht 436 zahf = ( 1._wp - MIN( 1._wp , ABS( ff_f(ji,jj) * z1_f20 ) ) ) * zDaht 437 ahtu(ji,jj,1) = ( MAX( zaht_min, ahtu(ji,jj,1) ) + zaht ) ! min value zaht_min 438 ahtv(ji,jj,1) = ( MAX( zaht_min, ahtv(ji,jj,1) ) + zahf ) ! increase within 20S-20N 439 END DO 440 END DO 441 DO jk = 1, jpkm1 ! deeper value = surface value + mask for all levels 333 442 ahtu(:,:,jk) = ahtu(:,:,1) * umask(:,:,jk) 334 443 ahtv(:,:,jk) = ahtv(:,:,1) * vmask(:,:,jk) … … 338 447 IF( ln_traldf_lap ) THEN ! laplacian operator |u| e /12 339 448 DO jk = 1, jpkm1 340 ahtu(:,:,jk) = ABS( ub(:,:,jk) ) * e1u(:,:) * r1_12 449 ahtu(:,:,jk) = ABS( ub(:,:,jk) ) * e1u(:,:) * r1_12 ! n.b. ub,vb are masked 341 450 ahtv(:,:,jk) = ABS( vb(:,:,jk) ) * e2v(:,:) * r1_12 342 451 END DO … … 355 464 CALL iom_put( "ahtv_3d", ahtv(:,:,:) ) ! 3D v-eddy diffusivity coeff. 356 465 ! 357 !!gm : THE IF below is to be checked (comes from Seb)358 466 IF( ln_ldfeiv ) THEN 359 467 CALL iom_put( "aeiu_2d", aeiu(:,:,1) ) ! surface u-EIV coeff. … … 372 480 !! ** Purpose : initialization of the eiv coeff. from namelist choices. 373 481 !! 374 !! ** Method : 375 !! 376 !! ** Action : aeiu , aeiv : EIV coeff. at u- & v-points 377 !! l_ldfeiv_time : =T if EIV coefficients vary with time 378 !!---------------------------------------------------------------------- 379 INTEGER :: jk ! dummy loop indices 380 INTEGER :: ierr, inum, ios ! local integer 381 ! 382 NAMELIST/namtra_ldfeiv/ ln_ldfeiv , ln_ldfeiv_dia, & ! eddy induced velocity (eiv) 383 & nn_aei_ijk_t, rn_aeiv_0 ! eiv coefficient 482 !! ** Method : the eiv diffusivity coef. specification depends on: 483 !! nn_aei_ijk_t = 0 => = constant 484 !! ! 485 !! = 10 => = F(z) : constant with a reduction of 1/4 with depth 486 !! ! 487 !! =-20 => = F(i,j) = shape read in 'eddy_diffusivity.nc' file 488 !! = 20 = F(i,j) = F(e1,e2) or F(e1^3,e2^3) (lap or bilap case) 489 !! = 21 = F(i,j,t) = F(growth rate of baroclinic instability) 490 !! ! 491 !! =-30 => = F(i,j,k) = shape read in 'eddy_diffusivity.nc' file 492 !! = 30 = F(i,j,k) = 2D (case 20) + decrease with depth (case 10) 493 !! 494 !! ** Action : aeiu , aeiv : initialized one for all or l_ldftra_time set to true 495 !! l_ldfeiv_time : =T if EIV coefficients vary with time 496 !!---------------------------------------------------------------------- 497 INTEGER :: jk ! dummy loop indices 498 INTEGER :: ierr, inum, ios, inn ! local integer 499 REAL(wp) :: zah_max, zUfac ! - scalar 500 !! 501 NAMELIST/namtra_eiv/ ln_ldfeiv , ln_ldfeiv_dia, & ! eddy induced velocity (eiv) 502 & nn_aei_ijk_t, rn_Ue, rn_Le ! eiv coefficient 384 503 !!---------------------------------------------------------------------- 385 504 ! … … 390 509 ENDIF 391 510 ! 392 REWIND( numnam_ref ) ! Namelist namtra_ ldfeiv in reference namelist : eddy induced velocity param.393 READ ( numnam_ref, namtra_ ldfeiv, IOSTAT = ios, ERR = 901)394 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namtra_ ldfeiv in reference namelist', lwp )395 ! 396 REWIND( numnam_cfg ) ! Namelist namtra_ ldfeiv in configuration namelist : eddy induced velocity param.397 READ ( numnam_cfg, namtra_ ldfeiv, IOSTAT = ios, ERR = 902 )398 902 IF( ios > 0 ) CALL ctl_nam ( ios , 'namtra_ ldfeiv in configuration namelist', lwp )399 IF(lwm) WRITE ( numond, namtra_ ldfeiv )511 REWIND( numnam_ref ) ! Namelist namtra_eiv in reference namelist : eddy induced velocity param. 512 READ ( numnam_ref, namtra_eiv, IOSTAT = ios, ERR = 901) 513 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namtra_eiv in reference namelist', lwp ) 514 ! 515 REWIND( numnam_cfg ) ! Namelist namtra_eiv in configuration namelist : eddy induced velocity param. 516 READ ( numnam_cfg, namtra_eiv, IOSTAT = ios, ERR = 902 ) 517 902 IF( ios > 0 ) CALL ctl_nam ( ios , 'namtra_eiv in configuration namelist', lwp ) 518 IF(lwm) WRITE ( numond, namtra_eiv ) 400 519 401 520 IF(lwp) THEN ! control print 402 WRITE(numout,*) ' Namelist namtra_ldfeiv : ' 403 WRITE(numout,*) ' Eddy Induced Velocity (eiv) param. ln_ldfeiv = ', ln_ldfeiv 404 WRITE(numout,*) ' eiv streamfunction & velocity diag. ln_ldfeiv_dia = ', ln_ldfeiv_dia 405 WRITE(numout,*) ' eddy induced velocity coef. rn_aeiv_0 = ', rn_aeiv_0 406 WRITE(numout,*) ' type of time-space variation nn_aei_ijk_t = ', nn_aei_ijk_t 521 WRITE(numout,*) ' Namelist namtra_eiv : ' 522 WRITE(numout,*) ' Eddy Induced Velocity (eiv) param. ln_ldfeiv = ', ln_ldfeiv 523 WRITE(numout,*) ' eiv streamfunction & velocity diag. ln_ldfeiv_dia = ', ln_ldfeiv_dia 524 WRITE(numout,*) ' coefficients :' 525 WRITE(numout,*) ' type of time-space variation nn_aei_ijk_t = ', nn_aht_ijk_t 526 WRITE(numout,*) ' lateral diffusive velocity (if cst) rn_Ue = ', rn_Ue, ' m/s' 527 WRITE(numout,*) ' lateral diffusive length (if cst) rn_Le = ', rn_Le, ' m' 407 528 WRITE(numout,*) 408 529 ENDIF 409 530 ! 410 IF( ln_ldfeiv .AND. ln_traldf_blp ) CALL ctl_stop( 'ldf_eiv_init: eddy induced velocity ONLY with laplacian diffusivity' ) 411 412 ! ! Parameter control 413 l_ldfeiv_time = .FALSE. 414 ! 415 IF( ln_ldfeiv ) THEN ! allocate the aei arrays 531 l_ldfeiv_time = .FALSE. ! no time variation except in case defined below 532 ! 533 ! 534 IF( .NOT.ln_ldfeiv ) THEN !== Parametrization not used ==! 535 ! 536 IF(lwp) WRITE(numout,*) ' ==>>> eddy induced velocity param is NOT used' 537 ln_ldfeiv_dia = .FALSE. 538 ! 539 ELSE !== use the parametrization ==! 540 ! 541 IF(lwp) WRITE(numout,*) ' ==>>> use eddy induced velocity parametrization' 542 IF(lwp) WRITE(numout,*) 543 ! 544 IF( ln_traldf_blp ) CALL ctl_stop( 'ldf_eiv_init: eddy induced velocity ONLY with laplacian diffusivity' ) 545 ! 546 ! != allocate the aei arrays 416 547 ALLOCATE( aeiu(jpi,jpj,jpk), aeiv(jpi,jpj,jpk), STAT=ierr ) 417 548 IF( ierr /= 0 ) CALL ctl_stop('STOP', 'ldf_eiv: failed to allocate arrays') 418 549 ! 419 SELECT CASE( nn_aei_ijk_t ) ! Specification of space time variations of eaiu, aeiv 420 ! 421 CASE( 0 ) !== constant ==! 422 IF(lwp) WRITE(numout,*) ' ==>>> eddy induced velocity coef. = constant = ', rn_aeiv_0 423 aeiu(:,:,:) = rn_aeiv_0 424 aeiv(:,:,:) = rn_aeiv_0 425 ! 426 CASE( 10 ) !== fixed profile ==! 550 ! != Specification of space-time variations of eaiu, aeiv 551 ! 552 aeiu(:,:,jpk) = 0._wp ! last level always 0 553 aeiv(:,:,jpk) = 0._wp 554 ! ! value of EIV coef. (laplacian operator) 555 zUfac = r1_2 *rn_Ue ! velocity factor 556 inn = 1 ! L-exponent 557 aei0 = zUfac * rn_Le**inn ! mixing coefficient 558 zah_max = zUfac * (ra*rad)**inn ! maximum reachable coefficient (value at the Equator) 559 560 SELECT CASE( nn_aei_ijk_t ) !* Specification of space-time variations 561 ! 562 CASE( 0 ) !-- constant --! 563 IF(lwp) WRITE(numout,*) ' ==>>> eddy induced velocity coef. = constant = ', aei0, ' m2/s' 564 aeiu(:,:,1:jpkm1) = aei0 565 aeiv(:,:,1:jpkm1) = aei0 566 ! 567 CASE( 10 ) !-- fixed profile --! 427 568 IF(lwp) WRITE(numout,*) ' ==>>> eddy induced velocity coef. = F( depth )' 428 aeiu(:,:,1) = rn_aeiv_0 ! constant surface value 429 aeiv(:,:,1) = rn_aeiv_0 430 CALL ldf_c1d( 'TRA', r1_4, aeiu(:,:,1), aeiv(:,:,1), aeiu, aeiv ) 431 ! 432 CASE ( -20 ) !== fixed horizontal shape read in file ==! 433 IF(lwp) WRITE(numout,*) ' ==>>> tracer mixing coef. = F(i,j) read in eddy_diffusivity_2D.nc file' 569 IF(lwp) WRITE(numout,*) ' surface eddy diffusivity = constant = ', aht0, ' m2/s' 570 aeiu(:,:,1) = aei0 ! constant surface value 571 aeiv(:,:,1) = aei0 572 CALL ldf_c1d( 'TRA', aeiu(:,:,1), aeiv(:,:,1), aeiu, aeiv ) 573 ! 574 CASE ( -20 ) !-- fixed horizontal shape read in file --! 575 IF(lwp) WRITE(numout,*) ' ==>>> eddy induced velocity coef. = F(i,j) read in eddy_diffusivity_2D.nc file' 434 576 CALL iom_open ( 'eddy_induced_velocity_2D.nc', inum ) 435 577 CALL iom_get ( inum, jpdom_data, 'aeiu', aeiu(:,:,1) ) 436 578 CALL iom_get ( inum, jpdom_data, 'aeiv', aeiv(:,:,1) ) 437 579 CALL iom_close( inum ) 438 DO jk = 2, jpk 580 DO jk = 2, jpkm1 439 581 aeiu(:,:,jk) = aeiu(:,:,1) 440 582 aeiv(:,:,jk) = aeiv(:,:,1) 441 583 END DO 442 584 ! 443 CASE( 20 ) !== fixed horizontal shape ==! 444 IF(lwp) WRITE(numout,*) ' ==>>> tracer mixing coef. = F( e1, e2 ) or F( e1^3, e2^3 ) (lap or bilap case)' 445 CALL ldf_c2d( 'TRA', 'LAP', rn_aeiv_0, aeiu, aeiv ) ! surface value proportional to scale factor 446 ! 447 CASE( 21 ) !== time varying 2D field ==! 448 IF(lwp) WRITE(numout,*) ' ==>>> tracer mixing coef. = F( latitude, longitude, time )' 449 IF(lwp) WRITE(numout,*) ' = F( growth rate of baroclinic instability )' 585 CASE( 20 ) !-- fixed horizontal shape --! 586 IF(lwp) WRITE(numout,*) ' ==>>> eddy induced velocity coef. = F( e1, e2 )' 587 IF(lwp) WRITE(numout,*) ' using a fixed diffusive velocity = ', rn_Ue, ' m/s and Le = Max(e1,e2)' 588 IF(lwp) WRITE(numout,*) ' maximum reachable coefficient (at the Equator) = ', zah_max, ' m2/s for e1=1°)' 589 CALL ldf_c2d( 'TRA', zUfac , inn , aeiu, aeiv ) ! value proportional to scale factor^inn 590 ! 591 CASE( 21 ) !-- time varying 2D field --! 592 IF(lwp) WRITE(numout,*) ' ==>>> eddy induced velocity coef. = F( latitude, longitude, time )' 593 IF(lwp) WRITE(numout,*) ' = F( growth rate of baroclinic instability )' 594 IF(lwp) WRITE(numout,*) ' maximum allowed value: aei0 = ', aei0, ' m2/s' 450 595 ! 451 596 l_ldfeiv_time = .TRUE. ! will be calculated by call to ldf_tra routine in step.F90 452 597 ! 453 CASE( -30 ) !== fixed 3D shape read in file ==!454 IF(lwp) WRITE(numout,*) ' ==>>> tracer mixingcoef. = F(i,j,k) read in eddy_diffusivity_3D.nc file'598 CASE( -30 ) !-- fixed 3D shape read in file --! 599 IF(lwp) WRITE(numout,*) ' ==>>> eddy induced velocity coef. = F(i,j,k) read in eddy_diffusivity_3D.nc file' 455 600 CALL iom_open ( 'eddy_induced_velocity_3D.nc', inum ) 456 601 CALL iom_get ( inum, jpdom_data, 'aeiu', aeiu ) … … 458 603 CALL iom_close( inum ) 459 604 ! 460 CASE( 30 ) !== fixed 3D shape ==! 461 IF(lwp) WRITE(numout,*) ' ==>>> tracer mixing coef. = F( latitude, longitude, depth )' 462 CALL ldf_c2d( 'TRA', 'LAP', rn_aeiv_0, aeiu, aeiv ) ! surface value proportional to scale factor 463 ! ! reduction with depth 464 CALL ldf_c1d( 'TRA', r1_4, aeiu(:,:,1), aeiv(:,:,1), aeiu, aeiv ) 605 CASE( 30 ) !-- fixed 3D shape --! 606 IF(lwp) WRITE(numout,*) ' ==>>> eddy induced velocity coef. = F( latitude, longitude, depth )' 607 CALL ldf_c2d( 'TRA', zUfac , inn , aeiu, aeiv ) ! surface value proportional to scale factor^inn 608 CALL ldf_c1d( 'TRA', aeiu(:,:,1), aeiv(:,:,1), aeiu, aeiv ) ! reduction with depth 465 609 ! 466 610 CASE DEFAULT … … 468 612 END SELECT 469 613 ! 470 ELSE 471 IF(lwp) WRITE(numout,*) ' ==>>> eddy induced velocity param is NOT used neither diagnosed' 472 ln_ldfeiv_dia = .FALSE. 614 IF( .NOT.l_ldfeiv_time ) THEN !* mask if No time variation 615 DO jk = 1, jpkm1 616 aeiu(:,:,jk) = aeiu(:,:,jk) * umask(:,:,jk) 617 ahtv(:,:,jk) = ahtv(:,:,jk) * vmask(:,:,jk) 618 END DO 619 ENDIF 620 ! 473 621 ENDIF 474 622 ! … … 493 641 INTEGER :: ji, jj, jk ! dummy loop indices 494 642 REAL(wp) :: zfw, ze3w, zn2, z1_f20, zaht, zaht_min, zzaei ! local scalars 495 REAL(wp), DIMENSION(jpi,jpj) :: zn, zah, zhw, z ross, zaeiw ! 2D workspace496 !!---------------------------------------------------------------------- 497 ! 498 zn (:,:) = 0._wp! Local initialization499 zhw 500 zah 501 z ross(:,:) = 0._wp643 REAL(wp), DIMENSION(jpi,jpj) :: zn, zah, zhw, zRo, zaeiw ! 2D workspace 644 !!---------------------------------------------------------------------- 645 ! 646 zn (:,:) = 0._wp ! Local initialization 647 zhw(:,:) = 5._wp 648 zah(:,:) = 0._wp 649 zRo(:,:) = 0._wp 502 650 ! ! Compute lateral diffusive coefficient at T-point 503 651 IF( ln_traldf_triad ) THEN … … 538 686 END DO 539 687 END DO 540 END 688 ENDIF 541 689 542 690 DO jj = 2, jpjm1 543 691 DO ji = fs_2, fs_jpim1 ! vector opt. 544 692 zfw = MAX( ABS( 2. * omega * SIN( rad * gphit(ji,jj) ) ) , 1.e-10 ) 545 ! Rossby radius at w-point taken < 40km and > 2km546 z ross(ji,jj) = MAX( MIN( .4 * zn(ji,jj) / zfw, 40.e3 ), 2.e3)693 ! Rossby radius at w-point taken betwenn 2 km and 40km 694 zRo(ji,jj) = MAX( 2.e3 , MIN( .4 * zn(ji,jj) / zfw, 40.e3 ) ) 547 695 ! Compute aeiw by multiplying Ro^2 and T^-1 548 zaeiw(ji,jj) = z ross(ji,jj) * zross(ji,jj) * SQRT( zah(ji,jj) / zhw(ji,jj) ) * tmask(ji,jj,1)696 zaeiw(ji,jj) = zRo(ji,jj) * zRo(ji,jj) * SQRT( zah(ji,jj) / zhw(ji,jj) ) * tmask(ji,jj,1) 549 697 END DO 550 698 END DO … … 554 702 DO jj = 2, jpjm1 555 703 DO ji = fs_2, fs_jpim1 ! vector opt. 556 zzaei = MIN( 1._wp, ABS( ff_t(ji,jj) * z1_f20 ) ) * zaeiw(ji,jj) 704 zzaei = MIN( 1._wp, ABS( ff_t(ji,jj) * z1_f20 ) ) * zaeiw(ji,jj) ! tropical decrease 557 705 zaeiw(ji,jj) = MIN( zzaei , paei0 ) ! Max value = paei0 558 706 END DO … … 620 768 DO jj = 1, jpjm1 621 769 DO ji = 1, fs_jpim1 ! vector opt. 622 zpsi_uw(ji,jj,jk) = - 0.25_wp* e2u(ji,jj) * ( wslpi(ji,jj,jk ) + wslpi(ji+1,jj,jk) ) &623 & 624 zpsi_vw(ji,jj,jk) = - 0.25_wp* e1v(ji,jj) * ( wslpj(ji,jj,jk ) + wslpj(ji,jj+1,jk) ) &625 & 770 zpsi_uw(ji,jj,jk) = - r1_4 * e2u(ji,jj) * ( wslpi(ji,jj,jk ) + wslpi(ji+1,jj,jk) ) & 771 & * ( aeiu (ji,jj,jk-1) + aeiu (ji ,jj,jk) ) * umask(ji,jj,jk) 772 zpsi_vw(ji,jj,jk) = - r1_4 * e1v(ji,jj) * ( wslpj(ji,jj,jk ) + wslpj(ji,jj+1,jk) ) & 773 & * ( aeiv (ji,jj,jk-1) + aeiv (ji,jj ,jk) ) * vmask(ji,jj,jk) 626 774 END DO 627 775 END DO
Note: See TracChangeset
for help on using the changeset viewer.
