Changeset 2030 for branches/DEV_r2006_merge_TRA_TRC
- Timestamp:
- 2010-07-29T14:19:46+02:00 (14 years ago)
- Location:
- branches/DEV_r2006_merge_TRA_TRC/NEMO/TOP_SRC/TRP
- Files:
-
- 6 added
- 21 deleted
- 8 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/DEV_r2006_merge_TRA_TRC/NEMO/TOP_SRC/TRP/trcbbl.F90
r1606 r2030 1 1 MODULE trcbbl 2 2 !!====================================================================== 3 3 !! *** MODULE trcbbl *** 4 4 !! Ocean passive tracers physics : advective and/or diffusive bottom boundary 5 5 !! layer scheme 6 6 !!====================================================================== 7 !! History : 8.0 ! 96-06 (L. Mortier) Original code 8 !! 8.0 ! 97-11 (G. Madec) Optimization 9 !! 8.5 ! 02-08 (G. Madec) free form + modules 10 !! 9.0 ! 04-03 (C. Ethe) Adaptation for passive tracers 11 !! ! 07-02 (C. Deltel) Diagnose ML trends for passive tracers 7 !!============================================================================== 8 !! History : OPA ! 1996-06 (L. Mortier) Original code 9 !! 8.0 ! 1997-11 (G. Madec) Optimization 10 !! NEMO 1.0 ! 2002-08 (G. Madec) free form + modules 11 !! - ! 2004-01 (A. de Miranda, G. Madec, J.M. Molines ) add advective bbl 12 !! 3.3 ! 2009-11 (G. Madec) merge trabbl and trabbl_adv + style + optimization 13 !! - ! 2010-04 (G. Madec) Campin & Goosse advective bbl 14 !! - ! 2010-06 (C. Ethe, G. Madec) merge TRA-TRC 12 15 !!---------------------------------------------------------------------- 13 #if defined key_top && ( defined key_trcbbl_dif || defined key_trcbbl_adv ) && ! defined key_c1d16 #if defined key_top && defined key_trabbl 14 17 !!---------------------------------------------------------------------- 15 !! 'key_trcbbl_dif' or diffusive bottom boundary layer 16 !! 'key_trcbbl_adv' advective bottom boundary layer 18 !! 'key_trabbl diffusive or/and adevective bottom boundary layer 17 19 !!---------------------------------------------------------------------- 18 !! trc_bbl_dif : update the passive tracer trends due to the bottom 19 !! boundary layer (diffusive only) 20 !! trc_bbl_adv : update the passive tracer trends due to the bottom 21 !! boundary layer (advective and/or diffusive) 20 !! trc_bbl : update the tracer trends due to the bottom boundary layer (advective and/or diffusive) 22 21 !!---------------------------------------------------------------------- 23 22 USE oce_trc ! ocean dynamics and active tracers variables 24 23 USE trc ! ocean passive tracers variables 25 USE trctrp_lec ! passive tracers transport 24 USE trcnam_trp ! passive tracers transport namelist variables 25 USE trabbl ! 26 26 USE prtctl_trc ! Print control for debbuging 27 USE eosbn2 28 USE lbclnk 29 USE trdmld_trc 30 USE trdmld_trc_oce 27 USE trdmod_oce 28 USE trdtra 31 29 32 IMPLICIT NONE 33 PRIVATE 30 PUBLIC trc_bbl ! routine called by step.F90 34 31 35 PUBLIC trc_bbl_dif ! routine called by step.F9036 PUBLIC trc_bbl_adv ! routine called by step.F9037 38 # if defined key_trcbbl_dif39 LOGICAL, PUBLIC, PARAMETER :: lk_trcbbl_dif = .TRUE. !: diffusive bottom boundary layer flag40 # else41 LOGICAL, PUBLIC, PARAMETER :: lk_trcbbl_dif = .FALSE. !: diffusive bottom boundary layer flag42 # endif43 44 # if defined key_trcbbl_adv45 LOGICAL, PUBLIC, PARAMETER :: lk_trcbbl_adv = .TRUE. !: advective bottom boundary layer flag46 REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) :: u_trc_bbl !: veloc. involved in the advective BBL47 REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) :: v_trc_bbl !: veloc. involved in the advective BBL48 REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) :: w_trc_bbl !: vertic. increment of veloc. due to adv. BBL49 ! ! only affect tracer vertical advection50 # else51 LOGICAL, PUBLIC, PARAMETER :: lk_trcbbl_adv = .FALSE. !: advective bottom boundary layer flag52 # endif53 54 INTEGER, DIMENSION(jpi,jpj) :: mbkt, mbku, mbkv55 32 56 33 !! * Substitutions 57 34 # include "top_substitute.h90" 58 35 !!---------------------------------------------------------------------- 59 !! TOP 1.0 , LOCEAN-IPSL (2005)60 !! $ Header: /home/opalod/NEMOCVSROOT/NEMO/TOP_SRC/TRP/trcbbl.F90,v 1.12 2006/09/12 11:10:13 opalod Exp$36 !! NEMO/OPA 3.2 , LOCEAN-IPSL (2009) 37 !! $Id$ 61 38 !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 62 39 !!---------------------------------------------------------------------- … … 64 41 CONTAINS 65 42 66 SUBROUTINE trc_bbl_dif( kt ) 43 44 SUBROUTINE trc_bbl( kt ) 67 45 !!---------------------------------------------------------------------- 68 !! *** ROUTINE trc_bbl_dif *** 46 !! *** ROUTINE bbl *** 47 !! 48 !! ** Purpose : Compute the before tracer (t & s) trend associated 49 !! with the bottom boundary layer and add it to the general trend 50 !! of tracer equations. 69 51 !! 70 !! ** Purpose : Compute the before tracer trend associated 71 !! with the bottom boundary layer and add it to the general trend 72 !! of tracer equations. The bottom boundary layer is supposed to be 73 !! a purely diffusive bottom boundary layer. 74 !! 75 !! ** Method : When the product grad( rho) * grad(h) < 0 (where grad 76 !! is an along bottom slope gradient) an additional lateral diffu- 77 !! sive trend along the bottom slope is added to the general tracer 78 !! trend, otherwise nothing is done. 79 !! Second order operator (laplacian type) with variable coefficient 80 !! computed as follow for temperature (idem on s): 81 !! difft = 1/(e1t*e2t*e3t) { di-1[ ahbt e2u*e3u/e1u di[ztb] ] 82 !! + dj-1[ ahbt e1v*e3v/e2v dj[ztb] ] } 83 !! where ztb is a 2D array: the bottom ocean temperature and ahtb 84 !! is a time and space varying diffusive coefficient defined by: 85 !! ahbt = zahbp if grad(rho).grad(h) < 0 86 !! = 0. otherwise. 87 !! Note that grad(.) is the along bottom slope gradient. grad(rho) 88 !! is evaluated using the local density (i.e. referenced at the 89 !! local depth). Typical value of ahbt is 2000 m2/s (equivalent to 90 !! a downslope velocity of 20 cm/s if the condition for slope 91 !! convection is satified) 92 !! Add this before trend to the general trend tra of the 93 !! botton ocean tracer point: 94 !! tra = tra + difft 95 !! 96 !! ** Action : - update tra at the bottom level with the bottom 97 !! boundary layer trend 98 !! 99 !! References : 100 !! Beckmann, A., and R. Doscher, 1997, J. Phys.Oceanogr., 581-591. 101 !!---------------------------------------------------------------------- 102 USE oce, ONLY : ztrtrd => ua ! use ua as 3D workspace 103 !! 104 INTEGER, INTENT( in ) :: kt ! ocean time-step 105 INTEGER :: ji, jj, jn ! dummy loop indices 106 INTEGER :: ik, iku, ikv 107 INTEGER :: ii0, ii1, ij0, ij1 ! temporary integers 108 INTEGER :: iku1, iku2, ikv1, ikv2 ! temporary intergers 109 REAL(wp) :: ze3u, ze3v ! temporary scalars 110 REAL(wp) :: zbtr, ztra 111 #if ! defined key_off_tra 112 REAL(wp) :: zgdrho, zalbet, zsign, zt, zs, zh 113 REAL(wp), DIMENSION(jpi,jpj) :: zki, zkj 114 #endif 115 REAL(wp), DIMENSION(jpi,jpj) :: zkx, zky ! temporary workspace arrays 116 REAL(wp), DIMENSION(jpi,jpj) :: ztnb, zsnb, zdep 117 REAL(wp), DIMENSION(jpi,jpj) :: ztrb, zahu, zahv 118 52 !!---------------------------------------------------------------------- 53 INTEGER, INTENT( in ) :: kt ! ocean time-step 119 54 CHARACTER (len=22) :: charout 120 REAL(wp) :: fsalbt, pft, pfs, pfh ! statement function 121 !!---------------------------------------------------------------------- 122 ! ratio alpha/beta 123 ! ================ 124 ! fsalbt: ratio of thermal over saline expension coefficients 125 ! pft : potential temperature in degrees celcius 126 ! pfs : salinity anomaly (s-35) in psu 127 ! pfh : depth in meters 128 129 fsalbt( pft, pfs, pfh ) = & 130 ( ( ( -0.255019e-07 * pft + 0.298357e-05 ) * pft & 131 - 0.203814e-03 ) * pft & 132 + 0.170907e-01 ) * pft & 133 + 0.665157e-01 & 134 +(-0.678662e-05 * pfs - 0.846960e-04 * pft + 0.378110e-02 ) * pfs & 135 + ( ( - 0.302285e-13 * pfh & 136 - 0.251520e-11 * pfs & 137 + 0.512857e-12 * pft * pft ) * pfh & 138 - 0.164759e-06 * pfs & 139 +( 0.791325e-08 * pft - 0.933746e-06 ) * pft & 140 + 0.380374e-04 ) * pfh 55 REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: ztrtrd 141 56 !!---------------------------------------------------------------------- 142 57 58 #if ! defined key_offline 59 ! Online coupling with dynamics : Computation of bbl coef and bbl transport 60 ! Offline coupling with dynamics : Read bbl coef and bbl transport from input files 61 CALL bbl( kt, 'TRC' ) 62 l_bbl = .FALSE. 63 #endif 143 64 144 IF( kt == nittrc000 ) CALL trc_bbl_init 145 146 147 ! 0. 2D fields of bottom temperature and salinity, and bottom slope 148 ! ----------------------------------------------------------------- 149 ! mbathy= number of w-level, minimum value=1 (cf dommsk.F) 150 151 # if defined key_vectopt_loop 152 jj = 1 153 DO ji = 1, jpij ! vector opt. (forced unrolling) 154 # else 155 DO jj = 1, jpj 156 DO ji = 1, jpi 157 # endif 158 ik = mbkt(ji,jj) ! index of the bottom ocean T-level 159 ztnb(ji,jj) = tn(ji,jj,ik) * tmask(ji,jj,1) ! masked now T and S at ocean bottom 160 zsnb(ji,jj) = sn(ji,jj,ik) * tmask(ji,jj,1) 161 zdep(ji,jj) = fsdept(ji,jj,ik) ! depth of the ocean bottom T-level 162 # if ! defined key_vectopt_loop 163 END DO 164 # endif 165 END DO 166 167 IF( ln_zps ) THEN ! partial steps correction 168 169 # if defined key_vectopt_loop 170 jj = 1 171 DO ji = 1, jpij-jpi ! vector opt. (forced unrolling) 172 # else 173 DO jj = 1, jpjm1 174 DO ji = 1, jpim1 175 # endif 176 iku1 = MAX( mbathy(ji+1,jj )-1, 1 ) 177 iku2 = MAX( mbathy(ji ,jj )-1, 1 ) 178 ikv1 = MAX( mbathy(ji ,jj+1)-1, 1 ) 179 ikv2 = MAX( mbathy(ji ,jj )-1, 1 ) 180 ze3u = MIN( fse3u(ji,jj,iku1), fse3u(ji,jj,iku2) ) 181 ze3v = MIN( fse3v(ji,jj,ikv1), fse3v(ji,jj,ikv2) ) 182 zahu(ji,jj) = atrcbbl * e2u(ji,jj) * ze3u / e1u(ji,jj) * umask(ji,jj,1) 183 zahv(ji,jj) = atrcbbl * e1v(ji,jj) * ze3v / e2v(ji,jj) * vmask(ji,jj,1) 184 # if ! defined key_vectopt_loop 185 END DO 186 # endif 187 END DO 188 ELSE ! z-coordinate - full steps or s-coordinate 189 # if defined key_vectopt_loop 190 jj = 1 191 DO ji = 1, jpij-jpi ! vector opt. (forced unrolling) 192 # else 193 DO jj = 1, jpjm1 194 DO ji = 1, jpim1 195 # endif 196 iku = mbku(ji,jj) 197 ikv = mbkv(ji,jj) 198 zahu(ji,jj) = atrcbbl * e2u(ji,jj) * fse3u(ji,jj,iku) / e1u(ji,jj) * umask(ji,jj,1) 199 zahv(ji,jj) = atrcbbl * e1v(ji,jj) * fse3v(ji,jj,ikv) / e2v(ji,jj) * vmask(ji,jj,1) 200 # if ! defined key_vectopt_loop 201 END DO 202 # endif 203 END DO 65 IF( l_trdtrc ) THEN 66 ALLOCATE( ztrtrd(jpi,jpj,jpk,jptra) ) ! temporary save of trends 67 ztrtrd(:,:,:,:) = tra(:,:,:,:) 204 68 ENDIF 205 69 206 #if defined key_off_tra 207 !!===================================================================== 208 !! I. OFFLINE VERSION OF DIFFUSIVE BBL 209 !!===================================================================== 210 211 ! 1. Criteria of additional bottom diffusivity : grad(rho).grad(h) < 0 212 ! -------------------------------------------------------------------- 213 214 ! Only used in the online version of diffusive BBL (see below) 215 216 ! 2. Additional second order diffusive trends 217 ! ------------------------------------------- 218 ! ! =========== 219 DO jn = 1, jptra ! tracer loop 220 ! ! =========== 70 !* Diffusive bbl : 71 IF( nn_bbl_ldf == 1 ) THEN 72 ! 73 CALL tra_bbl_dif( trb, tra, jptra ) 74 IF( ln_ctl ) THEN 75 WRITE(charout, FMT="(' bbl_dif')") ; CALL prt_ctl_trc_info(charout) 76 CALL prt_ctl_trc( tab4d=tra, mask=tmask, clinfo=ctrcnm, clinfo2='trd' ) 77 ENDIF 78 ! 79 END IF 221 80 222 IF( l_trdtrc ) ztrtrd(:,:,:) = tra(:,:,:,jn) 223 224 ! first derivative (gradient) 225 # if defined key_vectopt_loop 226 jj = 1 227 DO ji = 1, jpij ! vector opt. (forced unrolling) 228 # else 229 DO jj = 1, jpj 230 DO ji = 1, jpi 231 # endif 232 ik = mbkt(ji,jj) 233 ztrb(ji,jj) = trb(ji,jj,ik,jn) * tmask(ji,jj,1) 234 # if ! defined key_vectopt_loop 235 END DO 236 # endif 237 END DO 81 !* Advective bbl : bbl upstream advective trends added to the tracer trends 82 IF( nn_bbl_adv /= 0 ) THEN 83 ! 84 CALL tra_bbl_adv( trb, tra, jptra ) 85 IF( ln_ctl ) THEN 86 WRITE(charout, FMT="(' bbl_adv')") ; CALL prt_ctl_trc_info(charout) 87 CALL prt_ctl_trc( tab4d=tra, mask=tmask, clinfo=ctrcnm, clinfo2='trd' ) 88 ENDIF 89 ! 90 END IF 238 91 239 # if defined key_vectopt_loop 240 jj = 1 241 DO ji = 1, jpij-jpi ! vector opt. (forced unrolling) 242 # else 243 DO jj = 1, jpjm1 244 DO ji = 1, jpim1 245 # endif 246 zkx(ji,jj) = bblx(ji,jj) * zahu(ji,jj) * ( ztrb(ji+1,jj) - ztrb(ji,jj) ) 247 zky(ji,jj) = bbly(ji,jj) * zahv(ji,jj) * ( ztrb(ji,jj+1) - ztrb(ji,jj) ) 248 # if ! defined key_vectopt_loop 249 END DO 250 # endif 251 END DO 252 253 #else 254 !!===================================================================== 255 !! II. ONLINE VERSION OF DIFFUSIVE BBL 256 !!===================================================================== 257 258 ! 1. Criteria of additional bottom diffusivity : grad(rho).grad(h) < 0 259 ! -------------------------------------------------------------------- 260 ! Sign of the local density gradient along the i- and j-slopes 261 ! multiplied by the slope of the ocean bottom 262 SELECT CASE ( nn_eos ) 263 264 CASE ( 0 ) ! Jackett and McDougall (1994) formulation 265 266 # if defined key_vectopt_loop 267 jj = 1 268 DO ji = 1, jpij-jpi ! vector opt. (forced unrolling) 269 # else 270 DO jj = 1, jpjm1 271 DO ji = 1, jpim1 272 # endif 273 ! temperature, salinity anomalie and depth 274 zt = 0.5 * ( ztnb(ji,jj) + ztnb(ji+1,jj) ) 275 zs = 0.5 * ( zsnb(ji,jj) + zsnb(ji+1,jj) ) - 35.0 276 zh = 0.5 * ( zdep(ji,jj) + zdep(ji+1,jj) ) 277 ! masked ratio alpha/beta 278 zalbet = fsalbt( zt, zs, zh )*umask(ji,jj,1) 279 ! local density gradient along i-bathymetric slope 280 zgdrho = zalbet * ( ztnb(ji+1,jj) - ztnb(ji,jj) ) & 281 - ( zsnb(ji+1,jj) - zsnb(ji,jj) ) 282 ! sign of local i-gradient of density multiplied by the i-slope 283 zsign = SIGN( 0.5, - zgdrho * ( zdep(ji+1,jj) - zdep(ji,jj) ) ) 284 zki(ji,jj) = ( 0.5 - zsign ) * zahu(ji,jj) 285 # if ! defined key_vectopt_loop 286 END DO 287 # endif 288 END DO 289 290 # if defined key_vectopt_loop 291 jj = 1 292 DO ji = 1, jpij-jpi ! vector opt. (forced unrolling) 293 # else 294 DO jj = 1, jpjm1 295 DO ji = 1, jpim1 296 # endif 297 ! temperature, salinity anomalie and depth 298 zt = 0.5 * ( ztnb(ji,jj+1) + ztnb(ji,jj) ) 299 zs = 0.5 * ( zsnb(ji,jj+1) + zsnb(ji,jj) ) - 35.0 300 zh = 0.5 * ( zdep(ji,jj+1) + zdep(ji,jj) ) 301 ! masked ratio alpha/beta 302 zalbet = fsalbt( zt, zs, zh )*vmask(ji,jj,1) 303 ! local density gradient along j-bathymetric slope 304 zgdrho = zalbet * ( ztnb(ji,jj+1) - ztnb(ji,jj) ) & 305 - ( zsnb(ji,jj+1) - zsnb(ji,jj) ) 306 ! sign of local j-gradient of density multiplied by the j-slope 307 zsign = SIGN( 0.5, -zgdrho * ( zdep(ji,jj+1) - zdep(ji,jj) ) ) 308 zkj(ji,jj) = ( 0.5 - zsign ) * zahv(ji,jj) 309 # if ! defined key_vectopt_loop 310 END DO 311 # endif 312 END DO 313 314 CASE ( 1 ) ! Linear formulation function of temperature only 315 316 # if defined key_vectopt_loop 317 jj = 1 318 DO ji = 1, jpij-jpi ! vector opt. (forced unrolling) 319 # else 320 DO jj = 1, jpjm1 321 DO ji = 1, jpim1 322 # endif 323 ! local density gradient along i-bathymetric slope 324 zgdrho = ( ztnb(ji+1,jj) - ztnb(ji,jj) ) 325 ! sign of local i-gradient of density multiplied by the i-slope 326 zsign = SIGN( 0.5, - zgdrho * ( zdep(ji+1,jj) - zdep(ji,jj) ) ) 327 zki(ji,jj) = ( 0.5 - zsign ) * zahu(ji,jj) 328 # if ! defined key_vectopt_loop 329 END DO 330 # endif 331 END DO 332 333 # if defined key_vectopt_loop 334 jj = 1 335 DO ji = 1, jpij-jpi ! vector opt. (forced unrolling) 336 # else 337 DO jj = 1, jpjm1 338 DO ji = 1, jpim1 339 # endif 340 ! local density gradient along j-bathymetric slope 341 zgdrho = ( ztnb(ji,jj+1) - ztnb(ji,jj) ) 342 ! sign of local j-gradient of density multiplied by the j-slope 343 zsign = SIGN( 0.5, -zgdrho * ( zdep(ji,jj+1) - zdep(ji,jj) ) ) 344 zkj(ji,jj) = ( 0.5 - zsign ) * zahv(ji,jj) 345 346 # if ! defined key_vectopt_loop 347 END DO 348 # endif 349 END DO 350 351 CASE ( 2 ) ! Linear formulation function of temperature and salinity 352 353 DO jj = 1, jpjm1 354 DO ji = 1, fs_jpim1 ! vector opt. 355 ! local density gradient along i-bathymetric slope 356 zgdrho = - ( rn_beta*( zsnb(ji+1,jj) - zsnb(ji,jj) ) & 357 - rn_alpha*( ztnb(ji+1,jj) - ztnb(ji,jj) ) ) 358 ! sign of local i-gradient of density multiplied by the i-slope 359 zsign = SIGN( 0.5, - zgdrho * ( zdep(ji+1,jj) - zdep(ji,jj) ) ) 360 zki(ji,jj) = ( 0.5 - zsign ) * zahu(ji,jj) 92 IF( l_trdtrc ) THEN ! save the horizontal diffusive trends for further diagnostics 93 DO jn = 1, jptra 94 ztrtrd(:,:,:,jn) = tra(:,:,:,jn) - ztrtrd(:,:,:,jn) 95 CALL trd_tra( kt, 'TRC', jn, jptra_trd_ldf, ztrtrd(:,:,:,jn) ) 361 96 END DO 362 END DO 363 364 DO jj = 1, jpjm1 365 DO ji = 1, fs_jpim1 ! vector opt. 366 ! local density gradient along j-bathymetric slope 367 zgdrho = - ( rn_beta*( zsnb(ji,jj+1) - zsnb(ji,jj) ) & 368 - rn_alpha*( ztnb(ji,jj+1) - ztnb(ji,jj) ) ) 369 ! sign of local j-gradient of density multiplied by the j-slope 370 zsign = sign( 0.5, -zgdrho * ( zdep(ji,jj+1) - zdep(ji,jj) ) ) 371 zkj(ji,jj) = ( 0.5 - zsign ) * zahv(ji,jj) 372 END DO 373 END DO 374 375 CASE DEFAULT 376 377 WRITE(ctmp1,*) ' bad flag value for nn_eos = ', nn_eos 378 CALL ctl_stop( ctmp1 ) 379 380 END SELECT 381 382 ! 2. Additional second order diffusive trends 383 ! ------------------------------------------- 384 ! ! =========== 385 DO jn = 1, jptra ! tracer loop 386 ! ! =========== 387 IF( l_trdtrc ) ztrtrd(:,:,:) = tra(:,:,:,jn) 388 389 ! first derivative (gradient) 390 # if defined key_vectopt_loop 391 jj = 1 392 DO ji = 1, jpij ! vector opt. (forced unrolling) 393 # else 394 DO jj = 1, jpj 395 DO ji = 1, jpi 396 # endif 397 ik = mbkt(ji,jj) 398 ztrb(ji,jj) = trb(ji,jj,ik,jn) * tmask(ji,jj,1) 399 # if ! defined key_vectopt_loop 400 END DO 401 # endif 402 END DO 403 # if defined key_vectopt_loop 404 jj = 1 405 DO ji = 1, jpij-jpi ! vector opt. (forced unrolling) 406 # else 407 DO jj = 1, jpjm1 408 DO ji = 1, jpim1 409 # endif 410 zkx(ji,jj) = zki(ji,jj) * ( ztrb(ji+1,jj) - ztrb(ji,jj) ) 411 zky(ji,jj) = zkj(ji,jj) * ( ztrb(ji,jj+1) - ztrb(ji,jj) ) 412 # if ! defined key_vectopt_loop 413 END DO 414 # endif 415 END DO 416 #endif 417 418 !!===================================================================== 419 !! III. COMMON CODE FOR OFFLINE/ONLINE VERSIONS OF DIFFUSIVE BBL 420 !!===================================================================== 421 422 IF( cp_cfg == "orca" ) THEN 423 424 SELECT CASE ( jp_cfg ) 425 ! ! ======================= 426 CASE ( 2 ) ! ORCA_R2 configuration 427 ! ! ======================= 428 ! Gibraltar enhancement of BBL 429 ij0 = 102 ; ij1 = 102 430 ii0 = 139 ; ii1 = 140 431 zkx( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 4.e0 * zkx( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) 432 zky( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 4.e0 * zky( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) 433 434 ! Red Sea enhancement of BBL 435 ij0 = 88 ; ij1 = 88 436 ii0 = 161 ; ii1 = 162 437 zkx( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 10.e0 * zkx( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) 438 zky( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 10.e0 * zky( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) 439 440 ! ! ======================= 441 CASE ( 4 ) ! ORCA_R4 configuration 442 ! ! ======================= 443 ! Gibraltar enhancement of BBL 444 ij0 = 52 ; ij1 = 52 445 ii0 = 70 ; ii1 = 71 446 zkx( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 4.e0 * zkx( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) 447 zky( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 4.e0 * zky( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) 448 449 END SELECT 450 451 ENDIF 452 453 ! second derivative (divergence) and add to the general tracer trend 454 # if defined key_vectopt_loop 455 jj = 1 456 DO ji = jpi+2, jpij-jpi-1 ! vector opt. (forced unrolling) 457 # else 458 DO jj = 2, jpjm1 459 DO ji = 2, jpim1 460 # endif 461 ik = MAX( mbathy(ji,jj)-1, 1 ) 462 zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,ik) ) 463 ztra = ( zkx(ji,jj) - zkx(ji-1,jj ) & 464 & + zky(ji,jj) - zky(ji ,jj-1) ) * zbtr 465 tra(ji,jj,ik,jn) = tra(ji,jj,ik,jn) + ztra 466 # if ! defined key_vectopt_loop 467 END DO 468 # endif 469 END DO 470 471 ! save the trends for diagnostic 472 IF( l_trdtrc ) THEN 473 ztrtrd(:,:,:) = tra(:,:,:,jn) - ztrtrd(:,:,:) 474 IF (luttrd(jn)) CALL trd_mod_trc( ztrtrd, jn, jptrc_trd_bbl, kt ) 475 END IF 476 ! ! =========== 477 END DO ! tracer loop 478 ! ! =========== 479 480 IF( ln_ctl ) THEN ! print mean trends (used for debugging) 481 WRITE(charout, FMT="('bbl - dif')") 482 CALL prt_ctl_trc_info(charout) 483 CALL prt_ctl_trc(tab4d=tra, mask=tmask, clinfo=ctrcnm,clinfo2='trd') 97 DEALLOCATE( ztrtrd ) 484 98 ENDIF 485 486 END SUBROUTINE trc_bbl_dif 487 488 # if defined key_trcbbl_adv 489 !!---------------------------------------------------------------------- 490 !! 'key_trcbbl_adv' advective bottom boundary layer 491 !!---------------------------------------------------------------------- 492 # include "trcbbl_adv.h90" 493 # else 494 !!---------------------------------------------------------------------- 495 !! Default option : NO advective bottom boundary layer 496 !!---------------------------------------------------------------------- 497 SUBROUTINE trc_bbl_adv (kt ) ! Empty routine 498 INTEGER, INTENT(in) :: kt 499 WRITE(*,*) 'trc_bbl_adv: You should not have seen this print! error?', kt 500 END SUBROUTINE trc_bbl_adv 501 # endif 502 503 SUBROUTINE trc_bbl_init 504 !!---------------------------------------------------------------------- 505 !! *** ROUTINE trc_bbl_init *** 506 !! 507 !! ** Purpose : Initialization for the bottom boundary layer scheme. 508 !!---------------------------------------------------------------------- 509 INTEGER :: ji, jj 510 REAL(wp), DIMENSION(jpi,jpj) :: zmbk 511 !!---------------------------------------------------------------------- 512 513 DO jj = 1, jpj 514 DO ji = 1, jpi 515 mbkt(ji,jj) = MAX( mbathy(ji,jj) - 1, 1 ) ! vertical index of the bottom ocean T-level 516 END DO 517 END DO 518 519 DO jj = 1, jpjm1 520 DO ji = 1, jpim1 521 mbku(ji,jj) = MAX( MIN( mbathy(ji+1,jj ), mbathy(ji,jj) ) - 1, 1 ) 522 mbkv(ji,jj) = MAX( MIN( mbathy(ji ,jj+1), mbathy(ji,jj) ) - 1, 1 ) 523 END DO 524 END DO 525 526 zmbk(:,:) = FLOAT( mbku (:,:) ) 527 CALL lbc_lnk(zmbk,'U',1.) 528 mbku (:,:) = MAX( INT( zmbk(:,:) ), 1 ) 529 530 zmbk(:,:) = FLOAT( mbkv (:,:) ) 531 CALL lbc_lnk(zmbk,'V',1.) 532 mbkv (:,:) = MAX( INT( zmbk(:,:) ), 1 ) 533 534 # if defined key_trcbbl_adv 535 w_trc_bbl(:,:,:) = 0.e0 ! initialisation of w_trc_bbl to zero 536 # endif 537 538 END SUBROUTINE trc_bbl_init 99 ! 100 END SUBROUTINE trc_bbl 539 101 540 102 #else … … 542 104 !! Dummy module : No bottom boundary layer scheme 543 105 !!---------------------------------------------------------------------- 544 LOGICAL, PUBLIC, PARAMETER :: lk_trcbbl_dif = .FALSE. !: diff bbl flag545 LOGICAL, PUBLIC, PARAMETER :: lk_trcbbl_adv = .FALSE. !: adv bbl flag546 106 CONTAINS 547 SUBROUTINE trc_bbl_dif (kt ) ! Empty routine 548 INTEGER, INTENT(in) :: kt 549 WRITE(*,*) 'trc_bbl_dif: You should not have seen this print! error?', kt 550 END SUBROUTINE trc_bbl_dif 551 SUBROUTINE trc_bbl_adv (kt ) ! Empty routine 552 INTEGER, INTENT(in) :: kt 553 WRITE(*,*) 'trc_bbl_adv: You should not have seen this print! error?', kt 554 END SUBROUTINE trc_bbl_adv 107 SUBROUTINE trc_bbl( kt ) ! Empty routine 108 WRITE(*,*) 'tra_bbl: You should not have seen this print! error?', kt 109 END SUBROUTINE trc_bbl 555 110 #endif 556 111 -
branches/DEV_r2006_merge_TRA_TRC/NEMO/TOP_SRC/TRP/trcdmp.F90
r1175 r2030 4 4 !! Ocean physics: internal restoring trend on passive tracers 5 5 !!====================================================================== 6 !! History : 7.0 ! (G. Madec) Original code7 !! ! 96-01 (G. Madec)8 !! ! 97-05 (H. Loukos) adapted for passive tracers9 !! 8.5 ! 02-08 (G. Madec )free form + modules10 !! 9.0 ! 04-03 (C. Ethe) free form + modules11 !! ! 07-02 (C. Deltel) Diagnose ML trends for passive tracers6 !! History : OPA ! 1991-03 (O. Marti, G. Madec) Original code 7 !! ! 1996-01 (G. Madec) statement function for e3 8 !! ! 1997-05 (H. Loukos) adapted for passive tracers 9 !! NEMO 9.0 ! 2004-03 (C. Ethe) free form + modules 10 !! 3.2 ! 2007-02 (C. Deltel) Diagnose ML trends for passive tracers 11 !! 3.3 ! 2010-06 (C. Ethe, G. Madec) merge TRA-TRC 12 12 !!---------------------------------------------------------------------- 13 13 #if defined key_top && defined key_trcdmp … … 17 17 !! trc_dmp : update the tracer trend with the internal damping 18 18 !! trc_dmp_init : initialization, namlist read, parameters control 19 !! trccof_zoom : restoring coefficient for zoom domain20 !! trccof : restoring coefficient for global domain21 !! cofdis : compute the distance to the coastline22 19 !!---------------------------------------------------------------------- 23 20 USE oce_trc ! ocean dynamics and tracers variables 24 21 USE trc ! ocean passive tracers variables 25 USE trc trp_lec ! passive tracers transport22 USE trcnam_trp ! passive tracers transport namelist variables 26 23 USE trcdta 24 USE tradmp 27 25 USE prtctl_trc ! Print control for debbuging 28 USE trdmld_trc 29 USE trdmld_trc_oce 26 USE trdtra 30 27 31 28 IMPLICIT NONE … … 35 32 36 33 LOGICAL , PUBLIC, PARAMETER :: lk_trcdmp = .TRUE. !: internal damping flag 37 REAL(wp), DIMENSION(jpi,jpj,jpk,jptra) :: restotr ! restoring coeff. on tracers (s-1) 34 ! !!* Namelist namtrc_dmp : passive tracer newtonian damping * 35 INTEGER :: nn_hdmp_tr = -1 ! = 0/-1/'latitude' for damping over passive tracer 36 INTEGER :: nn_zdmp_tr = 0 ! = 0/1/2 flag for damping in the mixed layer 37 REAL(wp) :: rn_surf_tr = 50. ! surface time scale for internal damping [days] 38 REAL(wp) :: rn_bot_tr = 360. ! bottom time scale for internal damping [days] 39 REAL(wp) :: rn_dep_tr = 800. ! depth of transition between rn_surf and rn_bot [meters] 40 INTEGER :: nn_file_tr = 2 ! = 1 create a damping.coeff NetCDF file 41 42 REAL(wp), DIMENSION(jpi,jpj,jpk) :: restotr ! restoring coeff. on tracers (s-1) 38 43 39 44 !! * Substitutions … … 66 71 !! - save the trends ('key_trdmld_trc') 67 72 !!---------------------------------------------------------------------- 68 USE oce, ONLY : ztrtrd => ua ! use ua as 3D workspace69 73 !! 70 74 INTEGER, INTENT( in ) :: kt ! ocean time-step index 75 !! 71 76 INTEGER :: ji, jj, jk, jn ! dummy loop indices 72 REAL(wp) :: zt est, ztra !!!, zdt! temporary scalars77 REAL(wp) :: ztra ! temporary scalars 73 78 CHARACTER (len=22) :: charout 79 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ztrtrd 74 80 !!---------------------------------------------------------------------- 75 81 … … 78 84 IF( kt == nittrc000 ) CALL trc_dmp_init 79 85 86 IF( l_trdtrc ) ALLOCATE( ztrtrd(jpi,jpj,jpk) ) ! temporary save of trends 80 87 81 88 ! 1. Newtonian damping trends on tracer fields 82 89 ! -------------------------------------------- 83 ! compute the newtonian damping trends depending on nmldmptr84 85 !!! zdt = rdt * FLOAT( ndttrc )86 87 90 ! Initialize the input fields for newtonian damping 88 CALL dta_trc( kt ) 89 91 CALL trc_dta( kt ) 90 92 ! ! =========== 91 93 DO jn = 1, jptra ! tracer loop … … 94 96 95 97 IF( lutini(jn) ) THEN 96 97 SELECT CASE ( nmldmptr ) 98 99 CASE( 0 ) ! newtonian damping throughout the water column 100 98 ! 99 SELECT CASE ( nn_zdmp_trc ) 100 ! 101 CASE( 0 ) !== newtonian damping throughout the water column ==! 101 102 DO jk = 1, jpkm1 102 103 DO jj = 2, jpjm1 103 104 DO ji = fs_2, fs_jpim1 ! vector opt. 104 ztra = restotr(ji,jj,jk,jn) * ( trdta(ji,jj,jk,jn) - trb(ji,jj,jk,jn) ) 105 ! add the trends to the general tracer trends 106 !! trn(ji,jj,jk,jn) = trn(ji,jj,jk,jn) + ztra * zdt 105 ztra = restotr(ji,jj,jk) * ( trdta(ji,jj,jk,jn) - trb(ji,jj,jk,jn) ) 107 106 tra(ji,jj,jk,jn) = tra(ji,jj,jk,jn) + ztra 108 107 END DO 109 108 END DO 110 109 END DO 111 112 CASE ( 1 ) ! no damping in the turbocline (avt > 5 cm2/s)110 ! 111 CASE ( 1 ) !== no damping in the turbocline (avt > 5 cm2/s) ==! 113 112 DO jk = 1, jpkm1 114 113 DO jj = 2, jpjm1 115 114 DO ji = fs_2, fs_jpim1 ! vector opt. 116 ztest = avt(ji,jj,jk) - 5.e-4 117 IF( ztest < 0. ) THEN 118 ztra = restotr(ji,jj,jk,jn) * ( trdta(ji,jj,jk,jn) - trb(ji,jj,jk,jn) ) 119 ELSE 120 ztra = 0.e0 115 IF( avt(ji,jj,jk) <= 5.e-4 ) THEN 116 ztra = restotr(ji,jj,jk) * ( trdta(ji,jj,jk,jn) - trb(ji,jj,jk,jn) ) 117 tra(ji,jj,jk,jn) = tra(ji,jj,jk,jn) + ztra 121 118 ENDIF 122 ! add the trends to the general tracer trends123 !! trn(ji,jj,jk,jn) = trn(ji,jj,jk,jn) + ztra * zdt124 tra(ji,jj,jk,jn) = tra(ji,jj,jk,jn) + ztra125 # if defined key_trc_diatrd126 ! save the trends for diagnostics127 IF (luttrd(jn)) trtrd(ji,jj,jk,ikeep(jn),jpdiatrc) = ztra128 # endif129 130 119 END DO 131 120 END DO 132 121 END DO 133 134 CASE ( 2 ) ! no damping in the mixed layer122 ! 123 CASE ( 2 ) !== no damping in the mixed layer ==! 135 124 DO jk = 1, jpkm1 136 125 DO jj = 2, jpjm1 … … 138 127 IF( fsdept(ji,jj,jk) >= hmlp (ji,jj) ) THEN 139 128 ztra = restotr(ji,jj,jk,jn) * ( trdta(ji,jj,jk,jn) - trb(ji,jj,jk,jn) ) 140 ELSE 141 ztra = 0.e0 142 ENDIF 143 ! add the trends to the general tracer trends 144 !! trn(ji,jj,jk,jn) = trn(ji,jj,jk,jn) + ztra * zdt 145 tra(ji,jj,jk,jn) = tra(ji,jj,jk,jn) + ztra 146 # if defined key_trc_diatrd 147 ! save the trends for diagnostics 148 IF (luttrd(jn)) trtrd(ji,jj,jk,ikeep(jn),jpdiatrc) = ztra 149 # endif 150 129 tra(ji,jj,jk,jn) = tra(ji,jj,jk,jn) + ztra 130 END IF 151 131 END DO 152 132 END DO 153 133 END DO 154 134 ! 155 135 END SELECT 156 136 ! 157 137 ENDIF 158 138 ! 159 139 IF( l_trdtrc ) THEN 160 140 ztrtrd(:,:,:) = tra(:,:,:,jn) - ztrtrd(:,:,:) 161 IF (luttrd(jn)) CALL trd_mod_trc( ztrtrd, jn, jptrc_trd_dmp, kt ) ! trends diags.141 CALL trd_tra( kt, 'TRC', jn, jptra_trd_dmp, ztrtrd ) 162 142 END IF 163 143 ! ! =========== 164 144 END DO ! tracer loop 165 145 ! ! =========== 166 167 IF( ln_ctl ) THEN! print mean trends (used for debugging)168 WRITE(charout, FMT="('dmp')")169 CALL prt_ctl_trc_info( charout)170 CALL prt_ctl_trc( tab4d=tra, mask=tmask, clinfo=ctrcnm,clinfo2='trd' )146 IF( l_trdtrc ) DEALLOCATE( ztrtrd ) 147 ! ! print mean trends (used for debugging) 148 IF( ln_ctl ) THEN 149 WRITE(charout, FMT="('dmp ')") ; CALL prt_ctl_trc_info(charout) 150 CALL prt_ctl_trc( tab4d=tra, mask=tmask, clinfo=ctrcnm, clinfo2='trd' ) 171 151 ENDIF 172 173 trb(:,:,:,:) = trn(:,:,:,:) 174 152 ! 175 153 END SUBROUTINE trc_dmp 176 154 … … 186 164 !!---------------------------------------------------------------------- 187 165 188 SELECT CASE ( ndmptr ) 189 190 CASE ( -1 ) ! ORCA: damping in Red & Med Seas only 191 IF(lwp) WRITE(numout,*) ' tracer damping in the Med & Red seas only' 192 193 CASE ( 1:90 ) ! Damping poleward of 'ndmptr' degrees 194 IF(lwp) WRITE(numout,*) ' tracer damping poleward of', ndmptr, ' degrees' 195 166 SELECT CASE ( nn_hdmp_tr ) 167 CASE ( -1 ) ; IF(lwp) WRITE(numout,*) ' tracer damping in the Med & Red seas only' 168 CASE ( 1:90 ) ; IF(lwp) WRITE(numout,*) ' tracer damping poleward of', nn_hdmp_tr, ' degrees' 196 169 CASE DEFAULT 197 WRITE(ctmp1,*) ' bad flag value for n dmptr = ', ndmptr170 WRITE(ctmp1,*) ' bad flag value for nn_hdmp_tr = ', nn_hdmp_tr 198 171 CALL ctl_stop(ctmp1) 199 200 172 END SELECT 201 173 202 203 SELECT CASE ( nmldmptr ) 204 205 CASE ( 0 ) ! newtonian damping throughout the water column 206 IF(lwp) WRITE(numout,*) ' tracer damping throughout the water column' 207 208 CASE ( 1 ) ! no damping in the turbocline (avt > 5 cm2/s) 209 IF(lwp) WRITE(numout,*) ' no tracer damping in the turbocline' 210 211 CASE ( 2 ) ! no damping in the mixed layer 212 IF(lwp) WRITE(numout,*) ' no tracer damping in the mixed layer' 213 174 SELECT CASE ( nn_zdmp_tr ) 175 CASE ( 0 ) ; IF(lwp) WRITE(numout,*) ' tracer damping throughout the water column' 176 CASE ( 1 ) ; IF(lwp) WRITE(numout,*) ' no tracer damping in the turbocline (avt > 5 cm2/s)' 177 CASE ( 2 ) ; IF(lwp) WRITE(numout,*) ' no tracer damping in the mixed layer' 214 178 CASE DEFAULT 215 WRITE(ctmp1,*) ' bad flag value for nmldmptr = ', nmldmptr179 WRITE(ctmp1,*) 'bad flag value for nn_zdmp_tr = ', nn_zdmp_tr 216 180 CALL ctl_stop(ctmp1) 217 218 181 END SELECT 219 182 220 ! Damping coefficients initialization 221 ! ----------------------------------- 222 IF( lzoom ) THEN 223 CALL trccof_zoom 224 ELSE 225 CALL trccof 183 IF( .NOT. lk_dtatrc ) & 184 & CALL ctl_stop( 'no passive tracer data define key_dtatrc' ) 185 186 IF( .NOT. lk_tradmp ) & 187 & CALL ctl_stop( 'passive trace damping need key_tradmp to compute damping coef.' ) 188 ! 189 ! ! Damping coefficients initialization 190 IF( lzoom ) THEN ; CALL dtacof_zoom( restotr ) 191 ELSE ; CALL dtacof( nn_hdmp_tr, rn_surf_tr, rn_bot_tr, rn_dep_tr, & 192 & nn_file_tr, 'TRC' , restotr ) 226 193 ENDIF 227 194 ! 228 195 END SUBROUTINE trc_dmp_init 229 230 231 SUBROUTINE trccof_zoom232 !!----------------------------------------------------------------------233 !! *** ROUTINE trccof_zoom ***234 !!235 !! ** Purpose : Compute the damping coefficient for zoom domain236 !!237 !! ** Method : - set along closed boundary due to zoom a damping over238 !! 6 points with a max time scale of 5 days.239 !! - ORCA arctic/antarctic zoom: set the damping along240 !! south/north boundary over a latitude strip.241 !!242 !! ** Action : - restotr, the damping coeff. passive tracers243 !!244 !! History :245 !! 9.0 ! 03-09 (G. Madec) Original code246 !! 9.0 ! 04-03 (C. Ethe) adapted for passive tracers247 !!----------------------------------------------------------------------248 !! * Local declarations249 INTEGER :: ji, jj, jk, jn ! dummy loop indices250 REAL(wp) :: &251 zlat, zlat0, zlat1, zlat2 ! temporary scalar252 REAL(wp), DIMENSION(6) :: &253 zfact ! temporary workspace254 !!----------------------------------------------------------------------255 256 zfact(1) = 1.257 zfact(2) = 1.258 zfact(3) = 11./12.259 zfact(4) = 8./12.260 zfact(5) = 4./12.261 zfact(6) = 1./12.262 zfact(:) = zfact(:) / ( 5. * rday ) ! 5 days max restoring time scale263 264 restotr(:,:,:,:) = 0.e0265 266 ! damping along the forced closed boundary over 6 grid-points267 DO jn = 1, 6268 IF( lzoom_w ) restotr( mi0(jn+jpizoom):mi1(jn+jpizoom), : , : , : ) = zfact(jn) ! west closed269 IF( lzoom_s ) restotr( : , mj0(jn+jpjzoom):mj1(jn+jpjzoom), : , : ) = zfact(jn) ! south closed270 IF( lzoom_e ) restotr( mi0(jpiglo+jpizoom-1-jn):mi1(jpiglo+jpizoom-1-jn) , : , : , : ) &271 & = zfact(jn) ! east closed272 IF( lzoom_n ) restotr( : , mj0(jpjglo+jpjzoom-1-jn):mj1(jpjglo+jpjzoom-1-jn) , : , : ) &273 & = zfact(jn) ! north closed274 END DO275 276 277 IF( lzoom_arct .AND. lzoom_anta ) THEN278 279 ! ====================================================280 ! ORCA configuration : arctic zoom or antarctic zoom281 ! ====================================================282 283 IF(lwp) WRITE(numout,*)284 IF(lwp .AND. lzoom_arct ) WRITE(numout,*) ' trccof_zoom : ORCA Arctic zoom'285 IF(lwp .AND. lzoom_arct ) WRITE(numout,*) ' trccof_zoom : ORCA Antarctic zoom'286 IF(lwp) WRITE(numout,*)287 288 ! ... Initialization :289 ! zlat0 : latitude strip where resto decreases290 ! zlat1 : resto = 1 before zlat1291 ! zlat2 : resto decreases from 1 to 0 between zlat1 and zlat2292 restotr(:,:,:,:) = 0.e0293 zlat0 = 10.294 zlat1 = 30.295 zlat2 = zlat1 + zlat0296 297 ! ... Compute arrays resto ; value for internal damping : 5 days298 DO jn = 1, jptra299 DO jk = 2, jpkm1300 DO jj = 1, jpj301 DO ji = 1, jpi302 zlat = ABS( gphit(ji,jj) )303 IF ( zlat1 <= zlat .AND. zlat <= zlat2 ) THEN304 restotr(ji,jj,jk,jn) = 0.5 * ( 1./(5.*rday) ) * &305 ( 1. - COS(rpi*(zlat2-zlat)/zlat0) )306 ELSE IF ( zlat < zlat1 ) THEN307 restotr(ji,jj,jk,jn) = 1./(5.*rday)308 ENDIF309 END DO310 END DO311 END DO312 END DO313 314 ENDIF315 316 ! ... Mask resto array317 DO jn = 1, jptra318 restotr(:,:,:,jn) = restotr(:,:,:,jn) * tmask(:,:,:)319 END DO320 321 322 END SUBROUTINE trccof_zoom323 324 SUBROUTINE trccof325 !!----------------------------------------------------------------------326 !! *** ROUTINE trccof ***327 !!328 !! ** Purpose : Compute the damping coefficient329 !!330 !! ** Method : Arrays defining the damping are computed for each grid331 !! point passive tracers (restotr)332 !! Damping depends on distance to coast, depth and latitude333 !!334 !! ** Action : - restotr, the damping coeff. for passive tracers335 !!336 !! History :337 !! 5.0 ! 91-03 (O. Marti, G. Madec) Original code338 !! ! 92-06 (M. Imbard) doctor norme339 !! ! 96-01 (G. Madec) statement function for e3340 !! ! 98-07 (M. Imbard, G. Madec) ORCA version341 !! ! 00-08 (G. Madec, D. Ludicone)342 !! 8.2 ! 04-03 (H. Loukos) adapted for passive tracers343 !! ! 04-02 (O. Aumont, C. Ethe) rewritten for debuging and update344 !!----------------------------------------------------------------------345 !! * Modules used346 USE iom347 USE ioipsl348 349 !! * Local declarations350 INTEGER :: ji, jj, jk, jn ! dummy loop indices351 INTEGER :: itime352 INTEGER :: ii0, ii1, ij0, ij1 ! " "353 INTEGER :: &354 idmp, & ! logical unit for file restoring damping term355 icot ! logical unit for file distance to the coast356 357 CHARACTER (len=32) :: clname3358 REAL(wp) :: &359 zdate0, zinfl, zlon, & ! temporary scalars360 zlat, zlat0, zlat1, zlat2, & ! " "361 zsdmp, zbdmp ! " "362 REAL(wp), DIMENSION(jpk) :: &363 gdept, zhfac364 REAL(wp), DIMENSION(jpi,jpj) :: &365 zmrs366 REAL(wp), DIMENSION(jpi,jpj,jpk) :: &367 zdct368 !!----------------------------------------------------------------------369 370 ! ====================================371 ! ORCA configuration : global domain372 ! ====================================373 374 IF(lwp) WRITE(numout,*)375 IF(lwp) WRITE(numout,*) ' trccof : Global domain of ORCA'376 IF(lwp) WRITE(numout,*) ' ------------------------------'377 378 379 ! ... Initialization :380 ! zdct() : distant to the coastline381 ! resto() : array of restoring coeff.382 383 zdct (:,:,:) = 0.e0384 restotr(:,:,:,:) = 0.e0385 386 387 IF ( ndmptr > 0 ) THEN388 389 ! ------------------------------------390 ! Damping poleward of 'ndmptr' degrees391 ! ------------------------------------392 393 IF(lwp) WRITE(numout,*)394 IF(lwp) WRITE(numout,*) ' Damping poleward of ', ndmptr,' deg.'395 IF(lwp) WRITE(numout,*)396 397 ! ... Distance to coast (zdct)398 399 IF(lwp) WRITE(numout,*)400 IF(lwp) WRITE(numout,*) ' dtacof : distance to coast file'401 CALL iom_open ( 'dist.coast.trc.nc', icot )402 IF( icot > 0 ) THEN403 CALL iom_get ( icot, jpdom_data, 'Tcoast', zdct )404 CALL iom_close (icot)405 ELSE406 ! ... Compute and save the distance-to-coast array (output in zdct)407 CALL cofdis( zdct )408 ENDIF409 410 411 ! ... Compute arrays resto412 ! zinfl : distance of influence for damping term413 ! zlat0 : latitude strip where resto decreases414 ! zlat1 : resto = 0 between -zlat1 and zlat1415 ! zlat2 : resto increases from 0 to 1 between |zlat1| and |zlat2|416 ! and resto = 1 between |zlat2| and 90 deg.417 zinfl = 1000.e3418 zlat0 = 10419 zlat1 = ndmptr420 zlat2 = zlat1 + zlat0421 422 DO jn = 1, jptra423 DO jj = 1, jpj424 DO ji = 1, jpi425 zlat = ABS( gphit(ji,jj) )426 IF ( zlat1 <= zlat .AND. zlat <= zlat2 ) THEN427 restotr(ji,jj,1,jn) = 0.5 * ( 1. - COS(rpi*(zlat-zlat1)/zlat0 ) )428 ELSEIF ( zlat > zlat2 ) THEN429 restotr(ji,jj,1,jn) = 1.430 ENDIF431 END DO432 END DO433 END DO434 435 ! ... North Indian ocean (20N/30N x 45E/100E) : resto=0436 IF ( ndmptr == 20 ) THEN437 DO jn = 1, jptra438 DO jj = 1, jpj439 DO ji = 1, jpi440 zlat = gphit(ji,jj)441 zlon = MOD( glamt(ji,jj), 360. )442 IF ( zlat1 < zlat .AND. zlat < zlat2 .AND. &443 45. < zlon .AND. zlon < 100. ) THEN444 restotr(ji,jj,1,jn) = 0.445 ENDIF446 END DO447 END DO448 END DO449 ENDIF450 451 zsdmp = 1./(sdmptr * rday)452 zbdmp = 1./(bdmptr * rday)453 DO jn = 1, jptra454 DO jk = 2, jpkm1455 DO jj = 1, jpj456 DO ji = 1, jpi457 zdct(ji,jj,jk) = MIN( zinfl, zdct(ji,jj,jk) )458 459 ! ... Decrease the value in the vicinity of the coast460 restotr(ji,jj,jk,jn) = restotr(ji,jj,1,jn)*0.5 &461 & * ( 1. - COS( rpi*zdct(ji,jj,jk)/zinfl) )462 463 ! ... Vertical variation from zsdmp (sea surface) to zbdmp (bottom)464 restotr(ji,jj,jk,jn) = restotr(ji,jj,jk,jn) &465 & * ( zbdmp + (zsdmp-zbdmp)*EXP(-fsdept(ji,jj,jk)/hdmptr) )466 END DO467 END DO468 END DO469 END DO470 471 ENDIF472 473 474 IF( cp_cfg == "orca" .AND. ( ndmptr > 0 .OR. ndmptr == -1 ) ) THEN475 476 ! ! =========================477 ! ! Med and Red Sea damping478 ! ! =========================479 IF(lwp)WRITE(numout,*)480 IF(lwp)WRITE(numout,*) ' ORCA configuration: Damping in Med and Red Seas'481 482 483 zmrs(:,:) = 0.e0 ! damping term on the Med or Red Sea484 485 SELECT CASE ( jp_cfg )486 ! ! =======================487 CASE ( 4 ) ! ORCA_R4 configuration488 ! ! =======================489 490 ! Mediterranean Sea491 ij0 = 50 ; ij1 = 56492 ii0 = 81 ; ii1 = 91 ; zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1.e0493 ij0 = 50 ; ij1 = 55494 ii0 = 75 ; ii1 = 80 ; zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1.e0495 ij0 = 52 ; ij1 = 53496 ii0 = 70 ; ii1 = 74 ; zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1.e0497 ! Smooth transition from 0 at surface to 1./rday at the 18th level in Med and Red Sea498 DO jk = 1, 17499 zhfac (jk) = 0.5*( 1.- COS( rpi*(jk-1)/16. ) ) / rday500 END DO501 DO jk = 18, jpkm1502 zhfac (jk) = 1./rday503 END DO504 505 ! ! =======================506 CASE ( 2 ) ! ORCA_R2 configuration507 ! ! =======================508 509 ! Mediterranean Sea510 ij0 = 96 ; ij1 = 110511 ii0 = 157 ; ii1 = 181 ; zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1.e0512 ij0 = 100 ; ij1 = 110513 ii0 = 144 ; ii1 = 156 ; zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1.e0514 ij0 = 100 ; ij1 = 103515 ii0 = 139 ; ii1 = 143 ; zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1.e0516 ! Decrease before Gibraltar Strait517 ij0 = 101 ; ij1 = 102518 ii0 = 139 ; ii1 = 141 ; zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.e0519 ii0 = 142 ; ii1 = 142 ; zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1.e0 / 90.e0520 ii0 = 143 ; ii1 = 143 ; zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.40e0521 ii0 = 144 ; ii1 = 144 ; zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.75e0522 ! Red Sea523 ij0 = 87 ; ij1 = 96524 ii0 = 147 ; ii1 = 163 ; zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1.e0525 ! Decrease before Bab el Mandeb Strait526 ij0 = 91 ; ij1 = 91527 ii0 = 153 ; ii1 = 160 ; zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.80e0528 ij0 = 90 ; ij1 = 90529 ii0 = 153 ; ii1 = 160 ; zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.40e0530 ij0 = 89 ; ij1 = 89531 ii0 = 158 ; ii1 = 160 ; zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1.e0 / 90.e0532 ij0 = 88 ; ij1 = 88533 ii0 = 160 ; ii1 = 163 ; zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.e0534 ! Smooth transition from 0 at surface to 1./rday at the 18th level in Med and Red Sea535 DO jk = 1, 17536 zhfac (jk) = 0.5*( 1.- COS( rpi*(jk-1)/16. ) ) / rday537 END DO538 DO jk = 18, jpkm1539 zhfac (jk) = 1./rday540 END DO541 542 ! ! =======================543 CASE ( 05 ) ! ORCA_R05 configuration544 ! ! =======================545 546 ! Mediterranean Sea547 ii0 = 568 ; ii1 = 574548 ij0 = 324 ; ij1 = 333 ; zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1.e0549 ii0 = 575 ; ii1 = 658550 ij0 = 314 ; ij1 = 366 ; zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1.e0551 ! Black Sea (remaining part552 ii0 = 641 ; ii1 = 651553 ij0 = 367 ; ij1 = 372 ; zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1.e0554 ! Decrease before Gibraltar Strait555 ii0 = 324 ; ii1 = 333556 ij0 = 565 ; ij1 = 565 ; zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1.e0 / 90.e0557 ij0 = 566 ; ij1 = 566 ; zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.40558 ij0 = 567 ; ij1 = 567 ; zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.75559 ! Red Sea560 ii0 = 641 ; ii1 = 665561 ij0 = 270 ; ij1 = 310 ; zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1.e0562 ! Decrease before Bab el Mandeb Strait563 ii0 = 666 ; ii1 = 675564 ij0 = 270 ; ij1 = 290565 DO ji = mi0(ii0), mi1(ii1)566 zmrs( ji , mj0(ij0):mj1(ij1) ) = 0.1 * ABS( FLOAT(ji - mi1(ii1)) )567 END DO568 zsdmp = 1./(sdmptr * rday)569 zbdmp = 1./(bdmptr * rday)570 DO jk = 1, jpk571 zhfac (jk) = ( zbdmp + (zsdmp-zbdmp) * EXP(-fsdept(1,1,jk)/hdmptr) )572 END DO573 574 ! ! ========================575 CASE ( 025 ) ! ORCA_R025 configuration576 577 CALL ctl_stop( ' Not yet implemented in ORCA_R025' )578 579 END SELECT580 581 DO jn = 1, jptra582 DO jk = 1, jpkm1583 restotr(:,:,jk,jn) = zmrs(:,:) * zhfac(jk) + ( 1. - zmrs(:,:) ) * restotr(:,:,jk,jn)584 END DO585 586 ! Mask resto array and set to 0 first and last levels587 restotr(:,:, : ,jn) = restotr(:,:,:,jn) * tmask(:,:,:)588 restotr(:,:, 1 ,jn) = 0.e0589 restotr(:,:,jpk,jn) = 0.e0590 END DO591 592 ELSE593 ! ------------594 ! No damping595 ! ------------596 CALL ctl_stop( 'Choose a correct value of ndmp or DO NOT defined key_tradmp' )597 598 ENDIF599 600 ! ----------------------------601 ! Create Print damping array602 ! ----------------------------603 604 ! ndmpftr : = 1 create a damping.coeff NetCDF file605 606 IF( ndmpftr == 1 ) THEN607 DO jn = 1, jptra608 IF(lwp) WRITE(numout,*) ' create damping.coeff.nc file ',jn609 itime = 0610 clname3 = 'damping.coeff'//ctrcnm(jn)611 CALL ymds2ju( 0 , 1 , 1 , 0.e0 , zdate0 )612 CALL restini( 'NONE', jpi , jpj , glamt, gphit, &613 & jpk , gdept , clname3, itime, zdate0, &614 & rdt , idmp , domain_id=nidom)615 CALL restput( idmp, 'Resto', jpi, jpj, jpk, 0 , restotr(:,:,:,jn) )616 CALL restclo( idmp )617 END DO618 ENDIF619 620 621 END SUBROUTINE trccof622 623 624 SUBROUTINE cofdis ( pdct )625 !!----------------------------------------------------------------------626 !! *** ROUTINE cofdis ***627 !!628 !! ** Purpose : Compute the distance between ocean T-points and the629 !! ocean model coastlines. Save the distance in a NetCDF file.630 !!631 !! ** Method : For each model level, the distance-to-coast is632 !! computed as follows :633 !! - The coastline is defined as the serie of U-,V-,F-points634 !! that are at the ocean-land bound.635 !! - For each ocean T-point, the distance-to-coast is then636 !! computed as the smallest distance (on the sphere) between the637 !! T-point and all the coastline points.638 !! - For land T-points, the distance-to-coast is set to zero.639 !! C A U T I O N : Computation not yet implemented in mpp case.640 !!641 !! ** Action : - pdct, distance to the coastline (argument)642 !! - NetCDF file 'trc.dist.coast.nc'643 !!644 !! History :645 !! 7.0 ! 01-02 (M. Imbard) Original code646 !! 8.1 ! 01-02 (G. Madec, E. Durand)647 !! 8.5 ! 02-08 (G. Madec, E. Durand) Free form, F90648 !!----------------------------------------------------------------------649 !! * Modules used650 USE ioipsl651 652 !! * Arguments653 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( out ) :: &654 pdct ! distance to the coastline655 656 !! * local declarations657 INTEGER :: ji, jj, jk, jl ! dummy loop indices658 INTEGER :: iju, ijt ! temporary integers659 INTEGER :: icoast, itime660 INTEGER :: &661 icot ! logical unit for file distance to the coast662 LOGICAL, DIMENSION(jpi,jpj) :: &663 llcotu, llcotv, llcotf ! ???664 CHARACTER (len=32) :: clname665 REAL(wp) :: zdate0666 REAL(wp), DIMENSION(jpi,jpj) :: &667 zxt, zyt, zzt, & ! cartesian coordinates for T-points668 zmask669 REAL(wp), DIMENSION(3*jpi*jpj) :: &670 zxc, zyc, zzc, zdis ! temporary workspace671 !!----------------------------------------------------------------------672 673 ! 0. Initialization674 ! -----------------675 IF(lwp) WRITE(numout,*)676 IF(lwp) WRITE(numout,*) 'cofdis : compute the distance to coastline'677 IF(lwp) WRITE(numout,*) '~~~~~~'678 IF(lwp) WRITE(numout,*)679 IF( lk_mpp ) &680 & CALL ctl_stop(' Computation not yet implemented with key_mpp_...', &681 & ' Rerun the code on another computer or ', &682 & ' create the "dist.coast.nc" file using IDL' )683 684 685 pdct(:,:,:) = 0.e0686 zxt(:,:) = cos( rad * gphit(:,:) ) * cos( rad * glamt(:,:) )687 zyt(:,:) = cos( rad * gphit(:,:) ) * sin( rad * glamt(:,:) )688 zzt(:,:) = sin( rad * gphit(:,:) )689 690 691 ! 1. Loop on vertical levels692 ! --------------------------693 ! ! ===============694 DO jk = 1, jpkm1 ! Horizontal slab695 ! ! ===============696 ! Define the coastline points (U, V and F)697 DO jj = 2, jpjm1698 DO ji = 2, jpim1699 zmask(ji,jj) = ( tmask(ji,jj+1,jk) + tmask(ji+1,jj+1,jk) &700 & + tmask(ji,jj ,jk) + tmask(ji+1,jj ,jk) )701 llcotu(ji,jj) = ( tmask(ji,jj, jk) + tmask(ji+1,jj ,jk) == 1. )702 llcotv(ji,jj) = ( tmask(ji,jj ,jk) + tmask(ji ,jj+1,jk) == 1. )703 llcotf(ji,jj) = ( zmask(ji,jj) > 0. ) .AND. ( zmask(ji,jj) < 4. )704 END DO705 END DO706 707 ! Lateral boundaries conditions708 llcotu(:, 1 ) = umask(:, 2 ,jk) == 1709 llcotu(:,jpj) = umask(:,jpjm1,jk) == 1710 llcotv(:, 1 ) = vmask(:, 2 ,jk) == 1711 llcotv(:,jpj) = vmask(:,jpjm1,jk) == 1712 llcotf(:, 1 ) = fmask(:, 2 ,jk) == 1713 llcotf(:,jpj) = fmask(:,jpjm1,jk) == 1714 715 IF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6 ) THEN716 llcotu( 1 ,:) = llcotu(jpim1,:)717 llcotu(jpi,:) = llcotu( 2 ,:)718 llcotv( 1 ,:) = llcotv(jpim1,:)719 llcotv(jpi,:) = llcotv( 2 ,:)720 llcotf( 1 ,:) = llcotf(jpim1,:)721 llcotf(jpi,:) = llcotf( 2 ,:)722 ELSE723 llcotu( 1 ,:) = umask( 2 ,:,jk) == 1724 llcotu(jpi,:) = umask(jpim1,:,jk) == 1725 llcotv( 1 ,:) = vmask( 2 ,:,jk) == 1726 llcotv(jpi,:) = vmask(jpim1,:,jk) == 1727 llcotf( 1 ,:) = fmask( 2 ,:,jk) == 1728 llcotf(jpi,:) = fmask(jpim1,:,jk) == 1729 ENDIF730 IF( nperio == 3 .OR. nperio == 4 ) THEN731 DO ji = 1, jpim1732 iju = jpi - ji + 1733 llcotu(ji,jpj ) = llcotu(iju,jpj-2)734 llcotf(ji,jpj-1) = llcotf(iju,jpj-2)735 llcotf(ji,jpj ) = llcotf(iju,jpj-3)736 END DO737 DO ji = jpi/2, jpi-1738 iju = jpi - ji + 1739 llcotu(ji,jpjm1) = llcotu(iju,jpjm1)740 END DO741 DO ji = 2, jpi742 ijt = jpi - ji + 2743 llcotv(ji,jpj-1) = llcotv(ijt,jpj-2)744 llcotv(ji,jpj ) = llcotv(ijt,jpj-3)745 END DO746 ENDIF747 IF( nperio == 5 .OR. nperio == 6 ) THEN748 DO ji = 1, jpim1749 iju = jpi - ji750 llcotu(ji,jpj ) = llcotu(iju,jpj-1)751 llcotf(ji,jpj ) = llcotf(iju,jpj-2)752 END DO753 DO ji = jpi/2, jpi-1754 iju = jpi - ji755 llcotf(ji,jpjm1) = llcotf(iju,jpjm1)756 END DO757 DO ji = 1, jpi758 ijt = jpi - ji + 1759 llcotv(ji,jpj ) = llcotv(ijt,jpj-1)760 END DO761 DO ji = jpi/2+1, jpi762 ijt = jpi - ji + 1763 llcotv(ji,jpjm1) = llcotv(ijt,jpjm1)764 END DO765 ENDIF766 767 ! Compute cartesian coordinates of coastline points768 ! and the number of coastline points769 770 icoast = 0771 DO jj = 1, jpj772 DO ji = 1, jpi773 IF( llcotf(ji,jj) ) THEN774 icoast = icoast + 1775 zxc(icoast) = COS( rad*gphif(ji,jj) ) * COS( rad*glamf(ji,jj) )776 zyc(icoast) = COS( rad*gphif(ji,jj) ) * SIN( rad*glamf(ji,jj) )777 zzc(icoast) = SIN( rad*gphif(ji,jj) )778 ENDIF779 IF( llcotu(ji,jj) ) THEN780 icoast = icoast+1781 zxc(icoast) = COS( rad*gphiu(ji,jj) ) * COS( rad*glamu(ji,jj) )782 zyc(icoast) = COS( rad*gphiu(ji,jj) ) * SIN( rad*glamu(ji,jj) )783 zzc(icoast) = SIN( rad*gphiu(ji,jj) )784 ENDIF785 IF( llcotv(ji,jj) ) THEN786 icoast = icoast+1787 zxc(icoast) = COS( rad*gphiv(ji,jj) ) * COS( rad*glamv(ji,jj) )788 zyc(icoast) = COS( rad*gphiv(ji,jj) ) * SIN( rad*glamv(ji,jj) )789 zzc(icoast) = SIN( rad*gphiv(ji,jj) )790 ENDIF791 END DO792 END DO793 794 ! Distance for the T-points795 796 DO jj = 1, jpj797 DO ji = 1, jpi798 IF( tmask(ji,jj,jk) == 0. ) THEN799 pdct(ji,jj,jk) = 0.800 ELSE801 DO jl = 1, icoast802 zdis(jl) = ( zxt(ji,jj) - zxc(jl) )**2 &803 + ( zyt(ji,jj) - zyc(jl) )**2 &804 + ( zzt(ji,jj) - zzc(jl) )**2805 END DO806 pdct(ji,jj,jk) = ra * SQRT( MINVAL( zdis(1:icoast) ) )807 ENDIF808 END DO809 END DO810 ! ! ===============811 END DO ! End of slab812 ! ! ===============813 814 815 ! 2. Create the distance to the coast file in NetCDF format816 ! ----------------------------------------------------------817 clname = 'trc.dist.coast'818 itime = 0819 CALL ymds2ju( 0 , 1 , 1 , 0.e0 , zdate0 )820 CALL restini( 'NONE', jpi , jpj , glamt, gphit , &821 jpk , gdept , clname, itime, zdate0, &822 rdt , icot , domain_id=nidom )823 CALL restput( icot, 'Tcoast', jpi, jpj, jpk, 0, pdct )824 CALL restclo( icot )825 826 END SUBROUTINE cofdis827 828 196 #else 829 197 !!---------------------------------------------------------------------- … … 837 205 END SUBROUTINE trc_dmp 838 206 #endif 839 840 207 !!====================================================================== 841 208 END MODULE trcdmp -
branches/DEV_r2006_merge_TRA_TRC/NEMO/TOP_SRC/TRP/trcnxt.F90
r1271 r2030 5 5 !!====================================================================== 6 6 !!====================================================================== 7 !! History : 7.0 ! 91-11 (G. Madec) Original code 8 !! ! 93-03 (M. Guyon) symetrical conditions 9 !! ! 95-02 (M. Levy) passive tracers 10 !! ! 96-02 (G. Madec & M. Imbard) opa release 8.0 11 !! 8.0 ! 96-04 (A. Weaver) Euler forward step 12 !! 8.2 ! 99-02 (G. Madec, N. Grima) semi-implicit pressure grad. 13 !! 8.5 ! 02-08 (G. Madec) F90: Free form and module 14 !! ! 02-11 (C. Talandier, A-M Treguier) Open boundaries 15 !! 9.0 ! 04-03 (C. Ethe) passive tracers 16 !! ! 07-02 (C. Deltel) Diagnose ML trends for passive tracers 7 !! History : 7.0 ! 1991-11 (G. Madec) Original code 8 !! ! 1993-03 (M. Guyon) symetrical conditions 9 !! ! 1995-02 (M. Levy) passive tracers 10 !! ! 1996-02 (G. Madec & M. Imbard) opa release 8.0 11 !! 8.0 ! 1996-04 (A. Weaver) Euler forward step 12 !! 8.2 ! 1999-02 (G. Madec, N. Grima) semi-implicit pressure grad. 13 !! NEMO 1.0 ! 2002-08 (G. Madec) F90: Free form and module 14 !! ! 2002-08 (G. Madec) F90: Free form and module 15 !! ! 2002-11 (C. Talandier, A-M Treguier) Open boundaries 16 !! ! 2004-03 (C. Ethe) passive tracers 17 !! ! 2007-02 (C. Deltel) Diagnose ML trends for passive tracers 18 !! 2.0 ! 2006-02 (L. Debreu, C. Mazauric) Agrif implementation 19 !! 3.0 ! 2008-06 (G. Madec) time stepping always done in trazdf 20 !! 3.1 ! 2009-02 (G. Madec, R. Benshila) re-introduce the vvl option 21 !! 3.3 ! 2010-06 (C. Ethe, G. Madec) Merge TRA-TRC 17 22 !!---------------------------------------------------------------------- 18 23 #if defined key_top … … 24 29 !! * Modules used 25 30 USE oce_trc ! ocean dynamics and tracers variables 26 USE tr p_trc ! ocean passive tracers variables31 USE trc ! ocean passive tracers variables 27 32 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 28 USE trctrp_lec ! pasive tracers transport29 33 USE prtctl_trc ! Print control for debbuging 30 USE trdmld_trc 31 USE trdmld_trc_oce 34 USE trdmod_oce 35 USE trdtra 36 USE tranxt 32 37 USE agrif_top_update 33 38 USE agrif_top_interp … … 38 43 !! * Routine accessibility 39 44 PUBLIC trc_nxt ! routine called by step.F90 45 46 REAL(wp), DIMENSION(jpk) :: r2dt_t 40 47 !!---------------------------------------------------------------------- 41 48 !! TOP 1.0 , LOCEAN-IPSL (2005) … … 70 77 !! ** Action : - update trb, trn 71 78 !!---------------------------------------------------------------------- 72 USE oce, ONLY : ztrtrd => ua ! use ua as 3D workspace73 79 !! * Arguments 74 80 INTEGER, INTENT( in ) :: kt ! ocean time-step index 75 81 !! * Local declarations 76 INTEGER :: j i, jj, jk, jn ! dummy loop indices82 INTEGER :: jk, jn ! dummy loop indices 77 83 REAL(wp) :: zfact ! temporary scalar 78 84 CHARACTER (len=22) :: charout 85 REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: ztrdt 79 86 !!---------------------------------------------------------------------- 80 87 … … 84 91 ENDIF 85 92 93 ! Update after tracer on domain lateral boundaries 86 94 DO jn = 1, jptra 95 CALL lbc_lnk( tra(:,:,:,jn), 'T', 1. ) 96 END DO 87 97 88 ! 0. Lateral boundary conditions on tra (T-point, unchanged sign)89 ! ---------------------------------============90 CALL lbc_lnk( tra(:,:,:,jn), 'T', 1. )91 92 ! ! ===============93 DO jk = 1, jpk ! Horizontal slab94 ! ! ===============95 ! 1. Leap-frog scheme (only in explicit case, otherwise the96 ! ------------------- time stepping is already done in trczdf)97 IF( l_trczdf_exp .AND. ( ln_trcadv_cen2 .OR. ln_trcadv_tvd) ) THEN98 zfact = 2. * rdttra(jk) * FLOAT(ndttrc)99 IF( neuler == 0 .AND. kt == nittrc000 ) zfact = rdttra(jk) * FLOAT(ndttrc)100 tra(:,:,jk,jn) = ( trb(:,:,jk,jn) + zfact * tra(:,:,jk,jn) ) * tmask(:,:,jk)101 ENDIF102 103 END DO104 98 105 99 #if defined key_obc 106 CALL ctl_stop( ' Passive tracers and Open Boundary condition can not be used together ' & 107 & ' Check in trc_nxt routine' ) 100 !! CALL obc_trc( kt ) ! OBC open boundaries 101 #endif 102 #if defined key_bdy 103 !! CALL bdy_trc( kt ) ! BDY open boundaries 104 #endif 105 #if defined key_agrif 106 CALL Agrif_trc ! AGRIF zoom boundaries 108 107 #endif 109 108 109 110 ! set time step size (Euler/Leapfrog) 111 IF( neuler == 0 .AND. kt == nittrc000) THEN ; r2dt_t(:) = rdttra(:) * FLOAT( nn_dttrc ) ! at nit000 (Euler) 112 ELSEIF( kt <= nittrc000 + 1 ) THEN ; r2dt_t(:) = 2.* rdttra(:) * FLOAT( nn_dttrc ) ! at nit000 or nit000+1 (Leapfrog) 113 ENDIF 114 115 ! trends computation initialisation 116 IF( l_trdtrc ) THEN 117 ALLOCATE( ztrdt(jpi,jpj,jpk,jptra) ) !* store now fields before applying the Asselin filter 118 ztrdt(:,:,:,:) = trn(:,:,:,:) 119 ENDIF 120 121 ! Leap-Frog + Asselin filter time stepping 122 IF( lk_vvl ) THEN ; CALL tra_nxt_vvl( kt , nittrc000, & 123 & trb, trn, tra, jptra ) ! variable volume level (vvl) 124 ELSE ; CALL tra_nxt_fix( kt , nittrc000, & 125 & trb, trn, tra, jptra ) ! fixed volume level 126 ENDIF 127 110 128 #if defined key_agrif 111 ! ! =============== 112 END DO ! End of slab 113 ! ! =============== 114 ! Interp tracers on boundaries (coarse => fine) 115 CALL Agrif_trc 116 ! ! =============== 117 DO jn = 1, jptra ! Horizontal slab 118 ! ! =============== 119 #endif 129 ! Update tracer at AGRIF zoom boundaries 130 IF( .NOT.Agrif_Root() ) CALL Agrif_Update_Trc( kt ) ! children only 131 #endif 120 132 121 DO jk = 1, jpk 122 123 ! 2. Time filter and swap of arrays 124 ! --------------------------------- 125 IF ( ln_trcadv_cen2 .OR. ln_trcadv_tvd ) THEN 126 127 IF( neuler == 0 .AND. kt == nittrc000 ) THEN 128 DO jj = 1, jpj 129 DO ji = 1, jpi 130 trb(ji,jj,jk,jn) = trn(ji,jj,jk,jn) 131 trn(ji,jj,jk,jn) = tra(ji,jj,jk,jn) 132 tra(ji,jj,jk,jn) = 0. 133 END DO 134 END DO 135 IF( l_trdtrc ) ztrtrd(:,:,:) = 0.e0 ! no trend 136 ELSE 137 IF( l_trdtrc ) THEN ! Asselin trend 138 DO jj = 1, jpj 139 DO ji = 1, jpi 140 ztrtrd(ji,jj,jk) = atfp * ( trb(ji,jj,jk,jn) - 2*trn(ji,jj,jk,jn) + tra(ji,jj,jk,jn) ) 141 END DO 142 END DO 143 ENDIF 144 145 DO jj = 1, jpj 146 DO ji = 1, jpi 147 trb(ji,jj,jk,jn) = atfp * ( trb(ji,jj,jk,jn) + tra(ji,jj,jk,jn) ) + atfp1 * trn(ji,jj,jk,jn) 148 trn(ji,jj,jk,jn) = tra(ji,jj,jk,jn) 149 tra(ji,jj,jk,jn) = 0. 150 END DO 151 END DO 152 ENDIF 153 ELSE ! >> EULER-FORWARD schemes (SMOLAR, MUSCL) 154 IF( l_trdtrc ) ztrtrd(:,:,:) = 0.e0 ! no trend 155 156 DO jj = 1, jpj 157 DO ji = 1, jpi 158 trb(ji,jj,jk,jn) = tra(ji,jj,jk,jn) 159 trn(ji,jj,jk,jn) = tra(ji,jj,jk,jn) 160 tra(ji,jj,jk,jn) = 0. 161 END DO 162 END DO 163 164 ENDIF 165 ! ! =============== 166 END DO ! End of slab 167 ! ! =============== 168 169 IF( l_trdtrc ) THEN ! trends 170 DO jk = 1, jpk 171 zfact = 2. * rdttra(jk) * FLOAT(ndttrc) 172 ztrtrd(:,:,jk) = ztrtrd(:,:,jk) / zfact ! n.b. ztrtrd=0 in Euler-forward case 133 ! trends computation 134 IF( l_trdtrc ) THEN ! trends 135 DO jn = 1, jptra 136 DO jk = 1, jpkm1 137 zfact = 1.e0 / r2dt_t(jk) 138 ztrdt(:,:,jk,jn) = ( trb(:,:,jk,jn) - ztrdt(:,:,jk,jn) ) * zfact 139 CALL trd_tra( kt, 'TRC', jn, jptra_trd_atf, ztrdt ) 173 140 END DO 174 IF (luttrd(jn)) CALL trd_mod_trc( ztrtrd, jn, jptrc_trd_atf, kt ) 175 ENDIF 176 ! ! =========== 177 END DO ! tracer loop 178 ! ! =========== 179 141 END DO 142 DEALLOCATE( ztrdt ) 143 END IF 144 ! 180 145 IF(ln_ctl) THEN ! print mean trends (used for debugging) 181 146 WRITE(charout, FMT="('nxt')") … … 183 148 CALL prt_ctl_trc(tab4d=trn, mask=tmask, clinfo=ctrcnm) 184 149 ENDIF 185 186 #if defined key_agrif 187 IF (.NOT.Agrif_Root()) CALL Agrif_Update_Trc( kt ) 188 #endif 189 190 150 ! 191 151 END SUBROUTINE trc_nxt 192 152 -
branches/DEV_r2006_merge_TRA_TRC/NEMO/TOP_SRC/TRP/trcrad.F90
r1257 r2030 14 14 !!---------------------------------------------------------------------- 15 15 USE oce_trc ! ocean dynamics and tracers variables 16 USE tr p_trc! ocean passive tracers variables17 USE trdm ld_trc18 USE trd mld_trc_oce16 USE trc ! ocean passive tracers variables 17 USE trdmod_oce 18 USE trdtra 19 19 USE lib_mpp 20 20 USE prtctl_trc ! Print control for debbuging … … 139 139 DO ji = 1, jpi 140 140 zvolk = cvol(ji,jj,jk) 141 # if defined key_ off_degrad141 # if defined key_degrad 142 142 zvolk = zvolk * facvol(ji,jj,jk) 143 143 # endif … … 180 180 ztrtrdb(:,:,:) = ( ptrb(:,:,:,jn) - ztrtrdb(:,:,:) ) * zs2rdt 181 181 ztrtrdn(:,:,:) = ( ptrn(:,:,:,jn) - ztrtrdn(:,:,:) ) * zs2rdt 182 IF (luttrd(jn)) CALL trd_mod_trc( ztrtrdb, jn, jptrc_trd_radb, kt) ! Asselin-like trend handling183 IF (luttrd(jn)) CALL trd_mod_trc( ztrtrdn, jn, jptrc_trd_radn, kt) ! standard trend handling182 CALL trd_tra( kt, 'TRC', jn, jptra_trd_radb, ztrtrdb ) ! Asselin-like trend handling 183 CALL trd_tra( kt, 'TRC', jn, jptra_trd_radn, ztrtrdn ) ! standard trend handling 184 184 ! 185 185 ENDIF … … 208 208 IF( l_trdtrc ) THEN 209 209 ! 210 zs2rdt = 1. / ( 2. * rdt * FLOAT(n dttrc) )210 zs2rdt = 1. / ( 2. * rdt * FLOAT(nn_dttrc) ) 211 211 ztrtrdb(:,:,:) = ( ptrb(:,:,:,jn) - ztrtrdb(:,:,:) ) * zs2rdt 212 212 ztrtrdn(:,:,:) = ( ptrn(:,:,:,jn) - ztrtrdn(:,:,:) ) * zs2rdt 213 IF (luttrd(jn)) CALL trd_mod_trc( ztrtrdb, jn, jptrc_trd_radb, kt) ! Asselin-like trend handling214 IF (luttrd(jn)) CALL trd_mod_trc( ztrtrdn, jn, jptrc_trd_radn, kt) ! standard trend handling213 CALL trd_tra( kt, 'TRC', jn, jptra_trd_radb, ztrtrdb ) ! Asselin-like trend handling 214 CALL trd_tra( kt, 'TRC', jn, jptra_trd_radn, ztrtrdn ) ! standard trend handling 215 215 ! 216 216 ENDIF -
branches/DEV_r2006_merge_TRA_TRC/NEMO/TOP_SRC/TRP/trcsbc.F90
r1739 r2030 4 4 !! Ocean passive tracers: surface boundary condition 5 5 !!====================================================================== 6 !! History : 8.2 ! 98-10 (G. Madec, G. Roullet, M. Imbard) Original code7 !! 8.2 ! 01-02 (D. Ludicone) sea ice and free surface8 !! 8.5 ! 02-06 (G. Madec) F90: Free form and module9 !! 9.0 ! 04-03 (C. Ethe) adapted for passive tracers10 !! ! 06-08 (C. Deltel) Diagnose ML trends for passive tracers6 !! History : 8.2 ! 1998-10 (G. Madec, G. Roullet, M. Imbard) Original code 7 !! 8.2 ! 2001-02 (D. Ludicone) sea ice and free surface 8 !! 8.5 ! 2002-06 (G. Madec) F90: Free form and module 9 !! 9.0 ! 2004-03 (C. Ethe) adapted for passive tracers 10 !! ! 2006-08 (C. Deltel) Diagnose ML trends for passive tracers 11 11 !!============================================================================== 12 12 #if defined key_top … … 18 18 !! * Modules used 19 19 USE oce_trc ! ocean dynamics and active tracers variables 20 USE tr p_trc ! ocean passive tracers variables20 USE trc ! ocean passive tracers variables 21 21 USE prtctl_trc ! Print control for debbuging 22 USE trdm ld_trc23 USE trd mld_trc_oce22 USE trdmod_oce 23 USE trdtra 24 24 25 25 IMPLICIT NONE … … 80 80 ! 0. initialization 81 81 zsrau = 1. / rau0 82 IF( .NOT. ln_sco ) zse3t = 1. / fse3t(1,1,1) 82 #if defined key_zco 83 zse3t = 1. / e3t_0(1) 84 #endif 83 85 84 86 DO jn = 1, jptra … … 88 90 DO jj = 2, jpj 89 91 DO ji = fs_2, fs_jpim1 ! vector opt. 90 IF( ln_sco ) zse3t = 1. / fse3t(ji,jj,1) 91 ! concent./dilut. effect 92 ztra = emps(ji,jj) * zsrau * trn(ji,jj,1,jn) * zse3t * tmask(ji,jj,1) 93 ! add the trend to the general tracer trend 92 #if ! defined key_zco 93 zse3t = 1. / fse3t(ji,jj,1) 94 #endif 95 IF( lk_vvl ) THEN ; ztra = 0.e0 ! No concent./dilut. effect 96 ELSE ; ztra = emps(ji,jj) * zsrau * trn(ji,jj,1,jn) * zse3t * tmask(ji,jj,1) 97 ENDIF 94 98 tra(ji,jj,1,jn) = tra(ji,jj,1,jn) + ztra 95 #if defined key_trc_diatrd96 IF( luttrd(jn) ) &97 & trtrd(ji,jj,1,ikeep(jn),jpdiatrc) = trtrd(ji,jj,1,ikeep(jn),jpdiatrc) + ztra98 #endif99 100 99 END DO 101 100 END DO … … 103 102 IF( l_trdtrc ) THEN 104 103 ztrtrd(:,:,:) = tra(:,:,:,jn) - ztrtrd(:,:,:) 105 IF (luttrd(jn)) CALL trd_mod_trc(ztrtrd, jn, jptrc_trd_sbc, kt)104 CALL trd_tra( kt, 'TRC', jn, jptra_trd_nsr, ztrtrd ) 106 105 END IF 107 108 106 ! ! =========== 109 107 END DO ! tracer loop … … 111 109 112 110 IF( l_trdtrc ) DEALLOCATE( ztrtrd ) 113 114 115 IF(ln_ctl) THEN ! print mean trends (used for debugging) 116 WRITE(charout, FMT="('sbc')") 117 CALL prt_ctl_trc_info(charout) 118 CALL prt_ctl_trc(tab4d=tra, mask=tmask, clinfo=ctrcnm,clinfo2='trd') 111 IF( ln_ctl ) THEN 112 WRITE(charout, FMT="('sbc ')") ; CALL prt_ctl_trc_info(charout) 113 CALL prt_ctl_trc( tab4d=tra, mask=tmask, clinfo=ctrcnm, clinfo2='trd' ) 119 114 ENDIF 120 115 -
branches/DEV_r2006_merge_TRA_TRC/NEMO/TOP_SRC/TRP/trctrp.F90
r1953 r2030 5 5 !!====================================================================== 6 6 !! History : 1.0 ! 2004-03 (C. Ethe) Original code 7 !! 3.3 ! 2010-07 (C. Ethe) Merge TRA-TRC 7 8 !!---------------------------------------------------------------------- 8 9 #if defined key_top … … 13 14 !!---------------------------------------------------------------------- 14 15 USE oce_trc ! ocean dynamics and active tracers variables 15 USE trp_trc ! ocean passive tracers variables 16 USE trctrp_lec ! passive tracers transport parameters 17 USE prtctl_trc ! Print control for debbuging 18 16 USE trc ! ocean passive tracers variables 17 USE trcnam_trp ! passive tracers transport namelist variables 18 USE trabbl ! bottom boundary layer (trc_bbl routine) 19 19 USE trcbbl ! bottom boundary layer (trc_bbl routine) 20 #if ! defined key_offline 21 USE zdfkpp ! KPP non-local tracer fluxes (trc_kpp routine) 22 #endif 20 23 USE trcdmp ! internal damping (trc_dmp routine) 21 22 USE trcldf_bilapg ! lateral mixing (trc_ldf_bilapg routine) 23 USE trcldf_bilap ! lateral mixing (trc_ldf_bilap routine) 24 USE trcldf_iso ! lateral mixing (trc_ldf_iso routine) 25 USE trcldf_iso_zps ! lateral mixing (trc_ldf_iso_zps routine) 26 USE trcldf_lap ! lateral mixing (trc_ldf_lap routine) 27 24 USE trcldf ! lateral mixing (trc_ldf routine) 25 USE trcadv ! advection (trc_adv routine) 26 USE trczdf ! vertical diffusion (trc_zdf routine) 28 27 USE trcnxt ! time-stepping (trc_nxt routine) 29 28 USE trcrad ! positivity (trc_rad routine) 30 31 USE trcadv_cen2 ! 2nd order centered advection (trc_adv_cen2 routine)32 USE trcadv_muscl ! MUSCL advection (trc_adv_muscl routine)33 USE trcadv_muscl2 ! MUSCL2 advection (trc_adv_muscl2 routine)34 USE trcadv_tvd ! TVD advection (trc_adv_tvd routine)35 USE trcadv_smolar ! SMOLAR advection (trc_adv_smolar routine)36 37 USE trczdf_exp ! vertical diffusion (trc_zdf_exp routine)38 USE trczdf_imp ! vertical diffusion (trc_zdf_exp routine)39 USE trczdf_iso ! vertical diffusion (trc_zdf_exp routine)40 USE trczdf_iso_vopt ! vertical diffusion (trc_zdf_exp routine)41 29 USE trcsbc ! surface boundary condition (trc_sbc routine) 42 43 USE zpshde_trc ! partial step: hor. derivative (zps_hde_trc routine) 30 USE zpshde ! partial step: hor. derivative (zps_hde_trc routine) 44 31 45 32 #if defined key_agrif … … 62 49 CONTAINS 63 50 64 SUBROUTINE trc_trp( k t)51 SUBROUTINE trc_trp( kstp ) 65 52 !!---------------------------------------------------------------------- 66 53 !! *** ROUTINE trc_trp *** … … 71 58 !! - Update the passive tracers 72 59 !!---------------------------------------------------------------------- 73 INTEGER, INTENT( in ) :: kt ! ocean time-step index 74 !! 75 CHARACTER (len=25) :: charout 60 INTEGER, INTENT( in ) :: kstp ! ocean time-step index 76 61 !! --------------------------------------------------------------------- 77 78 CALL trc_sbc( kt ) ! surface boundary condition 79 # if defined key_trcbbc 80 !!gm bug : this should be control during the initialisation phase, not here! 81 CALL ctl_stop( ' Bottom heat flux not yet implemented with passive tracer ' & 82 & ' Check in trc_trp routine ' ) 83 # endif 84 ! ! bottom boundary condition 85 IF( lk_trcbbl_dif ) CALL trc_bbl_dif( kt ) ! diffusive bottom boundary layer scheme 86 IF( lk_trcbbl_adv ) CALL trc_bbl_adv( kt ) ! advective (and/or diffusive) bottom boundary layer scheme 87 88 IF( lk_trcdmp ) CALL trc_dmp( kt ) ! internal damping trends 89 90 ! ! horizontal & vertical advection 91 IF( ln_trcadv_cen2 ) CALL trc_adv_cen2 ( kt ) ! 2nd order centered scheme 92 IF( ln_trcadv_muscl ) CALL trc_adv_muscl ( kt ) ! MUSCL scheme 93 IF( ln_trcadv_muscl2 ) CALL trc_adv_muscl2( kt ) ! MUSCL2 scheme 94 IF( ln_trcadv_tvd ) CALL trc_adv_tvd ( kt ) ! TVD scheme 95 IF( ln_trcadv_smolar ) CALL trc_adv_smolar( kt ) ! SMOLARKIEWICZ scheme 96 97 98 IF( n_cla == 1 ) THEN 99 !!gm bug : this should be control during the initialisation phase, not here! 100 WRITE(ctmp1,*) ' Cross Land Advection not yet implemented with passive tracer n_cla = ',n_cla 101 CALL ctl_stop( ctmp1 ) 102 ENDIF 103 104 ! ! lateral mixing 105 IF( l_trcldf_bilapg ) CALL trc_ldf_bilapg ( kt ) ! s-coord. horizontal bilaplacian 106 IF( l_trcldf_bilap ) CALL trc_ldf_bilap ( kt ) ! iso-level bilaplacian 107 IF( l_trcldf_iso ) CALL trc_ldf_iso ( kt ) ! iso-neutral laplacian 108 IF( l_trcldf_iso_zps ) CALL trc_ldf_iso_zps( kt ) ! partial step iso-neutral laplacian 109 IF( l_trcldf_lap ) CALL trc_ldf_lap ( kt ) ! iso-level laplacian 110 62 IF( .NOT. lk_trc_c1d ) THEN 63 ! 64 CALL trc_sbc( kstp ) ! surface boundary condition 65 IF( lk_trabbl ) CALL trc_bbl( kstp ) ! advective (and/or diffusive) bottom boundary layer scheme 66 IF( lk_trcdmp ) CALL trc_dmp( kstp ) ! internal damping trends 67 CALL trc_adv( kstp ) ! horizontal & vertical advection 68 CALL trc_ldf( kstp ) ! lateral mixing 69 #if ! defined key_offline 70 IF( lk_zdfkpp ) CALL trc_kpp( kstp ) ! KPP non-local tracer fluxes 71 #endif 111 72 #if defined key_agrif 112 IF(.NOT. Agrif_Root()) CALL Agrif_Sponge_trc! tracers sponge73 IF(.NOT. Agrif_Root()) CALL Agrif_Sponge_trc ! tracers sponge 113 74 #endif 114 115 ! ! vertical diffusion116 IF( l_trczdf_exp ) CALL trc_zdf_exp ( kt ) ! explicit time stepping (time splitting scheme)117 IF( l_trczdf_imp ) CALL trc_zdf_imp ( kt ) ! implicit time stepping (euler backward)118 IF( l_trczdf_iso ) CALL trc_zdf_iso ( kt ) ! isopycnal119 IF( l_trczdf_iso_vo ) CALL trc_zdf_iso_vopt( kt ) ! vector opt. isopycnal120 121 CALL trc_nxt( kt ) ! tracer fields at next time step 122 123 IF( ln_trcrad ) CALL trc_rad( kt ) ! Correct artificial negative concentrations 124 ! ! especially useful when isopycnal mixing is used125 !126 127 IF( ln_zps .AND. .NOT. lk_trc_c1d ) & ! Partial steps: now horizontal gradient of passive128 & CALL zps_hde_trc( kt, trb, gtru, gtrv ) ! tracers at the bottom ocean level75 CALL trc_zdf( kstp ) ! vertical mixing and after tracer fields 76 CALL trc_nxt( kstp ) ! tracer fields at next time step 77 IF( ln_trcrad ) CALL trc_rad( kstp ) ! Correct artificial negative concentrations 78 IF( ln_zps ) CALL zps_hde_trc( kstp, jptra, trb, gtru, gtrv ) ! Partial steps: now horizontal gradient of passive 79 ! tracers at the bottom ocean level 80 ELSE 81 CALL trc_sbc( kstp ) ! surface boundary condition 82 #if ! defined key_offline 83 IF( lk_zdfkpp ) CALL trc_kpp( kstp ) ! KPP non-local tracer fluxes 84 #endif 85 CALL trc_zdf( kstp ) ! vertical mixing and after tracer fields 86 CALL trc_nxt( kstp ) ! tracer fields at next time step 87 IF( ln_trcrad ) CALL trc_rad( kstp ) ! Correct artificial negative concentrations 88 ! 89 END IF 129 90 ! 130 91 END SUBROUTINE trc_trp -
branches/DEV_r2006_merge_TRA_TRC/NEMO/TOP_SRC/TRP/trdmld_trc.F90
r1685 r2030 16 16 !! trd_mld_trc_init : initialization step 17 17 !!---------------------------------------------------------------------- 18 USE trp_trc ! tracer definitions (trn, trb, tra, etc.) 19 USE oce_trc ! needed for namelist logicals, and euphotic layer arrays 20 USE trctrp_lec 21 USE trdmld_trc_oce ! definition of main arrays used for trends computations 18 USE trc ! tracer definitions (trn, trb, tra, etc.) 19 USE dom_oce ! domain definition 20 USE zdfmxl , ONLY : nmln !: number of level in the mixed layer 21 USE zdf_oce , ONLY : avt !: vert. diffusivity coef. at w-point for temp 22 # if defined key_zdfddm 23 USE zdfddm , ONLY : avs !: salinity vertical diffusivity coeff. at w-point 24 # endif 25 USE trcnam_trp ! passive tracers transport namelist variables 26 USE trdmod_trc_oce ! definition of main arrays used for trends computations 22 27 USE in_out_manager ! I/O manager 23 28 USE dianam ! build the name of file (routine) … … 29 34 USE sms_pisces 30 35 USE sms_lobster 31 USE trc32 36 33 37 IMPLICIT NONE 34 38 PRIVATE 35 39 36 INTERFACE trd_mod_trc37 MODULE PROCEDURE trd_mod_trc_trp, trd_mod_trc_bio38 END INTERFACE39 40 PUBLIC trd_mod_trc ! routine called by step.F9041 40 PUBLIC trd_mld_trc 42 41 PUBLIC trd_mld_bio 43 42 PUBLIC trd_mld_trc_init 43 PUBLIC trd_mld_trc_zint 44 PUBLIC trd_mld_bio_zint 44 45 45 46 CHARACTER (LEN=40) :: clhstnam ! name of the trends NetCDF file … … 66 67 67 68 CONTAINS 68 69 SUBROUTINE trd_mod_trc_trp( ptrtrd, kjn, ktrd, kt )70 !!----------------------------------------------------------------------71 !! *** ROUTINE trd_mod_trc ***72 !!----------------------------------------------------------------------73 #if defined key_trcbbl_adv74 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zun, zvn ! temporary arrays75 #else76 USE oce_trc, zun => un ! When no bbl, zun == un77 USE oce_trc, zvn => vn ! When no bbl, zvn == vn78 #endif79 INTEGER, INTENT( in ) :: kt ! time step80 INTEGER, INTENT( in ) :: kjn ! tracer index81 INTEGER, INTENT( in ) :: ktrd ! tracer trend index82 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( inout ) :: ptrtrd ! Temperature or U trend83 !!----------------------------------------------------------------------84 85 IF( kt == nittrc000 ) THEN86 ! IF(lwp)WRITE(numout,*)87 ! IF(lwp)WRITE(numout,*) 'trd_mod_trc:'88 ! IF(lwp)WRITE(numout,*) '~~~~~~~~~~~~'89 ENDIF90 91 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>92 ! Mixed layer trends for passive tracers93 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<94 95 SELECT CASE ( ktrd )96 CASE ( jptrc_trd_xad ) ; CALL trd_mld_trc_zint( ptrtrd, jpmld_trc_xad , '3D', kjn )97 CASE ( jptrc_trd_yad ) ; CALL trd_mld_trc_zint( ptrtrd, jpmld_trc_yad , '3D', kjn )98 CASE ( jptrc_trd_zad ) ; CALL trd_mld_trc_zint( ptrtrd, jpmld_trc_zad , '3D', kjn )99 CASE ( jptrc_trd_ldf ) ; CALL trd_mld_trc_zint( ptrtrd, jpmld_trc_ldf , '3D', kjn )100 CASE ( jptrc_trd_xei ) ; CALL trd_mld_trc_zint( ptrtrd, jpmld_trc_xei , '3D', kjn )101 CASE ( jptrc_trd_yei ) ; CALL trd_mld_trc_zint( ptrtrd, jpmld_trc_yei , '3D', kjn )102 CASE ( jptrc_trd_bbl ) ; CALL trd_mld_trc_zint( ptrtrd, jpmld_trc_bbl , '3D', kjn )103 CASE ( jptrc_trd_zdf )104 IF( ln_trcldf_iso ) THEN105 CALL trd_mld_trc_zint( ptrtrd, jpmld_trc_ldf, '3D', kjn )106 ELSE107 CALL trd_mld_trc_zint( ptrtrd, jpmld_trc_zdf, '3D', kjn )108 ENDIF109 CASE ( jptrc_trd_zei ) ; CALL trd_mld_trc_zint( ptrtrd, jpmld_trc_zei , '3D', kjn )110 CASE ( jptrc_trd_dmp ) ; CALL trd_mld_trc_zint( ptrtrd, jpmld_trc_dmp , '3D', kjn )111 CASE ( jptrc_trd_sbc ) ; CALL trd_mld_trc_zint( ptrtrd, jpmld_trc_sbc , '2D', kjn )112 CASE ( jptrc_trd_sms ) ; CALL trd_mld_trc_zint( ptrtrd, jpmld_trc_sms , '3D', kjn )113 CASE ( jptrc_trd_bbc ) ; CALL trd_mld_trc_zint( ptrtrd, jpmld_trc_bbc , '3D', kjn )114 CASE ( jptrc_trd_radb ) ; CALL trd_mld_trc_zint( ptrtrd, jpmld_trc_radb , '3D', kjn )115 CASE ( jptrc_trd_radn ) ; CALL trd_mld_trc_zint( ptrtrd, jpmld_trc_radn , '3D', kjn )116 CASE ( jptrc_trd_atf ) ; CALL trd_mld_trc_zint( ptrtrd, jpmld_trc_atf , '3D', kjn )117 END SELECT118 119 120 END SUBROUTINE trd_mod_trc_trp121 122 SUBROUTINE trd_mod_trc_bio( ptrbio, ktrd, kt )123 !!----------------------------------------------------------------------124 !! *** ROUTINE trd_mod_bio ***125 !!----------------------------------------------------------------------126 127 INTEGER, INTENT( in ) :: kt ! time step128 INTEGER, INTENT( in ) :: ktrd ! bio trend index129 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( inout ) :: ptrbio ! Bio trend130 !!----------------------------------------------------------------------131 132 CALL trd_mld_bio_zint( ptrbio, ktrd ) ! Verticaly integrated biological trends133 134 END SUBROUTINE trd_mod_trc_bio135 136 69 137 70 SUBROUTINE trd_mld_trc_zint( ptrc_trdmld, ktrd, ctype, kjn ) … … 170 103 171 104 ! ... Set nmld(ji,jj) = index of first T point below control surf. or outside mixed-layer 172 SELECT CASE ( n ctls_trc ) ! choice of the control surface105 SELECT CASE ( nn_ctls_trc ) ! choice of the control surface 173 106 CASE ( -2 ) ; STOP 'trdmld_trc : not ready ' ! -> isopycnal surface (see ???) 174 107 #if defined key_pisces || defined key_lobster … … 177 110 CASE ( 0 ) ; nmld_trc(:,:) = nmln(:,:) ! -> ML with density criterion (see zdfmxl) 178 111 CASE ( 1 ) ; nmld_trc(:,:) = nbol_trc(:,:) ! -> read index from file 179 CASE ( 2: ) ; n ctls_trc = MIN( nctls_trc, jpktrd_trc - 1 )180 nmld_trc(:,:) = n ctls_trc + 1 ! -> model level112 CASE ( 2: ) ; nn_ctls_trc = MIN( nn_ctls_trc, jpktrd_trc - 1 ) 113 nmld_trc(:,:) = nn_ctls_trc + 1 ! -> model level 181 114 END SELECT 182 115 … … 281 214 tmltrd_bio(:,:,:) = 0.e0 ! <<< reset trend arrays to zero 282 215 ! ... Set nmld(ji,jj) = index of first T point below control surf. or outside mixed-layer 283 SELECT CASE ( n ctls_trc ) ! choice of the control surface216 SELECT CASE ( nn_ctls_trc ) ! choice of the control surface 284 217 CASE ( -2 ) ; STOP 'trdmld_trc : not ready ' ! -> isopycnal surface (see ???) 285 218 CASE ( -1 ) ; nmld_trc(:,:) = neln(:,:) ! -> euphotic layer with light criterion 286 219 CASE ( 0 ) ; nmld_trc(:,:) = nmln(:,:) ! -> ML with density criterion (see zdfmxl) 287 220 CASE ( 1 ) ; nmld_trc(:,:) = nbol_trc(:,:) ! -> read index from file 288 CASE ( 2: ) ; n ctls_trc = MIN( nctls_trc, jpktrd_trc - 1 )289 nmld_trc(:,:) = n ctls_trc + 1 ! -> model level221 CASE ( 2: ) ; nn_ctls_trc = MIN( nn_ctls_trc, jpktrd_trc - 1 ) 222 nmld_trc(:,:) = nn_ctls_trc + 1 ! -> model level 290 223 END SELECT 291 224 … … 380 313 !! of time-averaged mixed-layer fields, this residual WILL NOT BE ZERO 381 314 !! over the first two analysis windows (except if restart). 382 !! N.B. For ORCA2_LIM, use e.g. ntrc_trc=5, ucf_trc=1., nctls_trc=8315 !! N.B. For ORCA2_LIM, use e.g. ntrc_trc=5, rn_ucf_trc=1., nctls_trc=8 383 316 !! for checking residuals. 384 317 !! On a NEC-SX5 computer, this typically leads to: … … 421 354 REAL(wp), DIMENSION(jpi,jpj,jpltrd_trc,jptra) :: ztmltrd2 ! -+ 422 355 !! 423 REAL(wp), DIMENSION(jpi,jpj) :: z2d ! temporary array, used for eiv arrays424 356 CHARACTER (LEN= 5) :: clvar 425 357 #if defined key_dimgout … … 431 363 IF( llwarn ) THEN ! warnings 432 364 IF( ( nittrc000 /= nit000 ) & 433 .OR.( n dttrc /= 1 ) ) THEN365 .OR.( nn_dttrc /= 1 ) ) THEN 434 366 435 367 WRITE(numout,*) 'Be careful, trends diags never validated' … … 450 382 DO ji = 1,jpi 451 383 ik = nmld_trc(ji,jj) 452 zavt = avt(ji,jj,ik)384 zavt = fstravs(ji,jj,ik) 453 385 DO jn = 1, jptra 454 IF( l uttrd(jn) ) &386 IF( ln_trdtrc(jn) ) & 455 387 tmltrd_trc(ji,jj,jpmld_trc_zdf,jn) = - zavt / fse3w(ji,jj,ik) * tmask(ji,jj,ik) & 456 388 & * ( trn(ji,jj,ik-1,jn) - trn(ji,jj,ik,jn) ) & … … 462 394 DO jn = 1, jptra 463 395 ! ... Remove this K_z trend from the iso-neutral diffusion term (if any) 464 IF( l uttrd(jn) ) &396 IF( ln_trdtrc(jn) ) & 465 397 tmltrd_trc(:,:,jpmld_trc_ldf,jn) = tmltrd_trc(:,:,jpmld_trc_ldf,jn) - tmltrd_trc(:,:,jpmld_trc_zdf,jn) 466 398 … … 473 405 ! therefore we do not call lbc_lnk in GYRE config. (closed basin, no cyclic B.C.) 474 406 DO jn = 1, jptra 475 IF( l uttrd(jn) ) THEN407 IF( ln_trdtrc(jn) ) THEN 476 408 DO jl = 1, jpltrd_trc 477 409 CALL lbc_lnk( tmltrd_trc(:,:,jl,jn), 'T', 1. ) ! lateral boundary conditions … … 492 424 IF( kt > nittrc000 ) THEN 493 425 DO jn = 1, jptra 494 IF( l uttrd(jn) ) THEN426 IF( ln_trdtrc(jn) ) THEN 495 427 tmlb_trc (:,:,jn) = tml_trc (:,:,jn) 496 428 tmlatfn_trc(:,:,jn) = tmltrd_trc(:,:,jpmld_trc_atf,jn) … … 505 437 DO jk = 1, jpktrd_trc ! - 1 ??? 506 438 DO jn = 1, jptra 507 IF( l uttrd(jn) ) &439 IF( ln_trdtrc(jn) ) & 508 440 tml_trc(:,:,jn) = tml_trc(:,:,jn) + wkx_trc(:,:,jk) * trn(:,:,jk,jn) 509 441 END DO … … 515 447 ! 516 448 DO jn = 1, jptra 517 IF( l uttrd(jn) ) THEN449 IF( ln_trdtrc(jn) ) THEN 518 450 tmlbb_trc (:,:,jn) = tmlb_trc (:,:,jn) ; tmlbn_trc (:,:,jn) = tml_trc (:,:,jn) 519 451 tmlatfb_trc(:,:,jn) = tmlatfn_trc(:,:,jn) ; tmlradb_trc(:,:,jn) = tmlradn_trc(:,:,jn) … … 544 476 ! ... Cumulate over BOTH physical contributions AND over time steps 545 477 DO jn = 1, jptra 546 IF( l uttrd(jn) ) THEN478 IF( ln_trdtrc(jn) ) THEN 547 479 DO jl = 1, jpltrd_trc 548 480 tmltrdm_trc(:,:,jn) = tmltrdm_trc(:,:,jn) + tmltrd_trc(:,:,jl,jn) … … 552 484 553 485 DO jn = 1, jptra 554 IF( l uttrd(jn) ) THEN486 IF( ln_trdtrc(jn) ) THEN 555 487 ! ... Special handling of the Asselin trend 556 488 tmlatfm_trc(:,:,jn) = tmlatfm_trc(:,:,jn) + tmlatfn_trc(:,:,jn) … … 573 505 574 506 ! Convert to appropriate physical units 575 tmltrd_trc(:,:,:,:) = tmltrd_trc(:,:,:,:) * ucf_trc507 tmltrd_trc(:,:,:,:) = tmltrd_trc(:,:,:,:) * rn_ucf_trc 576 508 577 509 itmod = kt - nittrc000 + 1 578 510 it = kt 579 511 580 MODULO_NTRD : IF( MOD( itmod, n trd_trc ) == 0 ) THEN ! nitend MUST be multiple of ntrd_trc512 MODULO_NTRD : IF( MOD( itmod, nn_trd_trc ) == 0 ) THEN ! nitend MUST be multiple of nn_trd_trc 581 513 ! 582 514 ztmltot (:,:,:) = 0.e0 ! reset arrays to zero … … 591 523 592 524 DO jn = 1, jptra 593 IF( l uttrd(jn) ) THEN525 IF( ln_trdtrc(jn) ) THEN 594 526 !-- Compute total trends (use rdttrc instead of rdt ???) 595 IF ( ln_trcadv_ smolar .OR. ln_trcadv_muscl .OR. ln_trcadv_muscl2 ) THEN ! EULER-FORWARD schemes527 IF ( ln_trcadv_muscl .OR. ln_trcadv_muscl2 ) THEN ! EULER-FORWARD schemes 596 528 ztmltot(:,:,jn) = ( tml_trc(:,:,jn) - tmlbn_trc(:,:,jn) )/rdt 597 529 ELSE ! LEAP-FROG schemes … … 629 561 !-- Compute passive tracer total trends 630 562 DO jn = 1, jptra 631 IF( l uttrd(jn) ) THEN563 IF( ln_trdtrc(jn) ) THEN 632 564 tml_sum_trc(:,:,jn) = tmlbn_trc(:,:,jn) + 2 * ( tml_sum_trc(:,:,jn) - tml_trc(:,:,jn) ) + tml_trc(:,:,jn) 633 565 ztmltot2 (:,:,jn) = ( tml_sum_trc(:,:,jn) - tml_sumb_trc(:,:,jn) ) / ( 2.*rdt ) ! now tracer unit is /sec … … 637 569 !-- Compute passive tracer residuals 638 570 DO jn = 1, jptra 639 IF( l uttrd(jn) ) THEN571 IF( ln_trdtrc(jn) ) THEN 640 572 ! 641 573 DO jl = 1, jpltrd_trc … … 680 612 DO jn = 1, jptra 681 613 682 IF( l uttrd(jn) ) THEN614 IF( ln_trdtrc(jn) ) THEN 683 615 WRITE(numout, *) 684 616 WRITE(numout, *) '>>>>>>>>>>>>>>>>>> TRC TRACER jn =', jn, ' <<<<<<<<<<<<<<<<<<' … … 777 709 rmld_sum_trc(:,:) = rmld_sum_trc(:,:) / (2*zfn) ! similar to tml_sum and sml_sum 778 710 DO jn = 1, jptra 779 IF( l uttrd(jn) ) THEN711 IF( ln_trdtrc(jn) ) THEN 780 712 ! For passive tracer instantaneous diagnostics 781 713 tmlbb_trc (:,:,jn) = tmlb_trc (:,:,jn) ; tmlbn_trc (:,:,jn) = tml_trc (:,:,jn) … … 791 723 ! III.4 Convert to appropriate physical units 792 724 ! ------------------------------------------- 793 ztmltot (:,:,jn) = ztmltot (:,:,jn) * ucf_trc/zfn ! instant diags794 ztmlres (:,:,jn) = ztmlres (:,:,jn) * ucf_trc/zfn795 ztmlatf (:,:,jn) = ztmlatf (:,:,jn) * ucf_trc/zfn796 ztmlrad (:,:,jn) = ztmlrad (:,:,jn) * ucf_trc/zfn725 ztmltot (:,:,jn) = ztmltot (:,:,jn) * rn_ucf_trc/zfn ! instant diags 726 ztmlres (:,:,jn) = ztmlres (:,:,jn) * rn_ucf_trc/zfn 727 ztmlatf (:,:,jn) = ztmlatf (:,:,jn) * rn_ucf_trc/zfn 728 ztmlrad (:,:,jn) = ztmlrad (:,:,jn) * rn_ucf_trc/zfn 797 729 tml_sum_trc (:,:,jn) = tml_sum_trc (:,:,jn) / (2*zfn) ! mean diags 798 ztmltot2 (:,:,jn) = ztmltot2 (:,:,jn) * ucf_trc/zfn2799 ztmltrd2 (:,:,:,jn) = ztmltrd2 (:,:,:,jn) * ucf_trc/zfn2800 ztmlatf2 (:,:,jn) = ztmlatf2 (:,:,jn) * ucf_trc/zfn2801 ztmlrad2 (:,:,jn) = ztmlrad2 (:,:,jn) * ucf_trc/zfn2802 ztmlres2 (:,:,jn) = ztmlres2 (:,:,jn) * ucf_trc/zfn2730 ztmltot2 (:,:,jn) = ztmltot2 (:,:,jn) * rn_ucf_trc/zfn2 731 ztmltrd2 (:,:,:,jn) = ztmltrd2 (:,:,:,jn) * rn_ucf_trc/zfn2 732 ztmlatf2 (:,:,jn) = ztmlatf2 (:,:,jn) * rn_ucf_trc/zfn2 733 ztmlrad2 (:,:,jn) = ztmlrad2 (:,:,jn) * rn_ucf_trc/zfn2 734 ztmlres2 (:,:,jn) = ztmlres2 (:,:,jn) * rn_ucf_trc/zfn2 803 735 ENDIF 804 736 END DO … … 820 752 ! ---------------------------------- 821 753 822 IF( lwp .AND. MOD( itmod , n trd_trc ) == 0 ) THEN754 IF( lwp .AND. MOD( itmod , nn_trd_trc ) == 0 ) THEN 823 755 WRITE(numout,*) ' ' 824 756 WRITE(numout,*) 'trd_mld_trc : write passive tracer trends in the NetCDF file :' … … 834 766 DO jn = 1, jptra 835 767 ! 836 IF( luttrd(jn) ) THEN 837 !-- Specific treatment for EIV trends 838 ! WARNING : When eiv is switched on but key_diaeiv is not, we do NOT diagnose 839 ! u_eiv, v_eiv, and w_eiv : the exact eiv advective trends thus cannot be computed, 840 ! only their sum makes sense => mask directional contrib. to avoid confusion 841 z2d(:,:) = tmltrd_trc(:,:,jpmld_trc_xei,jn) + tmltrd_trc(:,:,jpmld_trc_yei,jn) & 842 & + tmltrd_trc(:,:,jpmld_trc_zei,jn) 843 #if ( defined key_trcldf_eiv && defined key_diaeiv ) 844 tmltrd_trc(:,:,jpmld_trc_xei,jn) = -999. 845 tmltrd_trc(:,:,jpmld_trc_yei,jn) = -999. 846 tmltrd_trc(:,:,jpmld_trc_zei,jn) = -999. 847 #endif 768 IF( ln_trdtrc(jn) ) THEN 848 769 CALL histwrite( nidtrd(jn), "mxl_depth", it, rmld_trc(:,:), ndimtrd1, ndextrd1 ) 849 770 !-- Output the fields … … 864 785 & it, ztmlatf(:,:,jn), ndimtrd1, ndextrd1 ) 865 786 866 CALL histwrite( nidtrd(jn), trim(clvar//ctrd_trc( jpltrd_trc+1,2)), & ! now total EIV : jpltrd_trc + 1867 & it, z2d(:,:), ndimtrd1, ndextrd1 )868 !869 787 ENDIF 870 788 END DO … … 872 790 IF( kt == nitend ) THEN 873 791 DO jn = 1, jptra 874 IF( l uttrd(jn) ) CALL histclo( nidtrd(jn) )792 IF( ln_trdtrc(jn) ) CALL histclo( nidtrd(jn) ) 875 793 END DO 876 794 ENDIF … … 881 799 DO jn = 1, jptra 882 800 ! 883 IF( luttrd(jn) ) THEN 884 !-- Specific treatment for EIV trends 885 ! WARNING : see above 886 z2d(:,:) = ztmltrd2(:,:,jpmld_trc_xei,jn) + ztmltrd2(:,:,jpmld_trc_yei,jn) & 887 & + ztmltrd2(:,:,jpmld_trc_zei,jn) 888 889 #if ( defined key_trcldf_eiv && defined key_diaeiv ) 890 ztmltrd2(:,:,jpmld_trc_xei,jn) = -999. 891 ztmltrd2(:,:,jpmld_trc_yei,jn) = -999. 892 ztmltrd2(:,:,jpmld_trc_zei,jn) = -999. 893 #endif 801 IF( ln_trdtrc(jn) ) THEN 894 802 CALL histwrite( nidtrd(jn), "mxl_depth", it, rmld_sum_trc(:,:), ndimtrd1, ndextrd1 ) 895 803 !-- Output the fields … … 911 819 & it, ztmlatf2(:,:,jn), ndimtrd1, ndextrd1 ) 912 820 913 CALL histwrite( nidtrd(jn), trim(clvar//ctrd_trc( jpltrd_trc+1,2)), & ! now total EIV : jpltrd_trc + 1914 & it, z2d(:,:), ndimtrd1, ndextrd1 )915 916 821 ENDIF 917 822 ! … … 919 824 IF( kt == nitend ) THEN 920 825 DO jn = 1, jptra 921 IF( l uttrd(jn) ) CALL histclo( nidtrd(jn) )826 IF( ln_trdtrc(jn) ) CALL histclo( nidtrd(jn) ) 922 827 END DO 923 828 ENDIF … … 931 836 # endif /* key_dimgout */ 932 837 933 IF( MOD( itmod, n trd_trc ) == 0 ) THEN838 IF( MOD( itmod, nn_trd_trc ) == 0 ) THEN 934 839 ! 935 840 ! Reset cumulative arrays to zero … … 1012 917 IF( llwarn ) THEN 1013 918 IF( ( nittrc000 /= nit000 ) & 1014 .OR.( n dttrc /= 1 ) ) THEN919 .OR.( nn_dttrc /= 1 ) ) THEN 1015 920 1016 921 WRITE(numout,*) 'Be careful, trends diags never validated' … … 1058 963 1059 964 ! Convert to appropriate physical units 1060 tmltrd_bio(:,:,:) = tmltrd_bio(:,:,:) * ucf_trc1061 1062 MODULO_NTRD : IF( MOD( kt, n trd_trc ) == 0 ) THEN ! nitend MUST be multiple of ntrd965 tmltrd_bio(:,:,:) = tmltrd_bio(:,:,:) * rn_ucf_trc 966 967 MODULO_NTRD : IF( MOD( kt, nn_trd_trc ) == 0 ) THEN ! nitend MUST be multiple of ntrd 1063 968 ! 1064 969 zfn = float(nmoymltrdbio) ; zfn2 = zfn * zfn … … 1114 1019 ! III.4 Convert to appropriate physical units 1115 1020 ! ------------------------------------------- 1116 ztmltrdbio2 (:,:,:) = ztmltrdbio2 (:,:,:) * ucf_trc/zfn21021 ztmltrdbio2 (:,:,:) = ztmltrdbio2 (:,:,:) * rn_ucf_trc/zfn2 1117 1022 1118 1023 END IF MODULO_NTRD … … 1136 1041 it = kt 1137 1042 1138 IF( lwp .AND. MOD( itmod , n trd_trc ) == 0 ) THEN1043 IF( lwp .AND. MOD( itmod , nn_trd_trc ) == 0 ) THEN 1139 1044 WRITE(numout,*) ' ' 1140 1045 WRITE(numout,*) 'trd_mld_bio : write ML bio trends in the NetCDF file :' … … 1176 1081 # endif /* key_dimgout */ 1177 1082 1178 IF( MOD( itmod, n trd_trc ) == 0 ) THEN1083 IF( MOD( itmod, nn_trd_trc ) == 0 ) THEN 1179 1084 ! 1180 1085 ! III.5 Reset cumulative arrays to zero … … 1216 1121 INTEGER :: ilseq, jl, jn 1217 1122 REAL(wp) :: zjulian, zsto, zout 1218 CHARACTER (LEN=40) :: clop , cleiv1123 CHARACTER (LEN=40) :: clop 1219 1124 CHARACTER (LEN=15) :: csuff 1220 1125 CHARACTER (LEN=12) :: clmxl 1221 1126 CHARACTER (LEN=16) :: cltrcu 1222 1127 CHARACTER (LEN= 5) :: clvar 1223 1224 NAMELIST/namtoptrd/ ntrd_trc, nctls_trc, ucf_trc, &1225 ln_trdmld_trc_restart, ln_trdmld_trc_instant, luttrd1226 1128 1227 1129 !!---------------------------------------------------------------------- … … 1241 1143 ! I.1 Check consistency of user defined preferences 1242 1144 ! ------------------------------------------------- 1243 #if defined key_trcldf_eiv 1244 IF( lk_trdmld_trc .AND. ln_trcldf_iso ) THEN 1245 WRITE(numout,cform_war) 1246 WRITE(numout,*) ' You asked for ML diagnostics with iso-neutral diffusion ' 1247 WRITE(numout,*) ' and eiv physics. ' 1248 WRITE(numout,*) ' Yet, key_diaeiv is NOT switched on, so the eddy induced ' 1249 WRITE(numout,*) ' velocity is not diagnosed. ' 1250 WRITE(numout,*) ' Therefore, we cannot deduce the eiv advective trends. ' 1251 WRITE(numout,*) ' Only THE SUM of the i,j,k directional contributions then ' 1252 WRITE(numout,*) ' makes sense => To avoid any confusion, we choosed to mask ' 1253 WRITE(numout,*) ' these i,j,k directional contributions (with -999.) ' 1254 nwarn = nwarn + 1 1255 ENDIF 1256 # endif 1257 1258 IF( ( lk_trdmld_trc ) .AND. ( MOD( nitend, ntrd_trc ) /= 0 ) ) THEN 1145 1146 IF( ( lk_trdmld_trc ) .AND. ( MOD( nitend, nn_trd_trc ) /= 0 ) ) THEN 1259 1147 WRITE(numout,cform_err) 1260 1148 WRITE(numout,*) ' Your nitend parameter, nitend = ', nitend 1261 1149 WRITE(numout,*) ' is no multiple of the trends diagnostics frequency ' 1262 WRITE(numout,*) ' you defined, n trd_trc = ', ntrd_trc1150 WRITE(numout,*) ' you defined, nn_trd_trc = ', nn_trd_trc 1263 1151 WRITE(numout,*) ' This will not allow you to restart from this simulation. ' 1264 1152 WRITE(numout,*) ' You should reconsider this choice. ' … … 1280 1168 IF( lldebug ) THEN 1281 1169 WRITE(numout,*) ' ln_trcadv_muscl = ' , ln_trcadv_muscl 1282 WRITE(numout,*) ' ln_trcadv_smolar = ' , ln_trcadv_smolar1283 1170 WRITE(numout,*) ' ln_trdmld_trc_instant = ', ln_trdmld_trc_instant 1284 ENDIF1285 1286 IF( ln_trcadv_smolar .AND. .NOT. ln_trdmld_trc_instant ) THEN1287 WRITE(numout,cform_err)1288 WRITE(numout,*) ' Currently, you can NOT use simultaneously tracer Smolark. '1289 WRITE(numout,*) ' advection and window averaged diagnostics of ML trends. '1290 WRITE(numout,*) ' WHY? Everything in trdmld_trc is coded for leap-frog, and '1291 WRITE(numout,*) ' Smolarkiewicz scheme is Euler forward. '1292 WRITE(numout,*) ' In particuliar, entrainment trend would be FALSE. However '1293 WRITE(numout,*) ' this residual is correct for instantaneous ML diagnostics.'1294 WRITE(numout,*)1295 nstop = nstop + 11296 1171 ENDIF 1297 1172 … … 1364 1239 ! I.3 Read control surface from file ctlsurf_idx 1365 1240 ! ---------------------------------------------- 1366 IF( n ctls_trc == 1 ) THEN1241 IF( nn_ctls_trc == 1 ) THEN 1367 1242 CALL ctl_opn( inum, 'ctlsurf_idx', 'OLD', 'UNFORMATTED', 'SEQUENTIAL', -1, numout, lwp ) 1368 1243 READ ( inum ) nbol_trc … … 1378 1253 #else 1379 1254 ! clmxl = legend root for netCDF output 1380 IF( n ctls_trc == 0 ) THEN ! control surface = mixed-layer with density criterion1255 IF( nn_ctls_trc == 0 ) THEN ! control surface = mixed-layer with density criterion 1381 1256 clmxl = 'Mixed Layer ' 1382 ELSE IF( n ctls_trc == 1 ) THEN ! control surface = read index from file1257 ELSE IF( nn_ctls_trc == 1 ) THEN ! control surface = read index from file 1383 1258 clmxl = ' Bowl ' 1384 ELSE IF( n ctls_trc >= 2 ) THEN ! control surface = model level1385 WRITE(clmxl,'(A10,I2,1X)') 'Levels 1 -', n ctls_trc1259 ELSE IF( nn_ctls_trc >= 2 ) THEN ! control surface = model level 1260 WRITE(clmxl,'(A10,I2,1X)') 'Levels 1 -', nn_ctls_trc 1386 1261 ENDIF 1387 1262 … … 1395 1270 STOP 'trd_mld_trc : this was never checked. Comment this line to proceed...' 1396 1271 ENDIF 1397 zsto = n trd_trc * rdt1272 zsto = nn_trd_trc * rdt 1398 1273 clop = "inst("//TRIM(clop)//")" 1399 1274 # else … … 1401 1276 zsto = rdt ! inst. diags : we use IOIPSL time averaging 1402 1277 ELSE 1403 zsto = n trd_trc * rdt ! mean diags : we DO NOT use any IOIPSL time averaging1278 zsto = nn_trd_trc * rdt ! mean diags : we DO NOT use any IOIPSL time averaging 1404 1279 ENDIF 1405 1280 clop = "ave("//TRIM(clop)//")" 1406 1281 # endif 1407 zout = n trd_trc * rdt1282 zout = nn_trd_trc * rdt 1408 1283 1409 1284 IF(lwp) WRITE (numout,*) ' netCDF initialization' … … 1424 1299 ! ==> choose them according to trdmld_trc_oce.F90 <== 1425 1300 1426 #if defined key_diaeiv1427 cleiv = " (*** only total EIV is meaningful ***)" ! eiv advec. trends require u_eiv, v_eiv1428 #else1429 cleiv = " "1430 #endif1431 1301 ctrd_trc(jpmld_trc_xad ,1) = " Zonal advection" ; ctrd_trc(jpmld_trc_xad ,2) = "_xad" 1432 1302 ctrd_trc(jpmld_trc_yad ,1) = " Meridional advection" ; ctrd_trc(jpmld_trc_yad ,2) = "_yad" … … 1434 1304 ctrd_trc(jpmld_trc_ldf ,1) = " Lateral diffusion" ; ctrd_trc(jpmld_trc_ldf ,2) = "_ldf" 1435 1305 ctrd_trc(jpmld_trc_zdf ,1) = " Vertical diff. (Kz)" ; ctrd_trc(jpmld_trc_zdf ,2) = "_zdf" 1436 ctrd_trc(jpmld_trc_xei ,1) = " Zonal EIV advection"//cleiv ; ctrd_trc(jpmld_trc_xei ,2) = "_xei"1437 ctrd_trc(jpmld_trc_yei ,1) = " Merid. EIV advection"//cleiv ; ctrd_trc(jpmld_trc_yei ,2) = "_yei"1438 ctrd_trc(jpmld_trc_zei ,1) = " Vertical EIV advection"//cleiv ; ctrd_trc(jpmld_trc_zei ,2) = "_zei"1439 ctrd_trc(jpmld_trc_bbc ,1) = " Geothermal flux" ; ctrd_trc(jpmld_trc_bbc ,2) = "_bbc"1440 1306 ctrd_trc(jpmld_trc_bbl ,1) = " Adv/diff. Bottom boundary layer" ; ctrd_trc(jpmld_trc_bbl ,2) = "_bbl" 1441 1307 ctrd_trc(jpmld_trc_dmp ,1) = " Tracer damping" ; ctrd_trc(jpmld_trc_dmp ,2) = "_dmp" … … 1445 1311 ctrd_trc(jpmld_trc_radn ,1) = " Correct negative concentrations" ; ctrd_trc(jpmld_trc_radn ,2) = "_radn" 1446 1312 ctrd_trc(jpmld_trc_atf ,1) = " Asselin time filter" ; ctrd_trc(jpmld_trc_atf ,2) = "_atf" 1447 ctrd_trc(jpltrd_trc+1 ,1) = " Total EIV"//cleiv ; ctrd_trc(jpltrd_trc+1 ,2) = "_tei"1448 1313 1449 1314 DO jn = 1, jptra 1450 1315 !-- Create a NetCDF file and enter the define mode 1451 IF( l uttrd(jn) ) THEN1316 IF( ln_trdtrc(jn) ) THEN 1452 1317 csuff="ML_"//ctrcnm(jn) 1453 CALL dia_nam( clhstnam, n trd_trc, csuff )1318 CALL dia_nam( clhstnam, nn_trd_trc, csuff ) 1454 1319 CALL histbeg( clhstnam, jpi, glamt, jpj, gphit, & 1455 & 1, jpi, 1, jpj, nittrc000-n dttrc, zjulian, rdt, nh_t(jn), nidtrd(jn), domain_id=nidom )1320 & 1, jpi, 1, jpj, nittrc000-nn_dttrc, zjulian, rdt, nh_t(jn), nidtrd(jn), domain_id=nidom ) 1456 1321 1457 1322 !-- Define the ML depth variable … … 1464 1329 #if defined key_lobster 1465 1330 !-- Create a NetCDF file and enter the define mode 1466 CALL dia_nam( clhstnam, n trd_trc, 'trdbio' )1331 CALL dia_nam( clhstnam, nn_trd_trc, 'trdbio' ) 1467 1332 CALL histbeg( clhstnam, jpi, glamt, jpj, gphit, & 1468 & 1, jpi, 1, jpj, nittrc000-n dttrc, zjulian, rdt, nh_tb, nidtrdbio, domain_id=nidom )1333 & 1, jpi, 1, jpj, nittrc000-nn_dttrc, zjulian, rdt, nh_tb, nidtrdbio, domain_id=nidom ) 1469 1334 #endif 1470 1335 1471 1336 !-- Define physical units 1472 IF( ucf_trc == 1. ) THEN1337 IF( rn_ucf_trc == 1. ) THEN 1473 1338 cltrcu = "(mmole-N/m3)/sec" ! all passive tracers have the same unit 1474 ELSEIF ( ucf_trc == 3600.*24.) THEN ! ??? trop long : seulement (mmole-N/m3)1339 ELSEIF ( rn_ucf_trc == 3600.*24.) THEN ! ??? trop long : seulement (mmole-N/m3) 1475 1340 cltrcu = "(mmole-N/m3)/day" ! ??? apparait dans les sorties netcdf 1476 1341 ELSE … … 1485 1350 DO jn = 1, jptra 1486 1351 ! 1487 IF( l uttrd(jn) ) THEN1352 IF( ln_trdtrc(jn) ) THEN 1488 1353 clvar = trim(ctrcnm(jn))//"ml" ! e.g. detml, zooml, no3ml, etc. 1489 1354 CALL histdef(nidtrd(jn), clvar, clmxl//" "//trim(ctrcnm(jn))//" Mixed Layer ", & … … 1504 1369 CALL histdef(nidtrd(jn), trim(clvar//ctrd_trc(jpmld_trc_atf,2)), clmxl//" "//clvar//ctrd_trc(jpmld_trc_atf,1), & 1505 1370 & cltrcu, jpi, jpj, nh_t(jn), 1 , 1, 1 , -99 , 32, clop, zout, zout ) ! IOIPSL: NO time mean 1506 1507 CALL histdef(nidtrd(jn), trim(clvar//ctrd_trc(jpltrd_trc+1,2)), clmxl//" "//clvar//ctrd_trc(jpltrd_trc+1 ,1), &1508 & cltrcu, jpi, jpj, nh_t(jn), 1 , 1, 1 , -99 , 32, clop, zsto, zout ) ! Total EIV1509 1371 ! 1510 1372 ENDIF … … 1520 1382 !-- Leave IOIPSL/NetCDF define mode 1521 1383 DO jn = 1, jptra 1522 IF( l uttrd(jn) ) CALL histend( nidtrd(jn) )1384 IF( ln_trdtrc(jn) ) CALL histend( nidtrd(jn) ) 1523 1385 END DO 1524 1386 … … 1539 1401 !!---------------------------------------------------------------------- 1540 1402 1541 INTERFACE trd_mod_trc1542 MODULE PROCEDURE trd_mod_trc_trp, trd_mod_trc_bio1543 END INTERFACE1544 1545 1403 CONTAINS 1546 1404 … … 1554 1412 WRITE(*,*) 'trd_mld_bio: You should not have seen this print! error?', kt 1555 1413 END SUBROUTINE trd_mld_bio 1556 1557 SUBROUTINE trd_mod_trc_bio( ptrbio, ktrd, kt )1558 INTEGER , INTENT( in ) :: kt ! time step1559 INTEGER , INTENT( in ) :: ktrd ! bio trend index1560 REAL, DIMENSION(:,:,:), INTENT( inout ) :: ptrbio ! Bio trend1561 WRITE(*,*) 'trd_mod_trc_bio : You should not have seen this print! error?', ptrbio(1,1,1)1562 WRITE(*,*) ' " " : You should not have seen this print! error?', ktrd1563 WRITE(*,*) ' " " : You should not have seen this print! error?', kt1564 END SUBROUTINE trd_mod_trc_bio1565 1566 SUBROUTINE trd_mod_trc_trp( ptrtrd, kjn, ktrd, kt )1567 INTEGER , INTENT( in ) :: kt ! time step1568 INTEGER , INTENT( in ) :: kjn ! tracer index1569 INTEGER , INTENT( in ) :: ktrd ! tracer trend index1570 REAL, DIMENSION(:,:,:), INTENT( inout ) :: ptrtrd ! Temperature or U trend1571 WRITE(*,*) 'trd_mod_trc_trp : You should not have seen this print! error?', ptrtrd(1,1,1)1572 WRITE(*,*) ' " " : You should not have seen this print! error?', kjn1573 WRITE(*,*) ' " " : You should not have seen this print! error?', ktrd1574 WRITE(*,*) ' " " : You should not have seen this print! error?', kt1575 END SUBROUTINE trd_mod_trc_trp1576 1414 1577 1415 SUBROUTINE trd_mld_trc_zint( ptrc_trdmld, ktrd, ctype, kjn ) -
branches/DEV_r2006_merge_TRA_TRC/NEMO/TOP_SRC/TRP/trdmld_trc_rst.F90
r1473 r2030 12 12 USE in_out_manager ! I/O manager 13 13 USE iom ! I/O module 14 USE trc ! for n dttrc ctrcnm15 USE trdm ld_trc_oce ! for lk_trdmld_trc14 USE trc ! for nn_dttrc ctrcnm 15 USE trdmod_trc_oce ! for lk_trdmld_trc 16 16 17 17 IMPLICIT NONE … … 45 45 !!-------------------------------------------------------------------------------- 46 46 47 IF( kt == nitrst - n dttrc .OR. nitend - nit000 + 1 < 2 * ndttrc ) THEN ! idem trcrst.F9047 IF( kt == nitrst - nn_dttrc .OR. nitend - nit000 + 1 < 2 * nn_dttrc ) THEN ! idem trcrst.F90 48 48 IF( nitrst > 1.0e9 ) THEN 49 49 WRITE(clkt,*) nitrst
Note: See TracChangeset
for help on using the changeset viewer.