- Timestamp:
- 2017-12-13T15:58:53+01:00 (6 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90
r7831 r9019 16 16 !! 3.7 ! 2015-11 (J. Chanut) free surface simplification 17 17 !! - ! 2016-12 (G. Madec, E. Clementi) update for Stoke-Drift divergence 18 !! 4.0 ! 2017-05 (G. Madec) drag coef. defined at t-point (zdfdrg.F90) 18 19 !!--------------------------------------------------------------------- 19 20 … … 27 28 USE dom_oce ! ocean space and time domain 28 29 USE sbc_oce ! surface boundary condition: ocean 29 USE zdf_oce ! Bottom friction coefts 30 USE zdf_oce ! vertical physics: variables 31 USE zdfdrg ! vertical physics: top/bottom drag coef. 30 32 USE sbcisf ! ice shelf variable (fwfisf) 31 33 USE sbcapr ! surface boundary condition: atmospheric pressure … … 40 42 USE updtide ! tide potential 41 43 USE sbcwave ! surface wave 44 USE diatmb ! Top,middle,bottom output 45 #if defined key_agrif 46 USE agrif_opa_interp ! agrif 47 #endif 48 #if defined key_asminc 49 USE asminc ! Assimilation increment 50 #endif 42 51 ! 43 52 USE in_out_manager ! I/O manager … … 47 56 USE iom ! IOM library 48 57 USE restart ! only for lrst_oce 49 USE wrk_nemo ! Memory Allocation50 58 USE timing ! Timing 51 USE diatmb ! Top,middle,bottom output52 #if defined key_agrif53 USE agrif_opa_interp ! agrif54 #endif55 #if defined key_asminc56 USE asminc ! Assimilation increment57 #endif58 59 59 60 60 IMPLICIT NONE … … 66 66 PUBLIC ts_rst ! " " " " 67 67 68 INTEGER, SAVE :: icycle ! Number of barotropic sub-steps for each internal step nn_baro <= 2.5 nn_baro69 REAL(wp),SAVE :: rdtbt ! Barotropic time step70 71 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:) :: wgtbtp1, wgtbtp2 !: 1st & 2nd weights used in time filtering of barotropic fields72 73 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: zwz !: ff_f/h at F points74 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: ftnw, ftne !: triad of coriolis parameter75 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: ftsw, ftse !: (only used with een vorticity scheme)76 77 68 !! Time filtered arrays at baroclinic time step: 78 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: un_adv , vn_adv !: Advection vel. at "now" barocl. step 69 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: un_adv , vn_adv !: Advection vel. at "now" barocl. step 70 71 INTEGER , SAVE :: icycle ! Number of barotropic sub-steps for each internal step nn_baro <= 2.5 nn_baro 72 REAL(wp), SAVE :: rdtbt ! Barotropic time step 73 ! 74 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:) :: wgtbtp1, wgtbtp2 ! 1st & 2nd weights used in time filtering of barotropic fields 75 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: zwz ! ff_f/h at F points 76 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: ftnw, ftne ! triad of coriolis parameter 77 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: ftsw, ftse ! (only used with een vorticity scheme) 78 79 REAL(wp) :: r1_12 = 1._wp / 12._wp ! local ratios 80 REAL(wp) :: r1_8 = 0.125_wp ! 81 REAL(wp) :: r1_4 = 0.25_wp ! 82 REAL(wp) :: r1_2 = 0.5_wp ! 79 83 80 84 !! * Substitutions 81 85 # include "vectopt_loop_substitute.h90" 82 86 !!---------------------------------------------------------------------- 83 !! NEMO/OPA 3.5 , NEMO Consortium (2013)87 !! NEMO/OPA 4.0 , NEMO Consortium (2017) 84 88 !! $Id$ 85 89 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) … … 137 141 INTEGER, INTENT(in) :: kt ! ocean time-step index 138 142 ! 139 LOGICAL :: ll_fw_start ! if true, forward integration140 LOGICAL :: ll_init ! if true, special startup of 2d equations141 LOGICAL :: ll_tmp1, ll_tmp2 ! local logical variables used in W/D142 143 INTEGER :: ji, jj, jk, jn ! dummy loop indices 143 INTEGER :: ikbu, ikbv, noffset ! local integers 144 INTEGER :: iktu, iktv ! local integers 145 REAL(wp) :: zmdi 146 REAL(wp) :: zraur, z1_2dt_b, z2dt_bf ! local scalars 147 REAL(wp) :: zx1, zy1, zx2, zy2 ! - - 148 REAL(wp) :: z1_12, z1_8, z1_4, z1_2 ! - - 149 REAL(wp) :: zu_spg, zv_spg ! - - 150 REAL(wp) :: zhura, zhvra ! - - 151 REAL(wp) :: za0, za1, za2, za3 ! - - 152 ! 153 REAL(wp), POINTER, DIMENSION(:,:) :: zsshp2_e 154 REAL(wp), POINTER, DIMENSION(:,:) :: zu_trd, zv_trd, zu_frc, zv_frc, zssh_frc 155 REAL(wp), POINTER, DIMENSION(:,:) :: zwx, zwy, zhdiv 156 REAL(wp), POINTER, DIMENSION(:,:) :: zhup2_e, zhvp2_e, zhust_e, zhvst_e 157 REAL(wp), POINTER, DIMENSION(:,:) :: zsshu_a, zsshv_a 158 REAL(wp), POINTER, DIMENSION(:,:) :: zhf 159 REAL(wp), POINTER, DIMENSION(:,:) :: zcpx, zcpy ! Wetting/Dying gravity filter coef. 144 LOGICAL :: ll_fw_start ! =T : forward integration 145 LOGICAL :: ll_init ! =T : special startup of 2d equations 146 LOGICAL :: ll_tmp1, ll_tmp2 ! local logical variables used in W/D 147 INTEGER :: ikbu, iktu, noffset ! local integers 148 INTEGER :: ikbv, iktv ! - - 149 REAL(wp) :: z1_2dt_b, z2dt_bf ! local scalars 150 REAL(wp) :: zx1, zx2, zu_spg, zhura ! - - 151 REAL(wp) :: zy1, zy2, zv_spg, zhvra ! - - 152 REAL(wp) :: za0, za1, za2, za3 ! - - 153 REAL(wp) :: zmdi, zztmp ! - - 154 REAL(wp), DIMENSION(jpi,jpj) :: zsshp2_e, zhf 155 REAL(wp), DIMENSION(jpi,jpj) :: zwx, zu_trd, zu_frc, zssh_frc 156 REAL(wp), DIMENSION(jpi,jpj) :: zwy, zv_trd, zv_frc, zhdiv 157 REAL(wp), DIMENSION(jpi,jpj) :: zsshu_a, zhup2_e, zhust_e 158 REAL(wp), DIMENSION(jpi,jpj) :: zsshv_a, zhvp2_e, zhvst_e 159 REAL(wp), DIMENSION(jpi,jpj) :: zCdU_u, zCdU_v ! top/bottom stress at u- & v-points 160 ! 161 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zcpx, zcpy ! Wetting/Dying gravity filter coef. 160 162 !!---------------------------------------------------------------------- 161 163 ! 162 IF( nn_timing == 1 ) CALL timing_start('dyn_spg_ts') 163 ! 164 ! !* Allocate temporary arrays 165 CALL wrk_alloc( jpi,jpj, zsshp2_e, zhdiv ) 166 CALL wrk_alloc( jpi,jpj, zu_trd, zv_trd) 167 CALL wrk_alloc( jpi,jpj, zwx, zwy, zssh_frc, zu_frc, zv_frc) 168 CALL wrk_alloc( jpi,jpj, zhup2_e, zhvp2_e, zhust_e, zhvst_e) 169 CALL wrk_alloc( jpi,jpj, zsshu_a, zsshv_a ) 170 CALL wrk_alloc( jpi,jpj, zhf ) 171 IF( ln_wd ) CALL wrk_alloc( jpi, jpj, zcpx, zcpy ) 164 IF( ln_timing ) CALL timing_start('dyn_spg_ts') 165 ! 166 IF( ln_wd ) ALLOCATE( zcpx(jpi,jpj), zcpy(jpi,jpj) ) 172 167 ! 173 168 zmdi=1.e+20 ! missing data indicator for masking 174 ! !* Local constant initialization 175 z1_12 = 1._wp / 12._wp 176 z1_8 = 0.125_wp 177 z1_4 = 0.25_wp 178 z1_2 = 0.5_wp 179 zraur = 1._wp / rau0 169 ! 180 170 ! ! reciprocal of baroclinic time step 181 171 IF( kt == nit000 .AND. neuler == 0 ) THEN ; z2dt_bf = rdt … … 210 200 CALL ts_wgt( ln_bt_av, ll_fw_start, icycle, wgtbtp1, wgtbtp2 ) 211 201 ! 202 ENDIF 203 ! 204 IF( ln_isfcav ) THEN ! top+bottom friction (ocean cavities) 205 DO jj = 2, jpjm1 206 DO ji = fs_2, fs_jpim1 ! vector opt. 207 zCdU_u(ji,jj) = 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) + rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) 208 zCdU_v(ji,jj) = 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) + rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) 209 END DO 210 END DO 211 ELSE ! bottom friction only 212 DO jj = 2, jpjm1 213 DO ji = fs_2, fs_jpim1 ! vector opt. 214 zCdU_u(ji,jj) = 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) 215 zCdU_v(ji,jj) = 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) 216 END DO 217 END DO 212 218 ENDIF 213 219 ! … … 263 269 !!gm 264 270 !! 265 IF ( .not.ln_sco ) THEN271 IF( .NOT.ln_sco ) THEN 266 272 267 273 !!gm agree the JC comment : this should be done in a much clear way … … 314 320 IF (.NOT.ln_bt_fw .AND.( neuler==0 .AND. kt==nit000+1 ) ) THEN 315 321 ll_fw_start=.FALSE. 316 CALL ts_wgt( ln_bt_av, ll_fw_start, icycle, wgtbtp1, wgtbtp2)322 CALL ts_wgt( ln_bt_av, ll_fw_start, icycle, wgtbtp1, wgtbtp2 ) 317 323 ENDIF 318 324 … … 363 369 zx2 = ( zwx(ji ,jj) + zwx(ji ,jj+1) ) * r1_e2v(ji,jj) 364 370 ! energy conserving formulation for planetary vorticity term 365 zu_trd(ji,jj) = z1_4 * ( zwz(ji ,jj-1) * zy1 + zwz(ji,jj) * zy2 )366 zv_trd(ji,jj) = -z1_4 * ( zwz(ji-1,jj ) * zx1 + zwz(ji,jj) * zx2 )371 zu_trd(ji,jj) = r1_4 * ( zwz(ji ,jj-1) * zy1 + zwz(ji,jj) * zy2 ) 372 zv_trd(ji,jj) = - r1_4 * ( zwz(ji-1,jj ) * zx1 + zwz(ji,jj) * zx2 ) 367 373 END DO 368 374 END DO … … 371 377 DO jj = 2, jpjm1 372 378 DO ji = fs_2, fs_jpim1 ! vector opt. 373 zy1 = z1_8 * ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) &379 zy1 = r1_8 * ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) & 374 380 & + zwy(ji ,jj ) + zwy(ji+1,jj ) ) * r1_e1u(ji,jj) 375 zx1 = - z1_8 * ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) &381 zx1 = - r1_8 * ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) & 376 382 & + zwx(ji ,jj ) + zwx(ji ,jj+1) ) * r1_e2v(ji,jj) 377 383 zu_trd(ji,jj) = zy1 * ( zwz(ji ,jj-1) + zwz(ji,jj) ) … … 383 389 DO jj = 2, jpjm1 384 390 DO ji = fs_2, fs_jpim1 ! vector opt. 385 zu_trd(ji,jj) = + z1_12 * r1_e1u(ji,jj) * ( ftne(ji,jj ) * zwy(ji ,jj ) &391 zu_trd(ji,jj) = + r1_12 * r1_e1u(ji,jj) * ( ftne(ji,jj ) * zwy(ji ,jj ) & 386 392 & + ftnw(ji+1,jj) * zwy(ji+1,jj ) & 387 393 & + ftse(ji,jj ) * zwy(ji ,jj-1) & 388 394 & + ftsw(ji+1,jj) * zwy(ji+1,jj-1) ) 389 zv_trd(ji,jj) = - z1_12 * r1_e2v(ji,jj) * ( ftsw(ji,jj+1) * zwx(ji-1,jj+1) &395 zv_trd(ji,jj) = - r1_12 * r1_e2v(ji,jj) * ( ftsw(ji,jj+1) * zwx(ji-1,jj+1) & 390 396 & + ftse(ji,jj+1) * zwx(ji ,jj+1) & 391 397 & + ftnw(ji,jj ) * zwx(ji-1,jj ) & … … 399 405 ! ! ---------------------------------------------------- 400 406 IF( .NOT.ln_linssh ) THEN ! Variable volume : remove surface pressure gradient 401 IF( ln_wd ) THEN ! Calculating and applying W/D gravity filters402 DO jj = 2, jpjm1403 DO ji = 2, jpim1404 ll_tmp1 = MIN( sshn(ji,jj) , sshn(ji+1,jj) ) > &405 & MAX( -ht_wd(ji,jj) , -ht_wd(ji+1,jj) ) .AND. &406 & MAX( sshn(ji,jj) + ht_wd(ji,jj), sshn(ji+1,jj) + ht_wd(ji+1,jj) ) &407 IF( ln_wd ) THEN ! Calculating and applying W/D gravity filters 408 DO jj = 2, jpjm1 409 DO ji = 2, jpim1 410 ll_tmp1 = MIN( sshn(ji,jj) , sshn(ji+1,jj) ) > & 411 & MAX( -ht_wd(ji,jj) , -ht_wd(ji+1,jj) ) .AND. & 412 & MAX( sshn(ji,jj) + ht_wd(ji,jj), sshn(ji+1,jj) + ht_wd(ji+1,jj) ) & 407 413 & > rn_wdmin1 + rn_wdmin2 408 ll_tmp2 = ( ABS( sshn(ji+1,jj) - sshn(ji ,jj)) > 1.E-12 ).AND.( & 409 & MAX( sshn(ji,jj) , sshn(ji+1,jj) ) > & 410 & MAX( -ht_wd(ji,jj) , -ht_wd(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 ) 411 412 IF(ll_tmp1) THEN 413 zcpx(ji,jj) = 1.0_wp 414 ELSE IF(ll_tmp2) THEN 415 ! no worries about sshn(ji+1,jj) - sshn(ji ,jj) = 0, it won't happen ! here 416 zcpx(ji,jj) = ABS( (sshn(ji+1,jj) + ht_wd(ji+1,jj) - sshn(ji,jj) - ht_wd(ji,jj)) & 417 & / (sshn(ji+1,jj) - sshn(ji ,jj)) ) 418 ELSE 419 zcpx(ji,jj) = 0._wp 420 END IF 421 422 ll_tmp1 = MIN( sshn(ji,jj) , sshn(ji,jj+1) ) > & 414 ll_tmp2 = ( ABS( sshn(ji+1,jj) - sshn(ji ,jj)) > 1.E-12 ).AND.( & 415 & MAX( sshn(ji,jj) , sshn(ji+1,jj) ) > & 416 & MAX( -ht_wd(ji,jj) , -ht_wd(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 ) 417 ! 418 IF(ll_tmp1) THEN 419 zcpx(ji,jj) = 1.0_wp 420 ELSE IF(ll_tmp2) THEN ! no worries about sshn(ji+1,jj) - sshn(ji ,jj) = 0, it won't happen ! here 421 zcpx(ji,jj) = ABS( (sshn(ji+1,jj) + ht_wd(ji+1,jj) - sshn(ji,jj) - ht_wd(ji,jj)) & 422 & / (sshn(ji+1,jj) - sshn(ji ,jj)) ) 423 ELSE 424 zcpx(ji,jj) = 0._wp 425 ENDIF 426 ! 427 ll_tmp1 = MIN( sshn(ji,jj) , sshn(ji,jj+1) ) > & 423 428 & MAX( -ht_wd(ji,jj) , -ht_wd(ji,jj+1) ) .AND. & 424 429 & MAX( sshn(ji,jj) + ht_wd(ji,jj), sshn(ji,jj+1) + ht_wd(ji,jj+1) ) & 425 430 & > rn_wdmin1 + rn_wdmin2 426 ll_tmp2 = ( ABS( sshn(ji,jj) - sshn(ji,jj+1)) > 1.E-12 ).AND.( &431 ll_tmp2 = ( ABS( sshn(ji,jj) - sshn(ji,jj+1)) > 1.E-12 ).AND.( & 427 432 & MAX( sshn(ji,jj) , sshn(ji,jj+1) ) > & 428 433 & MAX( -ht_wd(ji,jj) , -ht_wd(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 ) 429 430 IF(ll_tmp1) THEN431 zcpy(ji,jj) = 1.0_wp432 ELSE IF(ll_tmp2) THEN433 ! no worries about sshn(ji,jj+1) - sshn(ji,jj ) = 0, it won't happen ! here434 zcpy(ji,jj) = ABS( (sshn(ji,jj+1) + ht_wd(ji,jj+1) - sshn(ji,jj) - ht_wd(ji,jj)) &435 &/ (sshn(ji,jj+1) - sshn(ji,jj )) )436 ELSE437 zcpy(ji,jj) = 0._wp438 ENDIF439 END DO440 END DO441 442 DO jj = 2, jpjm1443 DO ji = 2, jpim1444 zu_trd(ji,jj) = zu_trd(ji,jj) - grav * ( sshn(ji+1,jj ) - sshn(ji ,jj ) ) &445 &* r1_e1u(ji,jj) * zcpx(ji,jj)446 zv_trd(ji,jj) = zv_trd(ji,jj) - grav * ( sshn(ji ,jj+1) - sshn(ji ,jj ) ) &447 &* r1_e2v(ji,jj) * zcpy(ji,jj)448 END DO449 END DO450 434 ! 435 IF(ll_tmp1) THEN 436 zcpy(ji,jj) = 1.0_wp 437 ELSE IF(ll_tmp2) THEN 438 ! no worries about sshn(ji,jj+1) - sshn(ji,jj ) = 0, it won't happen ! here 439 zcpy(ji,jj) = ABS( (sshn(ji,jj+1) + ht_wd(ji,jj+1) - sshn(ji,jj) - ht_wd(ji,jj)) & 440 & / (sshn(ji,jj+1) - sshn(ji,jj )) ) 441 ELSE 442 zcpy(ji,jj) = 0._wp 443 ENDIF 444 END DO 445 END DO 446 ! 447 DO jj = 2, jpjm1 448 DO ji = 2, jpim1 449 zu_trd(ji,jj) = zu_trd(ji,jj) - grav * ( sshn(ji+1,jj ) - sshn(ji ,jj ) ) & 450 & * r1_e1u(ji,jj) * zcpx(ji,jj) 451 zv_trd(ji,jj) = zv_trd(ji,jj) - grav * ( sshn(ji ,jj+1) - sshn(ji ,jj ) ) & 452 & * r1_e2v(ji,jj) * zcpy(ji,jj) 453 END DO 454 END DO 455 ! 451 456 ELSE 452 453 DO jj = 2, jpjm1454 DO ji = fs_2, fs_jpim1 ! vector opt.455 zu_trd(ji,jj) = zu_trd(ji,jj) - grav * ( sshn(ji+1,jj ) - sshn(ji ,jj ) ) * r1_e1u(ji,jj)456 zv_trd(ji,jj) = zv_trd(ji,jj) - grav * ( sshn(ji ,jj+1) - sshn(ji ,jj ) ) * r1_e2v(ji,jj)457 END DO458 END DO459 ENDIF460 461 ENDIF 462 457 ! 458 DO jj = 2, jpjm1 459 DO ji = fs_2, fs_jpim1 ! vector opt. 460 zu_trd(ji,jj) = zu_trd(ji,jj) - grav * ( sshn(ji+1,jj ) - sshn(ji ,jj ) ) * r1_e1u(ji,jj) 461 zv_trd(ji,jj) = zv_trd(ji,jj) - grav * ( sshn(ji ,jj+1) - sshn(ji ,jj ) ) * r1_e2v(ji,jj) 462 END DO 463 END DO 464 ENDIF 465 ! 466 ENDIF 467 ! 463 468 DO jj = 2, jpjm1 ! Remove coriolis term (and possibly spg) from barotropic trend 464 469 DO ji = fs_2, fs_jpim1 … … 468 473 END DO 469 474 ! 470 ! ! Add bottomstress contribution from baroclinic velocities:471 IF (ln_bt_fw) THEN475 ! ! Add BOTTOM stress contribution from baroclinic velocities: 476 IF( ln_bt_fw ) THEN 472 477 DO jj = 2, jpjm1 473 478 DO ji = fs_2, fs_jpim1 ! vector opt. … … 491 496 ! Note that the "unclipped" bottom friction parameter is used even with explicit drag 492 497 IF( ln_wd ) THEN 493 zu_frc(:,:) = zu_frc(:,:) + MAX(r1_hu_n(:,:) * bfrua(:,:),-1._wp / rdtbt) * zwx(:,:) 494 zv_frc(:,:) = zv_frc(:,:) + MAX(r1_hv_n(:,:) * bfrva(:,:),-1._wp / rdtbt) * zwy(:,:) 495 ELSE 496 zu_frc(:,:) = zu_frc(:,:) + r1_hu_n(:,:) * bfrua(:,:) * zwx(:,:) 497 zv_frc(:,:) = zv_frc(:,:) + r1_hv_n(:,:) * bfrva(:,:) * zwy(:,:) 498 END IF 499 ! 500 ! ! Add top stress contribution from baroclinic velocities: 501 IF( ln_bt_fw ) THEN 502 DO jj = 2, jpjm1 503 DO ji = fs_2, fs_jpim1 ! vector opt. 504 iktu = miku(ji,jj) 505 iktv = mikv(ji,jj) 506 zwx(ji,jj) = un(ji,jj,iktu) - un_b(ji,jj) ! NOW top baroclinic velocities 507 zwy(ji,jj) = vn(ji,jj,iktv) - vn_b(ji,jj) 498 zztmp = - 1._wp / rdtbt 499 DO jj = 2, jpjm1 500 DO ji = fs_2, fs_jpim1 ! vector opt. 501 zu_frc(ji,jj) = zu_frc(ji,jj) + MAX( r1_hu_n(ji,jj) * 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) , zztmp ) * zwx(ji,jj) 502 zv_frc(ji,jj) = zv_frc(ji,jj) + MAX( r1_hv_n(ji,jj) * 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) , zztmp ) * zwy(ji,jj) 508 503 END DO 509 504 END DO … … 511 506 DO jj = 2, jpjm1 512 507 DO ji = fs_2, fs_jpim1 ! vector opt. 513 iktu = miku(ji,jj) 514 iktv = mikv(ji,jj) 515 zwx(ji,jj) = ub(ji,jj,iktu) - ub_b(ji,jj) ! BEFORE top baroclinic velocities 516 zwy(ji,jj) = vb(ji,jj,iktv) - vb_b(ji,jj) 517 END DO 518 END DO 519 ENDIF 520 ! 521 ! Note that the "unclipped" top friction parameter is used even with explicit drag 522 zu_frc(:,:) = zu_frc(:,:) + r1_hu_n(:,:) * tfrua(:,:) * zwx(:,:) 523 zv_frc(:,:) = zv_frc(:,:) + r1_hv_n(:,:) * tfrva(:,:) * zwy(:,:) 508 zu_frc(ji,jj) = zu_frc(ji,jj) + r1_hu_n(ji,jj) * 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) * zwx(ji,jj) 509 zv_frc(ji,jj) = zv_frc(ji,jj) + r1_hv_n(ji,jj) * 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) * zwy(ji,jj) 510 END DO 511 END DO 512 ENDIF 513 ! 514 IF( ln_isfcav ) THEN ! Add TOP stress contribution from baroclinic velocities: 515 IF( ln_bt_fw ) THEN 516 DO jj = 2, jpjm1 517 DO ji = fs_2, fs_jpim1 ! vector opt. 518 iktu = miku(ji,jj) 519 iktv = mikv(ji,jj) 520 zwx(ji,jj) = un(ji,jj,iktu) - un_b(ji,jj) ! NOW top baroclinic velocities 521 zwy(ji,jj) = vn(ji,jj,iktv) - vn_b(ji,jj) 522 END DO 523 END DO 524 ELSE 525 DO jj = 2, jpjm1 526 DO ji = fs_2, fs_jpim1 ! vector opt. 527 iktu = miku(ji,jj) 528 iktv = mikv(ji,jj) 529 zwx(ji,jj) = ub(ji,jj,iktu) - ub_b(ji,jj) ! BEFORE top baroclinic velocities 530 zwy(ji,jj) = vb(ji,jj,iktv) - vb_b(ji,jj) 531 END DO 532 END DO 533 ENDIF 534 ! 535 ! Note that the "unclipped" top friction parameter is used even with explicit drag 536 DO jj = 2, jpjm1 537 DO ji = fs_2, fs_jpim1 ! vector opt. 538 zu_frc(ji,jj) = zu_frc(ji,jj) + r1_hu_n(ji,jj) * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * zwx(ji,jj) 539 zv_frc(ji,jj) = zv_frc(ji,jj) + r1_hv_n(ji,jj) * 0.5*( rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) * zwy(ji,jj) 540 END DO 541 END DO 542 ENDIF 524 543 ! 525 IF (ln_bt_fw) THEN ! Add wind forcing 526 zu_frc(:,:) = zu_frc(:,:) + zraur * utau(:,:) * r1_hu_n(:,:) 527 zv_frc(:,:) = zv_frc(:,:) + zraur * vtau(:,:) * r1_hv_n(:,:) 544 IF( ln_bt_fw ) THEN ! Add wind forcing 545 DO jj = 2, jpjm1 546 DO ji = fs_2, fs_jpim1 ! vector opt. 547 zu_frc(ji,jj) = zu_frc(ji,jj) + r1_rau0 * utau(ji,jj) * r1_hu_n(ji,jj) 548 zv_frc(ji,jj) = zv_frc(ji,jj) + r1_rau0 * vtau(ji,jj) * r1_hv_n(ji,jj) 549 END DO 550 END DO 528 551 ELSE 529 zu_frc(:,:) = zu_frc(:,:) + zraur * z1_2 * ( utau_b(:,:) + utau(:,:) ) * r1_hu_n(:,:) 530 zv_frc(:,:) = zv_frc(:,:) + zraur * z1_2 * ( vtau_b(:,:) + vtau(:,:) ) * r1_hv_n(:,:) 552 zztmp = r1_rau0 * r1_2 553 DO jj = 2, jpjm1 554 DO ji = fs_2, fs_jpim1 ! vector opt. 555 zu_frc(ji,jj) = zu_frc(ji,jj) + zztmp * ( utau_b(ji,jj) + utau(ji,jj) ) * r1_hu_n(ji,jj) 556 zv_frc(ji,jj) = zv_frc(ji,jj) + zztmp * ( vtau_b(ji,jj) + vtau(ji,jj) ) * r1_hv_n(ji,jj) 557 END DO 558 END DO 531 559 ENDIF 532 560 ! 533 IF ( ln_apr_dyn ) THEN! Add atm pressure forcing534 IF (ln_bt_fw) THEN561 IF( ln_apr_dyn ) THEN ! Add atm pressure forcing 562 IF( ln_bt_fw ) THEN 535 563 DO jj = 2, jpjm1 536 564 DO ji = fs_2, fs_jpim1 ! vector opt. … … 542 570 END DO 543 571 ELSE 572 zztmp = grav * r1_2 544 573 DO jj = 2, jpjm1 545 574 DO ji = fs_2, fs_jpim1 ! vector opt. 546 zu_spg = grav * z1_2* ( ssh_ib (ji+1,jj ) - ssh_ib (ji,jj) &547 & 548 zv_spg = grav * z1_2* ( ssh_ib (ji ,jj+1) - ssh_ib (ji,jj) &549 & 575 zu_spg = zztmp * ( ssh_ib (ji+1,jj ) - ssh_ib (ji,jj) & 576 & + ssh_ibb(ji+1,jj ) - ssh_ibb(ji,jj) ) * r1_e1u(ji,jj) 577 zv_spg = zztmp * ( ssh_ib (ji ,jj+1) - ssh_ib (ji,jj) & 578 & + ssh_ibb(ji ,jj+1) - ssh_ibb(ji,jj) ) * r1_e2v(ji,jj) 550 579 zu_frc(ji,jj) = zu_frc(ji,jj) + zu_spg 551 580 zv_frc(ji,jj) = zv_frc(ji,jj) + zv_spg … … 558 587 ! ! Surface net water flux and rivers 559 588 IF (ln_bt_fw) THEN 560 zssh_frc(:,:) = zraur* ( emp(:,:) - rnf(:,:) + fwfisf(:,:) )589 zssh_frc(:,:) = r1_rau0 * ( emp(:,:) - rnf(:,:) + fwfisf(:,:) ) 561 590 ELSE 562 zssh_frc(:,:) = zraur * z1_2 * ( emp(:,:) + emp_b(:,:) - rnf(:,:) - rnf_b(:,:) & 563 & + fwfisf(:,:) + fwfisf_b(:,:) ) 591 zztmp = r1_rau0 * r1_2 592 zssh_frc(:,:) = zztmp * ( emp(:,:) + emp_b(:,:) - rnf(:,:) - rnf_b(:,:) & 593 & + fwfisf(:,:) + fwfisf_b(:,:) ) 564 594 ENDIF 565 595 ! … … 657 687 DO jj = 2, jpjm1 ! Sea Surface Height at u- & v-points 658 688 DO ji = 2, fs_jpim1 ! Vector opt. 659 zwx(ji,jj) = z1_2 * ssumask(ji,jj) * r1_e1e2u(ji,jj) &689 zwx(ji,jj) = r1_2 * ssumask(ji,jj) * r1_e1e2u(ji,jj) & 660 690 & * ( e1e2t(ji ,jj) * zsshp2_e(ji ,jj) & 661 691 & + e1e2t(ji+1,jj) * zsshp2_e(ji+1,jj) ) 662 zwy(ji,jj) = z1_2 * ssvmask(ji,jj) * r1_e1e2v(ji,jj) &692 zwy(ji,jj) = r1_2 * ssvmask(ji,jj) * r1_e1e2v(ji,jj) & 663 693 & * ( e1e2t(ji,jj ) * zsshp2_e(ji,jj ) & 664 694 & + e1e2t(ji,jj+1) * zsshp2_e(ji,jj+1) ) … … 685 715 IF( .NOT.Agrif_Root() .AND. ln_bt_fw ) THEN 686 716 IF((nbondi == -1).OR.(nbondi == 2)) THEN 687 DO jj =1,jpj688 zwx(2 ,jj) = ubdy_w(jj) * e2u(2,jj)717 DO jj = 1, jpj 718 zwx(2:nbghostcells+1,jj) = ubdy_w(jj) * e2u(2:nbghostcells+1,jj) 689 719 END DO 690 720 ENDIF 691 721 IF((nbondi == 1).OR.(nbondi == 2)) THEN 692 722 DO jj=1,jpj 693 zwx(nlci- 2,jj) = ubdy_e(jj) * e2u(nlci-2,jj)723 zwx(nlci-nbghostcells-1:nlci-2,jj) = ubdy_e(jj) * e2u(nlci-nbghostcells-1:nlci-2,jj) 694 724 END DO 695 725 ENDIF 696 726 IF((nbondj == -1).OR.(nbondj == 2)) THEN 697 727 DO ji=1,jpi 698 zwy(ji,2 ) = vbdy_s(ji) * e1v(ji,2)728 zwy(ji,2:nbghostcells+1) = vbdy_s(ji) * e1v(ji,2:nbghostcells+1) 699 729 END DO 700 730 ENDIF 701 731 IF((nbondj == 1).OR.(nbondj == 2)) THEN 702 732 DO ji=1,jpi 703 zwy(ji,nlcj- 2) = vbdy_n(ji) * e1v(ji,nlcj-2)733 zwy(ji,nlcj-nbghostcells-1:nlcj-2) = vbdy_n(ji) * e1v(ji,nlcj-nbghostcells-1:nlcj-2) 704 734 END DO 705 735 ENDIF … … 734 764 DO jj = 2, jpjm1 735 765 DO ji = 2, jpim1 ! NO Vector Opt. 736 zsshu_a(ji,jj) = z1_2 * ssumask(ji,jj) * r1_e1e2u(ji,jj) &766 zsshu_a(ji,jj) = r1_2 * ssumask(ji,jj) * r1_e1e2u(ji,jj) & 737 767 & * ( e1e2t(ji ,jj ) * ssha_e(ji ,jj ) & 738 768 & + e1e2t(ji+1,jj ) * ssha_e(ji+1,jj ) ) 739 zsshv_a(ji,jj) = z1_2 * ssvmask(ji,jj) * r1_e1e2v(ji,jj) &769 zsshv_a(ji,jj) = r1_2 * ssvmask(ji,jj) * r1_e1e2v(ji,jj) & 740 770 & * ( e1e2t(ji ,jj ) * ssha_e(ji ,jj ) & 741 771 & + e1e2t(ji ,jj+1) * ssha_e(ji ,jj+1) ) … … 813 843 DO jj = 2, jpjm1 814 844 DO ji = 2, jpim1 815 zx1 = z1_2 * ssumask(ji ,jj) * r1_e1e2u(ji ,jj) &845 zx1 = r1_2 * ssumask(ji ,jj) * r1_e1e2u(ji ,jj) & 816 846 & * ( e1e2t(ji ,jj ) * zsshp2_e(ji ,jj) & 817 847 & + e1e2t(ji+1,jj ) * zsshp2_e(ji+1,jj ) ) 818 zy1 = z1_2 * ssvmask(ji ,jj) * r1_e1e2v(ji ,jj ) &848 zy1 = r1_2 * ssvmask(ji ,jj) * r1_e1e2v(ji ,jj ) & 819 849 & * ( e1e2t(ji ,jj ) * zsshp2_e(ji ,jj ) & 820 850 & + e1e2t(ji ,jj+1) * zsshp2_e(ji ,jj+1) ) … … 840 870 zx1 = ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) ) * r1_e2v(ji,jj) 841 871 zx2 = ( zwx(ji ,jj ) + zwx(ji ,jj+1) ) * r1_e2v(ji,jj) 842 zu_trd(ji,jj) = z1_4 * ( zwz(ji ,jj-1) * zy1 + zwz(ji,jj) * zy2 )843 zv_trd(ji,jj) =- z1_4 * ( zwz(ji-1,jj ) * zx1 + zwz(ji,jj) * zx2 )872 zu_trd(ji,jj) = r1_4 * ( zwz(ji ,jj-1) * zy1 + zwz(ji,jj) * zy2 ) 873 zv_trd(ji,jj) =-r1_4 * ( zwz(ji-1,jj ) * zx1 + zwz(ji,jj) * zx2 ) 844 874 END DO 845 875 END DO … … 848 878 DO jj = 2, jpjm1 849 879 DO ji = fs_2, fs_jpim1 ! vector opt. 850 zy1 = z1_8 * ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) &880 zy1 = r1_8 * ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) & 851 881 & + zwy(ji ,jj ) + zwy(ji+1,jj ) ) * r1_e1u(ji,jj) 852 zx1 = - z1_8 * ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) &882 zx1 = - r1_8 * ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) & 853 883 & + zwx(ji ,jj ) + zwx(ji ,jj+1) ) * r1_e2v(ji,jj) 854 884 zu_trd(ji,jj) = zy1 * ( zwz(ji ,jj-1) + zwz(ji,jj) ) … … 860 890 DO jj = 2, jpjm1 861 891 DO ji = fs_2, fs_jpim1 ! vector opt. 862 zu_trd(ji,jj) = + z1_12 * r1_e1u(ji,jj) * ( ftne(ji,jj ) * zwy(ji ,jj ) &892 zu_trd(ji,jj) = + r1_12 * r1_e1u(ji,jj) * ( ftne(ji,jj ) * zwy(ji ,jj ) & 863 893 & + ftnw(ji+1,jj) * zwy(ji+1,jj ) & 864 894 & + ftse(ji,jj ) * zwy(ji ,jj-1) & 865 895 & + ftsw(ji+1,jj) * zwy(ji+1,jj-1) ) 866 zv_trd(ji,jj) = - z1_12 * r1_e2v(ji,jj) * ( ftsw(ji,jj+1) * zwx(ji-1,jj+1) &896 zv_trd(ji,jj) = - r1_12 * r1_e2v(ji,jj) * ( ftsw(ji,jj+1) * zwx(ji-1,jj+1) & 867 897 & + ftse(ji,jj+1) * zwx(ji ,jj+1) & 868 898 & + ftnw(ji,jj ) * zwx(ji-1,jj ) & … … 885 915 ENDIF 886 916 ! 887 ! Add bottom stresses: 888 zu_trd(:,:) = zu_trd(:,:) + bfrua(:,:) * un_e(:,:) * hur_e(:,:) 889 zv_trd(:,:) = zv_trd(:,:) + bfrva(:,:) * vn_e(:,:) * hvr_e(:,:) 890 ! 891 ! Add top stresses: 892 zu_trd(:,:) = zu_trd(:,:) + tfrua(:,:) * un_e(:,:) * hur_e(:,:) 893 zv_trd(:,:) = zv_trd(:,:) + tfrva(:,:) * vn_e(:,:) * hvr_e(:,:) 917 DO jj = 2, jpjm1 ! Add top/bottom stresses: 918 DO ji = fs_2, fs_jpim1 ! vector opt. 919 zu_trd(ji,jj) = zu_trd(ji,jj) + zCdU_u(ji,jj) * un_e(ji,jj) * hur_e(ji,jj) 920 zv_trd(ji,jj) = zv_trd(ji,jj) + zCdU_v(ji,jj) * vn_e(ji,jj) * hvr_e(ji,jj) 921 END DO 922 END DO 894 923 ! 895 924 ! Surface pressure trend: 896 897 925 IF( ln_wd ) THEN 898 926 DO jj = 2, jpjm1 … … 1025 1053 vn_adv(:,:) = zwy(:,:) * r1_hv_n(:,:) 1026 1054 ELSE 1027 un_adv(:,:) = z1_2 * ( ub2_b(:,:) + zwx(:,:) ) * r1_hu_n(:,:)1028 vn_adv(:,:) = z1_2 * ( vb2_b(:,:) + zwy(:,:) ) * r1_hv_n(:,:)1055 un_adv(:,:) = r1_2 * ( ub2_b(:,:) + zwx(:,:) ) * r1_hu_n(:,:) 1056 vn_adv(:,:) = r1_2 * ( vb2_b(:,:) + zwy(:,:) ) * r1_hv_n(:,:) 1029 1057 END IF 1030 1058 … … 1044 1072 DO jj = 1, jpjm1 1045 1073 DO ji = 1, jpim1 ! NO Vector Opt. 1046 zsshu_a(ji,jj) = z1_2 * umask(ji,jj,1) * r1_e1e2u(ji,jj) &1074 zsshu_a(ji,jj) = r1_2 * umask(ji,jj,1) * r1_e1e2u(ji,jj) & 1047 1075 & * ( e1e2t(ji ,jj) * ssha(ji ,jj) & 1048 1076 & + e1e2t(ji+1,jj) * ssha(ji+1,jj) ) 1049 zsshv_a(ji,jj) = z1_2 * vmask(ji,jj,1) * r1_e1e2v(ji,jj) &1077 zsshv_a(ji,jj) = r1_2 * vmask(ji,jj,1) * r1_e1e2v(ji,jj) & 1050 1078 & * ( e1e2t(ji,jj ) * ssha(ji,jj ) & 1051 1079 & + e1e2t(ji,jj+1) * ssha(ji,jj+1) ) … … 1091 1119 IF( lrst_oce .AND.ln_bt_fw ) CALL ts_rst( kt, 'WRITE' ) 1092 1120 ! 1093 CALL wrk_dealloc( jpi,jpj, zsshp2_e, zhdiv ) 1094 CALL wrk_dealloc( jpi,jpj, zu_trd, zv_trd ) 1095 CALL wrk_dealloc( jpi,jpj, zwx, zwy, zssh_frc, zu_frc, zv_frc ) 1096 CALL wrk_dealloc( jpi,jpj, zhup2_e, zhvp2_e, zhust_e, zhvst_e ) 1097 CALL wrk_dealloc( jpi,jpj, zsshu_a, zsshv_a ) 1098 CALL wrk_dealloc( jpi,jpj, zhf ) 1099 IF( ln_wd ) CALL wrk_dealloc( jpi, jpj, zcpx, zcpy ) 1100 ! 1101 IF ( ln_diatmb ) THEN 1121 IF( ln_wd ) DEALLOCATE( zcpx, zcpy ) 1122 ! 1123 IF( ln_diatmb ) THEN 1102 1124 CALL iom_put( "baro_u" , un_b*umask(:,:,1)+zmdi*(1-umask(:,:,1 ) ) ) ! Barotropic U Velocity 1103 1125 CALL iom_put( "baro_v" , vn_b*vmask(:,:,1)+zmdi*(1-vmask(:,:,1 ) ) ) ! Barotropic V Velocity 1104 1126 ENDIF 1105 IF( nn_timing == 1 )CALL timing_stop('dyn_spg_ts')1127 IF( ln_timing ) CALL timing_stop('dyn_spg_ts') 1106 1128 ! 1107 1129 END SUBROUTINE dyn_spg_ts … … 1248 1270 INTEGER :: ji ,jj ! dummy loop indices 1249 1271 REAL(wp) :: zxr2, zyr2, zcmax ! local scalar 1250 REAL(wp), POINTER, DIMENSION(:,:) :: zcu1272 REAL(wp), DIMENSION(jpi,jpj) :: zcu 1251 1273 !!---------------------------------------------------------------------- 1252 1274 ! 1253 1275 ! Max courant number for ext. grav. waves 1254 !1255 CALL wrk_alloc( jpi,jpj, zcu )1256 1276 ! 1257 1277 DO jj = 1, jpj … … 1320 1340 ENDIF 1321 1341 ! 1322 CALL wrk_dealloc( jpi,jpj, zcu )1323 !1324 1342 END SUBROUTINE dyn_spg_ts_init 1325 1343
Note: See TracChangeset
for help on using the changeset viewer.
