- Timestamp:
- 2015-07-09T12:14:37+02:00 (9 years ago)
- Location:
- branches/UKMO/dev_r5107_hadgem3_cplseq/NEMOGCM/NEMO/OPA_SRC/TRA
- Files:
-
- 22 edited
-
eosbn2.F90 (modified) (7 diffs)
-
traadv.F90 (modified) (4 diffs)
-
traadv_cen2.F90 (modified) (1 diff)
-
traadv_muscl.F90 (modified) (2 diffs)
-
traadv_muscl2.F90 (modified) (1 diff)
-
traadv_qck.F90 (modified) (1 diff)
-
traadv_tvd.F90 (modified) (10 diffs)
-
traadv_ubs.F90 (modified) (1 diff)
-
trabbc.F90 (modified) (5 diffs)
-
traldf.F90 (modified) (1 diff)
-
traldf_bilap.F90 (modified) (4 diffs)
-
traldf_bilapg.F90 (modified) (1 diff)
-
traldf_iso.F90 (modified) (6 diffs)
-
traldf_iso_grif.F90 (modified) (3 diffs)
-
traldf_lap.F90 (modified) (3 diffs)
-
tranpc.F90 (modified) (12 diffs)
-
tranxt.F90 (modified) (6 diffs)
-
traqsr.F90 (modified) (12 diffs)
-
trasbc.F90 (modified) (3 diffs)
-
trazdf.F90 (modified) (1 diff)
-
trazdf_imp.F90 (modified) (4 diffs)
-
zpshde.F90 (modified) (14 diffs)
Legend:
- Unmodified
- Added
- Removed
-
branches/UKMO/dev_r5107_hadgem3_cplseq/NEMOGCM/NEMO/OPA_SRC/TRA/eosbn2.F90
r5477 r5572 47 47 USE lbclnk ! ocean lateral boundary conditions 48 48 USE timing ! Timing 49 USE stopar ! Stochastic T/S fluctuations 50 USE stopts ! Stochastic T/S fluctuations 49 51 50 52 IMPLICIT NONE … … 72 74 PUBLIC eos_init ! called by istate module 73 75 74 ! 75 INTEGER , PUBLIC :: nn_eos = 0 !:= 0/1/2 type of eq. of state and Brunt-Vaisala frequ.76 LOGICAL , PUBLIC :: ln_useCT = .FALSE.! determine if eos_pt_from_ct is used to compute sst_m76 ! !!* Namelist (nameos) * 77 INTEGER , PUBLIC :: nn_eos ! = 0/1/2 type of eq. of state and Brunt-Vaisala frequ. 78 LOGICAL , PUBLIC :: ln_useCT ! determine if eos_pt_from_ct is used to compute sst_m 77 79 78 80 ! !!! simplified eos coefficients … … 313 315 REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ) :: pdep ! depth [m] 314 316 ! 315 INTEGER :: ji, jj, jk ! dummy loop indices 316 REAL(wp) :: zt , zh , zs , ztm ! local scalars 317 REAL(wp) :: zn , zn0, zn1, zn2, zn3 ! - - 317 INTEGER :: ji, jj, jk, jsmp ! dummy loop indices 318 INTEGER :: jdof 319 REAL(wp) :: zt , zh , zstemp, zs , ztm ! local scalars 320 REAL(wp) :: zn , zn0, zn1, zn2, zn3 ! - - 321 REAL(wp), DIMENSION(:), ALLOCATABLE :: zn0_sto, zn_sto, zsign ! local vectors 318 322 !!---------------------------------------------------------------------- 319 323 ! … … 324 328 CASE( -1, 0 ) !== polynomial TEOS-10 / EOS-80 ==! 325 329 ! 326 DO jk = 1, jpkm1 327 DO jj = 1, jpj 328 DO ji = 1, jpi 329 ! 330 zh = pdep(ji,jj,jk) * r1_Z0 ! depth 331 zt = pts (ji,jj,jk,jp_tem) * r1_T0 ! temperature 332 zs = SQRT( ABS( pts(ji,jj,jk,jp_sal) + rdeltaS ) * r1_S0 ) ! square root salinity 333 ztm = tmask(ji,jj,jk) ! tmask 334 ! 335 zn3 = EOS013*zt & 336 & + EOS103*zs+EOS003 337 ! 338 zn2 = (EOS022*zt & 339 & + EOS112*zs+EOS012)*zt & 340 & + (EOS202*zs+EOS102)*zs+EOS002 341 ! 342 zn1 = (((EOS041*zt & 343 & + EOS131*zs+EOS031)*zt & 344 & + (EOS221*zs+EOS121)*zs+EOS021)*zt & 345 & + ((EOS311*zs+EOS211)*zs+EOS111)*zs+EOS011)*zt & 346 & + (((EOS401*zs+EOS301)*zs+EOS201)*zs+EOS101)*zs+EOS001 347 ! 348 zn0 = (((((EOS060*zt & 349 & + EOS150*zs+EOS050)*zt & 350 & + (EOS240*zs+EOS140)*zs+EOS040)*zt & 351 & + ((EOS330*zs+EOS230)*zs+EOS130)*zs+EOS030)*zt & 352 & + (((EOS420*zs+EOS320)*zs+EOS220)*zs+EOS120)*zs+EOS020)*zt & 353 & + ((((EOS510*zs+EOS410)*zs+EOS310)*zs+EOS210)*zs+EOS110)*zs+EOS010)*zt & 354 & + (((((EOS600*zs+EOS500)*zs+EOS400)*zs+EOS300)*zs+EOS200)*zs+EOS100)*zs+EOS000 355 ! 356 zn = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0 357 ! 358 prhop(ji,jj,jk) = zn0 * ztm ! potential density referenced at the surface 359 ! 360 prd(ji,jj,jk) = ( zn * r1_rau0 - 1._wp ) * ztm ! density anomaly (masked) 330 ! Stochastic equation of state 331 IF ( ln_sto_eos ) THEN 332 ALLOCATE(zn0_sto(1:2*nn_sto_eos)) 333 ALLOCATE(zn_sto(1:2*nn_sto_eos)) 334 ALLOCATE(zsign(1:2*nn_sto_eos)) 335 DO jsmp = 1, 2*nn_sto_eos, 2 336 zsign(jsmp) = 1._wp 337 zsign(jsmp+1) = -1._wp 338 END DO 339 ! 340 DO jk = 1, jpkm1 341 DO jj = 1, jpj 342 DO ji = 1, jpi 343 ! 344 ! compute density (2*nn_sto_eos) times: 345 ! (1) for t+dt, s+ds (with the random TS fluctutation computed in sto_pts) 346 ! (2) for t-dt, s-ds (with the opposite fluctuation) 347 DO jsmp = 1, nn_sto_eos*2 348 jdof = (jsmp + 1) / 2 349 zh = pdep(ji,jj,jk) * r1_Z0 ! depth 350 zt = (pts (ji,jj,jk,jp_tem) + pts_ran(ji,jj,jk,jp_tem,jdof) * zsign(jsmp)) * r1_T0 ! temperature 351 zstemp = pts (ji,jj,jk,jp_sal) + pts_ran(ji,jj,jk,jp_sal,jdof) * zsign(jsmp) 352 zs = SQRT( ABS( zstemp + rdeltaS ) * r1_S0 ) ! square root salinity 353 ztm = tmask(ji,jj,jk) ! tmask 354 ! 355 zn3 = EOS013*zt & 356 & + EOS103*zs+EOS003 357 ! 358 zn2 = (EOS022*zt & 359 & + EOS112*zs+EOS012)*zt & 360 & + (EOS202*zs+EOS102)*zs+EOS002 361 ! 362 zn1 = (((EOS041*zt & 363 & + EOS131*zs+EOS031)*zt & 364 & + (EOS221*zs+EOS121)*zs+EOS021)*zt & 365 & + ((EOS311*zs+EOS211)*zs+EOS111)*zs+EOS011)*zt & 366 & + (((EOS401*zs+EOS301)*zs+EOS201)*zs+EOS101)*zs+EOS001 367 ! 368 zn0_sto(jsmp) = (((((EOS060*zt & 369 & + EOS150*zs+EOS050)*zt & 370 & + (EOS240*zs+EOS140)*zs+EOS040)*zt & 371 & + ((EOS330*zs+EOS230)*zs+EOS130)*zs+EOS030)*zt & 372 & + (((EOS420*zs+EOS320)*zs+EOS220)*zs+EOS120)*zs+EOS020)*zt & 373 & + ((((EOS510*zs+EOS410)*zs+EOS310)*zs+EOS210)*zs+EOS110)*zs+EOS010)*zt & 374 & + (((((EOS600*zs+EOS500)*zs+EOS400)*zs+EOS300)*zs+EOS200)*zs+EOS100)*zs+EOS000 375 ! 376 zn_sto(jsmp) = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0_sto(jsmp) 377 END DO 378 ! 379 ! compute stochastic density as the mean of the (2*nn_sto_eos) densities 380 prhop(ji,jj,jk) = 0._wp ; prd(ji,jj,jk) = 0._wp 381 DO jsmp = 1, nn_sto_eos*2 382 prhop(ji,jj,jk) = prhop(ji,jj,jk) + zn0_sto(jsmp) ! potential density referenced at the surface 383 ! 384 prd(ji,jj,jk) = prd(ji,jj,jk) + ( zn_sto(jsmp) * r1_rau0 - 1._wp ) ! density anomaly (masked) 385 END DO 386 prhop(ji,jj,jk) = 0.5_wp * prhop(ji,jj,jk) * ztm / nn_sto_eos 387 prd (ji,jj,jk) = 0.5_wp * prd (ji,jj,jk) * ztm / nn_sto_eos 388 END DO 361 389 END DO 362 390 END DO 363 END DO 364 ! 391 DEALLOCATE(zn0_sto,zn_sto,zsign) 392 ! Non-stochastic equation of state 393 ELSE 394 DO jk = 1, jpkm1 395 DO jj = 1, jpj 396 DO ji = 1, jpi 397 ! 398 zh = pdep(ji,jj,jk) * r1_Z0 ! depth 399 zt = pts (ji,jj,jk,jp_tem) * r1_T0 ! temperature 400 zs = SQRT( ABS( pts(ji,jj,jk,jp_sal) + rdeltaS ) * r1_S0 ) ! square root salinity 401 ztm = tmask(ji,jj,jk) ! tmask 402 ! 403 zn3 = EOS013*zt & 404 & + EOS103*zs+EOS003 405 ! 406 zn2 = (EOS022*zt & 407 & + EOS112*zs+EOS012)*zt & 408 & + (EOS202*zs+EOS102)*zs+EOS002 409 ! 410 zn1 = (((EOS041*zt & 411 & + EOS131*zs+EOS031)*zt & 412 & + (EOS221*zs+EOS121)*zs+EOS021)*zt & 413 & + ((EOS311*zs+EOS211)*zs+EOS111)*zs+EOS011)*zt & 414 & + (((EOS401*zs+EOS301)*zs+EOS201)*zs+EOS101)*zs+EOS001 415 ! 416 zn0 = (((((EOS060*zt & 417 & + EOS150*zs+EOS050)*zt & 418 & + (EOS240*zs+EOS140)*zs+EOS040)*zt & 419 & + ((EOS330*zs+EOS230)*zs+EOS130)*zs+EOS030)*zt & 420 & + (((EOS420*zs+EOS320)*zs+EOS220)*zs+EOS120)*zs+EOS020)*zt & 421 & + ((((EOS510*zs+EOS410)*zs+EOS310)*zs+EOS210)*zs+EOS110)*zs+EOS010)*zt & 422 & + (((((EOS600*zs+EOS500)*zs+EOS400)*zs+EOS300)*zs+EOS200)*zs+EOS100)*zs+EOS000 423 ! 424 zn = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0 425 ! 426 prhop(ji,jj,jk) = zn0 * ztm ! potential density referenced at the surface 427 ! 428 prd(ji,jj,jk) = ( zn * r1_rau0 - 1._wp ) * ztm ! density anomaly (masked) 429 END DO 430 END DO 431 END DO 432 ENDIF 433 365 434 CASE( 1 ) !== simplified EOS ==! 366 435 ! … … 1183 1252 WRITE(numout,*) ' model uses Conservative Temperature' 1184 1253 WRITE(numout,*) ' Important: model must be initialized with CT and SA fields' 1254 ELSE 1255 WRITE(numout,*) ' model does not use Conservative Temperature' 1185 1256 ENDIF 1186 1257 ENDIF … … 1589 1660 END SELECT 1590 1661 ! 1662 rau0_rcp = rau0 * rcp 1591 1663 r1_rau0 = 1._wp / rau0 1592 1664 r1_rcp = 1._wp / rcp 1593 r1_rau0_rcp = 1._wp / ( rau0 * rcp )1665 r1_rau0_rcp = 1._wp / rau0_rcp 1594 1666 ! 1595 1667 IF(lwp) WRITE(numout,*) … … 1597 1669 IF(lwp) WRITE(numout,*) ' 1. / rau0 r1_rau0 = ', r1_rau0, ' m^3/kg' 1598 1670 IF(lwp) WRITE(numout,*) ' ocean specific heat rcp = ', rcp , ' J/Kelvin' 1671 IF(lwp) WRITE(numout,*) ' rau0 * rcp rau0_rcp = ', rau0_rcp 1599 1672 IF(lwp) WRITE(numout,*) ' 1. / ( rau0 * rcp ) r1_rau0_rcp = ', r1_rau0_rcp 1600 1673 ! -
branches/UKMO/dev_r5107_hadgem3_cplseq/NEMOGCM/NEMO/OPA_SRC/TRA/traadv.F90
r5477 r5572 26 26 USE cla ! cross land advection (cla_traadv routine) 27 27 USE ldftra_oce ! lateral diffusion coefficient on tracers 28 ! 28 29 USE in_out_manager ! I/O manager 29 30 USE iom ! I/O module … … 33 34 USE timing ! Timing 34 35 USE sbc_oce 36 USE diaptr ! Poleward heat transport 35 37 36 38 … … 111 113 ! 112 114 IF( ln_mle ) CALL tra_adv_mle( kt, nit000, zun, zvn, zwn, 'TRA' ) ! add the mle transport (if necessary) 115 ! 113 116 CALL iom_put( "uocetr_eff", zun ) ! output effective transport 114 117 CALL iom_put( "vocetr_eff", zvn ) 115 118 CALL iom_put( "wocetr_eff", zwn ) 116 119 ! 120 IF( ln_diaptr ) CALL dia_ptr( zvn ) ! diagnose the effective MSF 121 ! 122 117 123 SELECT CASE ( nadv ) !== compute advection trend and add it to general trend ==! 118 CASE ( 1 ) ; CALL tra_adv_cen2 ( kt, nit000, 'TRA', zun, zvn, zwn, tsb, tsn, tsa, jpts ) ! 2nd order centered119 CASE ( 2 ) ; CALL tra_adv_tvd ( kt, nit000, 'TRA', r2dtra, zun, zvn, zwn, tsb, tsn, tsa, jpts ) ! TVD120 CASE ( 3 ) ; CALL tra_adv_muscl ( kt, nit000, 'TRA', r2dtra, zun, zvn, zwn, tsb, tsa, jpts, ln_traadv_msc_ups ) ! MUSCL121 CASE ( 4 ) ; CALL tra_adv_muscl2 ( kt, nit000, 'TRA', r2dtra, zun, zvn, zwn, tsb, tsn, tsa, jpts ) ! MUSCL2122 CASE ( 5 ) ; CALL tra_adv_ubs ( kt, nit000, 'TRA', r2dtra, zun, zvn, zwn, tsb, tsn, tsa, jpts ) ! UBS123 CASE ( 6 ) ; CALL tra_adv_qck ( kt, nit000, 'TRA', r2dtra, zun, zvn, zwn, tsb, tsn, tsa, jpts ) ! QUICKEST124 CASE ( 7 ) ; CALL tra_adv_tvd_zts( kt, nit000, 'TRA', r2dtra, zun, zvn, zwn, tsb, tsn, tsa, jpts ) ! TVD ZTS124 CASE ( 1 ) ; CALL tra_adv_cen2 ( kt, nit000, 'TRA', zun, zvn, zwn, tsb, tsn, tsa, jpts ) ! 2nd order centered 125 CASE ( 2 ) ; CALL tra_adv_tvd ( kt, nit000, 'TRA', r2dtra, zun, zvn, zwn, tsb, tsn, tsa, jpts ) ! TVD 126 CASE ( 3 ) ; CALL tra_adv_muscl ( kt, nit000, 'TRA', r2dtra, zun, zvn, zwn, tsb, tsa, jpts, ln_traadv_msc_ups ) ! MUSCL 127 CASE ( 4 ) ; CALL tra_adv_muscl2 ( kt, nit000, 'TRA', r2dtra, zun, zvn, zwn, tsb, tsn, tsa, jpts ) ! MUSCL2 128 CASE ( 5 ) ; CALL tra_adv_ubs ( kt, nit000, 'TRA', r2dtra, zun, zvn, zwn, tsb, tsn, tsa, jpts ) ! UBS 129 CASE ( 6 ) ; CALL tra_adv_qck ( kt, nit000, 'TRA', r2dtra, zun, zvn, zwn, tsb, tsn, tsa, jpts ) ! QUICKEST 130 CASE ( 7 ) ; CALL tra_adv_tvd_zts( kt, nit000, 'TRA', r2dtra, zun, zvn, zwn, tsb, tsn, tsa, jpts ) ! TVD ZTS 125 131 ! 126 132 CASE (-1 ) !== esopa: test all possibility with control print ==! … … 206 212 IF( lk_esopa ) ioptio = 1 207 213 208 IF( ( ln_traadv_muscl .OR. ln_traadv_muscl2 .OR. ln_traadv_ubs .OR. ln_traadv_qck ) .AND. nn_isf .NE. 0 )&209 &CALL ctl_stop( 'Only traadv_cen2 and traadv_tvd is compatible with ice shelf cavity')214 IF( ( ln_traadv_muscl .OR. ln_traadv_muscl2 .OR. ln_traadv_ubs .OR. ln_traadv_qck .OR. ln_traadv_tvd_zts ) & 215 .AND. ln_isfcav ) CALL ctl_stop( 'Only traadv_cen2 and traadv_tvd is compatible with ice shelf cavity') 210 216 211 217 IF( ioptio /= 1 ) CALL ctl_stop( 'Choose ONE advection scheme in namelist namtra_adv' ) -
branches/UKMO/dev_r5107_hadgem3_cplseq/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_cen2.F90
r5477 r5572 279 279 END IF 280 280 ! ! "Poleward" heat and salt transports (contribution of upstream fluxes) 281 IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nn_fptr ) == 0 )) THEN282 IF( jn == jp_tem ) htr_adv(:) = ptr_ vj( zwy(:,:,:) )283 IF( jn == jp_sal ) str_adv(:) = ptr_ vj( zwy(:,:,:) )281 IF( cdtype == 'TRA' .AND. ln_diaptr ) THEN 282 IF( jn == jp_tem ) htr_adv(:) = ptr_sj( zwy(:,:,:) ) 283 IF( jn == jp_sal ) str_adv(:) = ptr_sj( zwy(:,:,:) ) 284 284 ENDIF 285 285 ! -
branches/UKMO/dev_r5107_hadgem3_cplseq/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_muscl.F90
r5477 r5572 21 21 USE trdtra ! tracers trends manager 22 22 USE dynspg_oce ! choice/control of key cpp for surface pressure gradient 23 USE sbcrnf 23 USE sbcrnf ! river runoffs 24 24 USE diaptr ! poleward transport diagnostics 25 25 ! … … 219 219 END IF 220 220 ! ! "Poleward" heat and salt transports (contribution of upstream fluxes) 221 IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nn_fptr ) == 0 )) THEN222 IF( jn == jp_tem ) htr_adv(:) = ptr_ vj( zwy(:,:,:) )223 IF( jn == jp_sal ) str_adv(:) = ptr_ vj( zwy(:,:,:) )221 IF( cdtype == 'TRA' .AND. ln_diaptr ) THEN 222 IF( jn == jp_tem ) htr_adv(:) = ptr_sj( zwy(:,:,:) ) 223 IF( jn == jp_sal ) str_adv(:) = ptr_sj( zwy(:,:,:) ) 224 224 ENDIF 225 225 -
branches/UKMO/dev_r5107_hadgem3_cplseq/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_muscl2.F90
r5477 r5572 200 200 201 201 ! ! "Poleward" heat and salt transports (contribution of upstream fluxes) 202 IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nn_fptr ) == 0 )) THEN203 IF( jn == jp_tem ) htr_adv(:) = ptr_ vj( zwy(:,:,:) )204 IF( jn == jp_sal ) str_adv(:) = ptr_ vj( zwy(:,:,:) )202 IF( cdtype == 'TRA' .AND. ln_diaptr ) THEN 203 IF( jn == jp_tem ) htr_adv(:) = ptr_sj( zwy(:,:,:) ) 204 IF( jn == jp_sal ) str_adv(:) = ptr_sj( zwy(:,:,:) ) 205 205 ENDIF 206 206 -
branches/UKMO/dev_r5107_hadgem3_cplseq/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_qck.F90
r5477 r5572 355 355 IF( l_trd ) CALL trd_tra( kt, cdtype, jn, jptra_yad, zwy, pvn, ptn(:,:,:,jn) ) 356 356 ! ! "Poleward" heat and salt transports (contribution of upstream fluxes) 357 IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nn_fptr ) == 0 )) THEN358 IF( jn == jp_tem ) htr_adv(:) = ptr_ vj( zwy(:,:,:) )359 IF( jn == jp_sal ) str_adv(:) = ptr_ vj( zwy(:,:,:) )357 IF( cdtype == 'TRA' .AND. ln_diaptr ) THEN 358 IF( jn == jp_tem ) htr_adv(:) = ptr_sj( zwy(:,:,:) ) 359 IF( jn == jp_sal ) str_adv(:) = ptr_sj( zwy(:,:,:) ) 360 360 ENDIF 361 361 ! -
branches/UKMO/dev_r5107_hadgem3_cplseq/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_tvd.F90
r5477 r5572 106 106 ENDIF 107 107 ! 108 zwi(:,:,:) = 0.e0 ; zwz(:,:,:) = 0.e0108 zwi(:,:,:) = 0.e0 ; 109 109 ! 110 110 ! ! =========== 111 111 DO jn = 1, kjpt ! tracer loop 112 112 ! ! =========== 113 ! 1. Bottom value : flux set to zero113 ! 1. Bottom and k=1 value : flux set to zero 114 114 ! ---------------------------------- 115 115 zwx(:,:,jpk) = 0.e0 ; zwz(:,:,jpk) = 0.e0 116 116 zwy(:,:,jpk) = 0.e0 ; zwi(:,:,jpk) = 0.e0 117 117 118 zwz(:,:,1 ) = 0._wp 118 119 ! 2. upstream advection with initial mass fluxes & intermediate update 119 120 ! -------------------------------------------------------------------- … … 134 135 135 136 ! upstream tracer flux in the k direction 137 ! Interior value 138 DO jk = 2, jpkm1 139 DO jj = 1, jpj 140 DO ji = 1, jpi 141 zfp_wk = pwn(ji,jj,jk) + ABS( pwn(ji,jj,jk) ) 142 zfm_wk = pwn(ji,jj,jk) - ABS( pwn(ji,jj,jk) ) 143 zwz(ji,jj,jk) = 0.5 * ( zfp_wk * ptb(ji,jj,jk,jn) + zfm_wk * ptb(ji,jj,jk-1,jn) ) * wmask(ji,jj,jk) 144 END DO 145 END DO 146 END DO 136 147 ! Surface value 137 148 IF( lk_vvl ) THEN 138 DO jj = 1, jpj 139 DO ji = 1, jpi 140 zwz(ji,jj, mikt(ji,jj) ) = 0.e0 ! volume variable 141 END DO 142 END DO 149 IF ( ln_isfcav ) THEN 150 DO jj = 1, jpj 151 DO ji = 1, jpi 152 zwz(ji,jj, mikt(ji,jj) ) = 0.e0 ! volume variable 153 END DO 154 END DO 155 ELSE 156 zwz(:,:,1) = 0.e0 ! volume variable 157 END IF 143 158 ELSE 144 DO jj = 1, jpj 145 DO ji = 1, jpi 146 zwz(ji,jj, mikt(ji,jj) ) = pwn(ji,jj,mikt(ji,jj)) * ptb(ji,jj,mikt(ji,jj),jn) ! linear free surface 147 END DO 148 END DO 159 IF ( ln_isfcav ) THEN 160 DO jj = 1, jpj 161 DO ji = 1, jpi 162 zwz(ji,jj, mikt(ji,jj) ) = pwn(ji,jj,mikt(ji,jj)) * ptb(ji,jj,mikt(ji,jj),jn) ! linear free surface 163 END DO 164 END DO 165 ELSE 166 zwz(:,:,1) = pwn(:,:,1) * ptb(:,:,1,jn) ! linear free surface 167 END IF 149 168 ENDIF 150 ! Interior value151 DO jj = 1, jpj152 DO ji = 1, jpi153 DO jk = mikt(ji,jj)+1, jpkm1154 zfp_wk = pwn(ji,jj,jk) + ABS( pwn(ji,jj,jk) )155 zfm_wk = pwn(ji,jj,jk) - ABS( pwn(ji,jj,jk) )156 zwz(ji,jj,jk) = 0.5 * ( zfp_wk * ptb(ji,jj,jk,jn) + zfm_wk * ptb(ji,jj,jk-1,jn) )157 END DO158 END DO159 END DO160 169 161 170 ! total advective trend … … 184 193 END IF 185 194 ! ! "Poleward" heat and salt transports (contribution of upstream fluxes) 186 IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nn_fptr ) == 0 )) THEN187 IF( jn == jp_tem ) htr_adv(:) = ptr_ vj( zwy(:,:,:) )188 IF( jn == jp_sal ) str_adv(:) = ptr_ vj( zwy(:,:,:) )195 IF( cdtype == 'TRA' .AND. ln_diaptr ) THEN 196 IF( jn == jp_tem ) htr_adv(:) = ptr_sj( zwy(:,:,:) ) 197 IF( jn == jp_sal ) str_adv(:) = ptr_sj( zwy(:,:,:) ) 189 198 ENDIF 190 199 … … 202 211 203 212 ! antidiffusive flux on k 204 zwz(:,:,1) = 0.e0 ! Surface value 205 ! 206 DO jj = 1, jpj 207 DO ji = 1, jpi 208 ik=mikt(ji,jj) 209 ! surface value 210 zwz(ji,jj,1:ik) = 0.e0 211 ! Interior value 212 DO jk = mikt(ji,jj)+1, jpkm1 213 ! Interior value 214 DO jk = 2, jpkm1 215 DO jj = 1, jpj 216 DO ji = 1, jpi 213 217 zwz(ji,jj,jk) = 0.5 * pwn(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji,jj,jk-1,jn) ) - zwz(ji,jj,jk) 214 218 END DO 215 219 END DO 216 220 END DO 221 ! surface value 222 IF ( ln_isfcav ) THEN 223 DO jj = 1, jpj 224 DO ji = 1, jpi 225 zwz(ji,jj,mikt(ji,jj)) = 0.e0 226 END DO 227 END DO 228 ELSE 229 zwz(:,:,1) = 0.e0 230 END IF 217 231 CALL lbc_lnk( zwx, 'U', -1. ) ; CALL lbc_lnk( zwy, 'V', -1. ) ! Lateral bondary conditions 218 232 CALL lbc_lnk( zwz, 'W', 1. ) … … 250 264 END IF 251 265 ! ! "Poleward" heat and salt transports (contribution of upstream fluxes) 252 IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nn_fptr ) == 0 )) THEN253 IF( jn == jp_tem ) htr_adv(:) = ptr_ vj( zwy(:,:,:) ) + htr_adv(:)254 IF( jn == jp_sal ) str_adv(:) = ptr_ vj( zwy(:,:,:) ) + str_adv(:)266 IF( cdtype == 'TRA' .AND. ln_diaptr ) THEN 267 IF( jn == jp_tem ) htr_adv(:) = ptr_sj( zwy(:,:,:) ) + htr_adv(:) 268 IF( jn == jp_sal ) str_adv(:) = ptr_sj( zwy(:,:,:) ) + str_adv(:) 255 269 ENDIF 256 270 ! … … 358 372 359 373 ! upstream tracer flux in the k direction 360 ! Surface value361 IF( lk_vvl ) THEN ; zwz(:,:, 1 ) = 0._wp ! volume variable362 ELSE ; zwz(:,:, 1 ) = pwn(:,:,1) * ptb(:,:,1,jn) ! linear free surface363 ENDIF364 374 ! Interior value 365 375 DO jk = 2, jpkm1 … … 372 382 END DO 373 383 END DO 384 ! Surface value 385 IF( lk_vvl ) THEN 386 IF ( ln_isfcav ) THEN 387 DO jj = 1, jpj 388 DO ji = 1, jpi 389 zwz(ji,jj, mikt(ji,jj) ) = 0.e0 ! volume variable + isf 390 END DO 391 END DO 392 ELSE 393 zwz(:,:,1) = 0.e0 ! volume variable + no isf 394 END IF 395 ELSE 396 IF ( ln_isfcav ) THEN 397 DO jj = 1, jpj 398 DO ji = 1, jpi 399 zwz(ji,jj, mikt(ji,jj) ) = pwn(ji,jj,mikt(ji,jj)) * ptb(ji,jj,mikt(ji,jj),jn) ! linear free surface + isf 400 END DO 401 END DO 402 ELSE 403 zwz(:,:,1) = pwn(:,:,1) * ptb(:,:,1,jn) ! linear free surface + no isf 404 END IF 405 ENDIF 374 406 375 407 ! total advective trend … … 398 430 END IF 399 431 ! ! "Poleward" heat and salt transports (contribution of upstream fluxes) 400 IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nn_fptr ) == 0 )) THEN401 IF( jn == jp_tem ) htr_adv(:) = ptr_ vj( zwy(:,:,:) )402 IF( jn == jp_sal ) str_adv(:) = ptr_ vj( zwy(:,:,:) )432 IF( cdtype == 'TRA' .AND. ln_diaptr ) THEN 433 IF( jn == jp_tem ) htr_adv(:) = ptr_sj( zwy(:,:,:) ) 434 IF( jn == jp_sal ) str_adv(:) = ptr_sj( zwy(:,:,:) ) 403 435 ENDIF 404 436 … … 524 556 END IF 525 557 ! ! "Poleward" heat and salt transports (contribution of upstream fluxes) 526 IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nn_fptr ) == 0 )) THEN527 IF( jn == jp_tem ) htr_adv(:) = ptr_ vj( zwy(:,:,:) ) + htr_adv(:)528 IF( jn == jp_sal ) str_adv(:) = ptr_ vj( zwy(:,:,:) ) + str_adv(:)558 IF( cdtype == 'TRA' .AND. ln_diaptr ) THEN 559 IF( jn == jp_tem ) htr_adv(:) = ptr_sj( zwy(:,:,:) ) + htr_adv(:) 560 IF( jn == jp_sal ) str_adv(:) = ptr_sj( zwy(:,:,:) ) + str_adv(:) 529 561 ENDIF 530 562 ! … … 580 612 & paft * tmask + zbig * ( 1._wp - tmask ) ) 581 613 582 DO j j = 2, jpjm1583 DO ji = fs_2, fs_jpim1 ! vector opt.584 DO jk = mikt(ji,jj), jpkm1585 ikm1 = MAX(jk-1,mikt(ji,jj))586 z2dtt = p2dt(jk)587 614 DO jk = 1, jpkm1 615 ikm1 = MAX(jk-1,1) 616 z2dtt = p2dt(jk) 617 DO jj = 2, jpjm1 618 DO ji = fs_2, fs_jpim1 ! vector opt. 619 588 620 ! search maximum in neighbourhood 589 621 zup = MAX( zbup(ji ,jj ,jk ), & -
branches/UKMO/dev_r5107_hadgem3_cplseq/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_ubs.F90
r5477 r5572 177 177 END IF 178 178 ! ! "Poleward" heat and salt transports (contribution of upstream fluxes) 179 IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nn_fptr ) == 0 )) THEN180 IF( jn == jp_tem ) htr_adv(:) = ptr_ vj( ztv(:,:,:) )181 IF( jn == jp_sal ) str_adv(:) = ptr_ vj( ztv(:,:,:) )179 IF( cdtype == 'TRA' .AND. ln_diaptr ) THEN 180 IF( jn == jp_tem ) htr_adv(:) = ptr_sj( ztv(:,:,:) ) 181 IF( jn == jp_sal ) str_adv(:) = ptr_sj( ztv(:,:,:) ) 182 182 ENDIF 183 183 -
branches/UKMO/dev_r5107_hadgem3_cplseq/NEMOGCM/NEMO/OPA_SRC/TRA/trabbc.F90
r5477 r5572 21 21 USE trdtra ! trends manager: tracers 22 22 USE in_out_manager ! I/O manager 23 USE iom ! I/O manager 24 USE fldread ! read input fields 25 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 26 USE lib_mpp ! distributed memory computing library 23 27 USE prtctl ! Print control 24 28 USE wrk_nemo ! Memory Allocation … … 37 41 38 42 REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: qgh_trd0 ! geothermal heating trend 43 TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_qgh ! structure of input qgh (file informations, fields read) 39 44 40 45 !! * Substitutions … … 92 97 END DO 93 98 ! 99 CALL lbc_lnk( tsa(:,:,:,jp_tem) , 'T', 1. ) 100 ! 94 101 IF( l_trdtra ) THEN ! Save the geothermal heat flux trend for diagnostics 95 102 ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:) … … 125 132 INTEGER :: inum ! temporary logical unit 126 133 INTEGER :: ios ! Local integer output status for namelist read 127 ! 128 NAMELIST/nambbc/ln_trabbc, nn_geoflx, rn_geoflx_cst 134 INTEGER :: ierror ! local integer 135 ! 136 TYPE(FLD_N) :: sn_qgh ! informations about the geotherm. field to be read 137 CHARACTER(len=256) :: cn_dir ! Root directory for location of ssr files 138 ! 139 NAMELIST/nambbc/ln_trabbc, nn_geoflx, rn_geoflx_cst, sn_qgh, cn_dir 129 140 !!---------------------------------------------------------------------- 130 141 … … 161 172 CASE ( 2 ) !* variable geothermal heat flux : read the geothermal fluxes in mW/m2 162 173 IF(lwp) WRITE(numout,*) ' *** variable geothermal heat flux' 163 CALL iom_open ( 'geothermal_heating.nc', inum ) 164 CALL iom_get ( inum, jpdom_data, 'heatflow', qgh_trd0 ) 165 CALL iom_close( inum ) 166 qgh_trd0(:,:) = r1_rau0_rcp * qgh_trd0(:,:) * 1.e-3 ! conversion in W/m2 174 ! 175 ALLOCATE( sf_qgh(1), STAT=ierror ) 176 IF( ierror > 0 ) THEN 177 CALL ctl_stop( 'tra_bbc_init: unable to allocate sf_qgh structure' ) ; 178 RETURN 179 ENDIF 180 ALLOCATE( sf_qgh(1)%fnow(jpi,jpj,1) ) 181 IF( sn_qgh%ln_tint )ALLOCATE( sf_qgh(1)%fdta(jpi,jpj,1,2) ) 182 ! fill sf_chl with sn_chl and control print 183 CALL fld_fill( sf_qgh, (/ sn_qgh /), cn_dir, 'tra_bbc_init', & 184 & 'bottom temperature boundary condition', 'nambbc' ) 185 186 CALL fld_read( nit000, 1, sf_qgh ) ! Read qgh data 187 qgh_trd0(:,:) = r1_rau0_rcp * sf_qgh(1)%fnow(:,:,1) * 1.e-3 ! conversion in W/m2 167 188 ! 168 189 CASE DEFAULT -
branches/UKMO/dev_r5107_hadgem3_cplseq/NEMOGCM/NEMO/OPA_SRC/TRA/traldf.F90
r5477 r5572 290 290 IF(lwp) WRITE(numout,*) ' homogeneous ocean T = ', zt0, ' S = ',zs0 291 291 292 ! Initialisation of gtui/gtvi in case of no cavity 293 IF ( .NOT. ln_isfcav ) THEN 294 gtui(:,:,:) = 0.0_wp 295 gtvi(:,:,:) = 0.0_wp 296 END IF 292 297 ! ! T & S profile (to be coded +namelist parameter 293 298 -
branches/UKMO/dev_r5107_hadgem3_cplseq/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_bilap.F90
r5477 r5572 116 116 END DO 117 117 END DO 118 119 118 ! !== Laplacian ==! 120 119 ! … … 125 124 END DO 126 125 END DO 126 ! 127 127 IF( ln_zps ) THEN ! set gradient at partial step level (last ocean level) 128 128 DO jj = 1, jpjm1 … … 130 130 IF( mbku(ji,jj) == jk ) ztu(ji,jj,jk) = zeeu(ji,jj) * pgu(ji,jj,jn) 131 131 IF( mbkv(ji,jj) == jk ) ztv(ji,jj,jk) = zeev(ji,jj) * pgv(ji,jj,jn) 132 ! (ISH)133 IF( miku(ji,jj) == jk ) ztu(ji,jj,jk) = zeeu(ji,jj) * pgui(ji,jj,jn)134 IF( mikv(ji,jj) == jk ) ztu(ji,jj,jk) = zeev(ji,jj) * pgvi(ji,jj,jn)135 132 END DO 136 133 END DO 137 134 ENDIF 135 ! (ISH) 136 IF( ln_zps .AND. ln_isfcav ) THEN ! set gradient at partial step level (first ocean level in a cavity) 137 DO jj = 1, jpjm1 138 DO ji = 1, jpim1 139 IF( miku(ji,jj) == MAX(jk,2) ) ztu(ji,jj,jk) = zeeu(ji,jj) * pgui(ji,jj,jn) 140 IF( mikv(ji,jj) == MAX(jk,2) ) ztu(ji,jj,jk) = zeev(ji,jj) * pgvi(ji,jj,jn) 141 END DO 142 END DO 143 ENDIF 144 ! 138 145 DO jj = 2, jpjm1 ! Second derivative (divergence) time the eddy diffusivity coefficient 139 146 DO ji = fs_2, fs_jpim1 ! vector opt. … … 166 173 ! 167 174 ! "zonal" mean lateral diffusive heat and salt transport 168 IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nn_fptr ) == 0 )) THEN169 IF( jn == jp_tem ) htr_ldf(:) = ptr_ vj( ztv(:,:,:) )170 IF( jn == jp_sal ) str_ldf(:) = ptr_ vj( ztv(:,:,:) )175 IF( cdtype == 'TRA' .AND. ln_diaptr ) THEN 176 IF( jn == jp_tem ) htr_ldf(:) = ptr_sj( ztv(:,:,:) ) 177 IF( jn == jp_sal ) str_ldf(:) = ptr_sj( ztv(:,:,:) ) 171 178 ENDIF 172 179 ! ! =========== -
branches/UKMO/dev_r5107_hadgem3_cplseq/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_bilapg.F90
r5477 r5572 247 247 ! ! =============== 248 248 ! "Poleward" diffusive heat or salt transport 249 IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( kaht == 2 ) .AND. ( MOD( kt, nn_fptr ) == 0 )) THEN249 IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( kaht == 2 ) ) THEN 250 250 ! note sign is reversed to give down-gradient diffusive transports (#1043) 251 IF( jn == jp_tem) htr_ldf(:) = ptr_ vj( -zftv(:,:,:) )252 IF( jn == jp_sal) str_ldf(:) = ptr_ vj( -zftv(:,:,:) )251 IF( jn == jp_tem) htr_ldf(:) = ptr_sj( -zftv(:,:,:) ) 252 IF( jn == jp_sal) str_ldf(:) = ptr_sj( -zftv(:,:,:) ) 253 253 ENDIF 254 254 -
branches/UKMO/dev_r5107_hadgem3_cplseq/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_iso.F90
r5477 r5572 28 28 USE in_out_manager ! I/O manager 29 29 USE iom ! I/O library 30 #if defined key_diaar531 30 USE phycst ! physical constants 32 31 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 33 #endif34 32 USE wrk_nemo ! Memory Allocation 35 33 USE timing ! Timing … … 106 104 ! 107 105 INTEGER :: ji, jj, jk, jn ! dummy loop indices 106 INTEGER :: ikt 108 107 REAL(wp) :: zmsku, zabe1, zcof1, zcoef3 ! local scalars 109 108 REAL(wp) :: zmskv, zabe2, zcof2, zcoef4 ! - - 110 109 REAL(wp) :: zcoef0, zbtr, ztra ! - - 111 #if defined key_diaar5112 REAL(wp) :: zztmp ! local scalar113 #endif114 110 REAL(wp), POINTER, DIMENSION(:,: ) :: z2d 115 111 REAL(wp), POINTER, DIMENSION(:,:,:) :: zdkt, zdk1t, zdit, zdjt, ztfw … … 149 145 END DO 150 146 END DO 147 148 ! partial cell correction 151 149 IF( ln_zps ) THEN ! partial steps correction at the last ocean level 152 150 DO jj = 1, jpjm1 153 151 DO ji = 1, fs_jpim1 ! vector opt. 154 152 ! IF useless if zpshde defines pgu everywhere 155 IF (mbku(ji,jj) > 1) zdit(ji,jj,mbku(ji,jj)) = pgu(ji,jj,jn) 156 IF (mbkv(ji,jj) > 1) zdjt(ji,jj,mbkv(ji,jj)) = pgv(ji,jj,jn) 157 ! (ISF) 153 zdit(ji,jj,mbku(ji,jj)) = pgu(ji,jj,jn) 154 zdjt(ji,jj,mbkv(ji,jj)) = pgv(ji,jj,jn) 155 END DO 156 END DO 157 ENDIF 158 IF( ln_zps .AND. ln_isfcav ) THEN ! partial steps correction at the first wet level beneath a cavity 159 DO jj = 1, jpjm1 160 DO ji = 1, fs_jpim1 ! vector opt. 158 161 IF (miku(ji,jj) > 1) zdit(ji,jj,miku(ji,jj)) = pgui(ji,jj,jn) 159 162 IF (mikv(ji,jj) > 1) zdjt(ji,jj,mikv(ji,jj)) = pgvi(ji,jj,jn) 160 163 END DO 161 164 END DO 162 END IF165 END IF 163 166 164 167 !!---------------------------------------------------------------------- 165 168 !! II - horizontal trend (full) 166 169 !!---------------------------------------------------------------------- 167 !CDIR PARALLEL DO PRIVATE( zdk1t ) 168 ! ! =============== 169 DO jj = 1, jpj ! Horizontal slab 170 ! ! =============== 171 DO ji = 1, jpi ! vector opt. 172 DO jk = mikt(ji,jj), jpkm1 173 ! 1. Vertical tracer gradient at level jk and jk+1 174 ! ------------------------------------------------ 175 ! surface boundary condition: zdkt(jk=1)=zdkt(jk=2) 176 zdk1t(ji,jj,jk) = ( ptb(ji,jj,jk,jn) - ptb(ji,jj,jk+1,jn) ) * tmask(ji,jj,jk+1) 177 ! 178 IF( jk == mikt(ji,jj) ) THEN ; zdkt(ji,jj,jk) = zdk1t(ji,jj,jk) 179 ELSE ; zdkt(ji,jj,jk) = ( ptb(ji,jj,jk-1,jn) - ptb(ji,jj,jk,jn) ) * tmask(ji,jj,jk) 180 ENDIF 170 !!!!!!!!!!CDIR PARALLEL DO PRIVATE( zdk1t ) 171 ! 1. Vertical tracer gradient at level jk and jk+1 172 ! ------------------------------------------------ 173 ! 174 ! interior value 175 DO jk = 2, jpkm1 176 DO jj = 1, jpj 177 DO ji = 1, jpi ! vector opt. 178 zdk1t(ji,jj,jk) = ( ptb(ji,jj,jk,jn ) - ptb(ji,jj,jk+1,jn) ) * wmask(ji,jj,jk+1) 179 ! 180 zdkt(ji,jj,jk) = ( ptb(ji,jj,jk-1,jn) - ptb(ji,jj,jk,jn ) ) * wmask(ji,jj,jk) 181 181 END DO 182 182 END DO 183 183 END DO 184 185 ! 2. Horizontal fluxes 186 ! -------------------- 187 DO jj = 1 , jpjm1 188 DO ji = 1, fs_jpim1 ! vector opt. 189 DO jk = mikt(ji,jj), jpkm1 184 ! surface boundary condition: zdkt(jk=1)=zdkt(jk=2) 185 zdk1t(:,:,1) = ( ptb(:,:,1,jn ) - ptb(:,:,2,jn) ) * wmask(:,:,2) 186 zdkt (:,:,1) = zdk1t(:,:,1) 187 IF ( ln_isfcav ) THEN 188 DO jj = 1, jpj 189 DO ji = 1, jpi ! vector opt. 190 ikt = mikt(ji,jj) ! surface level 191 zdk1t(ji,jj,ikt) = ( ptb(ji,jj,ikt,jn ) - ptb(ji,jj,ikt+1,jn) ) * wmask(ji,jj,ikt+1) 192 zdkt (ji,jj,ikt) = zdk1t(ji,jj,ikt) 193 END DO 194 END DO 195 END IF 196 197 ! 2. Horizontal fluxes 198 ! -------------------- 199 DO jk = 1, jpkm1 200 DO jj = 1 , jpjm1 201 DO ji = 1, fs_jpim1 ! vector opt. 190 202 zabe1 = ( fsahtu(ji,jj,jk) + pahtb0 ) * re2u_e1u(ji,jj) * fse3u_n(ji,jj,jk) 191 203 zabe2 = ( fsahtv(ji,jj,jk) + pahtb0 ) * re1v_e2v(ji,jj) * fse3v_n(ji,jj,jk) … … 208 220 END DO 209 221 END DO 210 END DO211 222 212 223 ! II.4 Second derivative (divergence) and add to the general trend 213 224 ! ---------------------------------------------------------------- 214 DO jj = 2 , jpjm1 215 DO ji = fs_2, fs_jpim1 ! vector opt. 216 DO jk = mikt(ji,jj), jpkm1 217 zbtr = 1.0 / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 225 DO jj = 2 , jpjm1 226 DO ji = fs_2, fs_jpim1 ! vector opt. 227 zbtr = 1.0 / ( e12t(ji,jj) * fse3t_n(ji,jj,jk) ) 218 228 ztra = zbtr * ( zftu(ji,jj,jk) - zftu(ji-1,jj,jk) + zftv(ji,jj,jk) - zftv(ji,jj-1,jk) ) 219 229 pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra … … 225 235 ! 226 236 ! "Poleward" diffusive heat or salt transports (T-S case only) 227 IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nn_fptr ) == 0 )) THEN237 IF( cdtype == 'TRA' .AND. ln_diaptr ) THEN 228 238 ! note sign is reversed to give down-gradient diffusive transports (#1043) 229 IF( jn == jp_tem) htr_ldf(:) = ptr_ vj( -zftv(:,:,:) )230 IF( jn == jp_sal) str_ldf(:) = ptr_ vj( -zftv(:,:,:) )239 IF( jn == jp_tem) htr_ldf(:) = ptr_sj( -zftv(:,:,:) ) 240 IF( jn == jp_sal) str_ldf(:) = ptr_sj( -zftv(:,:,:) ) 231 241 ENDIF 232 242 233 #if defined key_diaar5 234 IF( cdtype == 'TRA' .AND. jn == jp_tem ) THEN235 z2d(:,:) = 0._wp236 ! note sign is reversed to give down-gradient diffusive transports (#1043)237 zztmp = -1.0_wp * rau0 * rcp238 DO jk = 1, jpkm1239 DO jj = 2, jpjm1240 DO ji = fs_2, fs_jpim1 ! vector opt.241 z2d(ji,jj) = z2d(ji,jj) + zftu(ji,jj,jk)243 IF( iom_use("udiff_heattr") .OR. iom_use("vdiff_heattr") ) THEN 244 ! 245 IF( cdtype == 'TRA' .AND. jn == jp_tem ) THEN 246 z2d(:,:) = 0._wp 247 DO jk = 1, jpkm1 248 DO jj = 2, jpjm1 249 DO ji = fs_2, fs_jpim1 ! vector opt. 250 z2d(ji,jj) = z2d(ji,jj) + zftu(ji,jj,jk) 251 END DO 242 252 END DO 243 253 END DO 244 END DO 245 z2d(:,:) = zztmp * z2d(:,:) 246 CALL lbc_lnk( z2d, 'U', -1. ) 247 CALL iom_put( "udiff_heattr", z2d ) ! heat transport in i-direction 248 z2d(:,:) = 0._wp 249 DO jk = 1, jpkm1 250 DO jj = 2, jpjm1 251 DO ji = fs_2, fs_jpim1 ! vector opt. 252 z2d(ji,jj) = z2d(ji,jj) + zftv(ji,jj,jk) 254 z2d(:,:) = - rau0_rcp * z2d(:,:) ! note sign is reversed to give down-gradient diffusive transports (#1043) 255 CALL lbc_lnk( z2d, 'U', -1. ) 256 CALL iom_put( "udiff_heattr", z2d ) ! heat transport in i-direction 257 ! 258 z2d(:,:) = 0._wp 259 DO jk = 1, jpkm1 260 DO jj = 2, jpjm1 261 DO ji = fs_2, fs_jpim1 ! vector opt. 262 z2d(ji,jj) = z2d(ji,jj) + zftv(ji,jj,jk) 263 END DO 253 264 END DO 254 265 END DO 255 END DO256 z2d(:,:) = zztmp * z2d(:,:)257 CALL lbc_lnk( z2d, 'V', -1. )258 CALL iom_put( "vdiff_heattr", z2d ) ! heat transport in i-direction259 END IF260 #endif 266 z2d(:,:) = - rau0_rcp * z2d(:,:) ! note sign is reversed to give down-gradient diffusive transports (#1043) 267 CALL lbc_lnk( z2d, 'V', -1. ) 268 CALL iom_put( "vdiff_heattr", z2d ) ! heat transport in i-direction 269 END IF 270 ! 271 ENDIF 261 272 262 273 !!---------------------------------------------------------------------- … … 278 289 DO jj = 2, jpjm1 279 290 DO ji = fs_2, fs_jpim1 ! vector opt. 280 zcoef0 = - fsahtw(ji,jj,jk) * tmask(ji,jj,jk) * tmask(ji,jj,jk-1)291 zcoef0 = - fsahtw(ji,jj,jk) * wmask(ji,jj,jk) 281 292 ! 282 293 zmsku = 1./MAX( umask(ji ,jj,jk-1) + umask(ji-1,jj,jk) & -
branches/UKMO/dev_r5107_hadgem3_cplseq/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_iso_grif.F90
r5477 r5572 113 113 REAL(wp) :: ze1ur, zdxt, ze2vr, ze3wr, zdyt, zdzt 114 114 REAL(wp) :: zah, zah_slp, zaei_slp 115 #if defined key_diaar5116 REAL(wp) :: zztmp ! local scalar117 #endif118 115 REAL(wp), POINTER, DIMENSION(:,: ) :: z2d 119 116 REAL(wp), POINTER, DIMENSION(:,:,:) :: zdit, zdjt, ztfw … … 207 204 END DO 208 205 ! 209 #if defined key_iomput 210 IF( ln_traldf_gdia .AND. cdtype == 'TRA' ) THEN 211 CALL wrk_alloc( jpi , jpj , jpk , zw3d ) 212 DO jk=1,jpkm1 213 zw3d(:,:,jk) = (psix_eiv(:,:,jk+1) - psix_eiv(:,:,jk))/fse3u(:,:,jk) ! u_eiv = -dpsix/dz 214 END DO 215 zw3d(:,:,jpk) = 0._wp 216 CALL iom_put( "uoce_eiv", zw3d ) ! i-eiv current 217 218 DO jk=1,jpk-1 219 zw3d(:,:,jk) = (psiy_eiv(:,:,jk+1) - psiy_eiv(:,:,jk))/fse3v(:,:,jk) ! v_eiv = -dpsiy/dz 220 END DO 221 zw3d(:,:,jpk) = 0._wp 222 CALL iom_put( "voce_eiv", zw3d ) ! j-eiv current 223 224 DO jk=1,jpk-1 225 DO jj = 2, jpjm1 226 DO ji = fs_2, fs_jpim1 ! vector opt. 227 zw3d(ji,jj,jk) = (psiy_eiv(ji,jj,jk) - psiy_eiv(ji,jj-1,jk))/e2t(ji,jj) + & 228 & (psix_eiv(ji,jj,jk) - psix_eiv(ji-1,jj,jk))/e1t(ji,jj) ! w_eiv = dpsiy/dy + dpsiy/dx 229 END DO 230 END DO 231 END DO 232 zw3d(:,:,jpk) = 0._wp 233 CALL iom_put( "woce_eiv", zw3d ) ! vert. eiv current 234 CALL wrk_dealloc( jpi , jpj , jpk , zw3d ) 206 IF( iom_use("uoce_eiv") .OR. iom_use("voce_eiv") .OR. iom_use("woce_eiv") ) THEN 207 ! 208 IF( ln_traldf_gdia .AND. cdtype == 'TRA' ) THEN 209 CALL wrk_alloc( jpi , jpj , jpk , zw3d ) 210 DO jk=1,jpkm1 211 zw3d(:,:,jk) = (psix_eiv(:,:,jk+1) - psix_eiv(:,:,jk))/fse3u(:,:,jk) ! u_eiv = -dpsix/dz 212 END DO 213 zw3d(:,:,jpk) = 0._wp 214 CALL iom_put( "uoce_eiv", zw3d ) ! i-eiv current 215 216 DO jk=1,jpk-1 217 zw3d(:,:,jk) = (psiy_eiv(:,:,jk+1) - psiy_eiv(:,:,jk))/fse3v(:,:,jk) ! v_eiv = -dpsiy/dz 218 END DO 219 zw3d(:,:,jpk) = 0._wp 220 CALL iom_put( "voce_eiv", zw3d ) ! j-eiv current 221 222 DO jk=1,jpk-1 223 DO jj = 2, jpjm1 224 DO ji = fs_2, fs_jpim1 ! vector opt. 225 zw3d(ji,jj,jk) = (psiy_eiv(ji,jj,jk) - psiy_eiv(ji,jj-1,jk))/e2t(ji,jj) + & 226 & (psix_eiv(ji,jj,jk) - psix_eiv(ji-1,jj,jk))/e1t(ji,jj) ! w_eiv = dpsiy/dy + dpsiy/dx 227 END DO 228 END DO 229 END DO 230 zw3d(:,:,jpk) = 0._wp 231 CALL iom_put( "woce_eiv", zw3d ) ! vert. eiv current 232 CALL wrk_dealloc( jpi , jpj , jpk , zw3d ) 233 ENDIF 234 ! 235 235 ENDIF 236 #endif237 236 ! ! =========== 238 237 DO jn = 1, kjpt ! tracer loop … … 387 386 ! 388 387 ! ! "Poleward" diffusive heat or salt transports (T-S case only) 389 IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nn_fptr ) == 0 )) THEN390 IF( jn == jp_tem) htr_ldf(:) = ptr_ vj( zftv(:,:,:) ) ! 3.3 names391 IF( jn == jp_sal) str_ldf(:) = ptr_ vj( zftv(:,:,:) )388 IF( cdtype == 'TRA' .AND. ln_diaptr ) THEN 389 IF( jn == jp_tem) htr_ldf(:) = ptr_sj( zftv(:,:,:) ) ! 3.3 names 390 IF( jn == jp_sal) str_ldf(:) = ptr_sj( zftv(:,:,:) ) 392 391 ENDIF 393 392 394 #if defined key_diaar5 395 IF( cdtype == 'TRA' .AND. jn == jp_tem ) THEN 396 z2d(:,:) = 0._wp 397 zztmp = rau0 * rcp 398 DO jk = 1, jpkm1 399 DO jj = 2, jpjm1 400 DO ji = fs_2, fs_jpim1 ! vector opt. 401 z2d(ji,jj) = z2d(ji,jj) + zftu(ji,jj,jk) 402 END DO 403 END DO 404 END DO 405 z2d(:,:) = zztmp * z2d(:,:) 406 CALL lbc_lnk( z2d, 'U', -1. ) 407 CALL iom_put( "udiff_heattr", z2d ) ! heat transport in i-direction 408 z2d(:,:) = 0._wp 409 DO jk = 1, jpkm1 410 DO jj = 2, jpjm1 411 DO ji = fs_2, fs_jpim1 ! vector opt. 412 z2d(ji,jj) = z2d(ji,jj) + zftv(ji,jj,jk) 413 END DO 414 END DO 415 END DO 416 z2d(:,:) = zztmp * z2d(:,:) 417 CALL lbc_lnk( z2d, 'V', -1. ) 418 CALL iom_put( "vdiff_heattr", z2d ) ! heat transport in j-direction 419 END IF 420 #endif 393 IF( iom_use("udiff_heattr") .OR. iom_use("vdiff_heattr") ) THEN 394 ! 395 IF( cdtype == 'TRA' .AND. jn == jp_tem ) THEN 396 z2d(:,:) = 0._wp 397 DO jk = 1, jpkm1 398 DO jj = 2, jpjm1 399 DO ji = fs_2, fs_jpim1 ! vector opt. 400 z2d(ji,jj) = z2d(ji,jj) + zftu(ji,jj,jk) 401 END DO 402 END DO 403 END DO 404 z2d(:,:) = rau0_rcp * z2d(:,:) 405 CALL lbc_lnk( z2d, 'U', -1. ) 406 CALL iom_put( "udiff_heattr", z2d ) ! heat transport in i-direction 407 ! 408 z2d(:,:) = 0._wp 409 DO jk = 1, jpkm1 410 DO jj = 2, jpjm1 411 DO ji = fs_2, fs_jpim1 ! vector opt. 412 z2d(ji,jj) = z2d(ji,jj) + zftv(ji,jj,jk) 413 END DO 414 END DO 415 END DO 416 z2d(:,:) = rau0_rcp * z2d(:,:) 417 CALL lbc_lnk( z2d, 'V', -1. ) 418 CALL iom_put( "vdiff_heattr", z2d ) ! heat transport in i-direction 419 END IF 420 ! 421 ENDIF 421 422 ! 422 423 END DO -
branches/UKMO/dev_r5107_hadgem3_cplseq/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_lap.F90
r5477 r5572 102 102 END DO 103 103 END DO 104 IF( ln_zps ) THEN ! set gradient at partial step level 104 IF( ln_zps ) THEN ! set gradient at partial step level for the last ocean cell 105 105 DO jj = 1, jpjm1 106 106 DO ji = 1, fs_jpim1 ! vector opt. … … 116 116 ztv(ji,jj,jk) = zabe2 * pgv(ji,jj,jn) 117 117 ENDIF 118 119 ! (ISH) 118 END DO 119 END DO 120 ENDIF 121 ! (ISH) 122 IF( ln_zps .AND. ln_isfcav ) THEN ! set gradient at partial step level for the first ocean cell 123 ! into a cavity 124 DO jj = 1, jpjm1 125 DO ji = 1, fs_jpim1 ! vector opt. 120 126 ! ice shelf level level MAX(2,jk) => only where ice shelf 121 127 iku = miku(ji,jj) … … 148 154 ! 149 155 ! "Poleward" diffusive heat or salt transports 150 IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nn_fptr ) == 0 )) THEN151 IF( jn == jp_tem) htr_ldf(:) = ptr_ vj( ztv(:,:,:) )152 IF( jn == jp_sal) str_ldf(:) = ptr_ vj( ztv(:,:,:) )156 IF( cdtype == 'TRA' .AND. ln_diaptr ) THEN 157 IF( jn == jp_tem) htr_ldf(:) = ptr_sj( ztv(:,:,:) ) 158 IF( jn == jp_sal) str_ldf(:) = ptr_sj( ztv(:,:,:) ) 153 159 ENDIF 154 160 ! ! ================== -
branches/UKMO/dev_r5107_hadgem3_cplseq/NEMOGCM/NEMO/OPA_SRC/TRA/tranpc.F90
r5477 r5572 9 9 !! 3.0 ! 2008-06 (G. Madec) applied on ta, sa and called before tranxt in step.F90 10 10 !! 3.3 ! 2010-05 (C. Ethe, G. Madec) merge TRC-TRA 11 !! 3. 7 ! 2014-06(L. Brodeau) new algorithm based on local Brunt-Vaisala freq.11 !! 3.6 ! 2015-05 (L. Brodeau) new algorithm based on local Brunt-Vaisala freq. 12 12 !!---------------------------------------------------------------------- 13 13 … … 64 64 INTEGER :: ji, jj, jk ! dummy loop indices 65 65 INTEGER :: inpcc ! number of statically instable water column 66 INTEGER :: jiter, ikbot, ik , ikup, ikdown, ilayer, ikm! local integers66 INTEGER :: jiter, ikbot, ikp, ikup, ikdown, ilayer, ik_low ! local integers 67 67 LOGICAL :: l_bottom_reached, l_column_treated 68 68 REAL(wp) :: zta, zalfa, zsum_temp, zsum_alfa, zaw, zdz, zsum_z 69 69 REAL(wp) :: zsa, zbeta, zsum_sali, zsum_beta, zbw, zrw, z1_r2dt 70 REAL(wp), PARAMETER :: zn2_zero = 1.e-14_wp ! acceptance criteria for neutrality (N2==0) 70 71 REAL(wp), POINTER, DIMENSION(:) :: zvn2 ! vertical profile of N2 at 1 given point... 71 72 REAL(wp), POINTER, DIMENSION(:,:) :: zvts ! vertical profile of T and S at 1 given point... … … 75 76 REAL(wp), POINTER, DIMENSION(:,:,:) :: ztrdt, ztrds ! 3D workspace 76 77 ! 77 !!LB debug: 78 LOGICAL, PARAMETER :: l_LB_debug = .FALSE. 79 INTEGER :: ilc1, jlc1, klc1, nncpu 80 LOGICAL :: lp_monitor_point = .FALSE. 81 !!LB debug. 78 LOGICAL, PARAMETER :: l_LB_debug = .FALSE. ! set to true if you want to follow what is 79 INTEGER :: ilc1, jlc1, klc1, nncpu ! actually happening in a water column at point "ilc1, jlc1" 80 LOGICAL :: lp_monitor_point = .FALSE. ! in CPU domain "nncpu" 82 81 !!---------------------------------------------------------------------- 83 82 ! … … 97 96 ENDIF 98 97 99 !LB debug:100 IF( lwp .AND. l_LB_debug ) THEN101 WRITE(numout,*)102 WRITE(numout,*) 'LOLO: entering tra_npc, kt, narea =', kt, narea103 ENDIF104 !LBdebug: Monitoring of 1 column subject to convection...105 98 IF( l_LB_debug ) THEN 106 ! Location of 1 known convection spot to follow what's happening in the water column 107 ilc1 = 54 ; jlc1 = 15 ; ! Labrador ORCA1 4x4 cpus: 108 nncpu = 15 ; ! the CPU domain contains the convection spot 109 !ilc1 = 14 ; jlc1 = 13 ; ! Labrador ORCA1 8x8 cpus: 110 !nncpu = 54 ; ! the CPU domain contains the convection spot 99 ! Location of 1 known convection site to follow what's happening in the water column 100 ilc1 = 45 ; jlc1 = 3 ; ! ORCA2 4x4, Antarctic coast, more than 2 unstable portions in the water column... 101 nncpu = 1 ; ! the CPU domain contains the convection spot 111 102 klc1 = mbkt(ilc1,jlc1) ! bottom of the ocean for debug point... 112 103 ENDIF 113 !LBdebug. 114 115 CALL eos_rab( tsa, zab ) ! after alpha and beta 116 CALL bn2 ( tsa, zab, zn2 ) ! after Brunt-Vaisala 104 105 CALL eos_rab( tsa, zab ) ! after alpha and beta (given on T-points) 106 CALL bn2 ( tsa, zab, zn2 ) ! after Brunt-Vaisala (given on W-points) 117 107 118 108 inpcc = 0 … … 134 124 IF( ( ji == ilc1 ).AND.( jj == jlc1 ) ) lp_monitor_point = .TRUE. 135 125 ! writing only if on CPU domain where conv region is: 136 lp_monitor_point = (narea == nncpu).AND.lp_monitor_point 137 138 IF(lp_monitor_point) THEN 139 WRITE(numout,*) '' ;WRITE(numout,*) '' ; 140 WRITE(numout,'("Time step = ",i6.6," !!!")') kt 141 WRITE(numout,'(" *** BEFORE anything, N^2 for point ",i3,",",i3,":" )') ji,jj 142 DO jk = 1, klc1 143 WRITE(numout,*) jk, zvn2(jk) 144 END DO 145 WRITE(numout,*) ' ' 146 ENDIF 126 lp_monitor_point = (narea == nncpu).AND.lp_monitor_point 147 127 ENDIF !LB debug end 148 128 149 129 ikbot = mbkt(ji,jj) ! ikbot: ocean bottom T-level 150 ik = 1 ! because N2 is irrelevant at the surface level (will start at ik=2)130 ikp = 1 ! because N2 is irrelevant at the surface level (will start at ikp=2) 151 131 ilayer = 0 152 132 jiter = 0 … … 163 143 DO WHILE ( .NOT. l_bottom_reached ) 164 144 165 ik = ik+ 1145 ikp = ikp + 1 166 146 167 !! Checking level ikfor instability147 !! Testing level ikp for instability 168 148 !! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 169 170 IF( zvn2(ik) < 0. ) THEN ! Instability found! 171 172 ikm = ik ! first level whith negative N2 173 ilayer = ilayer + 1 ! yet another layer found.... 174 IF(jiter == 1) inpcc = inpcc + 1 175 176 IF(l_LB_debug .AND. lp_monitor_point) & 177 & WRITE(numout,*) 'Negative N2 at ik =', ikm, ' layer nb.', ilayer, & 178 & ' inpcc =', inpcc 179 180 !! Case we mix with upper regions where N2==0: 181 !! All the points above ikup where N2 == 0 must also be mixed => we go 182 !! upward to find a new ikup, where the layer doesn't have N2==0 183 ikup = ikm 184 DO jk = ikm, 2, -1 185 ikup = ikup - 1 186 IF( (zvn2(jk-1) > 0.).OR.(ikup == 1) ) EXIT 187 END DO 188 189 ! adjusting ikup if the upper part of the unstable column was neutral (N2=0) 190 IF((zvn2(ikup+1) == 0.).AND.(ikup /= 1)) ikup = ikup+1 ; 191 192 193 IF( lp_monitor_point ) WRITE(numout,*) ' => ikup is =', ikup, ' layer nb.', ilayer 194 149 IF( zvn2(ikp) < -zn2_zero ) THEN ! Instability found! 150 151 ilayer = ilayer + 1 ! yet another instable portion of the water column found.... 152 153 IF( lp_monitor_point ) THEN 154 WRITE(numout,*) 155 IF( ilayer == 1 .AND. jiter == 1 ) THEN ! first time a column is spoted with an instability 156 WRITE(numout,*) 157 WRITE(numout,*) 'Time step = ',kt,' !!!' 158 ENDIF 159 WRITE(numout,*) ' * Iteration #',jiter,': found instable portion #',ilayer, & 160 & ' in column! Starting at ikp =', ikp 161 WRITE(numout,*) ' *** N2 for point (i,j) = ',ji,' , ',jj 162 DO jk = 1, klc1 163 WRITE(numout,*) jk, zvn2(jk) 164 END DO 165 WRITE(numout,*) 166 ENDIF 167 168 169 IF( jiter == 1 ) inpcc = inpcc + 1 170 171 IF( lp_monitor_point ) WRITE(numout, *) 'Negative N2 at ikp =',ikp,' for layer #', ilayer 172 173 !! ikup is the uppermost point where mixing will start: 174 ikup = ikp - 1 ! ikup is always "at most at ikp-1", less if neutral levels overlying 175 176 !! If the points above ikp-1 have N2 == 0 they must also be mixed: 177 IF( ikp > 2 ) THEN 178 DO jk = ikp-1, 2, -1 179 IF( ABS(zvn2(jk)) < zn2_zero ) THEN 180 ikup = ikup - 1 ! 1 more upper level has N2=0 and must be added for the mixing 181 ELSE 182 EXIT 183 ENDIF 184 END DO 185 ENDIF 186 187 IF( ikup < 1 ) CALL ctl_stop( 'tra_npc : PROBLEM #1') 188 195 189 zsum_temp = 0._wp 196 190 zsum_sali = 0._wp … … 199 193 zsum_z = 0._wp 200 194 201 DO jk = ikup, ikbot+1 ! Inside the instable (and overlying neutral) portion of the column 202 ! 203 IF(l_LB_debug .AND. lp_monitor_point) WRITE(numout,*) ' -> summing for jk =', jk 195 DO jk = ikup, ikbot ! Inside the instable (and overlying neutral) portion of the column 204 196 ! 205 197 zdz = fse3t(ji,jj,jk) … … 209 201 zsum_beta = zsum_beta + zvab(jk,jp_sal)*zdz 210 202 zsum_z = zsum_z + zdz 211 ! 212 !! EXIT if we found the bottom of the unstable portion of the water column 213 IF( (zvn2(jk+1) > 0.).OR.(jk == ikbot ).OR.((jk==ikm).AND.(zvn2(jk+1) == 0.)) ) EXIT 203 ! 204 IF( jk == ikbot ) EXIT ! avoid array-index overshoot in case ikbot = jpk, cause we're calling jk+1 next line 205 !! EXIT when we have reached the last layer that is instable (N2<0) or neutral (N2=0): 206 IF( zvn2(jk+1) > zn2_zero ) EXIT 214 207 END DO 215 208 216 !ik = jk !LB remove? 217 ikdown = jk ! for the current unstable layer, ikdown is the deepest point with a negative N2 218 219 IF(l_LB_debug .AND. lp_monitor_point) & 220 & WRITE(numout,*) ' => ikdown =', ikdown, ' layer nb.', ilayer 221 222 ! Mixing Temperature and salinity between ikup and ikdown: 209 ikdown = jk ! for the current unstable layer, ikdown is the deepest point with a negative or neutral N2 210 IF( ikup == ikdown ) CALL ctl_stop( 'tra_npc : PROBLEM #2') 211 212 ! Mixing Temperature, salinity, alpha and beta from ikup to ikdown included: 223 213 zta = zsum_temp/zsum_z 224 214 zsa = zsum_sali/zsum_z … … 226 216 zbeta = zsum_beta/zsum_z 227 217 228 IF(l_LB_debug .AND. lp_monitor_point) THEN 218 IF( lp_monitor_point ) THEN 219 WRITE(numout,*) 'MIXED T, S, alfa and beta between ikup =',ikup, & 220 & ' and ikdown =',ikdown,', in layer #',ilayer 229 221 WRITE(numout,*) ' => Mean temp. in that portion =', zta 230 222 WRITE(numout,*) ' => Mean sali. in that portion =', zsa 231 WRITE(numout,*) ' => Mean Al phain that portion =', zalfa223 WRITE(numout,*) ' => Mean Alfa in that portion =', zalfa 232 224 WRITE(numout,*) ' => Mean Beta in that portion =', zbeta 233 225 ENDIF … … 240 232 zvab(jk,jp_sal) = zbeta 241 233 END DO 242 ! 243 !! Before updating N2, it is possible that another unstable 244 !! layer exists underneath the one we just homogeneized! 245 ik = ikdown 246 ! 247 ENDIF ! IF( zvn2(ik+1) < 0. ) THEN 248 ! 249 IF( ik == ikbot ) l_bottom_reached = .TRUE. 234 235 236 !! Updating N2 in the relvant portion of the water column 237 !! Temperature, Salinity, Alpha and Beta have been homogenized in the unstable portion 238 !! => Need to re-compute N2! will use Alpha and Beta! 239 240 ikup = MAX(2,ikup) ! ikup can never be 1 ! 241 ik_low = MIN(ikdown+1,ikbot) ! we must go 1 point deeper than ikdown! 242 243 DO jk = ikup, ik_low ! we must go 1 point deeper than ikdown! 244 245 !! Interpolating alfa and beta at W point: 246 zrw = (fsdepw(ji,jj,jk ) - fsdept(ji,jj,jk)) & 247 & / (fsdept(ji,jj,jk-1) - fsdept(ji,jj,jk)) 248 zaw = zvab(jk,jp_tem) * (1._wp - zrw) + zvab(jk-1,jp_tem) * zrw 249 zbw = zvab(jk,jp_sal) * (1._wp - zrw) + zvab(jk-1,jp_sal) * zrw 250 251 !! N2 at W point, doing exactly as in eosbn2.F90: 252 zvn2(jk) = grav*( zaw * ( zvts(jk-1,jp_tem) - zvts(jk,jp_tem) ) & 253 & - zbw * ( zvts(jk-1,jp_sal) - zvts(jk,jp_sal) ) ) & 254 & / fse3w(ji,jj,jk) * tmask(ji,jj,jk) 255 256 !! OR, faster => just considering the vertical gradient of density 257 !! as only the signa maters... 258 !zvn2(jk) = ( zaw * ( zvts(jk-1,jp_tem) - zvts(jk,jp_tem) ) & 259 ! & - zbw * ( zvts(jk-1,jp_sal) - zvts(jk,jp_sal) ) ) 260 261 END DO 262 263 ikp = MIN(ikdown+1,ikbot) 264 265 266 ENDIF !IF( zvn2(ikp) < 0. ) 267 268 269 IF( ikp == ikbot ) l_bottom_reached = .TRUE. 250 270 ! 251 271 END DO ! DO WHILE ( .NOT. l_bottom_reached ) 252 272 253 IF( ik /= ikbot ) STOP 'ERROR: tranpc.F90 => PROBLEM #1'273 IF( ikp /= ikbot ) CALL ctl_stop( 'tra_npc : PROBLEM #3') 254 274 255 ! ******* At this stage ik == ikbot ! *******275 ! ******* At this stage ikp == ikbot ! ******* 256 276 257 IF( ilayer > 0 ) THEN 258 !! least an unstable layer has been found 259 !! Temperature, Salinity, Alpha and Beta have been homogenized in the unstable portion 260 !! => Need to re-compute N2! will use Alpha and Beta! 277 IF( ilayer > 0 ) THEN !! least an unstable layer has been found 261 278 ! 262 DO jk = ikup+1, ikdown+1 ! we must go 1 point deeper than ikdown! 263 !! Doing exactly as in eosbn2.F90: 264 !! * Except that we only are interested in the sign of N2 !!! 265 !! => just considering the vertical gradient of density 266 zrw = (fsdepw(ji,jj,jk ) - fsdept(ji,jj,jk)) & 267 & / (fsdept(ji,jj,jk-1) - fsdept(ji,jj,jk)) 268 zaw = zvab(jk,jp_tem) * (1._wp - zrw) + zvab(jk-1,jp_tem) * zrw 269 zbw = zvab(jk,jp_sal) * (1._wp - zrw) + zvab(jk-1,jp_sal) * zrw 270 271 !zvn2(jk) = grav*( zaw * ( zvts(jk-1,jp_tem) - zvts(jk,jp_tem) ) & 272 ! & - zbw * ( zvts(jk-1,jp_sal) - zvts(jk,jp_sal) ) ) & 273 ! & / fse3w(ji,jj,jk) * tmask(ji,jj,jk) 274 zvn2(jk) = ( zaw * ( zvts(jk-1,jp_tem) - zvts(jk,jp_tem) ) & 275 & - zbw * ( zvts(jk-1,jp_sal) - zvts(jk,jp_sal) ) ) 276 END DO 277 278 IF(l_LB_debug .AND. lp_monitor_point) THEN 279 WRITE(numout, '(" *** After iteration #",i3.3,", N^2 for point ",i3,",",i3,":" )') & 280 & jiter, ji,jj 279 IF( lp_monitor_point ) THEN 280 WRITE(numout,*) 281 WRITE(numout,*) 'After ',jiter,' iteration(s), we neutralized ',ilayer,' instable layer(s)' 282 WRITE(numout,*) ' ==> N2 at i,j=',ji,',',jj,' now looks like this:' 281 283 DO jk = 1, klc1 282 284 WRITE(numout,*) jk, zvn2(jk) 283 285 END DO 284 WRITE(numout,*) ' '286 WRITE(numout,*) 285 287 ENDIF 286 287 ik = 1! starting again at the surface for the next iteration288 ! 289 ikp = 1 ! starting again at the surface for the next iteration 288 290 ilayer = 0 289 291 ENDIF 290 292 ! 291 IF( ik >= ikbot ) THEN 292 IF(l_LB_debug .AND. lp_monitor_point) WRITE(numout,*) ' --- exiting jiter loop ---' 293 l_column_treated = .TRUE. 294 ENDIF 293 IF( ikp >= ikbot ) l_column_treated = .TRUE. 295 294 ! 296 295 END DO ! DO WHILE ( .NOT. l_column_treated ) … … 300 299 tsa(ji,jj,:,jp_sal) = zvts(:,jp_sal) 301 300 302 !! lolo: Should we update something else????303 !! => like alpha and beta?304 305 IF( l_LB_debug .AND. lp_monitor_point) WRITE(numout,*) ''301 !! LB: Potentially some other global variable beside theta and S can be treated here 302 !! like BGC tracers. 303 304 IF( lp_monitor_point ) WRITE(numout,*) 306 305 307 306 ENDIF ! IF( tmask(ji,jj,3) == 1 ) THEN … … 321 320 CALL lbc_lnk( tsa(:,:,:,jp_tem), 'T', 1. ) ; CALL lbc_lnk( tsa(:,:,:,jp_sal), 'T', 1. ) 322 321 ! 323 IF(lwp) THEN 324 WRITE(numout,*) 'LOLO: exiting tra_npc, kt =', kt 325 WRITE(numout,*)' => number of statically instable water column : ',inpcc 326 WRITE(numout,*) '' ; WRITE(numout,*) '' 322 IF( lwp .AND. l_LB_debug ) THEN 323 WRITE(numout,*) 'Exiting tra_npc , kt = ',kt,', => numb. of statically instable water-columns: ', inpcc 324 WRITE(numout,*) 327 325 ENDIF 328 326 ! -
branches/UKMO/dev_r5107_hadgem3_cplseq/NEMOGCM/NEMO/OPA_SRC/TRA/tranxt.F90
r5477 r5572 27 27 USE dom_oce ! ocean space and time domain variables 28 28 USE sbc_oce ! surface boundary condition: ocean 29 USE sbcrnf ! river runoffs 29 30 USE zdf_oce ! ocean vertical mixing 30 31 USE domvvl ! variable volume … … 143 144 ELSE ! Leap-Frog + Asselin filter time stepping 144 145 ! 145 IF( lk_vvl ) THEN ; CALL tra_nxt_vvl( kt, nit000, 'TRA', tsb, tsn, tsa, jpts ) ! variable volume level (vvl) 146 ELSE ; CALL tra_nxt_fix( kt, nit000, 'TRA', tsb, tsn, tsa, jpts ) ! fixed volume level 146 IF( lk_vvl ) THEN ; CALL tra_nxt_vvl( kt, nit000, rdttra, 'TRA', tsb, tsn, tsa, & 147 & sbc_tsc, sbc_tsc_b, jpts ) ! variable volume level (vvl) 148 ELSE ; CALL tra_nxt_fix( kt, nit000, 'TRA', tsb, tsn, tsa, jpts ) ! fixed volume level 147 149 ENDIF 148 150 ENDIF … … 241 243 242 244 243 SUBROUTINE tra_nxt_vvl( kt, kit000, cdtype, ptb, ptn, pta, kjpt )245 SUBROUTINE tra_nxt_vvl( kt, kit000, p2dt, cdtype, ptb, ptn, pta, psbc_tc, psbc_tc_b, kjpt ) 244 246 !!---------------------------------------------------------------------- 245 247 !! *** ROUTINE tra_nxt_vvl *** … … 265 267 !! - (ta,sa) time averaged (t,s) (ln_dynhpg_imp = T) 266 268 !!---------------------------------------------------------------------- 267 INTEGER , INTENT(in ) :: kt ! ocean time-step index 268 INTEGER , INTENT(in ) :: kit000 ! first time step index 269 CHARACTER(len=3), INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator) 270 INTEGER , INTENT(in ) :: kjpt ! number of tracers 271 REAL(wp) , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt) :: ptb ! before tracer fields 272 REAL(wp) , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt) :: ptn ! now tracer fields 273 REAL(wp) , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt) :: pta ! tracer trend 269 INTEGER , INTENT(in ) :: kt ! ocean time-step index 270 INTEGER , INTENT(in ) :: kit000 ! first time step index 271 REAL(wp) , INTENT(in ), DIMENSION(jpk) :: p2dt ! time-step 272 CHARACTER(len=3), INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator) 273 INTEGER , INTENT(in ) :: kjpt ! number of tracers 274 REAL(wp) , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt) :: ptb ! before tracer fields 275 REAL(wp) , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt) :: ptn ! now tracer fields 276 REAL(wp) , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt) :: pta ! tracer trend 277 REAL(wp) , INTENT(in ), DIMENSION(jpi,jpj,kjpt) :: psbc_tc ! surface tracer content 278 REAL(wp) , INTENT(in ), DIMENSION(jpi,jpj,kjpt) :: psbc_tc_b ! before surface tracer content 279 274 280 !! 275 LOGICAL :: ll_tra , ll_tra_hpg, ll_traqsr! local logical281 LOGICAL :: ll_tra_hpg, ll_traqsr, ll_rnf ! local logical 276 282 INTEGER :: ji, jj, jk, jn ! dummy loop indices 277 283 REAL(wp) :: zfact1, ztc_a , ztc_n , ztc_b , ztc_f , ztc_d ! local scalar … … 286 292 ! 287 293 IF( cdtype == 'TRA' ) THEN 288 ll_tra = .TRUE. ! active tracers case289 294 ll_tra_hpg = ln_dynhpg_imp ! active tracers case and semi-implicit hpg 290 295 ll_traqsr = ln_traqsr ! active tracers case and solar penetration 296 ll_rnf = ln_rnf ! active tracers case and river runoffs 291 297 ELSE 292 ll_tra = .FALSE. ! passive tracers case293 298 ll_tra_hpg = .FALSE. ! passive tracers case or NO semi-implicit hpg 294 299 ll_traqsr = .FALSE. ! active tracers case and NO solar penetration 300 ll_rnf = .FALSE. ! passive tracers or NO river runoffs 295 301 ENDIF 296 302 ! 297 303 DO jn = 1, kjpt 298 304 DO jk = 1, jpkm1 299 zfact1 = atfp * rdttra(jk)305 zfact1 = atfp * p2dt(jk) 300 306 zfact2 = zfact1 / rau0 301 307 DO jj = 1, jpj … … 315 321 ztc_f = ztc_n + atfp * ztc_d 316 322 ! 317 IF( ll_tra .AND. jk == 1 ) THEN ! first level only for T & S318 ze3t_f = ze3t_f - zfact2 * ( emp_b(ji,jj) - emp(ji,jj) )319 ztc_f = ztc_f - zfact1 * ( sbc_tsc(ji,jj,jn) - sbc_tsc_b(ji,jj,jn) )323 IF( jk == 1 ) THEN ! first level 324 ze3t_f = ze3t_f - zfact2 * ( emp_b(ji,jj) - emp(ji,jj) + rnf(ji,jj) - rnf_b(ji,jj) ) 325 ztc_f = ztc_f - zfact1 * ( psbc_tc(ji,jj,jn) - psbc_tc_b(ji,jj,jn) ) 320 326 ENDIF 327 321 328 IF( ll_traqsr .AND. jn == jp_tem .AND. jk <= nksr ) & ! solar penetration (temperature only) 322 329 & ztc_f = ztc_f - zfact1 * ( qsr_hc(ji,jj,jk) - qsr_hc_b(ji,jj,jk) ) 323 330 324 ze3t_f = 1.e0 / ze3t_f 325 ptb(ji,jj,jk,jn) = ztc_f * ze3t_f ! ptb <-- ptn filtered 326 ptn(ji,jj,jk,jn) = pta(ji,jj,jk,jn) ! ptn <-- pta 327 ! 328 IF( ll_tra_hpg ) THEN ! semi-implicit hpg (T & S only) 329 ze3t_d = 1.e0 / ( ze3t_n + rbcp * ze3t_d ) 330 pta(ji,jj,jk,jn) = ze3t_d * ( ztc_n + rbcp * ztc_d ) ! ta <-- Brown & Campana average 331 ENDIF 331 IF( ll_rnf .AND. jk <= nk_rnf(ji,jj) ) & ! river runoffs 332 & ztc_f = ztc_f - zfact1 * ( rnf_tsc(ji,jj,jn) - rnf_tsc_b(ji,jj,jn) ) & 333 & * fse3t_n(ji,jj,jk) / h_rnf(ji,jj) 334 335 ze3t_f = 1.e0 / ze3t_f 336 ptb(ji,jj,jk,jn) = ztc_f * ze3t_f ! ptb <-- ptn filtered 337 ptn(ji,jj,jk,jn) = pta(ji,jj,jk,jn) ! ptn <-- pta 338 ! 339 IF( ll_tra_hpg ) THEN ! semi-implicit hpg (T & S only) 340 ze3t_d = 1.e0 / ( ze3t_n + rbcp * ze3t_d ) 341 pta(ji,jj,jk,jn) = ze3t_d * ( ztc_n + rbcp * ztc_d ) ! ta <-- Brown & Campana average 342 ENDIF 332 343 END DO 333 344 END DO -
branches/UKMO/dev_r5107_hadgem3_cplseq/NEMOGCM/NEMO/OPA_SRC/TRA/traqsr.F90
r5477 r5572 32 32 USE wrk_nemo ! Memory Allocation 33 33 USE timing ! Timing 34 USE sbc_ice, ONLY : lk_lim335 34 36 35 IMPLICIT NONE … … 38 37 39 38 PUBLIC tra_qsr ! routine called by step.F90 (ln_traqsr=T) 40 PUBLIC tra_qsr_init ! routine called by opa.F9039 PUBLIC tra_qsr_init ! routine called by nemogcm.F90 41 40 42 41 ! !!* Namelist namtra_qsr: penetrative solar radiation … … 50 49 REAL(wp), PUBLIC :: rn_si0 !: very near surface depth of extinction (RGB & 2 bands) 51 50 REAL(wp), PUBLIC :: rn_si1 !: deepest depth of extinction (water type I) (2 bands) 52 51 53 52 ! Module variables 54 53 REAL(wp) :: xsi0r !: inverse of rn_si0 … … 165 164 CALL iom_put( 'qsr3d', etot3 ) ! Shortwave Radiation 3D distribution 166 165 ! clem: store attenuation coefficient of the first ocean level 167 IF ( l k_lim3 .AND. ln_qsr_ice ) THEN166 IF ( ln_qsr_ice ) THEN 168 167 DO jj = 1, jpj 169 168 DO ji = 1, jpi 170 169 IF ( qsr(ji,jj) /= 0._wp ) THEN 171 170 fraqsr_1lev(ji,jj) = ( qsr_hc(ji,jj,1) / ( r1_rau0_rcp * qsr(ji,jj) ) ) 171 ELSE 172 fraqsr_1lev(ji,jj) = 1. 172 173 ENDIF 173 174 END DO … … 233 234 END DO 234 235 ! clem: store attenuation coefficient of the first ocean level 235 IF ( l k_lim3 .AND. ln_qsr_ice ) THEN236 IF ( ln_qsr_ice ) THEN 236 237 DO jj = 1, jpj 237 238 DO ji = 1, jpi … … 256 257 END DO 257 258 ! clem: store attenuation coefficient of the first ocean level 258 IF ( l k_lim3 .AND. ln_qsr_ice ) THEN259 IF ( ln_qsr_ice ) THEN 259 260 fraqsr_1lev(:,:) = etot3(:,:,1) / r1_rau0_rcp 260 261 ENDIF … … 279 280 END DO 280 281 ! clem: store attenuation coefficient of the first ocean level 281 IF ( l k_lim3 .AND. ln_qsr_ice ) THEN282 IF ( ln_qsr_ice ) THEN 282 283 DO jj = 1, jpj 283 284 DO ji = 1, jpi … … 298 299 END DO 299 300 ! clem: store attenuation coefficient of the first ocean level 300 IF ( l k_lim3 .AND. ln_qsr_ice ) THEN301 IF ( ln_qsr_ice ) THEN 301 302 fraqsr_1lev(:,:) = etot3(:,:,1) / r1_rau0_rcp 302 303 ENDIF … … 324 325 & 'at it= ', kt,' date= ', ndastp 325 326 IF(lwp) WRITE(numout,*) '~~~~' 326 CALL iom_rstput( kt, nitrst, numrow, 'qsr_hc_b', qsr_hc ) 327 CALL iom_rstput( kt, nitrst, numrow, 'qsr_hc_b' , qsr_hc ) 328 CALL iom_rstput( kt, nitrst, numrow, 'fraqsr_1lev', fraqsr_1lev ) ! default definition in sbcssm 327 329 ! 328 330 ENDIF … … 379 381 ! 380 382 IF( nn_timing == 1 ) CALL timing_start('tra_qsr_init') 381 !382 ! Default value for fraqsr_1lev383 IF( .NOT. ln_rstart ) THEN384 fraqsr_1lev(:,:) = 1._wp385 ENDIF386 383 ! 387 384 CALL wrk_alloc( jpi, jpj, zekb, zekg, zekr ) … … 412 409 WRITE(numout,*) ' RGB & 2 bands: shortess depth of extinction rn_si0 = ', rn_si0 413 410 WRITE(numout,*) ' 2 bands: longest depth of extinction rn_si1 = ', rn_si1 414 WRITE(numout,*) ' light penetration for ice-model LIM3 ln_qsr_ice = ', ln_qsr_ice415 411 ENDIF 416 412 … … 564 560 ENDIF 565 561 ! 562 ! initialisation of fraqsr_1lev used in sbcssm 563 IF( iom_varid( numror, 'fraqsr_1lev', ldstop = .FALSE. ) > 0 ) THEN 564 CALL iom_get( numror, jpdom_autoglo, 'fraqsr_1lev' , fraqsr_1lev ) 565 ELSE 566 fraqsr_1lev(:,:) = 1._wp ! default definition 567 ENDIF 568 ! 566 569 CALL wrk_dealloc( jpi, jpj, zekb, zekg, zekr ) 567 570 CALL wrk_dealloc( jpi, jpj, jpk, ze0, ze1, ze2, ze3, zea ) -
branches/UKMO/dev_r5107_hadgem3_cplseq/NEMOGCM/NEMO/OPA_SRC/TRA/trasbc.F90
r5477 r5572 9 9 !! 3.3 ! 2010-04 (M. Leclair, G. Madec) Forcing averaged over 2 time steps 10 10 !! - ! 2010-09 (C. Ethe, G. Madec) Merge TRA-TRC 11 !! 3.6 ! 2014-11 (P. Mathiot) isf melting forcing 11 12 !!---------------------------------------------------------------------- 12 13 … … 20 21 USE sbcmod ! ln_rnf 21 22 USE sbcrnf ! River runoff 23 USE sbcisf ! Ice shelf 22 24 USE traqsr ! solar radiation penetration 23 25 USE trd_oce ! trends: ocean variables … … 26 28 USE in_out_manager ! I/O manager 27 29 USE prtctl ! Print control 28 USE sbcrnf ! River runoff29 USE sbcisf ! Ice shelf30 USE sbcmod ! ln_rnf31 30 USE iom 32 31 USE lbclnk ! ocean lateral boundary conditions (or mpp link) -
branches/UKMO/dev_r5107_hadgem3_cplseq/NEMOGCM/NEMO/OPA_SRC/TRA/trazdf.F90
r5477 r5572 88 88 & tab3d_2=tsa(:,:,:,jp_sal), clinfo2= ' Sa: ', mask2=tmask, clinfo3='tra' ) 89 89 END SELECT 90 ! DRAKKAR SSS control { 91 ! JMM avoid negative salinities near river outlet ! Ugly fix 92 ! JMM : restore negative salinities to small salinities: 93 WHERE ( tsa(:,:,:,jp_sal) < 0._wp ) tsa(:,:,:,jp_sal) = 0.1_wp 90 94 91 95 IF( l_trdtra ) THEN ! save the vertical diffusive trends for further diagnostics -
branches/UKMO/dev_r5107_hadgem3_cplseq/NEMOGCM/NEMO/OPA_SRC/TRA/trazdf_imp.F90
r5477 r5572 122 122 DO jj=1, jpj 123 123 DO ji=1, jpi 124 zwt(ji,jj,1 :mikt(ji,jj)) = 0._wp124 zwt(ji,jj,1) = 0._wp 125 125 END DO 126 126 END DO … … 184 184 DO jj = 2, jpjm1 185 185 DO ji = fs_2, fs_jpim1 186 zwt(ji,jj,1:mikt(ji,jj)) = zwd(ji,jj,1:mikt(ji,jj)) 187 DO jk = mikt(ji,jj)+1, jpkm1 186 zwt(ji,jj,1) = zwd(ji,jj,1) 187 END DO 188 END DO 189 DO jk = 2, jpkm1 190 DO jj = 2, jpjm1 191 DO ji = fs_2, fs_jpim1 188 192 zwt(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwt(ji,jj,jk-1) 189 193 END DO … … 196 200 DO jj = 2, jpjm1 197 201 DO ji = fs_2, fs_jpim1 198 ze3tb = ( 1. - r_vvl ) + r_vvl * fse3t_b(ji,jj,mikt(ji,jj)) 199 ze3tn = ( 1. - r_vvl ) + r_vvl * fse3t(ji,jj,mikt(ji,jj)) 200 pta(ji,jj,mikt(ji,jj),jn) = ze3tb * ptb(ji,jj,mikt(ji,jj),jn) & 201 & + p2dt(mikt(ji,jj)) * ze3tn * pta(ji,jj,mikt(ji,jj),jn) 202 DO jk = mikt(ji,jj)+1, jpkm1 202 ze3tb = ( 1. - r_vvl ) + r_vvl * fse3t_b(ji,jj,1) 203 ze3tn = ( 1. - r_vvl ) + r_vvl * fse3t(ji,jj,1) 204 pta(ji,jj,1,jn) = ze3tb * ptb(ji,jj,1,jn) & 205 & + p2dt(1) * ze3tn * pta(ji,jj,1,jn) 206 END DO 207 END DO 208 DO jk = 2, jpkm1 209 DO jj = 2, jpjm1 210 DO ji = fs_2, fs_jpim1 203 211 ze3tb = ( 1. - r_vvl ) + r_vvl * fse3t_b(ji,jj,jk) 204 212 ze3tn = ( 1. - r_vvl ) + r_vvl * fse3t (ji,jj,jk) … … 213 221 DO ji = fs_2, fs_jpim1 214 222 pta(ji,jj,jpkm1,jn) = pta(ji,jj,jpkm1,jn) / zwt(ji,jj,jpkm1) * tmask(ji,jj,jpkm1) 215 DO jk = jpk-2, mikt(ji,jj), -1 223 END DO 224 END DO 225 DO jk = jpk-2, 1, -1 226 DO jj = 2, jpjm1 227 DO ji = fs_2, fs_jpim1 216 228 pta(ji,jj,jk,jn) = ( pta(ji,jj,jk,jn) - zws(ji,jj,jk) * pta(ji,jj,jk+1,jn) ) & 217 229 & / zwt(ji,jj,jk) * tmask(ji,jj,jk) -
branches/UKMO/dev_r5107_hadgem3_cplseq/NEMOGCM/NEMO/OPA_SRC/TRA/zpshde.F90
r5477 r5572 8 8 !! - ! 2004-03 (C. Ethe) adapted for passive tracers 9 9 !! 3.3 ! 2010-05 (C. Ethe, G. Madec) merge TRC-TRA 10 !! 3.6 ! 2014-11 (P. Mathiot) Add zps_hde_isf (needed to open a cavity) 10 11 !!====================================================================== 11 12 … … 27 28 PRIVATE 28 29 29 PUBLIC zps_hde ! routine called by step.F90 30 PUBLIC zps_hde ! routine called by step.F90 31 PUBLIC zps_hde_isf ! routine called by step.F90 30 32 31 33 !! * Substitutions … … 40 42 41 43 SUBROUTINE zps_hde( kt, kjpt, pta, pgtu, pgtv, & 44 & prd, pgru, pgrv ) 45 !!---------------------------------------------------------------------- 46 !! *** ROUTINE zps_hde *** 47 !! 48 !! ** Purpose : Compute the horizontal derivative of T, S and rho 49 !! at u- and v-points with a linear interpolation for z-coordinate 50 !! with partial steps. 51 !! 52 !! ** Method : In z-coord with partial steps, scale factors on last 53 !! levels are different for each grid point, so that T, S and rd 54 !! points are not at the same depth as in z-coord. To have horizontal 55 !! gradients again, we interpolate T and S at the good depth : 56 !! Linear interpolation of T, S 57 !! Computation of di(tb) and dj(tb) by vertical interpolation: 58 !! di(t) = t~ - t(i,j,k) or t(i+1,j,k) - t~ 59 !! dj(t) = t~ - t(i,j,k) or t(i,j+1,k) - t~ 60 !! This formulation computes the two cases: 61 !! CASE 1 CASE 2 62 !! k-1 ___ ___________ k-1 ___ ___________ 63 !! Ti T~ T~ Ti+1 64 !! _____ _____ 65 !! k | |Ti+1 k Ti | | 66 !! | |____ ____| | 67 !! ___ | | | ___ | | | 68 !! 69 !! case 1-> e3w(i+1) >= e3w(i) ( and e3w(j+1) >= e3w(j) ) then 70 !! t~ = t(i+1,j ,k) + (e3w(i+1) - e3w(i)) * dk(Ti+1)/e3w(i+1) 71 !! ( t~ = t(i ,j+1,k) + (e3w(j+1) - e3w(j)) * dk(Tj+1)/e3w(j+1) ) 72 !! or 73 !! case 2-> e3w(i+1) <= e3w(i) ( and e3w(j+1) <= e3w(j) ) then 74 !! t~ = t(i,j,k) + (e3w(i) - e3w(i+1)) * dk(Ti)/e3w(i ) 75 !! ( t~ = t(i,j,k) + (e3w(j) - e3w(j+1)) * dk(Tj)/e3w(j ) ) 76 !! Idem for di(s) and dj(s) 77 !! 78 !! For rho, we call eos which will compute rd~(t~,s~) at the right 79 !! depth zh from interpolated T and S for the different formulations 80 !! of the equation of state (eos). 81 !! Gradient formulation for rho : 82 !! di(rho) = rd~ - rd(i,j,k) or rd(i+1,j,k) - rd~ 83 !! 84 !! ** Action : compute for top interfaces 85 !! - pgtu, pgtv: horizontal gradient of tracer at u- & v-points 86 !! - pgru, pgrv: horizontal gradient of rho (if present) at u- & v-points 87 !!---------------------------------------------------------------------- 88 INTEGER , INTENT(in ) :: kt ! ocean time-step index 89 INTEGER , INTENT(in ) :: kjpt ! number of tracers 90 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: pta ! 4D tracers fields 91 REAL(wp), DIMENSION(jpi,jpj, kjpt), INTENT( out) :: pgtu, pgtv ! hor. grad. of ptra at u- & v-pts 92 REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ), OPTIONAL :: prd ! 3D density anomaly fields 93 REAL(wp), DIMENSION(jpi,jpj ), INTENT( out), OPTIONAL :: pgru, pgrv ! hor. grad of prd at u- & v-pts (bottom) 94 ! 95 INTEGER :: ji, jj, jn ! Dummy loop indices 96 INTEGER :: iku, ikv, ikum1, ikvm1 ! partial step level (ocean bottom level) at u- and v-points 97 REAL(wp) :: ze3wu, ze3wv, zmaxu, zmaxv ! temporary scalars 98 REAL(wp), DIMENSION(jpi,jpj) :: zri, zrj, zhi, zhj ! NB: 3rd dim=1 to use eos 99 REAL(wp), DIMENSION(jpi,jpj,kjpt) :: zti, ztj ! 100 !!---------------------------------------------------------------------- 101 ! 102 IF( nn_timing == 1 ) CALL timing_start( 'zps_hde') 103 ! 104 pgtu(:,:,:)=0.0_wp ; pgtv(:,:,:)=0.0_wp ; 105 zti (:,:,:)=0.0_wp ; ztj (:,:,:)=0.0_wp ; 106 zhi (:,: )=0.0_wp ; zhj (:,: )=0.0_wp ; 107 ! 108 DO jn = 1, kjpt !== Interpolation of tracers at the last ocean level ==! 109 ! 110 DO jj = 1, jpjm1 111 DO ji = 1, jpim1 112 iku = mbku(ji,jj) ; ikum1 = MAX( iku - 1 , 1 ) ! last and before last ocean level at u- & v-points 113 ikv = mbkv(ji,jj) ; ikvm1 = MAX( ikv - 1 , 1 ) ! if level first is a p-step, ik.m1=1 114 ze3wu = fse3w(ji+1,jj ,iku) - fse3w(ji,jj,iku) 115 ze3wv = fse3w(ji ,jj+1,ikv) - fse3w(ji,jj,ikv) 116 ! 117 ! i- direction 118 IF( ze3wu >= 0._wp ) THEN ! case 1 119 zmaxu = ze3wu / fse3w(ji+1,jj,iku) 120 ! interpolated values of tracers 121 zti (ji,jj,jn) = pta(ji+1,jj,iku,jn) + zmaxu * ( pta(ji+1,jj,ikum1,jn) - pta(ji+1,jj,iku,jn) ) 122 ! gradient of tracers 123 pgtu(ji,jj,jn) = umask(ji,jj,1) * ( zti(ji,jj,jn) - pta(ji,jj,iku,jn) ) 124 ELSE ! case 2 125 zmaxu = -ze3wu / fse3w(ji,jj,iku) 126 ! interpolated values of tracers 127 zti (ji,jj,jn) = pta(ji,jj,iku,jn) + zmaxu * ( pta(ji,jj,ikum1,jn) - pta(ji,jj,iku,jn) ) 128 ! gradient of tracers 129 pgtu(ji,jj,jn) = umask(ji,jj,1) * ( pta(ji+1,jj,iku,jn) - zti(ji,jj,jn) ) 130 ENDIF 131 ! 132 ! j- direction 133 IF( ze3wv >= 0._wp ) THEN ! case 1 134 zmaxv = ze3wv / fse3w(ji,jj+1,ikv) 135 ! interpolated values of tracers 136 ztj (ji,jj,jn) = pta(ji,jj+1,ikv,jn) + zmaxv * ( pta(ji,jj+1,ikvm1,jn) - pta(ji,jj+1,ikv,jn) ) 137 ! gradient of tracers 138 pgtv(ji,jj,jn) = vmask(ji,jj,1) * ( ztj(ji,jj,jn) - pta(ji,jj,ikv,jn) ) 139 ELSE ! case 2 140 zmaxv = -ze3wv / fse3w(ji,jj,ikv) 141 ! interpolated values of tracers 142 ztj (ji,jj,jn) = pta(ji,jj,ikv,jn) + zmaxv * ( pta(ji,jj,ikvm1,jn) - pta(ji,jj,ikv,jn) ) 143 ! gradient of tracers 144 pgtv(ji,jj,jn) = vmask(ji,jj,1) * ( pta(ji,jj+1,ikv,jn) - ztj(ji,jj,jn) ) 145 ENDIF 146 END DO 147 END DO 148 CALL lbc_lnk( pgtu(:,:,jn), 'U', -1. ) ; CALL lbc_lnk( pgtv(:,:,jn), 'V', -1. ) ! Lateral boundary cond. 149 ! 150 END DO 151 152 ! horizontal derivative of density anomalies (rd) 153 IF( PRESENT( prd ) ) THEN ! depth of the partial step level 154 pgru(:,:)=0.0_wp ; pgrv(:,:)=0.0_wp ; 155 DO jj = 1, jpjm1 156 DO ji = 1, jpim1 157 iku = mbku(ji,jj) 158 ikv = mbkv(ji,jj) 159 ze3wu = fse3w(ji+1,jj ,iku) - fse3w(ji,jj,iku) 160 ze3wv = fse3w(ji ,jj+1,ikv) - fse3w(ji,jj,ikv) 161 IF( ze3wu >= 0._wp ) THEN ; zhi(ji,jj) = fsdept(ji ,jj,iku) ! i-direction: case 1 162 ELSE ; zhi(ji,jj) = fsdept(ji+1,jj,iku) ! - - case 2 163 ENDIF 164 IF( ze3wv >= 0._wp ) THEN ; zhj(ji,jj) = fsdept(ji,jj ,ikv) ! j-direction: case 1 165 ELSE ; zhj(ji,jj) = fsdept(ji,jj+1,ikv) ! - - case 2 166 ENDIF 167 END DO 168 END DO 169 170 ! Compute interpolated rd from zti, ztj for the 2 cases at the depth of the partial 171 ! step and store it in zri, zrj for each case 172 CALL eos( zti, zhi, zri ) 173 CALL eos( ztj, zhj, zrj ) 174 175 ! Gradient of density at the last level 176 DO jj = 1, jpjm1 177 DO ji = 1, jpim1 178 iku = mbku(ji,jj) 179 ikv = mbkv(ji,jj) 180 ze3wu = fse3w(ji+1,jj ,iku) - fse3w(ji,jj,iku) 181 ze3wv = fse3w(ji ,jj+1,ikv) - fse3w(ji,jj,ikv) 182 IF( ze3wu >= 0._wp ) THEN ; pgru(ji,jj) = umask(ji,jj,1) * ( zri(ji ,jj ) - prd(ji,jj,iku) ) ! i: 1 183 ELSE ; pgru(ji,jj) = umask(ji,jj,1) * ( prd(ji+1,jj,iku) - zri(ji,jj ) ) ! i: 2 184 ENDIF 185 IF( ze3wv >= 0._wp ) THEN ; pgrv(ji,jj) = vmask(ji,jj,1) * ( zrj(ji,jj ) - prd(ji,jj,ikv) ) ! j: 1 186 ELSE ; pgrv(ji,jj) = vmask(ji,jj,1) * ( prd(ji,jj+1,ikv) - zrj(ji,jj ) ) ! j: 2 187 ENDIF 188 END DO 189 END DO 190 CALL lbc_lnk( pgru , 'U', -1. ) ; CALL lbc_lnk( pgrv , 'V', -1. ) ! Lateral boundary conditions 191 ! 192 END IF 193 ! 194 IF( nn_timing == 1 ) CALL timing_stop( 'zps_hde') 195 ! 196 END SUBROUTINE zps_hde 197 ! 198 SUBROUTINE zps_hde_isf( kt, kjpt, pta, pgtu, pgtv, & 42 199 & prd, pgru, pgrv, pmru, pmrv, pgzu, pgzv, pge3ru, pge3rv, & 43 & sgtu, sgtv, sgru, sgrv, smru, smrv, sgzu, sgzv, sge3ru, sge3rv)200 & pgtui, pgtvi, pgrui, pgrvi, pmrui, pmrvi, pgzui, pgzvi, pge3rui, pge3rvi ) 44 201 !!---------------------------------------------------------------------- 45 202 !! *** ROUTINE zps_hde *** … … 82 239 !! 83 240 !! ** Action : compute for top and bottom interfaces 84 !! - pgtu, pgtv, sgtu, sgtv: horizontal gradient of tracer at u- & v-points85 !! - pgru, pgrv, sgru, sgtv: horizontal gradient of rho (if present) at u- & v-points86 !! - pmru, pmrv, smru, smrv: horizontal sum of rho at u- & v- point (used in dynhpg with vvl)87 !! - pgzu, pgzv, sgzu, sgzv: horizontal gradient of z at u- and v- point (used in dynhpg with vvl)88 !! - pge3ru, pge3rv, sge3ru, sge3rv: horizontal gradient of rho weighted by local e3w at u- & v-points241 !! - pgtu, pgtv, pgtui, pgtvi: horizontal gradient of tracer at u- & v-points 242 !! - pgru, pgrv, pgrui, pgtvi: horizontal gradient of rho (if present) at u- & v-points 243 !! - pmru, pmrv, pmrui, pmrvi: horizontal sum of rho at u- & v- point (used in dynhpg with vvl) 244 !! - pgzu, pgzv, pgzui, pgzvi: horizontal gradient of z at u- and v- point (used in dynhpg with vvl) 245 !! - pge3ru, pge3rv, pge3rui, pge3rvi: horizontal gradient of rho weighted by local e3w at u- & v-points 89 246 !!---------------------------------------------------------------------- 90 247 INTEGER , INTENT(in ) :: kt ! ocean time-step index … … 92 249 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: pta ! 4D tracers fields 93 250 REAL(wp), DIMENSION(jpi,jpj, kjpt), INTENT( out) :: pgtu, pgtv ! hor. grad. of ptra at u- & v-pts 94 REAL(wp), DIMENSION(jpi,jpj, kjpt), INTENT( out) :: sgtu, sgtv! hor. grad. of stra at u- & v-pts (ISF)251 REAL(wp), DIMENSION(jpi,jpj, kjpt), INTENT( out) :: pgtui, pgtvi ! hor. grad. of stra at u- & v-pts (ISF) 95 252 REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ), OPTIONAL :: prd ! 3D density anomaly fields 96 253 REAL(wp), DIMENSION(jpi,jpj ), INTENT( out), OPTIONAL :: pgru, pgrv ! hor. grad of prd at u- & v-pts (bottom) … … 98 255 REAL(wp), DIMENSION(jpi,jpj ), INTENT( out), OPTIONAL :: pgzu, pgzv ! hor. grad of z at u- & v-pts (bottom) 99 256 REAL(wp), DIMENSION(jpi,jpj ), INTENT( out), OPTIONAL :: pge3ru, pge3rv ! hor. grad of prd weighted by local e3w at u- & v-pts (bottom) 100 REAL(wp), DIMENSION(jpi,jpj ), INTENT( out), OPTIONAL :: sgru, sgrv! hor. grad of prd at u- & v-pts (top)101 REAL(wp), DIMENSION(jpi,jpj ), INTENT( out), OPTIONAL :: smru, smrv! hor. sum of prd at u- & v-pts (top)102 REAL(wp), DIMENSION(jpi,jpj ), INTENT( out), OPTIONAL :: sgzu, sgzv! hor. grad of z at u- & v-pts (top)103 REAL(wp), DIMENSION(jpi,jpj ), INTENT( out), OPTIONAL :: sge3ru, sge3rv! hor. grad of prd weighted by local e3w at u- & v-pts (top)257 REAL(wp), DIMENSION(jpi,jpj ), INTENT( out), OPTIONAL :: pgrui, pgrvi ! hor. grad of prd at u- & v-pts (top) 258 REAL(wp), DIMENSION(jpi,jpj ), INTENT( out), OPTIONAL :: pmrui, pmrvi ! hor. sum of prd at u- & v-pts (top) 259 REAL(wp), DIMENSION(jpi,jpj ), INTENT( out), OPTIONAL :: pgzui, pgzvi ! hor. grad of z at u- & v-pts (top) 260 REAL(wp), DIMENSION(jpi,jpj ), INTENT( out), OPTIONAL :: pge3rui, pge3rvi ! hor. grad of prd weighted by local e3w at u- & v-pts (top) 104 261 ! 105 262 INTEGER :: ji, jj, jn ! Dummy loop indices … … 110 267 !!---------------------------------------------------------------------- 111 268 ! 112 IF( nn_timing == 1 ) CALL timing_start( 'zps_hde ')269 IF( nn_timing == 1 ) CALL timing_start( 'zps_hde_isf') 113 270 ! 114 271 pgtu(:,:,:)=0.0_wp ; pgtv(:,:,:)=0.0_wp ; 115 sgtu(:,:,:)=0.0_wp ; sgtv(:,:,:)=0.0_wp ;272 pgtui(:,:,:)=0.0_wp ; pgtvi(:,:,:)=0.0_wp ; 116 273 zti (:,:,:)=0.0_wp ; ztj (:,:,:)=0.0_wp ; 117 274 zhi (:,: )=0.0_wp ; zhj (:,: )=0.0_wp ; … … 256 413 zti(ji,jj,jn) = pta(ji+1,jj,iku,jn) + zmaxu * ( pta(ji+1,jj,iku+1,jn) - pta(ji+1,jj,iku,jn) ) 257 414 ! gradient of tracers 258 sgtu(ji,jj,jn) = umask(ji,jj,iku) * ( zti(ji,jj,jn) - pta(ji,jj,iku,jn) )415 pgtui(ji,jj,jn) = umask(ji,jj,iku) * ( zti(ji,jj,jn) - pta(ji,jj,iku,jn) ) 259 416 ELSE ! case 2 260 417 zmaxu = - ze3wu / fse3w(ji,jj,iku+1) … … 262 419 zti(ji,jj,jn) = pta(ji,jj,iku,jn) + zmaxu * ( pta(ji,jj,iku+1,jn) - pta(ji,jj,iku,jn) ) 263 420 ! gradient of tracers 264 sgtu(ji,jj,jn) = umask(ji,jj,iku) * ( pta(ji+1,jj,iku,jn) - zti(ji,jj,jn) )421 pgtui(ji,jj,jn) = umask(ji,jj,iku) * ( pta(ji+1,jj,iku,jn) - zti(ji,jj,jn) ) 265 422 ENDIF 266 423 ! … … 271 428 ztj(ji,jj,jn) = pta(ji,jj+1,ikv,jn) + zmaxv * ( pta(ji,jj+1,ikv+1,jn) - pta(ji,jj+1,ikv,jn) ) 272 429 ! gradient of tracers 273 sgtv(ji,jj,jn) = vmask(ji,jj,ikv) * ( ztj(ji,jj,jn) - pta(ji,jj,ikv,jn) )430 pgtvi(ji,jj,jn) = vmask(ji,jj,ikv) * ( ztj(ji,jj,jn) - pta(ji,jj,ikv,jn) ) 274 431 ELSE ! case 2 275 432 zmaxv = - ze3wv / fse3w(ji,jj,ikv+1) … … 277 434 ztj(ji,jj,jn) = pta(ji,jj,ikv,jn) + zmaxv * ( pta(ji,jj,ikv+1,jn) - pta(ji,jj,ikv,jn) ) 278 435 ! gradient of tracers 279 sgtv(ji,jj,jn) = vmask(ji,jj,ikv) * ( pta(ji,jj+1,ikv,jn) - ztj(ji,jj,jn) )436 pgtvi(ji,jj,jn) = vmask(ji,jj,ikv) * ( pta(ji,jj+1,ikv,jn) - ztj(ji,jj,jn) ) 280 437 ENDIF 281 438 END DO!! 282 439 END DO!! 283 CALL lbc_lnk( sgtu(:,:,jn), 'U', -1. ) ; CALL lbc_lnk( sgtv(:,:,jn), 'V', -1. ) ! Lateral boundary cond.440 CALL lbc_lnk( pgtui(:,:,jn), 'U', -1. ) ; CALL lbc_lnk( pgtvi(:,:,jn), 'V', -1. ) ! Lateral boundary cond. 284 441 ! 285 442 END DO … … 287 444 ! horizontal derivative of density anomalies (rd) 288 445 IF( PRESENT( prd ) ) THEN ! depth of the partial step level 289 sgru(:,:) =0.0_wp ; sgrv(:,:) =0.0_wp ;290 sgzu(:,:) =0.0_wp ; sgzv(:,:) =0.0_wp ;291 smru(:,:) =0.0_wp ; smru(:,:) =0.0_wp ;292 sge3ru(:,:)=0.0_wp ; sge3rv(:,:)=0.0_wp ;446 pgrui(:,:) =0.0_wp ; pgrvi(:,:) =0.0_wp ; 447 pgzui(:,:) =0.0_wp ; pgzvi(:,:) =0.0_wp ; 448 pmrui(:,:) =0.0_wp ; pmrui(:,:) =0.0_wp ; 449 pge3rui(:,:)=0.0_wp ; pge3rvi(:,:)=0.0_wp ; 293 450 294 451 DO jj = 1, jpjm1 … … 321 478 ze3wv = (gdepw_0(ji,jj+1,ikv+1) - gdept_0(ji,jj+1,ikv)) - (gdepw_0(ji,jj,ikv+1) - gdept_0(ji,jj,ikv)) 322 479 IF( ze3wu >= 0._wp ) THEN 323 sgzu(ji,jj) = (fsde3w(ji+1,jj,iku) + ze3wu) - fsde3w(ji,jj,iku)324 sgru(ji,jj) = umask(ji,jj,iku) * ( zri(ji,jj) - prd(ji,jj,iku) ) ! i: 1325 smru(ji,jj) = umask(ji,jj,iku) * ( zri(ji,jj) + prd(ji,jj,iku) ) ! i: 1326 sge3ru(ji,jj) = umask(ji,jj,iku+1) &480 pgzui (ji,jj) = (fsde3w(ji+1,jj,iku) + ze3wu) - fsde3w(ji,jj,iku) 481 pgrui (ji,jj) = umask(ji,jj,iku) * ( zri(ji,jj) - prd(ji,jj,iku) ) ! i: 1 482 pmrui (ji,jj) = umask(ji,jj,iku) * ( zri(ji,jj) + prd(ji,jj,iku) ) ! i: 1 483 pge3rui(ji,jj) = umask(ji,jj,iku+1) & 327 484 * ( (fse3w(ji+1,jj,iku+1) - ze3wu) * (zri(ji,jj ) + prd(ji+1,jj,iku+1) + 2._wp) & 328 485 - fse3w(ji ,jj,iku+1) * (prd(ji,jj,iku) + prd(ji ,jj,iku+1) + 2._wp) ) ! i: 1 329 486 ELSE 330 sgzu(ji,jj) = fsde3w(ji+1,jj,iku) - (fsde3w(ji,jj,iku) - ze3wu)331 sgru(ji,jj) = umask(ji,jj,iku) * ( prd(ji+1,jj,iku) - zri(ji,jj) ) ! i: 2332 smru(ji,jj) = umask(ji,jj,iku) * ( prd(ji+1,jj,iku) + zri(ji,jj) ) ! i: 2333 sge3ru(ji,jj) = umask(ji,jj,iku+1) &487 pgzui (ji,jj) = fsde3w(ji+1,jj,iku) - (fsde3w(ji,jj,iku) - ze3wu) 488 pgrui (ji,jj) = umask(ji,jj,iku) * ( prd(ji+1,jj,iku) - zri(ji,jj) ) ! i: 2 489 pmrui (ji,jj) = umask(ji,jj,iku) * ( prd(ji+1,jj,iku) + zri(ji,jj) ) ! i: 2 490 pge3rui(ji,jj) = umask(ji,jj,iku+1) & 334 491 * ( fse3w(ji+1,jj,iku+1) * (prd(ji+1,jj,iku) + prd(ji+1,jj,iku+1) + 2._wp) & 335 492 -(fse3w(ji ,jj,iku+1) + ze3wu) * (zri(ji,jj ) + prd(ji ,jj,iku+1) + 2._wp) ) ! i: 2 336 493 ENDIF 337 494 IF( ze3wv >= 0._wp ) THEN 338 sgzv(ji,jj) = (fsde3w(ji,jj+1,ikv) + ze3wv) - fsde3w(ji,jj,ikv)339 sgrv(ji,jj) = vmask(ji,jj,ikv) * ( zrj(ji,jj ) - prd(ji,jj,ikv) ) ! j: 1340 smrv(ji,jj) = vmask(ji,jj,ikv) * ( zrj(ji,jj ) + prd(ji,jj,ikv) ) ! j: 1341 sge3rv(ji,jj) = vmask(ji,jj,ikv+1) &495 pgzvi (ji,jj) = (fsde3w(ji,jj+1,ikv) + ze3wv) - fsde3w(ji,jj,ikv) 496 pgrvi (ji,jj) = vmask(ji,jj,ikv) * ( zrj(ji,jj ) - prd(ji,jj,ikv) ) ! j: 1 497 pmrvi (ji,jj) = vmask(ji,jj,ikv) * ( zrj(ji,jj ) + prd(ji,jj,ikv) ) ! j: 1 498 pge3rvi(ji,jj) = vmask(ji,jj,ikv+1) & 342 499 * ( (fse3w(ji,jj+1,ikv+1) - ze3wv) * ( zrj(ji,jj ) + prd(ji,jj+1,ikv+1) + 2._wp) & 343 500 - fse3w(ji,jj ,ikv+1) * ( prd(ji,jj,ikv) + prd(ji,jj ,ikv+1) + 2._wp) ) ! j: 1 344 501 ! + 2 due to the formulation in density and not in anomalie in hpg sco 345 502 ELSE 346 sgzv(ji,jj) = fsde3w(ji,jj+1,ikv) - (fsde3w(ji,jj,ikv) - ze3wv)347 sgrv(ji,jj) = vmask(ji,jj,ikv) * ( prd(ji,jj+1,ikv) - zrj(ji,jj) ) ! j: 2348 smrv(ji,jj) = vmask(ji,jj,ikv) * ( prd(ji,jj+1,ikv) + zrj(ji,jj) ) ! j: 2349 sge3rv(ji,jj) = vmask(ji,jj,ikv+1) &503 pgzvi (ji,jj) = fsde3w(ji,jj+1,ikv) - (fsde3w(ji,jj,ikv) - ze3wv) 504 pgrvi (ji,jj) = vmask(ji,jj,ikv) * ( prd(ji,jj+1,ikv) - zrj(ji,jj) ) ! j: 2 505 pmrvi (ji,jj) = vmask(ji,jj,ikv) * ( prd(ji,jj+1,ikv) + zrj(ji,jj) ) ! j: 2 506 pge3rvi(ji,jj) = vmask(ji,jj,ikv+1) & 350 507 * ( fse3w(ji,jj+1,ikv+1) * ( prd(ji,jj+1,ikv) + prd(ji,jj+1,ikv+1) + 2._wp) & 351 508 -(fse3w(ji,jj ,ikv+1) + ze3wv) * ( zrj(ji,jj ) + prd(ji,jj ,ikv+1) + 2._wp) ) ! j: 2 … … 353 510 END DO 354 511 END DO 355 CALL lbc_lnk( sgru , 'U', -1. ) ; CALL lbc_lnk( sgrv, 'V', -1. ) ! Lateral boundary conditions356 CALL lbc_lnk( smru , 'U', 1. ) ; CALL lbc_lnk( smrv, 'V', 1. ) ! Lateral boundary conditions357 CALL lbc_lnk( sgzu , 'U', -1. ) ; CALL lbc_lnk( sgzv, 'V', -1. ) ! Lateral boundary conditions358 CALL lbc_lnk( sge3ru , 'U', -1. ) ; CALL lbc_lnk( sge3rv, 'V', -1. ) ! Lateral boundary conditions512 CALL lbc_lnk( pgrui , 'U', -1. ) ; CALL lbc_lnk( pgrvi , 'V', -1. ) ! Lateral boundary conditions 513 CALL lbc_lnk( pmrui , 'U', 1. ) ; CALL lbc_lnk( pmrvi , 'V', 1. ) ! Lateral boundary conditions 514 CALL lbc_lnk( pgzui , 'U', -1. ) ; CALL lbc_lnk( pgzvi , 'V', -1. ) ! Lateral boundary conditions 515 CALL lbc_lnk( pge3rui , 'U', -1. ) ; CALL lbc_lnk( pge3rvi , 'V', -1. ) ! Lateral boundary conditions 359 516 ! 360 517 END IF 361 518 ! 362 IF( nn_timing == 1 ) CALL timing_stop( 'zps_hde') 363 ! 364 END SUBROUTINE zps_hde 365 519 IF( nn_timing == 1 ) CALL timing_stop( 'zps_hde_isf') 520 ! 521 END SUBROUTINE zps_hde_isf 366 522 !!====================================================================== 367 523 END MODULE zpshde
Note: See TracChangeset
for help on using the changeset viewer.
