Changeset 6004 for branches/2015/dev_r5836_NOC3_vvl_by_default/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90
- Timestamp:
- 2015-12-04T17:05:58+01:00 (8 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2015/dev_r5836_NOC3_vvl_by_default/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90
r5904 r6004 1 1 MODULE dynspg_ts 2 !!====================================================================== 3 !! *** MODULE dynspg_ts *** 4 !! Ocean dynamics: surface pressure gradient trend, split-explicit scheme 2 5 !!====================================================================== 3 6 !! History : 1.0 ! 2004-12 (L. Bessieres, G. Madec) Original code … … 11 14 !! 3.5 ! 2013-07 (J. Chanut) Switch to Forward-backward time stepping 12 15 !! 3.6 ! 2013-11 (A. Coward) Update for z-tilde compatibility 16 !! 3.7 ! 2015-11 (J. Chanut) free surface simplification 13 17 !!--------------------------------------------------------------------- 14 #if defined key_dynspg_ts 18 15 19 !!---------------------------------------------------------------------- 16 !! 'key_dynspg_ts' split explicit free surface17 !! ----------------------------------------------------------------------18 !! dyn_spg_ts : compute surface pressure gradient trend using a time-19 !! splitting scheme and add to the general trend20 !! dyn_spg_ts : compute surface pressure gradient trend using a time-splitting scheme 21 !! dyn_spg_ts_init: initialisation of the time-splitting scheme 22 !! ts_wgt : set time-splitting weights for temporal averaging (or not) 23 !! ts_rst : read/write time-splitting fields in restart file 20 24 !!---------------------------------------------------------------------- 21 25 USE oce ! ocean dynamics and tracers 22 26 USE dom_oce ! ocean space and time domain 23 27 USE sbc_oce ! surface boundary condition: ocean 28 USE zdf_oce ! Bottom friction coefts 24 29 USE sbcisf ! ice shelf variable (fwfisf) 25 USE dynspg_oce ! surface pressure gradient variables 30 USE sbcapr ! surface boundary condition: atmospheric pressure 31 USE dynadv , ONLY: ln_dynadv_vec 26 32 USE phycst ! physical constants 27 33 USE dynvor ! vorticity term 28 34 USE bdy_par ! for lk_bdy 29 USE bdytides ! open boundary condition data 35 USE bdytides ! open boundary condition data 30 36 USE bdydyn2d ! open boundary conditions on barotropic variables 31 37 USE sbctide ! tides 32 38 USE updtide ! tide potential 39 ! 40 USE in_out_manager ! I/O manager 33 41 USE lib_mpp ! distributed memory computing library 34 42 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 35 43 USE prtctl ! Print control 36 USE in_out_manager ! I/O manager37 44 USE iom ! IOM library 38 45 USE restart ! only for lrst_oce 39 USE zdf_oce ! Bottom friction coefts40 46 USE wrk_nemo ! Memory Allocation 41 47 USE timing ! Timing 42 USE sbcapr ! surface boundary condition: atmospheric pressure43 USE dynadv, ONLY: ln_dynadv_vec44 48 #if defined key_agrif 45 49 USE agrif_opa_interp ! agrif … … 60 64 REAL(wp),SAVE :: rdtbt ! Barotropic time step 61 65 62 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:) :: & 63 wgtbtp1, & ! Primary weights used for time filtering of barotropic variables 64 wgtbtp2 ! Secondary weights used for time filtering of barotropic variables 65 66 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: zwz ! ff/h at F points 67 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: ftnw, ftne ! triad of coriolis parameter 68 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: ftsw, ftse ! (only used with een vorticity scheme) 69 70 ! Arrays below are saved to allow testing of the "no time averaging" option 71 ! If this option is not retained, these could be replaced by temporary arrays 72 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: sshbb_e, sshb_e, & ! Instantaneous barotropic arrays 73 ubb_e, ub_e, & 74 vbb_e, vb_e 66 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:) :: wgtbtp1, wgtbtp2 !: 1st & 2nd weights used in time filtering of barotropic fields 67 68 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: zwz !: ff/h at F points 69 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: ftnw, ftne !: triad of coriolis parameter 70 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: ftsw, ftse !: (only used with een vorticity scheme) 71 72 !! Time filtered arrays at baroclinic time step: 73 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: un_adv , vn_adv !: Advection vel. at "now" barocl. step 74 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ub2_b , vb2_b !: Half step fluxes (ln_bt_fw=T) 75 #if defined key_agrif 76 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ub2_i_b, vb2_i_b !: Half step time integrated fluxes 77 #endif 78 79 !! Arrays at barotropic time step: ! bef before ! before ! now ! after ! 80 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ubb_e , ub_e , un_e , ua_e !: u-external velocity 81 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: vbb_e , vb_e , vn_e , va_e !: v-external velocity 82 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: sshbb_e , sshb_e , sshn_e , ssha_e !: external ssh 83 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hu_e !: external u-depth 84 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hv_e !: external v-depth 85 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hur_e !: inverse of u-depth 86 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hvr_e !: inverse of v-depth 75 87 76 88 !! * Substitutions … … 87 99 !! *** routine dyn_spg_ts_alloc *** 88 100 !!---------------------------------------------------------------------- 89 INTEGER :: ierr( 3)101 INTEGER :: ierr(5) 90 102 !!---------------------------------------------------------------------- 91 103 ierr(:) = 0 92 104 ! 93 ALLOCATE( sshb_e(jpi,jpj), sshbb_e(jpi,jpj), & 94 & ub_e(jpi,jpj) , vb_e(jpi,jpj) , & 95 & ubb_e(jpi,jpj) , vbb_e(jpi,jpj) , STAT= ierr(1) ) 96 ! 97 ALLOCATE( wgtbtp1(3*nn_baro), wgtbtp2(3*nn_baro), zwz(jpi,jpj), STAT= ierr(2) ) 98 ! 99 IF( ln_dynvor_een ) ALLOCATE( ftnw(jpi,jpj) , ftne(jpi,jpj) , & 100 & ftsw(jpi,jpj) , ftse(jpi,jpj) , STAT=ierr(3) ) 105 ALLOCATE( ssha_e(jpi,jpj), sshn_e(jpi,jpj), sshb_e(jpi,jpj), sshbb_e(jpi,jpj), & 106 & ua_e(jpi,jpj), un_e(jpi,jpj), ub_e(jpi,jpj), ubb_e(jpi,jpj), & 107 & va_e(jpi,jpj), vn_e(jpi,jpj), vb_e(jpi,jpj), vbb_e(jpi,jpj), & 108 & hu_e(jpi,jpj), hur_e(jpi,jpj), hv_e(jpi,jpj), hvr_e(jpi,jpj), STAT=ierr(1) ) 109 ! 110 ALLOCATE( wgtbtp1(3*nn_baro), wgtbtp2(3*nn_baro), zwz(jpi,jpj) , STAT=ierr(2) ) 111 ! 112 IF( ln_dynvor_een ) ALLOCATE( ftnw(jpi,jpj) , ftne(jpi,jpj) , & 113 & ftsw(jpi,jpj) , ftse(jpi,jpj) , STAT=ierr(3) ) 114 ! 115 ALLOCATE( ub2_b(jpi,jpj), vb2_b(jpi,jpj), un_adv(jpi,jpj), vn_adv(jpi,jpj) , STAT=ierr(4) ) 116 #if defined key_agrif 117 ALLOCATE( ub2_i_b(jpi,jpj), vb2_i_b(jpi,jpj) , STAT= ierr(5)) 118 #endif 101 119 ! 102 120 dyn_spg_ts_alloc = MAXVAL( ierr(:) ) 121 ! 103 122 IF( lk_mpp ) CALL mpp_sum( dyn_spg_ts_alloc ) 104 IF( dyn_spg_ts_alloc /= 0 ) CALL ctl_warn('dyn spg_oce_alloc: failed to allocate arrays')123 IF( dyn_spg_ts_alloc /= 0 ) CALL ctl_warn('dyn_spg_ts_alloc: failed to allocate arrays') 105 124 ! 106 125 END FUNCTION dyn_spg_ts_alloc … … 110 129 !!---------------------------------------------------------------------- 111 130 !! 112 !! ** Purpose : 113 !! -Compute the now trend due to the explicit time stepping114 !! of the quasi-linear barotropic system.131 !! ** Purpose : - Compute the now trend due to the explicit time stepping 132 !! of the quasi-linear barotropic system, and add it to the 133 !! general momentum trend. 115 134 !! 116 !! ** Method : 135 !! ** Method : - split-explicit schem (time splitting) : 117 136 !! Barotropic variables are advanced from internal time steps 118 137 !! "n" to "n+1" if ln_bt_fw=T … … 128 147 !! continuity equation taken at the baroclinic time steps. This 129 148 !! ensures tracers conservation. 130 !! - Update 3d trend (ua, va)with barotropic component.149 !! - (ua, va) momentum trend updated with barotropic component. 131 150 !! 132 !! References : Shchepetkin, A.F. and J.C. McWilliams, 2005: 133 !! The regional oceanic modeling system (ROMS): 134 !! a split-explicit, free-surface, 135 !! topography-following-coordinate oceanic model. 136 !! Ocean Modelling, 9, 347-404. 151 !! References : Shchepetkin and McWilliams, Ocean Modelling, 2005. 137 152 !!--------------------------------------------------------------------- 138 153 INTEGER, INTENT(in) :: kt ! ocean time-step index … … 149 164 REAL(wp) :: za0, za1, za2, za3 ! - - 150 165 ! 151 REAL(wp), POINTER, DIMENSION(:,:) :: zun_e, zvn_e,zsshp2_e152 REAL(wp), POINTER, DIMENSION(:,:) :: 153 REAL(wp), POINTER, DIMENSION(:,:) :: zu_sum, zv_sum,zwx, zwy, zhdiv154 REAL(wp), POINTER, DIMENSION(:,:) :: 155 REAL(wp), POINTER, DIMENSION(:,:) :: 156 REAL(wp), POINTER, DIMENSION(:,:) :: 166 REAL(wp), POINTER, DIMENSION(:,:) :: zsshp2_e 167 REAL(wp), POINTER, DIMENSION(:,:) :: zu_trd, zv_trd, zu_frc, zv_frc, zssh_frc 168 REAL(wp), POINTER, DIMENSION(:,:) :: zwx, zwy, zhdiv 169 REAL(wp), POINTER, DIMENSION(:,:) :: zhup2_e, zhvp2_e, zhust_e, zhvst_e 170 REAL(wp), POINTER, DIMENSION(:,:) :: zsshu_a, zsshv_a 171 REAL(wp), POINTER, DIMENSION(:,:) :: zhf 157 172 !!---------------------------------------------------------------------- 158 173 ! 159 174 IF( nn_timing == 1 ) CALL timing_start('dyn_spg_ts') 160 175 ! 176 ! !* Allocate temporary arrays 161 177 CALL wrk_alloc( jpi,jpj, zsshp2_e, zhdiv ) 162 CALL wrk_alloc( jpi,jpj, zu_trd, zv_trd , zun_e, zvn_e)163 CALL wrk_alloc( jpi,jpj, zwx, zwy, z u_sum, zv_sum, zssh_frc, zu_frc, zv_frc)178 CALL wrk_alloc( jpi,jpj, zu_trd, zv_trd) 179 CALL wrk_alloc( jpi,jpj, zwx, zwy, zssh_frc, zu_frc, zv_frc) 164 180 CALL wrk_alloc( jpi,jpj, zhup2_e, zhvp2_e, zhust_e, zhvst_e) 165 CALL wrk_alloc( jpi,jpj, zsshu_a, zsshv_a )181 CALL wrk_alloc( jpi,jpj, zsshu_a, zsshv_a ) 166 182 CALL wrk_alloc( jpi,jpj, zhf ) 167 183 ! … … 180 196 ll_fw_start = .FALSE. 181 197 ! ! time offset in steps for bdy data update 182 IF( .NOT.ln_bt_fw ) THEN ; noffset = -2*nn_baro183 ELSE ; noffset = 0198 IF( .NOT.ln_bt_fw ) THEN ; noffset = - nn_baro 199 ELSE ; noffset = 0 184 200 ENDIF 185 201 ! … … 194 210 ! 195 211 IF( ln_bt_fw .OR. neuler == 0 ) THEN 196 ll_fw_start =.TRUE.197 noffset = 0212 ll_fw_start =.TRUE. 213 noffset = 0 198 214 ELSE 199 ll_fw_start =.FALSE.215 ll_fw_start =.FALSE. 200 216 ENDIF 201 217 ! … … 212 228 ! and update depths at T-F points (ht and zhf resp.) at each barotropic time step 213 229 ! 214 IF 215 IF ( ln_dynvor_een ) THEN!== EEN scheme ==!230 IF( kt == nit000 .OR. .NOT.ln_linssh ) THEN 231 IF( ln_dynvor_een ) THEN !== EEN scheme ==! 216 232 SELECT CASE( nn_een_e3f ) !* ff/e3 at F-point 217 233 CASE ( 0 ) ! original formulation (masked averaging of e3t divided by 4) … … 219 235 DO ji = 1, jpim1 220 236 zwz(ji,jj) = ( ht_n(ji ,jj+1) + ht_n(ji+1,jj+1) + & 221 & ht_n(ji ,jj ) + ht_n(ji+1,jj ) ) / 4._wp237 & ht_n(ji ,jj ) + ht_n(ji+1,jj ) ) * 0.25_wp 222 238 IF( zwz(ji,jj) /= 0._wp ) zwz(ji,jj) = ff(ji,jj) / zwz(ji,jj) 223 239 END DO … … 407 423 zv_frc(:,:) = zv_frc(:,:) + r1_hv_n(:,:) * bfrva(:,:) * zwy(:,:) 408 424 ! 409 IF (ln_bt_fw) THEN ! Add wind forcing425 IF( ln_bt_fw ) THEN ! Add wind forcing 410 426 zu_frc(:,:) = zu_frc(:,:) + zraur * utau(:,:) * r1_hu_n(:,:) 411 427 zv_frc(:,:) = zv_frc(:,:) + zraur * vtau(:,:) * r1_hv_n(:,:) … … 477 493 ! 478 494 IF (ln_bt_fw) THEN ! FORWARD integration: start from NOW fields 479 sshn_e(:,:) = sshn(:,:)480 zun_e (:,:) =un_b(:,:)481 zvn_e (:,:) =vn_b(:,:)495 sshn_e(:,:) = sshn(:,:) 496 un_e (:,:) = un_b(:,:) 497 vn_e (:,:) = vn_b(:,:) 482 498 ! 483 499 hu_e (:,:) = hu_n(:,:) … … 486 502 hvr_e (:,:) = r1_hv_n(:,:) 487 503 ELSE ! CENTRED integration: start from BEFORE fields 488 sshn_e(:,:) = sshb(:,:)489 zun_e (:,:) =ub_b(:,:)490 zvn_e (:,:) =vb_b(:,:)504 sshn_e(:,:) = sshb(:,:) 505 un_e (:,:) = ub_b(:,:) 506 vn_e (:,:) = vb_b(:,:) 491 507 ! 492 508 hu_e (:,:) = hu_b(:,:) … … 502 518 va_b (:,:) = 0._wp 503 519 ssha (:,:) = 0._wp ! Sum for after averaged sea level 504 zu_sum(:,:) = 0._wp ! Sum for now transport issued from ts loop505 zv_sum(:,:) = 0._wp520 un_adv(:,:) = 0._wp ! Sum for now transport issued from ts loop 521 vn_adv(:,:) = 0._wp 506 522 ! ! ==================== ! 507 523 DO jn = 1, icycle ! sub-time-step loop ! … … 511 527 ! Update only tidal forcing at open boundaries 512 528 #if defined key_tide 513 IF ( lk_bdy .AND. lk_tide ) CALL bdy_dta_tides( kt, kit=jn, time_offset=(noffset+1))514 IF ( ln_tide_pot .AND. lk_tide ) CALL upd_tide ( kt, kit=jn, koffset=noffset)529 IF( lk_bdy .AND. lk_tide ) CALL bdy_dta_tides( kt, kit=jn, time_offset= noffset+1 ) 530 IF( ln_tide_pot .AND. lk_tide ) CALL upd_tide ( kt, kit=jn, time_offset= noffset ) 515 531 #endif 516 532 ! … … 527 543 528 544 ! Extrapolate barotropic velocities at step jit+0.5: 529 ua_e(:,:) = za1 * zun_e(:,:) + za2 * ub_e(:,:) + za3 * ubb_e(:,:)530 va_e(:,:) = za1 * zvn_e(:,:) + za2 * vb_e(:,:) + za3 * vbb_e(:,:)545 ua_e(:,:) = za1 * un_e(:,:) + za2 * ub_e(:,:) + za3 * ubb_e(:,:) 546 va_e(:,:) = za1 * vn_e(:,:) + za2 * vb_e(:,:) + za3 * vbb_e(:,:) 531 547 532 548 IF( .NOT.ln_linssh ) THEN !* Update ocean depth (variable volume case only) … … 589 605 ! Sum over sub-time-steps to compute advective velocities 590 606 za2 = wgtbtp2(jn) 591 zu_sum(:,:) = zu_sum(:,:) + za2 * zwx(:,:) * r1_e2u(:,:)592 zv_sum(:,:) = zv_sum(:,:) + za2 * zwy(:,:) * r1_e1v(:,:)607 un_adv(:,:) = un_adv(:,:) + za2 * zwx(:,:) * r1_e2u(:,:) 608 vn_adv(:,:) = vn_adv(:,:) + za2 * zwy(:,:) * r1_e1v(:,:) 593 609 ! 594 610 ! Set next sea level: … … 648 664 ! 649 665 ! Compute associated depths at U and V points: 650 IF( .NOT. ( ln_dynadv_vec .OR. ln_linssh )) THEN !* Vector form666 IF( .NOT.ln_linssh .AND. .NOT.ln_dynadv_vec ) THEN !* Vector form 651 667 ! 652 668 DO jj = 2, jpjm1 … … 671 687 ! zwy(:,:) = e1v(:,:) * va_e(:,:) * zhvp2_e(:,:) 672 688 ! 673 IF( ln_dynvor_ene .OR. ln_dynvor_mix ) THEN 689 IF( ln_dynvor_ene .OR. ln_dynvor_mix ) THEN !== energy conserving or mixed scheme ==! 674 690 DO jj = 2, jpjm1 675 691 DO ji = fs_2, fs_jpim1 ! vector opt. … … 683 699 END DO 684 700 ! 685 ELSEIF ( ln_dynvor_ens ) THEN 701 ELSEIF ( ln_dynvor_ens ) THEN !== enstrophy conserving scheme ==! 686 702 DO jj = 2, jpjm1 687 703 DO ji = fs_2, fs_jpim1 ! vector opt. … … 695 711 END DO 696 712 ! 697 ELSEIF ( ln_dynvor_een ) THEN !== energy and enstrophy conserving scheme ==!713 ELSEIF ( ln_dynvor_een ) THEN !== energy and enstrophy conserving scheme ==! 698 714 DO jj = 2, jpjm1 699 715 DO ji = fs_2, fs_jpim1 ! vector opt. … … 724 740 ! 725 741 ! Add bottom stresses: 726 zu_trd(:,:) = zu_trd(:,:) + bfrua(:,:) * zun_e(:,:) * hur_e(:,:)727 zv_trd(:,:) = zv_trd(:,:) + bfrva(:,:) * zvn_e(:,:) * hvr_e(:,:)742 zu_trd(:,:) = zu_trd(:,:) + bfrua(:,:) * un_e(:,:) * hur_e(:,:) 743 zv_trd(:,:) = zv_trd(:,:) + bfrva(:,:) * vn_e(:,:) * hvr_e(:,:) 728 744 ! 729 745 ! Surface pressure trend: … … 742 758 DO jj = 2, jpjm1 743 759 DO ji = fs_2, fs_jpim1 ! vector opt. 744 ua_e(ji,jj) = ( zun_e(ji,jj) & 745 & + rdtbt * ( zwx(ji,jj) & 746 & + zu_trd(ji,jj) & 747 & + zu_frc(ji,jj) ) & 748 & ) * umask(ji,jj,1) 749 750 va_e(ji,jj) = ( zvn_e(ji,jj) & 751 & + rdtbt * ( zwy(ji,jj) & 752 & + zv_trd(ji,jj) & 753 & + zv_frc(ji,jj) ) & 754 & ) * vmask(ji,jj,1) 755 END DO 756 END DO 757 760 ua_e(ji,jj) = ( un_e(ji,jj) & 761 & + rdtbt * ( zwx(ji,jj) & 762 & + zu_trd(ji,jj) & 763 & + zu_frc(ji,jj) ) ) * umask(ji,jj,1) 764 ! 765 va_e(ji,jj) = ( vn_e(ji,jj) & 766 & + rdtbt * ( zwy(ji,jj) & 767 & + zv_trd(ji,jj) & 768 & + zv_frc(ji,jj) ) ) * vmask(ji,jj,1) 769 END DO 770 END DO 771 ! 758 772 ELSE !* Flux form 759 773 DO jj = 2, jpjm1 760 774 DO ji = fs_2, fs_jpim1 ! vector opt. 761 762 zhura = umask(ji,jj,1)/(hu_0(ji,jj) + zsshu_a(ji,jj) + 1._wp - umask(ji,jj,1)) 763 zhvra = vmask(ji,jj,1)/(hv_0(ji,jj) + zsshv_a(ji,jj) + 1._wp - vmask(ji,jj,1)) 764 765 ua_e(ji,jj) = ( hu_e(ji,jj) * zun_e(ji,jj) & 766 & + rdtbt * ( zhust_e(ji,jj) * zwx(ji,jj) & 767 & + zhup2_e(ji,jj) * zu_trd(ji,jj) & 768 & + hu_n(ji,jj) * zu_frc(ji,jj) ) & 769 & ) * zhura 770 771 va_e(ji,jj) = ( hv_e(ji,jj) * zvn_e(ji,jj) & 772 & + rdtbt * ( zhvst_e(ji,jj) * zwy(ji,jj) & 773 & + zhvp2_e(ji,jj) * zv_trd(ji,jj) & 774 & + hv_n(ji,jj) * zv_frc(ji,jj) ) & 775 & ) * zhvra 775 zhura = umask(ji,jj,1) / ( hu_0(ji,jj) + zsshu_a(ji,jj) + 1._wp - umask(ji,jj,1) ) 776 zhvra = vmask(ji,jj,1) / ( hv_0(ji,jj) + zsshv_a(ji,jj) + 1._wp - vmask(ji,jj,1) ) 777 ! 778 ua_e(ji,jj) = ( hu_e(ji,jj) * un_e(ji,jj) & 779 & + rdtbt * ( zhust_e(ji,jj) * zwx(ji,jj) & 780 & + zhup2_e(ji,jj) * zu_trd(ji,jj) & 781 & + hu_n(ji,jj) * zu_frc(ji,jj) ) ) * zhura 782 ! 783 va_e(ji,jj) = ( hv_e(ji,jj) * vn_e(ji,jj) & 784 & + rdtbt * ( zhvst_e(ji,jj) * zwy(ji,jj) & 785 & + zhvp2_e(ji,jj) * zv_trd(ji,jj) & 786 & + hv_n(ji,jj) * zv_frc(ji,jj) ) ) * zhvra 776 787 END DO 777 788 END DO … … 779 790 ! 780 791 IF( .NOT.ln_linssh ) THEN !* Update ocean depth (variable volume case only) 781 ! ! ----------------------------------------------782 792 hu_e (:,:) = hu_0(:,:) + zsshu_a(:,:) 783 793 hv_e (:,:) = hv_0(:,:) + zsshv_a(:,:) … … 787 797 ENDIF 788 798 ! !* domain lateral boundary 789 ! ! -----------------------790 !791 799 CALL lbc_lnk_multi( ua_e, 'U', -1._wp, va_e , 'V', -1._wp ) 792 800 ! 793 801 #if defined key_bdy 794 802 ! ! open boundaries 795 IF( lk_bdy ) CALL bdy_dyn2d( jn, ua_e, va_e, zun_e, zvn_e, hur_e, hvr_e, ssha_e )803 IF( lk_bdy ) CALL bdy_dyn2d( jn, ua_e, va_e, un_e, vn_e, hur_e, hvr_e, ssha_e ) 796 804 #endif 797 805 #if defined key_agrif … … 801 809 ! ! ---- 802 810 ubb_e (:,:) = ub_e (:,:) 803 ub_e (:,:) = zun_e(:,:)804 zun_e(:,:) = ua_e (:,:)811 ub_e (:,:) = un_e (:,:) 812 un_e (:,:) = ua_e (:,:) 805 813 ! 806 814 vbb_e (:,:) = vb_e (:,:) 807 vb_e (:,:) = zvn_e(:,:)808 zvn_e(:,:) = va_e (:,:)815 vb_e (:,:) = vn_e (:,:) 816 vn_e (:,:) = va_e (:,:) 809 817 ! 810 818 sshbb_e(:,:) = sshb_e(:,:) … … 831 839 ! ----------------------------------------------------------------------------- 832 840 ! 833 ! At this stage ssha holds a time averaged value 834 ! ! Sea Surface Height at u-,v- and f-points 835 IF( .NOT.ln_linssh ) THEN ! (required only in non-linear free surface case) 841 ! Set advection velocity correction: 842 zwx(:,:) = un_adv(:,:) 843 zwy(:,:) = vn_adv(:,:) 844 IF( ( kt == nit000 .AND. neuler==0 ) .OR. .NOT.ln_bt_fw ) THEN 845 un_adv(:,:) = zwx(:,:) * r1_hu_n(:,:) 846 vn_adv(:,:) = zwy(:,:) * r1_hv_n(:,:) 847 ELSE 848 un_adv(:,:) = z1_2 * ( ub2_b(:,:) + zwx(:,:) ) * r1_hu_n(:,:) 849 vn_adv(:,:) = z1_2 * ( vb2_b(:,:) + zwy(:,:) ) * r1_hv_n(:,:) 850 END IF 851 852 IF( ln_bt_fw ) THEN ! Save integrated transport for next computation 853 ub2_b(:,:) = zwx(:,:) 854 vb2_b(:,:) = zwy(:,:) 855 ENDIF 856 ! 857 ! Update barotropic trend: 858 IF( ln_dynadv_vec .OR. ln_linssh ) THEN 859 DO jk=1,jpkm1 860 ua(:,:,jk) = ua(:,:,jk) + ( ua_b(:,:) - ub_b(:,:) ) * z1_2dt_b 861 va(:,:,jk) = va(:,:,jk) + ( va_b(:,:) - vb_b(:,:) ) * z1_2dt_b 862 END DO 863 ELSE 864 ! At this stage, ssha has been corrected: compute new depths at velocity points 836 865 DO jj = 1, jpjm1 837 866 DO ji = 1, jpim1 ! NO Vector Opt. … … 845 874 END DO 846 875 CALL lbc_lnk_multi( zsshu_a, 'U', 1._wp, zsshv_a, 'V', 1._wp ) ! Boundary conditions 847 ENDIF 848 ! 849 ! Set advection velocity correction: 850 IF( ( kt == nit000 .AND. neuler==0 ) .OR. .NOT.ln_bt_fw ) THEN 851 un_adv(:,:) = zu_sum(:,:) * r1_hu_n(:,:) 852 vn_adv(:,:) = zv_sum(:,:) * r1_hv_n(:,:) 853 ELSE 854 un_adv(:,:) = z1_2 * ( ub2_b(:,:) + zu_sum(:,:) ) * r1_hu_n(:,:) 855 vn_adv(:,:) = z1_2 * ( vb2_b(:,:) + zv_sum(:,:) ) * r1_hv_n(:,:) 856 END IF 857 858 IF( ln_bt_fw ) THEN ! Save integrated transport for next computation 859 ub2_b(:,:) = zu_sum(:,:) 860 vb2_b(:,:) = zv_sum(:,:) 861 ENDIF 862 ! 863 ! Update barotropic trend: 864 IF( ln_dynadv_vec .OR. ln_linssh ) THEN 865 DO jk=1,jpkm1 866 ua(:,:,jk) = ua(:,:,jk) + ( ua_b(:,:) - ub_b(:,:) ) * z1_2dt_b 867 va(:,:,jk) = va(:,:,jk) + ( va_b(:,:) - vb_b(:,:) ) * z1_2dt_b 868 END DO 869 ELSE 876 ! 870 877 DO jk=1,jpkm1 871 878 ua(:,:,jk) = ua(:,:,jk) + r1_hu_n(:,:) * ( ua_b(:,:) - ub_b(:,:) * hu_b(:,:) ) * z1_2dt_b … … 883 890 ! 884 891 END DO 892 ! 893 CALL iom_put( "ubar", un_adv(:,:) ) ! barotropic i-current 894 CALL iom_put( "vbar", vn_adv(:,:) ) ! barotropic i-current 885 895 ! 886 896 #if defined key_agrif … … 898 908 vb2_i_b(:,:) = vb2_i_b(:,:) + za1 * vb2_b(:,:) 899 909 ENDIF 900 !901 !902 910 #endif 903 !904 911 ! !* write time-spliting arrays in the restart 905 912 IF( lrst_oce .AND.ln_bt_fw ) CALL ts_rst( kt, 'WRITE' ) 906 913 ! 907 914 CALL wrk_dealloc( jpi,jpj, zsshp2_e, zhdiv ) 908 CALL wrk_dealloc( jpi,jpj, zu_trd, zv_trd , zun_e, zvn_e)909 CALL wrk_dealloc( jpi,jpj, zwx, zwy, z u_sum, zv_sum, zssh_frc, zu_frc, zv_frc )915 CALL wrk_dealloc( jpi,jpj, zu_trd, zv_trd ) 916 CALL wrk_dealloc( jpi,jpj, zwx, zwy, zssh_frc, zu_frc, zv_frc ) 910 917 CALL wrk_dealloc( jpi,jpj, zhup2_e, zhvp2_e, zhust_e, zhvst_e ) 911 CALL wrk_dealloc( jpi,jpj, zsshu_a, zsshv_a )918 CALL wrk_dealloc( jpi,jpj, zsshu_a, zsshv_a ) 912 919 CALL wrk_dealloc( jpi,jpj, zhf ) 913 920 ! … … 994 1001 END SUBROUTINE ts_wgt 995 1002 1003 996 1004 SUBROUTINE ts_rst( kt, cdrw ) 997 1005 !!--------------------------------------------------------------------- … … 1047 1055 END SUBROUTINE ts_rst 1048 1056 1049 SUBROUTINE dyn_spg_ts_init( kt ) 1057 1058 SUBROUTINE dyn_spg_ts_init 1050 1059 !!--------------------------------------------------------------------- 1051 1060 !! *** ROUTINE dyn_spg_ts_init *** … … 1053 1062 !! ** Purpose : Set time splitting options 1054 1063 !!---------------------------------------------------------------------- 1055 INTEGER , INTENT(in) :: kt ! ocean time-step 1056 ! 1057 INTEGER :: ji ,jj 1058 INTEGER :: ios ! Local integer output status for namelist read 1059 REAL(wp) :: zxr2, zyr2, zcmax 1060 REAL(wp), POINTER, DIMENSION(:,:) :: zcu 1061 !! 1062 NAMELIST/namsplit/ ln_bt_fw, ln_bt_av, ln_bt_nn_auto, & 1063 & nn_baro, rn_bt_cmax, nn_bt_flt 1064 INTEGER :: ji ,jj ! dummy loop indices 1065 REAL(wp) :: zxr2, zyr2, zcmax ! local scalar 1066 REAL(wp), POINTER, DIMENSION(:,:) :: zcu 1064 1067 !!---------------------------------------------------------------------- 1065 1068 ! 1066 REWIND( numnam_ref ) ! Namelist namsplit in reference namelist : time splitting parameters 1067 READ ( numnam_ref, namsplit, IOSTAT = ios, ERR = 901) 1068 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsplit in reference namelist', lwp ) 1069 1070 REWIND( numnam_cfg ) ! Namelist namsplit in configuration namelist : time splitting parameters 1071 READ ( numnam_cfg, namsplit, IOSTAT = ios, ERR = 902 ) 1072 902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsplit in configuration namelist', lwp ) 1073 IF(lwm) WRITE ( numond, namsplit ) 1074 ! 1075 ! ! Max courant number for ext. grav. waves 1076 ! 1077 CALL wrk_alloc( jpi, jpj, zcu ) 1078 ! 1079 IF( .NOT.ln_linssh ) THEN 1080 DO jj = 1, jpj 1081 DO ji =1, jpi 1082 zxr2 = r1_e1t(ji,jj) * r1_e1t(ji,jj) 1083 zyr2 = r1_e2t(ji,jj) * r1_e2t(ji,jj) 1084 zcu(ji,jj) = SQRT( grav * ht_0(ji,jj) * (zxr2 + zyr2) ) 1085 END DO 1086 END DO 1087 ELSE 1088 !!gm BUG ?? restartability issue if ssh changes are large.... 1089 !!gm We should just test this with ht_0 only, no? 1090 DO jj = 1, jpj 1091 DO ji =1, jpi 1092 zxr2 = r1_e1t(ji,jj) * r1_e1t(ji,jj) 1093 zyr2 = r1_e1t(ji,jj) * r1_e1t(ji,jj) 1094 zcu(ji,jj) = SQRT( grav * ht_n(ji,jj) * (zxr2 + zyr2) ) 1095 END DO 1096 END DO 1097 ENDIF 1098 1069 ! Max courant number for ext. grav. waves 1070 ! 1071 CALL wrk_alloc( jpi,jpj, zcu ) 1072 ! 1073 DO jj = 1, jpj 1074 DO ji =1, jpi 1075 zxr2 = r1_e1t(ji,jj) * r1_e1t(ji,jj) 1076 zyr2 = r1_e2t(ji,jj) * r1_e2t(ji,jj) 1077 zcu(ji,jj) = SQRT( grav * ht_0(ji,jj) * (zxr2 + zyr2) ) 1078 END DO 1079 END DO 1080 ! 1099 1081 zcmax = MAXVAL( zcu(:,:) ) 1100 1082 IF( lk_mpp ) CALL mpp_max( zcmax ) 1101 1083 1102 1084 ! Estimate number of iterations to satisfy a max courant number= rn_bt_cmax 1103 IF (ln_bt_nn_auto)nn_baro = CEILING( rdt / rn_bt_cmax * zcmax)1085 IF( ln_bt_auto ) nn_baro = CEILING( rdt / rn_bt_cmax * zcmax) 1104 1086 1105 1087 rdtbt = rdt / REAL( nn_baro , wp ) … … 1109 1091 IF(lwp) WRITE(numout,*) 'dyn_spg_ts : split-explicit free surface' 1110 1092 IF(lwp) WRITE(numout,*) '~~~~~~~~~~' 1111 IF( ln_bt_ nn_auto ) THEN1112 IF(lwp) WRITE(numout,*) ' ln_ts_ nn_auto=.true. Automatically set nn_baro '1093 IF( ln_bt_auto ) THEN 1094 IF(lwp) WRITE(numout,*) ' ln_ts_auto=.true. Automatically set nn_baro ' 1113 1095 IF(lwp) WRITE(numout,*) ' Max. courant number allowed: ', rn_bt_cmax 1114 1096 ELSE 1115 IF(lwp) WRITE(numout,*) ' ln_ts_ nn_auto=.false.: Use nn_baro in namelist '1097 IF(lwp) WRITE(numout,*) ' ln_ts_auto=.false.: Use nn_baro in namelist ' 1116 1098 ENDIF 1117 1099 … … 1131 1113 #if defined key_agrif 1132 1114 ! Restrict the use of Agrif to the forward case only 1133 IF ((.NOT.ln_bt_fw ).AND.(.NOT.Agrif_Root()))CALL ctl_stop( 'AGRIF not implemented if ln_bt_fw=.FALSE.' )1115 IF( .NOT.ln_bt_fw .AND. .NOT.Agrif_Root() ) CALL ctl_stop( 'AGRIF not implemented if ln_bt_fw=.FALSE.' ) 1134 1116 #endif 1135 1117 ! 1136 1118 IF(lwp) WRITE(numout,*) ' Time filter choice, nn_bt_flt: ', nn_bt_flt 1137 1119 SELECT CASE ( nn_bt_flt ) 1138 CASE( 0 ); IF(lwp) WRITE(numout,*) ' Dirac'1139 CASE( 1 ); IF(lwp) WRITE(numout,*) ' Boxcar: width = nn_baro'1140 CASE( 2 ); IF(lwp) WRITE(numout,*) ' Boxcar: width = 2*nn_baro'1141 1120 CASE( 0 ) ; IF(lwp) WRITE(numout,*) ' Dirac' 1121 CASE( 1 ) ; IF(lwp) WRITE(numout,*) ' Boxcar: width = nn_baro' 1122 CASE( 2 ) ; IF(lwp) WRITE(numout,*) ' Boxcar: width = 2*nn_baro' 1123 CASE DEFAULT ; CALL ctl_stop( 'unrecognised value for nn_bt_flt: should 0,1,2' ) 1142 1124 END SELECT 1143 1125 ! … … 1147 1129 IF(lwp) WRITE(numout,*) ' Maximum Courant number is :', zcmax 1148 1130 ! 1149 IF ((.NOT.ln_bt_av).AND.(.NOT.ln_bt_fw)) THEN1131 IF( .NOT.ln_bt_av .AND. .NOT.ln_bt_fw ) THEN 1150 1132 CALL ctl_stop( 'dynspg_ts ERROR: No time averaging => only forward integration is possible' ) 1151 1133 ENDIF 1152 IF 1134 IF( zcmax>0.9_wp ) THEN 1153 1135 CALL ctl_stop( 'dynspg_ts ERROR: Maximum Courant number is greater than 0.9: Inc. nn_baro !' ) 1154 1136 ENDIF 1155 1137 ! 1156 CALL wrk_dealloc( jpi, jpj,zcu )1138 CALL wrk_dealloc( jpi,jpj, zcu ) 1157 1139 ! 1158 1140 END SUBROUTINE dyn_spg_ts_init 1159 1141 1160 #else1161 !!---------------------------------------------------------------------------1162 !! Default case : Empty module No split explicit free surface1163 !!---------------------------------------------------------------------------1164 CONTAINS1165 INTEGER FUNCTION dyn_spg_ts_alloc() ! Dummy function1166 dyn_spg_ts_alloc = 01167 END FUNCTION dyn_spg_ts_alloc1168 SUBROUTINE dyn_spg_ts( kt ) ! Empty routine1169 INTEGER, INTENT(in) :: kt1170 WRITE(*,*) 'dyn_spg_ts: You should not have seen this print! error?', kt1171 END SUBROUTINE dyn_spg_ts1172 SUBROUTINE ts_rst( kt, cdrw ) ! Empty routine1173 INTEGER , INTENT(in) :: kt ! ocean time-step1174 CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag1175 WRITE(*,*) 'ts_rst : You should not have seen this print! error?', kt, cdrw1176 END SUBROUTINE ts_rst1177 SUBROUTINE dyn_spg_ts_init( kt ) ! Empty routine1178 INTEGER , INTENT(in) :: kt ! ocean time-step1179 WRITE(*,*) 'dyn_spg_ts_init : You should not have seen this print! error?', kt1180 END SUBROUTINE dyn_spg_ts_init1181 #endif1182 1183 1142 !!====================================================================== 1184 1143 END MODULE dynspg_ts
Note: See TracChangeset
for help on using the changeset viewer.
