Changeset 4315
- Timestamp:
- 2013-11-26T13:48:57+01:00 (10 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2013/dev_r3411_CNRS4_IOCRS/NEMOGCM/CONFIG/ORCA2_LIM_CRS_SAL/MY_SRC/trabbl_crs.F90
r4056 r4315 1 MODULE trabbl_crs 2 !!============================================================================== 3 !! *** MODULE trabbl *** 4 !! Ocean physics : advective and/or diffusive bottom boundary layer scheme 5 !!============================================================================== 6 !! History : OPA ! 1996-06 (L. Mortier) Original code 7 !! 8.0 ! 1997-11 (G. Madec) Optimization 8 !! NEMO 1.0 ! 2002-08 (G. Madec) free form + modules 9 !! - ! 2004-01 (A. de Miranda, G. Madec, J.M. Molines ) add advective bbl 10 !! 3.3 ! 2009-11 (G. Madec) merge trabbl and trabbl_adv + style + optimization 11 !! - ! 2010-04 (G. Madec) Campin & Goosse advective bbl 12 !! - ! 2010-06 (C. Ethe, G. Madec) merge TRA-TRC 13 !! - ! 2010-11 (G. Madec) add mbk. arrays associated to the deepest ocean level 14 !!---------------------------------------------------------------------- 15 #if defined key_trabbl || defined key_esopa 16 !!---------------------------------------------------------------------- 17 !! 'key_trabbl' or bottom boundary layer 18 !!---------------------------------------------------------------------- 19 !! tra_bbl_alloc : allocate trabbl arrays 20 !! tra_bbl : update the tracer trends due to the bottom boundary layer (advective and/or diffusive) 21 !! tra_bbl_dif : generic routine to compute bbl diffusive trend 22 !! tra_bbl_adv : generic routine to compute bbl advective trend 23 !! bbl : computation of bbl diffu. flux coef. & transport in bottom boundary layer 24 !! tra_bbl_init : initialization, namelist read, parameters control 25 !!---------------------------------------------------------------------- 26 USE oce ! ocean dynamics and active tracers 27 USE dom_oce ! ocean space and time domain 28 USE phycst ! physical constant 29 USE eosbn2 ! equation of state 30 USE trdmod_oce ! trends: ocean variables 31 USE trdtra ! trends: active tracers 32 USE iom ! IOM server 33 USE in_out_manager ! I/O manager 34 USE lbclnk ! ocean lateral boundary conditions 35 USE prtctl ! Print control 36 USE wrk_nemo ! Memory Allocation 37 USE timing ! Timing 38 USE crs 39 USE crslbclnk 40 USE crsiom 41 42 43 IMPLICIT NONE 44 PRIVATE 45 46 PUBLIC tra_bbl_crs ! routine called by step.F90 47 PUBLIC tra_bbl_init_crs ! routine called by opa.F90 48 PUBLIC tra_bbl_dif_crs ! routine called by trcbbl.F90 49 PUBLIC tra_bbl_adv_crs ! - - - - 50 PUBLIC bbl_crs ! routine called by trcbbl.F90 and dtadyn.F90 51 52 LOGICAL, PUBLIC, PARAMETER :: lk_trabbl_crs = .TRUE. !: bottom boundary layer flag 53 54 ! !!* Namelist nambbl * 55 INTEGER , PUBLIC :: nn_bbl_ldf = 0 !: =1 : diffusive bbl or not (=0) 56 INTEGER , PUBLIC :: nn_bbl_adv = 0 !: =1/2 : advective bbl or not (=0) 57 ! ! =1 : advective bbl using the bottom ocean velocity 58 ! ! =2 : - - using utr_bbl proportional to grad(rho) 59 REAL(wp), PUBLIC :: rn_ahtbbl_crs = 1.e4_wp !: along slope bbl diffusive coefficient [m2/s] 60 REAL(wp), PUBLIC :: rn_ahtbbl = 1.e3_wp 61 REAL(wp), PUBLIC :: rn_gambbl = 10.0_wp !: lateral coeff. for bottom boundary layer scheme [s] 62 63 LOGICAL , PUBLIC :: l_bbl !: flag to compute bbl diffu. flux coef and transport 64 65 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC :: utr_bbl_crs , vtr_bbl_crs ! u- (v-) transport in the bottom boundary layer 66 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC :: ahu_bbl_crs , ahv_bbl_crs ! masked diffusive bbl coeff. at u & v-pts 67 68 INTEGER , ALLOCATABLE, SAVE, DIMENSION(:,:) :: mbku_d , mbkv_d ! vertical index of the "lower" bottom ocean U/V-level 69 INTEGER , ALLOCATABLE, SAVE, DIMENSION(:,:) :: mgrhu , mgrhv ! = +/-1, sign of grad(H) in u-(v-)direction 70 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: ahu_bbl_0, ahv_bbl_0 ! diffusive bbl flux coefficients at u and v-points 71 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: e3u_bbl_0, e3v_bbl_0 ! thichness of the bbl (e3) at u and v-points 72 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: r1_e1e2t ! inverse of the cell surface at t-point [1/m2] 73 74 !! * Substitutions 75 # include "domzgr_substitute.h90" 76 # include "vectopt_loop_substitute.h90" 77 !!---------------------------------------------------------------------- 78 !! NEMO/OPA 3.3 , NEMO Consortium (2010) 79 !! $Id: trabbl.F90 3294 2012-01-28 16:44:18Z rblod $ 80 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 81 !!---------------------------------------------------------------------- 82 CONTAINS 83 84 INTEGER FUNCTION tra_bbl_alloc_crs() 85 !!---------------------------------------------------------------------- 86 !! *** FUNCTION tra_bbl_alloc *** 87 !!---------------------------------------------------------------------- 88 ALLOCATE( utr_bbl_crs (jpi,jpj) , ahu_bbl_crs (jpi,jpj) , mbku_d (jpi,jpj) , mgrhu(jpi,jpj) , & 89 & vtr_bbl_crs (jpi,jpj) , ahv_bbl_crs (jpi,jpj) , mbkv_d (jpi,jpj) , mgrhv(jpi,jpj) , & 90 & ahu_bbl_0(jpi,jpj) , ahv_bbl_0(jpi,jpj) , & 91 & e3u_bbl_0(jpi,jpj) , e3v_bbl_0(jpi,jpj) , r1_e1e2t(jpi,jpj) , STAT= tra_bbl_alloc_crs ) 92 ! 93 IF( lk_mpp ) CALL mpp_sum ( tra_bbl_alloc_crs ) 94 IF( tra_bbl_alloc_crs > 0 ) CALL ctl_warn('tra_bbl_alloc: allocation of arrays failed.') 95 END FUNCTION tra_bbl_alloc_crs 96 97 98 SUBROUTINE tra_bbl_crs( kt ) 99 !!---------------------------------------------------------------------- 100 !! *** ROUTINE bbl *** 101 !! 102 !! ** Purpose : Compute the before tracer (t & s) trend associated 103 !! with the bottom boundary layer and add it to the general 104 !! trend of tracer equations. 105 !! 106 !! ** Method : Depending on namtra_bbl namelist parameters the bbl 107 !! diffusive and/or advective contribution to the tracer trend 108 !! is added to the general tracer trend 109 !!---------------------------------------------------------------------- 110 INTEGER, INTENT( in ) :: kt ! ocean time-step 111 !! 112 REAL(wp), POINTER, DIMENSION(:,:,:) :: ztrdt, ztrds 113 !!---------------------------------------------------------------------- 114 ! 115 IF( nn_timing == 1 ) CALL timing_start( 'tra_bbl') 116 ! 117 IF( l_trdtra ) THEN !* Save ta and sa trends 118 CALL wrk_alloc( jpi, jpj, jpk, ztrdt, ztrds ) 119 ztrdt(:,:,:) = tsa_crs(:,:,:,jp_tem) 120 ztrds(:,:,:) = tsa_crs(:,:,:,jp_sal) 121 ENDIF 122 123 IF( l_bbl ) CALL bbl_crs( kt, nit000, 'TRA' ) !* bbl coef. and transport (only if not already done in trcbbl) 124 125 IF( nn_bbl_ldf == 1 ) THEN !* Diffusive bbl 126 ! 127 CALL tra_bbl_dif_crs( tsb_crs, tsa_crs, jpts ) 128 IF( ln_ctl ) & 129 CALL prt_ctl( tab3d_1=tsa_crs(:,:,:,jp_tem), clinfo1=' bbl_ldf - Ta: ', mask1=tmask_crs, & 130 & tab3d_2=tsa_crs(:,:,:,jp_sal), clinfo2= ' Sa: ', mask2=tmask_crs, clinfo3='tra' ) 131 ! lateral boundary conditions ; just need for outputs 132 CALL crs_lbc_lnk( ahu_bbl_crs, 'U', 1. ) ; CALL crs_lbc_lnk( ahv_bbl_crs, 'V', 1. ) 133 ! CALL crs_iom_put( "ahu_bbl_crs", ahu_bbl ) ! bbl diffusive flux i-coef 134 ! CALL crs_iom_put( "ahv_bbl_crs", ahv_bbl ) ! bbl diffusive flux j-coef 135 ! 136 END IF 137 138 IF( nn_bbl_adv /= 0 ) THEN !* Advective bbl 139 ! 140 CALL tra_bbl_adv_crs( tsb_crs, tsa_crs, jpts ) 141 IF(ln_ctl) & 142 CALL prt_ctl( tab3d_1=tsa_crs(:,:,:,jp_tem), clinfo1=' bbl_adv - Ta: ', mask1=tmask_crs, & 143 & tab3d_2=tsa_crs(:,:,:,jp_sal), clinfo2= ' Sa: ', mask2=tmask_crs, clinfo3='tra' ) 144 ! lateral boundary conditions ; just need for outputs 145 CALL crs_lbc_lnk( utr_bbl_crs, 'U', 1. ) ; CALL crs_lbc_lnk( vtr_bbl_crs, 'V', 1. ) 146 ! CALL crs_iom_put( "uoce_bbl", utr_bbl ) ! bbl i-transport 147 ! CALL dom_grid_crs !! cc return at the CRS grid 148 ! CALL crs_iom_put( "voce_bbl", vtr_bbl ) ! bbl j-transport 149 ! CALL dom_grid_crs !! cc return at the CRS grid 150 ! 151 END IF 152 153 IF( l_trdtra ) THEN ! save the horizontal diffusive trends for further diagnostics 154 ztrdt(:,:,:) = tsa_crs(:,:,:,jp_tem) - ztrdt(:,:,:) 155 ztrds(:,:,:) = tsa_crs(:,:,:,jp_sal) - ztrds(:,:,:) 156 CALL trd_tra( kt, 'TRA', jp_tem, jptra_trd_bbl, ztrdt ) 157 CALL trd_tra( kt, 'TRA', jp_sal, jptra_trd_bbl, ztrds ) 158 CALL wrk_dealloc( jpi, jpj, jpk, ztrdt, ztrds ) 159 ENDIF 160 ! 161 IF( nn_timing == 1 ) CALL timing_stop( 'tra_bbl') 162 ! 163 END SUBROUTINE tra_bbl_crs 164 165 166 SUBROUTINE tra_bbl_dif_crs( ptb, pta, kjpt ) 167 !!---------------------------------------------------------------------- 168 !! *** ROUTINE tra_bbl_dif *** 169 !! 170 !! ** Purpose : Computes the bottom boundary horizontal and vertical 171 !! advection terms. 172 !! 173 !! ** Method : 174 !! * diffusive bbl (nn_bbl_ldf=1) : 175 !! When the product grad( rho) * grad(h) < 0 (where grad is an 176 !! along bottom slope gradient) an additional lateral 2nd order 177 !! diffusion along the bottom slope is added to the general 178 !! tracer trend, otherwise the additional trend is set to 0. 179 !! A typical value of ahbt is 2000 m2/s (equivalent to 180 !! a downslope velocity of 20 cm/s if the condition for slope 181 !! convection is satified) 182 !! 183 !! ** Action : pta increased by the bbl diffusive trend 184 !! 185 !! References : Beckmann, A., and R. Doscher, 1997, J. Phys.Oceanogr., 581-591. 186 !! Campin, J.-M., and H. Goosse, 1999, Tellus, 412-430. 187 !!---------------------------------------------------------------------- 188 ! 189 INTEGER , INTENT(in ) :: kjpt ! number of tracers 190 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: ptb ! before and now tracer fields 191 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: pta ! tracer trend 192 ! 193 INTEGER :: ji, jj, jn ! dummy loop indices 194 INTEGER :: ik ! local integers 195 REAL(wp) :: zbtr ! local scalars 196 REAL(wp), POINTER, DIMENSION(:,:) :: zptb 197 !!---------------------------------------------------------------------- 198 ! 199 IF( nn_timing == 1 ) CALL timing_start('tra_bbl_dif') 200 ! 201 CALL wrk_alloc( jpi, jpj, zptb ) 202 ! 203 DO jn = 1, kjpt ! tracer loop 204 ! ! =========== 205 # if defined key_vectopt_loop 206 DO jj = 1, 1 ! vector opt. (forced unrolling) 207 DO ji = 1, jpij 208 #else 209 DO jj = 1, jpj 210 DO ji = 1, jpi 211 #endif 212 ik = mbkt_crs(ji,jj) ! bottom T-level index 213 zptb(ji,jj) = ptb(ji,jj,ik,jn) ! bottom before T and S 214 END DO 215 END DO 216 ! ! Compute the trend 217 # if defined key_vectopt_loop 218 DO jj = 1, 1 ! vector opt. (forced unrolling) 219 DO ji = jpi+1, jpij-jpi-1 220 # else 221 DO jj = 2, jpjm1 222 DO ji = 2, jpim1 223 # endif 224 ik = mbkt_crs(ji,jj) ! bottom T-level index 225 zbtr = r1_bt_crs(ji,jj,ik) 226 pta(ji,jj,ik,jn) = pta(ji,jj,ik,jn) & 227 & + ( ahu_bbl_crs(ji ,jj ) * ( zptb(ji+1,jj ) - zptb(ji ,jj ) ) & 228 & - ahu_bbl_crs(ji-1,jj ) * ( zptb(ji ,jj ) - zptb(ji-1,jj ) ) & 229 & + ahv_bbl_crs(ji ,jj ) * ( zptb(ji ,jj+1) - zptb(ji ,jj ) ) & 230 & - ahv_bbl_crs(ji ,jj-1) * ( zptb(ji ,jj ) - zptb(ji ,jj-1) ) ) * zbtr 231 END DO 232 END DO 233 ! ! =========== 234 END DO ! end tracer 235 ! ! =========== 236 CALL wrk_dealloc( jpi, jpj, zptb ) 237 ! 238 IF( nn_timing == 1 ) CALL timing_stop('tra_bbl_dif') 239 ! 240 END SUBROUTINE tra_bbl_dif_crs 241 242 243 SUBROUTINE tra_bbl_adv_crs( ptb, pta, kjpt ) 244 !!---------------------------------------------------------------------- 245 !! *** ROUTINE trc_bbl *** 246 !! 247 !! ** Purpose : Compute the before passive tracer trend associated 248 !! with the bottom boundary layer and add it to the general trend 249 !! of tracer equations. 250 !! ** Method : advective bbl (nn_bbl_adv = 1 or 2) : 251 !! nn_bbl_adv = 1 use of the ocean near bottom velocity as bbl velocity 252 !! nn_bbl_adv = 2 follow Campin and Goosse (1999) implentation i.e. 253 !! transport proportional to the along-slope density gradient 254 !! 255 !! References : Beckmann, A., and R. Doscher, 1997, J. Phys.Oceanogr., 581-591. 256 !! Campin, J.-M., and H. Goosse, 1999, Tellus, 412-430. 257 !!---------------------------------------------------------------------- 258 INTEGER , INTENT(in ) :: kjpt ! number of tracers 259 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: ptb ! before and now tracer fields 260 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: pta ! tracer trend 261 ! 262 INTEGER :: ji, jj, jk, jn ! dummy loop indices 263 INTEGER :: iis , iid , ijs , ijd ! local integers 264 INTEGER :: ikus, ikud, ikvs, ikvd ! - - 265 REAL(wp) :: zbtr, ztra ! local scalars 266 REAL(wp) :: zu_bbl, zv_bbl ! - - 267 !!---------------------------------------------------------------------- 268 ! 269 IF( nn_timing == 1 ) CALL timing_start( 'tra_bbl_adv') 270 ! ! =========== 271 DO jn = 1, kjpt ! tracer loop 272 ! ! =========== 273 # if defined key_vectopt_loop 274 DO jj = 1, 1 275 DO ji = 1, jpij-jpi-1 ! vector opt. (forced unrolling) 276 # else 277 DO jj = 1, jpjm1 278 DO ji = 1, jpim1 ! CAUTION start from i=1 to update i=2 when cyclic east-west 279 # endif 280 IF( utr_bbl_crs(ji,jj) /= 0.e0 ) THEN ! non-zero i-direction bbl advection 281 ! down-slope i/k-indices (deep) & up-slope i/k indices (shelf) 282 iid = ji + MAX( 0, mgrhu(ji,jj) ) ; iis = ji + 1 - MAX( 0, mgrhu(ji,jj) ) 283 ikud = mbku_d(ji,jj) ; ikus = mbku_crs(ji,jj) 284 zu_bbl = ABS( utr_bbl_crs(ji,jj) ) 285 ! 286 ! ! up -slope T-point (shelf bottom point) 287 zbtr = r1_bt_crs(iis,jj,ikus) 288 ztra = zu_bbl * ( ptb(iid,jj,ikus,jn) - ptb(iis,jj,ikus,jn) ) * zbtr 289 pta(iis,jj,ikus,jn) = pta(iis,jj,ikus,jn) + ztra 290 ! 291 DO jk = ikus, ikud-1 ! down-slope upper to down T-point (deep column) 292 zbtr = r1_bt_crs(iid,jj,jk) 293 ztra = zu_bbl * ( ptb(iid,jj,jk+1,jn) - ptb(iid,jj,jk,jn) ) * zbtr 294 pta(iid,jj,jk,jn) = pta(iid,jj,jk,jn) + ztra 295 END DO 296 ! 297 zbtr = r1_bt_crs(iid,jj,ikud) 298 ztra = zu_bbl * ( ptb(iis,jj,ikus,jn) - ptb(iid,jj,ikud,jn) ) * zbtr 299 pta(iid,jj,ikud,jn) = pta(iid,jj,ikud,jn) + ztra 300 ENDIF 301 ! 302 IF( vtr_bbl_crs(ji,jj) /= 0.e0 ) THEN ! non-zero j-direction bbl advection 303 ! down-slope j/k-indices (deep) & up-slope j/k indices (shelf) 304 ijd = jj + MAX( 0, mgrhv(ji,jj) ) ; ijs = jj + 1 - MAX( 0, mgrhv(ji,jj) ) 305 ikvd = mbkv_d(ji,jj) ; ikvs = mbkv_crs(ji,jj) 306 zv_bbl = ABS( vtr_bbl_crs(ji,jj) ) 307 ! 308 ! up -slope T-point (shelf bottom point) 309 zbtr = r1_bt_crs(ji,ijs,ikvs) 310 ztra = zv_bbl * ( ptb(ji,ijd,ikvs,jn) - ptb(ji,ijs,ikvs,jn) ) * zbtr 311 pta(ji,ijs,ikvs,jn) = pta(ji,ijs,ikvs,jn) + ztra 312 ! 313 DO jk = ikvs, ikvd-1 ! down-slope upper to down T-point (deep column) 314 zbtr = r1_bt_crs(ji,ijd,jk) 315 ztra = zv_bbl * ( ptb(ji,ijd,jk+1,jn) - ptb(ji,ijd,jk,jn) ) * zbtr 316 pta(ji,ijd,jk,jn) = pta(ji,ijd,jk,jn) + ztra 317 END DO 318 ! ! down-slope T-point (deep bottom point) 319 zbtr = r1_bt_crs(ji,ijd,ikvd) 320 ztra = zv_bbl * ( ptb(ji,ijs,ikvs,jn) - ptb(ji,ijd,ikvd,jn) ) * zbtr 321 pta(ji,ijd,ikvd,jn) = pta(ji,ijd,ikvd,jn) + ztra 322 ENDIF 323 END DO 324 ! 325 END DO 326 ! ! =========== 327 END DO ! end tracer 328 ! ! =========== 329 ! 330 IF( nn_timing == 1 ) CALL timing_stop( 'tra_bbl_adv') 331 ! 332 END SUBROUTINE tra_bbl_adv_crs 333 334 335 SUBROUTINE bbl_crs( kt, kit000, cdtype ) 336 !!---------------------------------------------------------------------- 337 !! *** ROUTINE bbl *** 338 !! 339 !! ** Purpose : Computes the bottom boundary horizontal and vertical 340 !! advection terms. 341 !! 342 !! ** Method : 343 !! * diffusive bbl (nn_bbl_ldf=1) : 344 !! When the product grad( rho) * grad(h) < 0 (where grad is an 345 !! along bottom slope gradient) an additional lateral 2nd order 346 !! diffusion along the bottom slope is added to the general 347 !! tracer trend, otherwise the additional trend is set to 0. 348 !! A typical value of ahbt is 2000 m2/s (equivalent to 349 !! a downslope velocity of 20 cm/s if the condition for slope 350 !! convection is satified) 351 !! * advective bbl (nn_bbl_adv=1 or 2) : 352 !! nn_bbl_adv = 1 use of the ocean velocity as bbl velocity 353 !! nn_bbl_adv = 2 follow Campin and Goosse (1999) implentation 354 !! i.e. transport proportional to the along-slope density gradient 355 !! 356 !! NB: the along slope density gradient is evaluated using the 357 !! local density (i.e. referenced at a common local depth). 358 !! 359 !! References : Beckmann, A., and R. Doscher, 1997, J. Phys.Oceanogr., 581-591. 360 !! Campin, J.-M., and H. Goosse, 1999, Tellus, 412-430. 361 !!---------------------------------------------------------------------- 362 ! 363 INTEGER , INTENT(in ) :: kt ! ocean time-step index 364 INTEGER , INTENT(in ) :: kit000 ! first time step index 365 CHARACTER(len=3), INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator) 366 !! 367 INTEGER :: ji, jj ! dummy loop indices 368 INTEGER :: ik ! local integers 369 INTEGER :: iis , iid , ijs , ijd ! - - 370 INTEGER :: ikus, ikud, ikvs, ikvd ! - - 371 REAL(wp) :: zsign, zsigna, zgbbl ! local scalars 372 REAL(wp) :: zgdrho, zt, zs, zh ! - - 373 !! 374 REAL(wp) :: fsalbt, fsbeta, pft, pfs, pfh ! statement function 375 REAL(wp), POINTER, DIMENSION(:,:) :: zub, zvb, ztb, zsb, zdep 376 !!----------------------- zv_bbl----------------------------------------------- 377 ! ratio alpha/beta = fsalbt : ratio of thermal over saline expension coefficients 378 ! ================ pft : potential temperature in degrees celcius 379 ! pfs : salinity anomaly (s-35) in psu 380 ! pfh : depth in meters 381 ! nn_eos = 0 (Jackett and McDougall 1994 formulation) 382 fsalbt( pft, pfs, pfh ) = & ! alpha/beta 383 ( ( ( -0.255019e-07 * pft + 0.298357e-05 ) * pft & 384 - 0.203814e-03 ) * pft & 385 + 0.170907e-01 ) * pft & 386 + 0.665157e-01 & 387 +(-0.678662e-05 * pfs - 0.846960e-04 * pft + 0.378110e-02 ) * pfs & 388 + ( ( - 0.302285e-13 * pfh & 389 - 0.251520e-11 * pfs & 390 + 0.512857e-12 * pft * pft ) * pfh & 391 - 0.164759e-06 * pfs & 392 +( 0.791325e-08 * pft - 0.933746e-06 ) * pft & 393 + 0.380374e-04 ) * pfh 394 fsbeta( pft, pfs, pfh ) = & ! beta 395 ( ( -0.415613e-09 * pft + 0.555579e-07 ) * pft & 396 - 0.301985e-05 ) * pft & 397 + 0.785567e-03 & 398 + ( 0.515032e-08 * pfs & 399 + 0.788212e-08 * pft - 0.356603e-06 ) * pfs & 400 +( ( 0.121551e-17 * pfh & 401 - 0.602281e-15 * pfs & 402 - 0.175379e-14 * pft + 0.176621e-12 ) * pfh & 403 + 0.408195e-10 * pfs & 404 + ( - 0.213127e-11 * pft + 0.192867e-09 ) * pft & 405 - 0.121555e-07 ) * pfh 406 !!---------------------------------------------------------------------- 407 408 ! 409 IF( nn_timing == 1 ) CALL timing_start( 'bbl') 410 ! 411 CALL wrk_alloc( jpi, jpj, zub, zvb, ztb, zsb, zdep ) 412 ! 413 414 IF( kt == kit000 ) THEN 415 IF(lwp) WRITE(numout,*) 416 IF(lwp) WRITE(numout,*) 'trabbl:bbl : Compute bbl velocities and diffusive coefficients in ', cdtype 417 IF(lwp) WRITE(numout,*) '~~~~~~~~~~' 418 ENDIF 419 420 ! !* bottom temperature, salinity, velocity and depth 421 #if defined key_vectopt_loop 422 DO jj = 1, 1 ! vector opt. (forced unrolling) 423 DO ji = 1, jpij 424 #else 425 DO jj = 1, jpj 426 DO ji = 1, jpi 427 #endif 428 ik = mbkt_crs(ji,jj) ! bottom T-level index 429 ztb (ji,jj) = tsb_crs(ji,jj,ik,jp_tem) * tmask_crs(ji,jj,1) ! bottom before T and S 430 zsb (ji,jj) = tsb_crs(ji,jj,ik,jp_sal) * tmask_crs(ji,jj,1) 431 zdep(ji,jj) = gdept_crs(ji,jj,ik) ! bottom T-level reference depth !! 432 433 zub(ji,jj) = un_crs(ji,jj,mbku_crs(ji,jj)) ! bottom velocity 434 zvb(ji,jj) = vn_crs(ji,jj,mbkv_crs(ji,jj)) 435 END DO 436 END DO 437 438 ! !-------------------! 439 IF( nn_bbl_ldf == 1 ) THEN ! diffusive bbl ! 440 ! !-------------------! 441 DO jj = 1, jpjm1 ! (criteria for non zero flux: grad(rho).grad(h) < 0 ) 442 DO ji = 1, jpim1 443 ! ! i-direction 444 zt = 0.5 * ( ztb (ji,jj) + ztb (ji+1,jj) ) ! T, S anomalie, and depth 445 zs = 0.5 * ( zsb (ji,jj) + zsb (ji+1,jj) ) - 35.0 446 zh = 0.5 * ( zdep(ji,jj) + zdep(ji+1,jj) ) 447 ! ! masked bbl i-gradient of density 448 zgdrho = ( fsalbt( zt, zs, zh ) * ( ztb(ji+1,jj) - ztb(ji,jj) ) & 449 & - ( zsb(ji+1,jj) - zsb(ji,jj) ) ) * umask_crs(ji,jj,1) 450 ! 451 zsign = SIGN( 0.5, - zgdrho * REAL( mgrhu(ji,jj) ) ) ! sign of ( i-gradient * i-slope ) 452 ahu_bbl_crs(ji,jj) = ( 0.5 - zsign ) * ahu_bbl_0(ji,jj) ! masked diffusive flux coeff. 453 ! 454 ! ! j-direction 455 zt = 0.5 * ( ztb (ji,jj+1) + ztb (ji,jj) ) ! T, S anomalie, and depth 456 zs = 0.5 * ( zsb (ji,jj+1) + zsb (ji,jj) ) - 35.0 457 zh = 0.5 * ( zdep(ji,jj+1) + zdep(ji,jj) ) 458 ! ! masked bbl j-gradient of density 459 zgdrho = ( fsalbt( zt, zs, zh ) * ( ztb(ji,jj+1) - ztb(ji,jj) ) & 460 & - ( zsb(ji,jj+1) - zsb(ji,jj) ) ) * vmask_crs(ji,jj,1) 461 ! 462 zsign = SIGN( 0.5, -zgdrho * REAL( mgrhv(ji,jj) ) ) ! sign of ( j-gradient * j-slope ) 463 ahv_bbl_crs(ji,jj) = ( 0.5 - zsign ) * ahv_bbl_0(ji,jj) 464 ! 465 END DO 466 END DO 467 ! 468 ENDIF 469 470 ! !-------------------! 471 IF( nn_bbl_adv /= 0 ) THEN ! advective bbl ! 472 ! !-------------------! 473 SELECT CASE ( nn_bbl_adv ) !* bbl transport type 474 ! 475 CASE( 1 ) != use of upper velocity 476 DO jj = 1, jpjm1 ! criteria: grad(rho).grad(h)<0 and grad(rho).grad(h)<0 477 DO ji = 1, fs_jpim1 ! vector opt. 478 ! ! i-direction 479 zt = 0.5 * ( ztb (ji,jj) + ztb (ji+1,jj) ) ! T, S anomalie, and depth 480 zs = 0.5 * ( zsb (ji,jj) + zsb (ji+1,jj) ) - 35.0 481 zh = 0.5 * ( zdep(ji,jj) + zdep(ji+1,jj) ) 482 ! ! masked bbl i-gradient of density 483 zgdrho = ( fsalbt( zt, zs, zh ) * ( ztb(ji+1,jj) - ztb(ji,jj) ) & 484 & - ( zsb(ji+1,jj) - zsb(ji,jj) ) ) * umask_crs(ji,jj,1) 485 ! 486 zsign = SIGN( 0.5, - zgdrho * REAL( mgrhu(ji,jj) ) ) ! sign of i-gradient * i-slope 487 zsigna= SIGN( 0.5, zub(ji,jj) * REAL( mgrhu(ji,jj) ) ) ! sign of u * i-slope 488 ! 489 ! ! bbl velocity 490 utr_bbl_crs(ji,jj) = ( 0.5 + zsigna ) * ( 0.5 - zsign ) * e2u_crs(ji,jj) * e3u_bbl_0(ji,jj) * zub(ji,jj) 491 ! 492 ! ! j-direction 493 zt = 0.5 * ( ztb (ji,jj+1) + ztb (ji,jj) ) ! T, S anomalie, and depth 494 zs = 0.5 * ( zsb (ji,jj+1) + zsb (ji,jj) ) - 35.0 495 zh = 0.5 * ( zdep(ji,jj+1) + zdep(ji,jj) ) 496 ! ! masked bbl j-gradient of density 497 zgdrho = ( fsalbt( zt, zs, zh ) * ( ztb(ji,jj+1) - ztb(ji,jj) ) & 498 & - ( zsb(ji,jj+1) - zsb(ji,jj) ) ) * vmask_crs(ji,jj,1) 499 zsign = SIGN( 0.5, - zgdrho * REAL( mgrhv(ji,jj) ) ) ! sign of j-gradient * j-slope 500 zsigna= SIGN( 0.5, zvb(ji,jj) * REAL( mgrhv(ji,jj) ) ) ! sign of u * i-slope 501 ! 502 ! ! bbl velocity 503 vtr_bbl_crs(ji,jj) = ( 0.5 + zsigna ) * ( 0.5 - zsign ) * e1v_crs(ji,jj) * e3v_bbl_0(ji,jj) * zvb(ji,jj) 504 END DO 505 END DO 506 ! 507 CASE( 2 ) != bbl velocity = F( delta rho ) 508 zgbbl = grav * rn_gambbl 509 DO jj = 1, jpjm1 ! criteria: rho_up > rho_down 510 DO ji = 1, fs_jpim1 ! vector opt. 511 ! ! i-direction 512 ! down-slope T-point i/k-index (deep) & up-slope T-point i/k-index (shelf) 513 iid = ji + MAX( 0, mgrhu(ji,jj) ) ; iis = ji + 1 - MAX( 0, mgrhu(ji,jj) ) 514 ikud = mbku_d(ji,jj) ; ikus = mbku_crs(ji,jj) 515 ! 516 ! ! mid-depth density anomalie (up-slope minus down-slope) 517 zt = 0.5 * ( ztb (ji,jj) + ztb (ji+1,jj) ) ! mid slope depth of T, S, and depth 518 zs = 0.5 * ( zsb (ji,jj) + zsb (ji+1,jj) ) - 35.0 519 zh = 0.5 * ( zdep(ji,jj) + zdep(ji+1,jj) ) 520 zgdrho = fsbeta( zt, zs, zh ) & 521 & * ( fsalbt( zt, zs, zh ) * ( ztb(iid,jj) - ztb(iis,jj) ) & 522 & - ( zsb(iid,jj) - zsb(iis,jj) ) ) * umask_crs(ji,jj,1) 523 zgdrho = MAX( 0.e0, zgdrho ) ! only if shelf is denser than deep 524 ! 525 ! ! bbl transport (down-slope direction) 526 utr_bbl_crs(ji,jj) = e2u_crs(ji,jj) * e3u_bbl_0(ji,jj) * zgbbl * zgdrho * REAL( mgrhu(ji,jj) ) 527 ! 528 ! ! j-direction 529 ! down-slope T-point j/k-index (deep) & of the up -slope T-point j/k-index (shelf) 530 ijd = jj + MAX( 0, mgrhv(ji,jj) ) ; ijs = jj + 1 - MAX( 0, mgrhv(ji,jj) ) 531 ikvd = mbkv_d(ji,jj) ; ikvs = mbkv_crs(ji,jj) 532 ! 533 ! ! mid-depth density anomalie (up-slope minus down-slope) 534 zt = 0.5 * ( ztb (ji,jj) + ztb (ji,jj+1) ) ! mid slope depth of T, S, and depth 535 zs = 0.5 * ( zsb (ji,jj) + zsb (ji,jj+1) ) - 35.0 536 zh = 0.5 * ( zdep(ji,jj) + zdep(ji,jj+1) ) 537 zgdrho = fsbeta( zt, zs, zh ) & 538 & * ( fsalbt( zt, zs, zh ) * ( ztb(ji,ijd) - ztb(ji,ijs) ) & 539 & - ( zsb(ji,ijd) - zsb(ji,ijs) ) ) * vmask_crs(ji,jj,1) 540 zgdrho = MAX( 0.e0, zgdrho ) ! only if shelf is denser than deep 541 ! 542 ! ! bbl transport (down-slope direction) 543 vtr_bbl_crs(ji,jj) = e1v_crs(ji,jj) * e3v_bbl_0(ji,jj) * zgbbl * zgdrho * REAL( mgrhv(ji,jj) ) 544 END DO 545 END DO 546 END SELECT 547 ! 548 ENDIF 549 ! 550 CALL wrk_dealloc( jpi, jpj, zub, zvb, ztb, zsb, zdep ) 551 ! 552 IF( nn_timing == 1 ) CALL timing_stop( 'bbl') 553 ! 554 END SUBROUTINE bbl_crs 555 556 557 SUBROUTINE tra_bbl_init_crs 558 !!---------------------------------------------------------------------- 559 !! *** ROUTINE tra_bbl_init *** 560 !! 561 !! ** Purpose : Initialization for the bottom boundary layer scheme. 562 !! 563 !! ** Method : Read the nambbl namelist and check the parameters 564 !! called by nemo_init at the first timestep (kit000) 565 !!---------------------------------------------------------------------- 566 INTEGER :: ji, jj ! dummy loop indices 567 INTEGER :: ii0, ii1, ij0, ij1 ! temporary integer 568 REAL(wp), POINTER, DIMENSION(:,:) :: zmbk 569 !! 570 NAMELIST/nambbl/ nn_bbl_ldf, nn_bbl_adv, rn_ahtbbl, rn_gambbl 571 !!---------------------------------------------------------------------- 572 ! 573 IF( nn_timing == 1 ) CALL timing_start( 'tra_bbl_init') 574 ! 575 CALL wrk_alloc( jpi, jpj, zmbk ) 576 ! 577 578 REWIND ( numnam ) !* Read Namelist nambbl : bottom boundary layer scheme 579 READ ( numnam, nambbl ) 580 ! 581 l_bbl = .TRUE. !* flag to compute bbl coef and transport 582 ! 583 IF(lwp) THEN !* Parameter control and print 584 WRITE(numout,*) 585 WRITE(numout,*) 'tra_bbl_init_crs : bottom boundary layer initialisation' 586 WRITE(numout,*) '~~~~~~~~~~~~' 587 WRITE(numout,*) ' Namelist nambbl : set bbl parameters' 588 WRITE(numout,*) ' diffusive bbl (=1) or not (=0) nn_bbl_ldf = ', nn_bbl_ldf 589 WRITE(numout,*) ' advective bbl (=1/2) or not (=0) nn_bbl_adv = ', nn_bbl_adv 590 WRITE(numout,*) ' diffusive bbl coefficient rn_ahtbbl = ', rn_ahtbbl, ' m2/s' 591 WRITE(numout,*) ' advective bbl coefficient rn_gambbl = ', rn_gambbl, ' s' 592 ENDIF 593 594 ! ! allocate trabbl arrays 595 IF( tra_bbl_alloc_crs() /= 0 ) CALL ctl_stop( 'STOP', 'tra_bbl_init : unable to allocate arrays' ) 596 597 IF( nn_bbl_adv == 1 ) WRITE(numout,*) ' * Advective BBL using upper velocity' 598 IF( nn_bbl_adv == 2 ) WRITE(numout,*) ' * Advective BBL using velocity = F( delta rho)' 599 600 IF( nn_eos /= 0 ) CALL ctl_stop ( ' bbl parameterisation requires eos = 0. We stop.' ) 601 602 603 ! !* inverse of surface of T-cells 604 r1_e1e2t(:,:) = 1._wp / e1e2w_msk(:,:,1) 605 606 ! !* vertical index of "deep" bottom u- and v-points 607 DO jj = 1, jpjm1 ! (the "shelf" bottom k-indices are mbku and mbkv) 608 DO ji = 1, jpim1 609 mbku_d(ji,jj) = MAX( mbkt_crs(ji+1,jj ) , mbkt_crs(ji,jj) ) ! >= 1 as mbkt=1 over land 610 mbkv_d(ji,jj) = MAX( mbkt_crs(ji ,jj+1) , mbkt_crs(ji,jj) ) 611 END DO 612 END DO 613 ! converte into REAL to use lbc_lnk ; impose a min value of 1 as a zero can be set in lbclnk 614 zmbk(:,:) = REAL( mbku_d(:,:), wp ) ; CALL crs_lbc_lnk(zmbk,'U',1.) ; mbku_d(:,:) = MAX( INT( zmbk(:,:) ), 1 ) 615 zmbk(:,:) = REAL( mbkv_d(:,:), wp ) ; CALL crs_lbc_lnk(zmbk,'V',1.) ; mbkv_d(:,:) = MAX( INT( zmbk(:,:) ), 1 ) 616 617 !* sign of grad(H) at u- and v-points 618 mgrhu(jpi,:) = 0. ; mgrhu(:,jpj) = 0. ; mgrhv(jpi,:) = 0. ; mgrhv(:,jpj) = 0. 619 DO jj = 1, jpjm1 620 DO ji = 1, jpim1 621 mgrhu(ji,jj) = INT( SIGN( 1.e0, gdept_crs(ji+1,jj,mbkt_crs(ji+1,jj)) - gdept_crs(ji,jj,mbkt_crs(ji,jj)) ) ) !! cc car non vvl 622 mgrhv(ji,jj) = INT( SIGN( 1.e0, gdept_crs(ji,jj+1,mbkt_crs(ji,jj+1)) - gdept_crs(ji,jj,mbkt_crs(ji,jj)) ) ) 623 END DO 624 END DO 625 626 DO jj = 1, jpjm1 !* bbl thickness at u- (v-) point 627 DO ji = 1, jpim1 ! minimum of top & bottom e3u_0 (e3v_0) 628 e3u_bbl_0(ji,jj) = MIN( e3u_crs(ji,jj,mbkt_crs(ji+1,jj )), e3u_crs(ji,jj,mbkt_crs(ji,jj)) ) 629 e3v_bbl_0(ji,jj) = MIN( e3v_crs(ji,jj,mbkt_crs(ji ,jj+1)), e3v_crs(ji,jj,mbkt_crs(ji,jj)) ) !! cc car non vvl 630 END DO 631 END DO 632 CALL crs_lbc_lnk( e3u_bbl_0, 'U', 1. ) ; CALL crs_lbc_lnk( e3v_bbl_0, 'V', 1. ) ! lateral boundary conditions 633 634 ! !* masked diffusive flux coefficients 635 ahu_bbl_0(:,:) = rn_ahtbbl_crs * e2u_crs(:,:) * e3u_bbl_0(:,:) / e1u_crs(:,:) * umask_crs(:,:,1) 636 ahv_bbl_0(:,:) = rn_ahtbbl_crs * e1v_crs(:,:) * e3v_bbl_0(:,:) / e2v_crs(:,:) * vmask_crs(:,:,1) 637 638 !! cc reflexion for the coarsening 639 640 ! IF( cp_cfg == "orca" ) THEN !* ORCA configuration : regional enhancement of ah_bbl 641 ! 642 ! SELECT CASE ( jp_cfg ) 643 ! CASE ( 2 ) ! ORCA_R2 644 ! ij0 = 102 ; ij1 = 102 ! Gibraltar enhancement of BBL 645 ! ii0 = 139 ; ii1 = 140 646 ! ahu_bbl_0(mi0_crs(ii0):mi1_crs(ii1),mj0_crs(ij0):mj1_crs(ij1)) & 647 ! & = 4.e0*ahu_bbl_0(mi0_crs(ii0):mi1_crs(ii1),mj0_crs(ij0):mj1_crs(ij1)) 648 ! ahv_bbl_0(mi0_crs(ii0):mi1_crs(ii1),mj0_crs(ij0):mj1_crs(ij1)) & 649 ! & = 4.e0*ahv_bbl_0(mi0_crs(ii0):mi1_crs(ii1),mj0_crs(ij0):mj1_crs(ij1)) 650 ! ! 651 ! ij0 = 88 ; ij1 = 88 ! Red Sea enhancement of BBL 652 ! ii0 = 161 ; ii1 = 162 653 ! ahu_bbl_0(mi0_crs(ii0):mi1_crs(ii1),mj0_crs(ij0):mj1(ij1)) & 654 ! & = 10.e0*ahu_bbl_0(mi0_crs(ii0):mi1_crs(ii1),mj0_crs(ij0):mj1_crs(ij1)) 655 ! ahv_bbl_0(mi0_crs(ii0):mi1_crs(ii1),mj0_crs(ij0):mj1(ij1)) & 656 ! & = 10.e0*ahv_bbl_0(mi0_crs(ii0):mi1_crs(ii1),mj0_crs(ij0):mj1_crs(ij1)) 657 ! ! 658 ! CASE ( 4 ) ! ORCA_R4 659 ! ij0 = 52 ; ij1 = 52 ! Gibraltar enhancement of BBL !!cc 660 ! ii0 = 70 ; ii1 = 71 661 ! ahu_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1)) = 4.e0*ahu_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1)) 662 ! ahv_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1)) = 4.e0*ahv_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1)) 663 ! END SELECT 664 ! ! 665 ! ENDIF 666 ! 667 CALL wrk_dealloc( jpi, jpj, zmbk ) 668 ! 669 IF( nn_timing == 1 ) CALL timing_stop( 'tra_bbl_init') 670 ! 671 END SUBROUTINE tra_bbl_init_crs 672 673 #else 674 !!---------------------------------------------------------------------- 675 !! Dummy module : No bottom boundary layer scheme 676 !!---------------------------------------------------------------------- 677 LOGICAL, PUBLIC, PARAMETER :: lk_trabbl_crs = .FALSE. !: bbl flag 678 CONTAINS 679 SUBROUTINE tra_bbl_init_crs ! Dummy routine 680 END SUBROUTINE tra_bbl_init_crs 681 SUBROUTINE tra_bbl_crs( kt ) ! Dummy routine 682 WRITE(*,*) 'tra_bbl: You should not have seen this print! error?', kt 683 END SUBROUTINE tra_bbl_crs 684 #endif 685 686 !!====================================================================== 687 END MODULE trabbl_crs
Note: See TracChangeset
for help on using the changeset viewer.
