Changeset 1502
- Timestamp:
- 2009-07-20T17:20:23+02:00 (15 years ago)
- Location:
- trunk/NEMO/OPA_SRC
- Files:
-
- 4 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMO/OPA_SRC/BDY/bdydyn.F90
r1146 r1502 7 7 !! - ! 2007-07 (D. Storkey) Move Flather implementation to separate routine. 8 8 !! 3.0 ! 2008-04 (NEMO team) add in the reference version 9 !! 3.2 ! 2008-04 (R. Benshila) consider velocity instead of transport 9 10 !!---------------------------------------------------------------------- 10 11 #if defined key_bdy … … 86 87 END DO 87 88 ! 88 CALL lbc_lnk( ua, 'U', 1. ) ! Boundary points should be updated89 CALL lbc_lnk( va, 'V', 1. ) !89 CALL lbc_lnk( ua, 'U', -1. ) ! Boundary points should be updated 90 CALL lbc_lnk( va, 'V', -1. ) ! 90 91 ! 91 92 ENDIF ! ln_bdy_dyn_frs … … 96 97 #if defined key_dynspg_exp || defined key_dynspg_ts 97 98 !! Option to use Flather with dynspg_flt not coded yet... 98 SUBROUTINE bdy_dyn_fla 99 SUBROUTINE bdy_dyn_fla( pssh ) 99 100 !!---------------------------------------------------------------------- 100 101 !! *** SUBROUTINE bdy_dyn_fla *** … … 116 117 !! continental shelf. Mem. Soc. R. Sci. Liege, Ser. 6,10, 141-164. 117 118 !!---------------------------------------------------------------------- 119 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pssh 120 118 121 INTEGER :: ib, igrd ! dummy loop indices 119 122 INTEGER :: ii, ij, iim1, iip1, ijm1, ijp1 ! 2D addresses 120 123 REAL(wp) :: zcorr ! Flather correction 124 REAL(wp) :: zforc ! temporary scalar 121 125 !!---------------------------------------------------------------------- 122 126 … … 127 131 IF(ln_bdy_dyn_fla .OR. ln_bdy_tides) THEN ! If these are both false, then this routine does nothing. 128 132 129 ! Fill temporary array with ssh data (here spgu): 130 igrd = 1 131 spgu(:,:) = 0.0 132 DO ib = 1, nblenrim(igrd) 133 ii = nbi(ib,igrd) 134 ij = nbj(ib,igrd) 135 IF( ln_bdy_dyn_fla ) spgu(ii, ij) = sshbdy(ib) 136 IF( ln_bdy_tides ) spgu(ii, ij) = spgu(ii, ij) + sshtide(ib) 137 END DO 133 ! Fill temporary array with ssh data (here spgu): 134 igrd = 1 135 spgu(:,:) = 0.0 136 DO ib = 1, nblenrim(igrd) 137 ii = nbi(ib,igrd) 138 ij = nbj(ib,igrd) 139 IF( ln_bdy_dyn_fla ) spgu(ii, ij) = sshbdy(ib) 140 IF( ln_bdy_tides ) spgu(ii, ij) = spgu(ii, ij) + sshtide(ib) 141 END DO 142 ! 143 igrd = 2 ! Flather bc on u-velocity; 144 ! ! remember that flagu=-1 if normal velocity direction is outward 145 ! ! I think we should rather use after ssh ? 146 DO ib = 1, nblenrim(igrd) 147 ii = nbi(ib,igrd) 148 ij = nbj(ib,igrd) 149 iim1 = ii + MAX( 0, INT( flagu(ib) ) ) ! T pts i-indice inside the boundary 150 iip1 = ii - MIN( 0, INT( flagu(ib) ) ) ! T pts i-indice outside the boundary 151 ! 152 zcorr = - flagu(ib) * SQRT( grav * hur_e(ii, ij) ) * ( pssh(iim1, ij) - spgu(iip1,ij) ) 153 zforc = ubtbdy(ib) + utide(ib) 154 ua_e(ii,ij) = zforc + zcorr * umask(ii,ij,1) 155 END DO 156 ! 157 igrd = 3 ! Flather bc on v-velocity 158 ! ! remember that flagv=-1 if normal velocity direction is outward 159 DO ib = 1, nblenrim(igrd) 160 ii = nbi(ib,igrd) 161 ij = nbj(ib,igrd) 162 ijm1 = ij + MAX( 0, INT( flagv(ib) ) ) ! T pts j-indice inside the boundary 163 ijp1 = ij - MIN( 0, INT( flagv(ib) ) ) ! T pts j-indice outside the boundary 164 ! 165 zcorr = - flagv(ib) * SQRT( grav * hvr_e(ii, ij) ) * ( pssh(ii, ijm1) - spgu(ii,ijp1) ) 166 zforc = vbtbdy(ib) + vtide(ib) 167 va_e(ii,ij) = zforc + zcorr * vmask(ii,ij,1) 168 END DO 169 CALL lbc_lnk( ua_e, 'U', 1. ) ! Boundary points should be updated 170 CALL lbc_lnk( va_e, 'V', 1. ) ! 171 ! 172 ENDIF ! ln_bdy_dyn_fla .or. ln_bdy_tides 138 173 ! 139 igrd = 2 ! Flather bc on u-velocity;140 ! ! remember that flagu=-1 if normal velocity direction is outward141 ! ! I think we should rather use after ssh ?142 DO ib = 1, nblenrim(igrd)143 ii = nbi(ib,igrd)144 ij = nbj(ib,igrd)145 iim1 = ii + MAX( 0, INT( flagu(ib) ) ) ! T pts i-indice inside the boundary146 iip1 = ii - MIN( 0, INT( flagu(ib) ) ) ! T pts i-indice outside the boundary147 !148 zcorr = - flagu(ib) * SQRT( grav / (hu_e(ii, ij) + 1.e-20) ) * ( sshn_e(iim1, ij) - spgu(iip1,ij) )149 ua_e(ii, ij) = ( ubtbdy(ib) + utide(ib) ) * hu_e(ii,ij)150 ua_e(ii,ij) = ua_e(ii,ij) + zcorr * umask(ii,ij,1) * hu_e(ii,ij)151 END DO152 !153 igrd = 3 ! Flather bc on v-velocity154 ! ! remember that flagv=-1 if normal velocity direction is outward155 DO ib = 1, nblenrim(igrd)156 ii = nbi(ib,igrd)157 ij = nbj(ib,igrd)158 ijm1 = ij + MAX( 0, INT( flagv(ib) ) ) ! T pts j-indice inside the boundary159 ijp1 = ij - MIN( 0, INT( flagv(ib) ) ) ! T pts j-indice outside the boundary160 !161 zcorr = - flagv(ib) * SQRT( grav / (hv_e(ii, ij) + 1.e-20) ) * ( sshn_e(ii, ijm1) - spgu(ii,ijp1) )162 va_e(ii, ij) = ( vbtbdy(ib) + vtide(ib) ) * hv_e(ii,ij)163 va_e(ii,ij) = va_e(ii,ij) + zcorr * vmask(ii,ij,1) * hv_e(ii,ij)164 END DO165 !166 CALL lbc_lnk( ua_e, 'U', 1. ) ! Boundary points should be updated167 CALL lbc_lnk( va_e, 'V', 1. ) !168 !169 ENDIF ! ln_bdy_dyn_fla .or. ln_bdy_tides170 171 174 END SUBROUTINE bdy_dyn_fla 172 175 #endif -
trunk/NEMO/OPA_SRC/DYN/dynnxt.F90
r1438 r1502 1 1 MODULE dynnxt 2 !!====================================================================== 2 !!========================================================================= 3 3 !! *** MODULE dynnxt *** 4 4 !! Ocean dynamics: time stepping 5 !!====================================================================== 6 !!====================================================================== 5 !!========================================================================= 7 6 !! History : OPA ! 1987-02 (P. Andrich, D. L Hostis) Original code 8 7 !! ! 1990-10 (C. Levy, G. Madec) … … 15 14 !! 2.0 ! 2005-11 (V. Garnier) Surface pressure gradient organization 16 15 !! 2.3 ! 2007-07 (D. Storkey) Calls to BDY routines. 17 !! 3.2 ! 2009-0 4 (G. Madec, R.Benshila)) re-introduce the vvl option18 !!---------------------------------------------------------------------- 16 !! 3.2 ! 2009-06 (G. Madec, R.Benshila) re-introduce the vvl option 17 !!------------------------------------------------------------------------- 19 18 20 !!---------------------------------------------------------------------- 21 !! dyn_nxt : update the horizontal velocity from the momentum trend22 !!---------------------------------------------------------------------- 19 !!------------------------------------------------------------------------- 20 !! dyn_nxt : obtain the next (after) horizontal velocity 21 !!------------------------------------------------------------------------- 23 22 USE oce ! ocean dynamics and tracers 24 23 USE dom_oce ! ocean space and time domain 25 USE in_out_manager ! I/O manager 24 USE dynspg_oce ! type of surface pressure gradient 25 USE dynadv ! dynamics: vector invariant versus flux form 26 USE domvvl ! variable volume 26 27 USE obc_oce ! ocean open boundary conditions 27 28 USE obcdyn ! open boundary condition for momentum (obc_dyn routine) … … 31 32 USE bdydta ! unstructured open boundary conditions 32 33 USE bdydyn ! unstructured open boundary conditions 33 USE dynspg_oce ! type of surface pressure gradient 34 USE agrif_opa_update 35 USE agrif_opa_interp 36 USE in_out_manager ! I/O manager 34 37 USE lbclnk ! lateral boundary condition (or mpp link) 35 38 USE prtctl ! Print control 36 USE agrif_opa_update37 USE agrif_opa_interp38 USE domvvl ! variable volume39 39 40 40 IMPLICIT NONE … … 57 57 !! *** ROUTINE dyn_nxt *** 58 58 !! 59 !! ** Purpose : Compute the after horizontal velocity from the 60 !! momentum trend. 61 !! 62 !! ** Method : Apply lateral boundary conditions on the trends (ua,va) 63 !! through calls to routine lbc_lnk. 64 !! After velocity is compute using a leap-frog scheme environment: 65 !! (ua,va) = (ub,vb) + 2 rdt (ua,va) 66 !! Note that if lk_dynspg_flt=T, the time stepping has already been 67 !! performed in dynspg module 68 !! Time filter applied on now horizontal velocity to avoid the 69 !! divergence of two consecutive time-steps and swap of dynamics 70 !! arrays to start the next time step: 71 !! (ub,vb) = (un,vn) + atfp [ (ub,vb) + (ua,va) - 2 (un,vn) ] 72 !! (un,vn) = (ua,va) 73 !! 74 !! ** Action : - Update ub,vb arrays, the before horizontal velocity 75 !! - Update un,vn arrays, the now horizontal velocity 76 !! 59 !! ** Purpose : Compute the after horizontal velocity. Apply the boundary 60 !! condition on the after velocity, achieved the time stepping 61 !! by applying the Asselin filter on now fields and swapping 62 !! the fields. 63 !! 64 !! ** Method : * After velocity is compute using a leap-frog scheme: 65 !! (ua,va) = (ub,vb) + 2 rdt (ua,va) 66 !! Note that with flux form advection and variable volume layer 67 !! (lk_vvl=T), the leap-frog is applied on thickness weighted 68 !! velocity. 69 !! Note also that in filtered free surface (lk_dynspg_flt=T), 70 !! the time stepping has already been done in dynspg module 71 !! 72 !! * Apply lateral boundary conditions on after velocity 73 !! at the local domain boundaries through lbc_lnk call, 74 !! at the radiative open boundaries (lk_obc=T), 75 !! at the relaxed open boundaries (lk_bdy=T), and 76 !! at the AGRIF zoom boundaries (lk_agrif=T) 77 !! 78 !! * Apply the time filter applied and swap of the dynamics 79 !! arrays to start the next time step: 80 !! (ub,vb) = (un,vn) + atfp [ (ub,vb) + (ua,va) - 2 (un,vn) ] 81 !! (un,vn) = (ua,va). 82 !! Note that with flux form advection and variable volume layer 83 !! (lk_vvl=T), the time filter is applied on thickness weighted 84 !! velocity. 85 !! 86 !! ** Action : ub,vb filtered before horizontal velocity of next time-step 87 !! un,vn now horizontal velocity of next time-step 77 88 !!---------------------------------------------------------------------- 78 89 INTEGER, INTENT( in ) :: kt ! ocean time-step index … … 85 96 REAL(wp) :: ze3v_b, ze3v_n, ze3v_a ! - - 86 97 REAL(wp) :: zuf , zvf ! - - 87 !!----------------------------------------------------------------------98 !!---------------------------------------------------------------------- 88 99 89 100 IF( kt == nit000 ) THEN … … 93 104 ENDIF 94 105 95 ! Local constant initialization 96 z2dt = 2. * rdt 106 #if defined key_dynspg_flt 107 ! 108 ! Next velocity : Leap-frog time stepping already done in dynspg_flt.F routine 109 ! ------------- 110 111 ! Update after velocity on domain lateral boundaries (only local domain required) 112 ! -------------------------------------------------- 113 CALL lbc_lnk( ua, 'U', -1. ) ! local domain boundaries 114 CALL lbc_lnk( va, 'V', -1. ) 115 ! 116 #else 117 ! Next velocity : Leap-frog time stepping 118 ! ------------- 119 z2dt = 2. * rdt ! Euler or leap-frog time step 97 120 IF( neuler == 0 .AND. kt == nit000 ) z2dt = rdt 98 99 !! Explicit physics with thickness weighted updates 100 101 ! Lateral boundary conditions on ( ua, va ) 102 CALL lbc_lnk( ua, 'U', -1. ) ; CALL lbc_lnk( va, 'V', -1. ) 103 104 ! Next velocity 105 ! ------------- 106 #if defined key_dynspg_flt 107 ! Leap-frog time stepping already done in dynspg_flt.F routine 108 #else 109 IF( lk_vvl ) THEN ! Varying levels 121 ! 122 IF( ln_dynadv_vec .OR. .NOT. lk_vvl ) THEN ! applied on velocity 110 123 DO jk = 1, jpkm1 111 ua(:,:,jk) = ( ub(:,:,jk) * fse3u_b(:,:,jk) &112 & + z2dt * ua(:,:,jk) * fse3u_n(:,:,jk) ) &113 & / fse3u_a(:,:,jk) * umask(:,:,jk)114 va(:,:,jk) = ( vb(:,:,jk) * fse3v_b(:,:,jk) &115 & + z2dt * va(:,:,jk) * fse3v_n(:,:,jk) ) &116 & / fse3v_a(:,:,jk) * vmask(:,:,jk)117 END DO118 ELSE119 DO jk = 1, jpkm1120 ! Leap-frog time stepping121 124 ua(:,:,jk) = ( ub(:,:,jk) + z2dt * ua(:,:,jk) ) * umask(:,:,jk) 122 125 va(:,:,jk) = ( vb(:,:,jk) + z2dt * va(:,:,jk) ) * vmask(:,:,jk) 123 126 END DO 124 ENDIF 125 127 ELSE ! applied on thickness weighted velocity 128 DO jk = 1, jpkm1 129 ua(:,:,jk) = ( ub(:,:,jk) * fse3u_b(:,:,jk) & 130 & + z2dt * ua(:,:,jk) * fse3u_n(:,:,jk) ) & 131 & / fse3u_a(:,:,jk) * umask(:,:,jk) 132 va(:,:,jk) = ( vb(:,:,jk) * fse3v_b(:,:,jk) & 133 & + z2dt * va(:,:,jk) * fse3v_n(:,:,jk) ) & 134 & / fse3v_a(:,:,jk) * vmask(:,:,jk) 135 END DO 136 ENDIF 137 138 139 ! Update after velocity on domain lateral boundaries 140 ! -------------------------------------------------- 141 CALL lbc_lnk( ua, 'U', -1. ) !* local domain boundaries 142 CALL lbc_lnk( va, 'V', -1. ) 143 ! 126 144 # if defined key_obc 127 ! Update (ua,va) along open boundaries (only in the rigid-lid case)145 ! !* OBC open boundaries 128 146 CALL obc_dyn( kt ) 129 147 ! 130 148 IF ( lk_dynspg_exp .OR. lk_dynspg_ts ) THEN 131 ! Flather boundary condition : 132 ! - Update sea surface height on each open boundary 133 ! sshn (= after ssh) for explicit case 134 ! sshn_b (= after ssha_b) for time-splitting case 135 ! - Correct the barotropic velocities 149 ! Flather boundary condition : - Update sea surface height on each open boundary 150 ! sshn (= after ssh ) for explicit case 151 ! sshn_b (= after ssha_b) for time-splitting case 152 ! - Correct the barotropic velocities 136 153 CALL obc_dyn_bt( kt ) 137 154 ! 138 ! Boundary conditions on sshn ( after ssh) 139 CALL lbc_lnk( sshn, 'T', 1. ) 155 !!gm ERROR - potential BUG: sshn should not be modified at this stage !! ssh_nxt not alrady called 156 CALL lbc_lnk( sshn, 'T', 1. ) ! Boundary conditions on sshn 140 157 ! 141 158 IF( ln_vol_cst ) CALL obc_vol( kt ) … … 143 160 IF(ln_ctl) CALL prt_ctl( tab2d_1=sshn, clinfo1=' ssh : ', mask1=tmask ) 144 161 ENDIF 145 162 ! 146 163 # elif defined key_bdy 147 ! Update (ua,va) along open boundaries (for exp or ts options).148 IF( lk_dynspg_exp .OR. lk_dynspg_ts ) THEN 164 ! !* BDY open boundaries 165 IF( lk_dynspg_exp .OR. lk_dynspg_ts ) THEN ! except for filtered option 149 166 ! 150 167 CALL bdy_dyn_frs( kt ) 151 168 ! 152 IF( ln_bdy_ fla ) THEN169 IF( ln_bdy_dyn_fla ) THEN 153 170 ua_e(:,:) = 0.e0 154 171 va_e(:,:) = 0.e0 155 172 ! Set these variables for use in bdy_dyn_fla 156 hu _e(:,:) = hu(:,:)157 hv _e(:,:) = hv(:,:)173 hur_e(:,:) = hur(:,:) 174 hvr_e(:,:) = hvr(:,:) 158 175 DO jk = 1, jpkm1 !! Vertically integrated momentum trends 159 176 ua_e(:,:) = ua_e(:,:) + fse3u(:,:,jk) * umask(:,:,jk) * ua(:,:,jk) 160 177 va_e(:,:) = va_e(:,:) + fse3v(:,:,jk) * vmask(:,:,jk) * va(:,:,jk) 161 178 END DO 179 ua_e(:,:) = ua_e(:,:) * hur(:,:) 180 va_e(:,:) = va_e(:,:) * hvr(:,:) 162 181 DO jk = 1 , jpkm1 163 ua(:,:,jk) = ua(:,:,jk) - ua_e(:,:) * hur(:,:)164 va(:,:,jk) = va(:,:,jk) - va_e(:,:) * hvr(:,:)182 ua(:,:,jk) = ua(:,:,jk) - ua_e(:,:) 183 va(:,:,jk) = va(:,:,jk) - va_e(:,:) 165 184 END DO 166 185 CALL bdy_dta_bt( kt+1, 0) 167 CALL bdy_dyn_fla 186 CALL bdy_dyn_fla( sshn_b ) 187 DO jk = 1 , jpkm1 188 ua(:,:,jk) = ( ua(:,:,jk) + ua_e(:,:) ) * umask(:,:,jk) 189 va(:,:,jk) = ( va(:,:,jk) + va_e(:,:) ) * vmask(:,:,jk) 190 END DO 168 191 ENDIF 169 192 ! 170 DO jk = 1 , jpkm1171 ua(:,:,jk) = ua(:,:,jk) + ua_e(:,:) * hur(:,:)172 va(:,:,jk) = va(:,:,jk) + va_e(:,:) * hvr(:,:)173 END DO174 193 ENDIF 175 194 # endif 176 195 ! 177 196 # if defined key_agrif 178 CALL Agrif_dyn( kt ) 197 CALL Agrif_dyn( kt ) !* AGRIF zoom boundaries 179 198 # endif 180 199 #endif … … 182 201 ! Time filter and swap of dynamics arrays 183 202 ! ------------------------------------------ 184 IF( neuler == 0 .AND. kt == nit000 ) THEN 185 DO jk = 1, jpkm1 186 un(:,:,jk) = ua(:,:,jk) 203 IF( neuler == 0 .AND. kt == nit000 ) THEN !* Euler at first time-step: only swap 204 DO jk = 1, jpkm1 205 un(:,:,jk) = ua(:,:,jk) ! un <-- ua 187 206 vn(:,:,jk) = va(:,:,jk) 188 207 END DO 189 ELSE 190 IF( lk_vvl ) THEN ! Varying levels 191 DO jk = 1, jpkm1 ! filter applied on thickness weighted velocities 208 ELSE !* Leap-Frog : Asselin filter and swap 209 IF( ln_dynadv_vec .OR. .NOT. lk_vvl ) THEN ! applied on velocity 210 DO jk = 1, jpkm1 211 DO jj = 1, jpj 212 DO ji = 1, jpi 213 zuf = atfp * ( ub(ji,jj,jk) + ua(ji,jj,jk) ) + atfp1 * un(ji,jj,jk) 214 zvf = atfp * ( vb(ji,jj,jk) + va(ji,jj,jk) ) + atfp1 * vn(ji,jj,jk) 215 ! 216 ub(ji,jj,jk) = zuf ! ub <-- filtered velocity 217 vb(ji,jj,jk) = zvf 218 un(ji,jj,jk) = ua(ji,jj,jk) ! un <-- ua 219 vn(ji,jj,jk) = va(ji,jj,jk) 220 END DO 221 END DO 222 END DO 223 ELSE ! applied on thickness weighted velocity 224 DO jk = 1, jpkm1 192 225 DO jj = 1, jpj 193 226 DO ji = 1, jpi … … 206 239 zve3b = vb(ji,jj,jk) * ze3v_b 207 240 ! 208 ub(ji,jj,jk) = ( atfp * ( zue3b + zue3a ) + atfp1 * zue3n ) & 209 & / ( atfp * ( ze3u_b + ze3u_a ) + atfp1 * ze3u_n ) * umask(ji,jj,jk) 210 vb(ji,jj,jk) = ( atfp * ( zve3b + zve3a ) + atfp1 * zve3n ) & 211 & / ( atfp * ( ze3v_b + ze3v_a ) + atfp1 * ze3v_n ) * vmask(ji,jj,jk) 212 un(ji,jj,jk) = ua(ji,jj,jk) * umask(ji,jj,jk) 213 vn(ji,jj,jk) = va(ji,jj,jk) * vmask(ji,jj,jk) 214 END DO 215 END DO 216 END DO 217 ELSE ! Fixed levels 218 DO jk = 1, jpkm1 ! filter applied on velocities 219 DO jj = 1, jpj 220 DO ji = 1, jpi 221 zuf = atfp * ( ub(ji,jj,jk) + ua(ji,jj,jk) ) + atfp1 * un(ji,jj,jk) 222 zvf = atfp * ( vb(ji,jj,jk) + va(ji,jj,jk) ) + atfp1 * vn(ji,jj,jk) 223 ub(ji,jj,jk) = zuf 241 zuf = ( atfp * ( zue3b + zue3a ) + atfp1 * zue3n ) & 242 & / ( atfp * ( ze3u_b + ze3u_a ) + atfp1 * ze3u_n ) * umask(ji,jj,jk) 243 zvf = ( atfp * ( zve3b + zve3a ) + atfp1 * zve3n ) & 244 & / ( atfp * ( ze3v_b + ze3v_a ) + atfp1 * ze3v_n ) * vmask(ji,jj,jk) 245 ! 246 ub(ji,jj,jk) = zuf ! ub <-- filtered velocity 224 247 vb(ji,jj,jk) = zvf 225 un(ji,jj,jk) = ua(ji,jj,jk) 248 un(ji,jj,jk) = ua(ji,jj,jk) ! un <-- ua 226 249 vn(ji,jj,jk) = va(ji,jj,jk) 227 250 END DO … … 232 255 233 256 #if defined key_agrif 257 ! Update velocity at AGRIF zoom boundaries 234 258 IF (.NOT.Agrif_Root()) CALL Agrif_Update_Dyn( kt ) 235 259 #endif … … 240 264 END SUBROUTINE dyn_nxt 241 265 242 !!====================================================================== 266 !!========================================================================= 243 267 END MODULE dynnxt -
trunk/NEMO/OPA_SRC/DYN/dynspg_oce.F90
r1438 r1502 40 40 41 41 #if defined key_dynspg_ts || defined key_vvl || defined key_esopa 42 !! Time splitting variables 43 !! ------------------------ 44 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: & ! variables averaged over the barotropic loop 45 un_e , vn_e, & ! vertically integrated horizontal velocities (now) 46 ua_e , va_e ! vertically integrated horizontal velocities (after) 47 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: & ! variables of the explicit barotropic loop 48 sshn_e, ssha_e, & ! sea surface heigth (now, after) 49 hu_e , hv_e ! depth arrays for the barotropic solution 42 !! Time splitting variables ! variables of the explicit barotropic loop 43 !! ------------------------ 44 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: ua_e , va_e ! barotropic velocities (after) 45 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: sshn_e, ssha_e, sshn_b ! sea surface heigth (now, after, average) 46 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: hu_e , hv_e ! depth arrays for the barotropic solution 47 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: hur_e , hvr_e ! inverse of depth arrays 50 48 #endif 51 49 -
trunk/NEMO/OPA_SRC/DYN/dynspg_ts.F90
r1438 r1502 1 1 MODULE dynspg_ts 2 2 !!====================================================================== 3 !! History : 1.0 ! 04-12 (L. Bessieres, G. Madec) Original code 4 !! - ! 05-11 (V. Garnier, G. Madec) optimization 5 !! - ! 06-08 (S. Masson) distributed restart using iom 6 !! 2.0 ! 07-07 (D. Storkey) calls to BDY routines 7 !! - ! 08-01 (R. Benshila) change averaging method 3 !! History : 1.0 ! 2004-12 (L. Bessieres, G. Madec) Original code 4 !! - ! 2005-11 (V. Garnier, G. Madec) optimization 5 !! - ! 2006-08 (S. Masson) distributed restart using iom 6 !! 2.0 ! 2007-07 (D. Storkey) calls to BDY routines 7 !! - ! 2008-01 (R. Benshila) change averaging method 8 !! 3.2 ! 2009-07 (R. Benshila, G. Madec) Complete revisit associated to vvl reactivation 8 9 !!--------------------------------------------------------------------- 9 10 #if defined key_dynspg_ts || defined key_esopa 10 11 !!---------------------------------------------------------------------- 11 12 !! 'key_dynspg_ts' free surface cst volume with time splitting 12 !!----------------------------------------------------------------------13 13 !!---------------------------------------------------------------------- 14 14 !! dyn_spg_ts : compute surface pressure gradient trend using a time- … … 45 45 PUBLIC ts_rst ! routine called by istate.F90 46 46 47 47 48 REAL(wp), DIMENSION(jpi,jpj) :: ftnw, ftne ! triad of coriolis parameter 48 49 REAL(wp), DIMENSION(jpi,jpj) :: ftsw, ftse ! (only used with een vorticity scheme) 50 51 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: un_b, vn_b ! averaged velocity 49 52 50 53 !! * Substitutions … … 66 69 !! gradient in case of free surface formulation with time-splitting. 67 70 !! Add it to the general trend of momentum equation. 68 !! Compute the free surface.69 71 !! 70 72 !! ** Method : Free surface formulation with time-splitting … … 75 77 !! the barotropic integration. 76 78 !! -2- Barotropic loop : updates of sea surface height (ssha_e) and 77 !! barotropic transports(ua_e and va_e) through barotropic79 !! barotropic velocity (ua_e and va_e) through barotropic 78 80 !! momentum and continuity integration. Barotropic former 79 81 !! variables are time averaging over the full barotropic cycle 80 !! (= 2 * baroclinic time step) and saved in z sshX_b, zuX_b82 !! (= 2 * baroclinic time step) and saved in zuX_b 81 83 !! and zvX_b (X specifying after, now or before). 82 84 !! -3- The new general trend becomes : 83 !! ua = ua - sum_k(ua)/H + ( zua_e - sum_k(ub) )/H85 !! ua = ua - sum_k(ua)/H + ( ua_e - sum_k(ub) ) 84 86 !! 85 87 !! ** Action : - Update (ua,va) with the surf. pressure gradient trend … … 87 89 !! References : Griffies et al., (2003): A technical guide to MOM4. NOAA/GFDL 88 90 !!--------------------------------------------------------------------- 89 INTEGER, INTENT( in ) :: kt ! ocean time-step index 90 !! 91 INTEGER :: ji, jj, jk, jit ! dummy loop indices 92 INTEGER :: icycle ! temporary scalar 93 REAL(wp) :: & 94 zraur, zcoef, z2dt_e, z2dt_b, zfac25, & ! temporary scalars 95 zfact1, zspgu, zcubt, zx1, zy1, & ! " " 96 zfact2, zspgv, zcvbt, zx2, zy2 ! " " 97 REAL(wp), DIMENSION(jpi,jpj) :: & 98 zcu, zcv, zwx, zwy, zhdiv, & ! temporary arrays 99 zua, zva, zub, zvb, & ! " " 100 zua_e, zva_e, zssha_e, & ! " " 101 zub_e, zvb_e, zsshb_e, & ! " " 102 zu_sum, zv_sum, zssh_sum 103 !! Variable volume 104 REAL(wp), DIMENSION(jpi,jpj) :: & 105 zspgu_1, zspgv_1, zsshun_e, zsshvn_e ! 2D workspace 91 INTEGER, INTENT(in) :: kt ! ocean time-step index 92 !! 93 INTEGER :: ji, jj, jk, jn ! dummy loop indices 94 INTEGER :: icycle ! temporary scalar 95 REAL(wp) :: zraur, zcoef, z2dt_e, z2dt_b, & ! temporary scalars 96 z1_8, zspgu, zcubt, zx1, zy1, & ! " " 97 z1_4, zspgv, zcvbt, zx2, zy2 ! " " 98 REAL(wp), DIMENSION(jpi,jpj) :: zhdiv, zsshb_e 99 !! 100 REAL(wp), DIMENSION(jpi,jpj) :: zsshun_e, zsshvn_e ! 2D workspace 101 !! 102 REAL(wp), DIMENSION(jpi,jpj) :: zcu, zwx, zua, zun, zub ! 2D workspace 103 REAL(wp), DIMENSION(jpi,jpj) :: zcv, zwy, zva, zvn, zvb ! - - 104 REAL(wp), DIMENSION(jpi,jpj) :: zun_e, zub_e, zu_sum ! 2D workspace 105 REAL(wp), DIMENSION(jpi,jpj) :: zvn_e, zvb_e, zv_sum ! - - 106 REAL(wp), DIMENSION(jpi,jpj) :: zssh_sum ! - - 106 107 !!---------------------------------------------------------------------- 107 108 108 ! Arrays initialization 109 ! --------------------- 110 zhdiv(:,:) = 0.e0 111 112 113 IF( kt == nit000 ) THEN 109 IF( kt == nit000 ) THEN !* initialisation 114 110 ! 115 111 IF(lwp) WRITE(numout,*) … … 117 113 IF(lwp) WRITE(numout,*) '~~~~~~~~~~ free surface with time splitting' 118 114 IF(lwp) WRITE(numout,*) ' Number of sub cycle in 1 time-step (2 rdt) : icycle = ', 2*nn_baro 119 115 ! 120 116 CALL ts_rst( nit000, 'READ' ) ! read or initialize the following fields: 121 ! ! sshb, sshn, sshb_b, sshn_b, un_b, vn_b 122 123 zssha_e(:,:) = sshn(:,:) 124 zua_e (:,:) = un_e(:,:) 125 zva_e (:,:) = vn_e(:,:) 126 hu_e (:,:) = hu (:,:) 127 hv_e (:,:) = hv (:,:) 117 ! ! un_b, vn_b 118 ! 119 ua_e (:,:) = un_b (:,:) 120 va_e (:,:) = vn_b (:,:) 121 hu_e (:,:) = hu (:,:) 122 hv_e (:,:) = hv (:,:) 123 hur_e (:,:) = hur (:,:) 124 hvr_e (:,:) = hvr (:,:) 128 125 IF( ln_dynvor_een ) THEN 129 126 ftne(1,:) = 0.e0 ; ftnw(1,:) = 0.e0 ; ftse(1,:) = 0.e0 ; ftsw(1,:) = 0.e0 … … 140 137 ENDIF 141 138 142 ! Substract the surface pressure gradient from the momentum 143 ! --------------------------------------------------------- 144 IF( lk_vvl ) THEN ! Variable volume 145 146 ! 0. Local initialization 147 ! ----------------------- 148 zspgu_1(:,:) = 0.e0 149 zspgv_1(:,:) = 0.e0 150 151 ! 1. Surface pressure gradient (now) 152 ! ---------------------------- 153 DO jj = 2, jpjm1 154 DO ji = fs_2, fs_jpim1 ! vector opt. 155 zspgu_1(ji,jj) = -grav * ( ( rhd(ji+1,jj ,1) + 1 ) * sshn(ji+1,jj ) & 156 & - ( rhd(ji ,jj ,1) + 1 ) * sshn(ji ,jj ) ) / e1u(ji,jj) 157 zspgv_1(ji,jj) = -grav * ( ( rhd(ji ,jj+1,1) + 1 ) * sshn(ji ,jj+1) & 158 & - ( rhd(ji ,jj ,1) + 1 ) * sshn(ji ,jj ) ) / e2v(ji,jj) 159 END DO 160 END DO 161 162 ! 2. Substract the surface pressure trend from the general trend 163 ! ------------------------------------------------------ 164 DO jk = 1, jpkm1 165 DO jj = 2, jpjm1 166 DO ji = fs_2, fs_jpim1 ! vector opt. 167 ua(ji,jj,jk) = ua(ji,jj,jk) - zspgu_1(ji,jj) 168 va(ji,jj,jk) = va(ji,jj,jk) - zspgv_1(ji,jj) 169 END DO 170 END DO 171 END DO 172 173 ENDIF 174 175 ! Local constant initialization 176 ! -------------------------------- 139 ! !* Local constant initialization 177 140 z2dt_b = 2.0 * rdt ! baroclinic time step 178 IF ( neuler == 0 .AND. kt == nit000 ) z2dt_b = rdt 179 zfact1 = 0.5 * 0.25 ! coefficient for vorticity estimates 180 zfact2 = 0.5 * 0.5 141 z1_8 = 0.5 * 0.25 ! coefficient for vorticity estimates 142 z1_4 = 0.5 * 0.5 181 143 zraur = 1. / rauw ! 1 / volumic mass of pure water 144 ! 145 zhdiv(:,:) = 0.e0 ! barotropic divergence 182 146 183 147 ! ----------------------------------------------------------------------------- 184 148 ! Phase 1 : Coupling between general trend and barotropic estimates (1st step) 185 149 ! ----------------------------------------------------------------------------- 186 187 ! Vertically integrated quantities 188 ! -------------------------------- 189 zua(:,:) = 0.e0 190 zva(:,:) = 0.e0 191 zub(:,:) = 0.e0 192 zvb(:,:) = 0.e0 193 zwx(:,:) = 0.e0 194 zwy(:,:) = 0.e0 195 196 ! vertical sum 197 IF( lk_vopt_loop ) THEN ! vector opt., forced unroll 198 DO jk = 1, jpkm1 150 ! 151 ! !* e3*d/dt(Ua), e3*Ub, e3*Vn (Vertically integrated) 152 ! ! -------------------------- 153 zua(:,:) = 0.e0 ; zun(:,:) = 0.e0 ; zub(:,:) = 0.e0 154 zva(:,:) = 0.e0 ; zvn(:,:) = 0.e0 ; zvb(:,:) = 0.e0 155 ! 156 DO jk = 1, jpkm1 157 #if defined key_vectopt_loop 158 DO jj = 1, 1 !Vector opt. => forced unrolling 199 159 DO ji = 1, jpij 200 ! ! Vertically integrated momentum trends 201 zua(ji,1) = zua(ji,1) + fse3u(ji,1,jk) * umask(ji,1,jk) * ua(ji,1,jk) 202 zva(ji,1) = zva(ji,1) + fse3v(ji,1,jk) * vmask(ji,1,jk) * va(ji,1,jk) 203 ! ! Vertically integrated transports (before) 204 zub(ji,1) = zub(ji,1) + fse3u_b(ji,1,jk) * ub(ji,1,jk) 205 zvb(ji,1) = zvb(ji,1) + fse3v_b(ji,1,jk) * vb(ji,1,jk) 206 ! ! Planetary vorticity transport fluxes (now) 207 zwx(ji,1) = zwx(ji,1) + e2u(ji,1) * fse3u(ji,1,jk) * un(ji,1,jk) 208 zwy(ji,1) = zwy(ji,1) + e1v(ji,1) * fse3v(ji,1,jk) * vn(ji,1,jk) 209 END DO 210 END DO 211 ELSE ! No vector opt. 212 DO jk = 1, jpkm1 213 ! ! Vertically integrated momentum trends 214 zua(:,:) = zua(:,:) + fse3u(:,:,jk) * umask(:,:,jk) * ua(:,:,jk) 215 zva(:,:) = zva(:,:) + fse3v(:,:,jk) * vmask(:,:,jk) * va(:,:,jk) 216 ! ! Vertically integrated transports (before) 217 zub(:,:) = zub(:,:) + fse3u_b(:,:,jk) * ub(:,:,jk) 218 zvb(:,:) = zvb(:,:) + fse3v_b(:,:,jk) * vb(:,:,jk) 219 ! ! Planetary vorticity (now) 220 zwx(:,:) = zwx(:,:) + e2u(:,:) * fse3u(:,:,jk) * un(:,:,jk) 221 zwy(:,:) = zwy(:,:) + e1v(:,:) * fse3v(:,:,jk) * vn(:,:,jk) 222 END DO 223 ENDIF 224 160 #else 161 DO jj = 1, jpj 162 DO ji = 1, jpi 163 #endif 164 ! ! now trend 165 zua(ji,jj) = zua(ji,jj) + fse3u (ji,jj,jk) * ua(ji,jj,jk) * umask(ji,jj,jk) 166 zva(ji,jj) = zva(ji,jj) + fse3v (ji,jj,jk) * va(ji,jj,jk) * vmask(ji,jj,jk) 167 ! ! now velocity 168 zun(ji,jj) = zun(ji,jj) + fse3u (ji,jj,jk) * un(ji,jj,jk) 169 zvn(ji,jj) = zvn(ji,jj) + fse3v (ji,jj,jk) * vn(ji,jj,jk) 170 ! ! before velocity 171 zub(ji,jj) = zub(ji,jj) + fse3u_b(ji,jj,jk) * ub(ji,jj,jk) 172 zvb(ji,jj) = zvb(ji,jj) + fse3v_b(ji,jj,jk) * vb(ji,jj,jk) 173 END DO 174 END DO 175 END DO 176 177 ! !* baroclinic momentum trend (remove the vertical mean trend) 178 DO jk = 1, jpkm1 ! -------------------------- 179 DO jj = 2, jpjm1 180 DO ji = fs_2, fs_jpim1 ! vector opt. 181 ua(ji,jj,jk) = ua(ji,jj,jk) - zua(ji,jj) * hur(ji,jj) 182 va(ji,jj,jk) = va(ji,jj,jk) - zva(ji,jj) * hvr(ji,jj) 183 END DO 184 END DO 185 END DO 186 187 ! !* barotropic Coriolis trends * H (vorticity scheme dependent) 188 ! ! ---------------------------==== 189 zwx(:,:) = zun(:,:) * e2u(:,:) ! now transport 190 zwy(:,:) = zvn(:,:) * e1v(:,:) 191 ! 225 192 IF( ln_dynvor_ene .OR. ln_dynvor_mix ) THEN ! energy conserving or mixed scheme 226 193 DO jj = 2, jpjm1 … … 231 198 zx2 = ( zwx(ji ,jj) + zwx(ji ,jj+1) ) / e2v(ji,jj) 232 199 ! energy conserving formulation for planetary vorticity term 233 zcu(ji,jj) = z fact2* ( ff(ji ,jj-1) * zy1 + ff(ji,jj) * zy2 )234 zcv(ji,jj) =-z fact2* ( ff(ji-1,jj ) * zx1 + ff(ji,jj) * zx2 )200 zcu(ji,jj) = z1_4 * ( ff(ji ,jj-1) * zy1 + ff(ji,jj) * zy2 ) 201 zcv(ji,jj) =-z1_4 * ( ff(ji-1,jj ) * zx1 + ff(ji,jj) * zx2 ) 235 202 END DO 236 203 END DO … … 239 206 DO jj = 2, jpjm1 240 207 DO ji = fs_2, fs_jpim1 ! vector opt. 241 zy1 = zfact1 * ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) & 242 + zwy(ji ,jj ) + zwy(ji+1,jj ) ) / e1u(ji,jj) 243 zx1 =-zfact1 * ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) & 244 + zwx(ji ,jj ) + zwx(ji ,jj+1) ) / e2v(ji,jj) 208 zy1 = z1_8 * ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) + zwy(ji,jj) + zwy(ji+1,jj ) ) / e1u(ji,jj) 209 zx1 = - z1_8 * ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) + zwx(ji,jj) + zwx(ji ,jj+1) ) / e2v(ji,jj) 245 210 zcu(ji,jj) = zy1 * ( ff(ji ,jj-1) + ff(ji,jj) ) 246 211 zcv(ji,jj) = zx1 * ( ff(ji-1,jj ) + ff(ji,jj) ) … … 249 214 ! 250 215 ELSEIF ( ln_dynvor_een ) THEN ! enstrophy and energy conserving scheme 251 zfac25 = 0.25252 216 DO jj = 2, jpjm1 253 217 DO ji = fs_2, fs_jpim1 ! vector opt. 254 zcu(ji,jj) = + zfac25 / e1u(ji,jj) & 255 & * ( ftne(ji,jj ) * zwy(ji ,jj ) + ftnw(ji+1,jj) * zwy(ji+1,jj ) & 256 & + ftse(ji,jj ) * zwy(ji ,jj-1) + ftsw(ji+1,jj) * zwy(ji+1,jj-1) ) 257 zcv(ji,jj) = - zfac25 / e2v(ji,jj) & 258 & * ( ftsw(ji,jj+1) * zwx(ji-1,jj+1) + ftse(ji,jj+1) * zwx(ji ,jj+1) & 259 & + ftnw(ji,jj ) * zwx(ji-1,jj ) + ftne(ji,jj ) * zwx(ji ,jj ) ) 218 zcu(ji,jj) = + z1_4 / e1u(ji,jj) * ( ftne(ji,jj ) * zwy(ji ,jj ) + ftnw(ji+1,jj) * zwy(ji+1,jj ) & 219 & + ftse(ji,jj ) * zwy(ji ,jj-1) + ftsw(ji+1,jj) * zwy(ji+1,jj-1) ) 220 zcv(ji,jj) = - z1_4 / e2v(ji,jj) * ( ftsw(ji,jj+1) * zwx(ji-1,jj+1) + ftse(ji,jj+1) * zwx(ji ,jj+1) & 221 & + ftnw(ji,jj ) * zwx(ji-1,jj ) + ftne(ji,jj ) * zwx(ji ,jj ) ) 260 222 END DO 261 223 END DO … … 263 225 ENDIF 264 226 265 266 ! Remove barotropic trend from general momentum trend 267 ! --------------------------------------------------- 268 DO jk = 1 , jpkm1 269 DO jj = 2, jpjm1 270 DO ji = fs_2, fs_jpim1 ! vector opt. 271 ua(ji,jj,jk) = ua(ji,jj,jk) - zua(ji,jj) * hur(ji,jj) 272 va(ji,jj,jk) = va(ji,jj,jk) - zva(ji,jj) * hvr(ji,jj) 273 END DO 274 END DO 275 END DO 276 277 ! Remove coriolis term from barotropic trend 278 ! ------------------------------------------ 279 DO jj = 2, jpjm1 227 ! !* Right-Hand-Side of the barotropic momentum equation 228 ! ! ---------------------------------------------------- 229 IF( lk_vvl ) THEN ! Variable volume : remove both Coriolis and Surface pressure gradient 230 DO jj = 2, jpjm1 231 DO ji = fs_2, fs_jpim1 ! vector opt. 232 zcu(ji,jj) = zcu(ji,jj) - grav * ( ( rhd(ji+1,jj ,1) + 1 ) * sshn(ji+1,jj ) & 233 & - ( rhd(ji ,jj ,1) + 1 ) * sshn(ji ,jj ) ) * hu(ji,jj) / e1u(ji,jj) 234 zcv(ji,jj) = zcv(ji,jj) - grav * ( ( rhd(ji ,jj+1,1) + 1 ) * sshn(ji ,jj+1) & 235 & - ( rhd(ji ,jj ,1) + 1 ) * sshn(ji ,jj ) ) * hv(ji,jj) / e2v(ji,jj) 236 END DO 237 END DO 238 ENDIF 239 240 DO jj = 2, jpjm1 ! Remove coriolis term (and possibly spg) from barotropic trend 280 241 DO ji = fs_2, fs_jpim1 281 242 zua(ji,jj) = zua(ji,jj) - zcu(ji,jj) … … 284 245 END DO 285 246 247 ! !* d/dt(Ua), Ub, Vn (Vertical mean velocity) 248 ! ! -------------------------- 249 zua(:,:) = zua(:,:) * hur(:,:) 250 zva(:,:) = zva(:,:) * hvr(:,:) 251 ! 252 IF( lk_vvl ) THEN 253 zub(:,:) = zub(:,:) * umask(:,:,1) / ( hu_0(:,:) + sshu_b(:,:) + 1.e0 - umask(:,:,1) ) 254 zvb(:,:) = zvb(:,:) * vmask(:,:,1) / ( hv_0(:,:) + sshv_b(:,:) + 1.e0 - vmask(:,:,1) ) 255 ELSE 256 zub(:,:) = zub(:,:) * hur(:,:) 257 zvb(:,:) = zvb(:,:) * hvr(:,:) 258 ENDIF 259 286 260 ! ----------------------------------------------------------------------- 287 261 ! Phase 2 : Integration of the barotropic equations with time splitting 288 262 ! ----------------------------------------------------------------------- 289 290 ! Initialisations 291 !---------------- 292 ! Number of iteration of the barotropic loop 293 icycle = 2 * nn_baro + 1 294 295 ! variables for the barotropic equations 296 zu_sum (:,:) = 0.e0 297 zv_sum (:,:) = 0.e0 298 zssh_sum(:,:) = 0.e0 299 hu_e (:,:) = hu (:,:) ! (barotropic) ocean depth at u-point 300 hv_e (:,:) = hv (:,:) ! (barotropic) ocean depth at v-point 301 zsshb_e (:,:) = sshn_e(:,:) ! (barotropic) sea surface height (before and now) 302 ! vertical sum 303 un_e (:,:) = 0.e0 304 vn_e (:,:) = 0.e0 305 IF( lk_vopt_loop ) THEN ! vector opt., forced unroll 306 DO jk = 1, jpkm1 307 DO ji = 1, jpij 308 un_e(ji,1) = un_e(ji,1) + fse3u(ji,1,jk) * un(ji,1,jk) 309 vn_e(ji,1) = vn_e(ji,1) + fse3v(ji,1,jk) * vn(ji,1,jk) 310 END DO 311 END DO 312 ELSE ! No vector opt. 313 DO jk = 1, jpkm1 314 un_e(:,:) = un_e(:,:) + fse3u(:,:,jk) * un(:,:,jk) 315 vn_e(:,:) = vn_e(:,:) + fse3v(:,:,jk) * vn(:,:,jk) 316 END DO 317 ENDIF 318 zub_e (:,:) = un_e(:,:) 319 zvb_e (:,:) = vn_e(:,:) 320 321 263 ! 264 ! ! ==================== ! 265 ! ! Initialisations ! 266 ! ! ==================== ! 267 icycle = 2 * nn_baro ! Number of barotropic sub time-step 268 269 ! ! Start from NOW field 270 hu_e (:,:) = hu (:,:) ! ocean depth at u- and v-points 271 hv_e (:,:) = hv (:,:) 272 hur_e (:,:) = hur (:,:) ! ocean depth inverted at u- and v-points 273 hvr_e (:,:) = hvr (:,:) 274 zsshb_e(:,:) = sshn (:,:) ! sea surface height (before and now) 275 sshn_e (:,:) = sshn (:,:) 276 277 zun_e (:,:) = un_b (:,:) ! barotropic velocity (external) 278 zvn_e (:,:) = vn_b (:,:) 279 zub_e (:,:) = un_b (:,:) 280 zvb_e (:,:) = vn_b (:,:) 281 282 zu_sum (:,:) = un_b (:,:) ! summation 283 zv_sum (:,:) = vn_b (:,:) 284 zssh_sum(:,:) = sshn (:,:) 285 286 #if defined key_obc 322 287 ! set ssh corrections to 0 323 288 ! ssh corrections are applied to normal velocities (Flather's algorithm) and averaged over the barotropic loop 324 #if defined key_obc325 289 IF( lp_obc_east ) sshfoe_b(:,:) = 0.e0 326 290 IF( lp_obc_west ) sshfow_b(:,:) = 0.e0 … … 329 293 #endif 330 294 331 ! Barotropic integration over 2 baroclinic time steps 332 ! --------------------------------------------------- 333 334 ! ! ==================== ! 335 DO jit = 1, icycle ! sub-time-step loop ! 336 ! ! ==================== ! 295 ! ! ==================== ! 296 DO jn = 1, icycle ! sub-time-step loop ! (from NOW to AFTER+1) 297 ! ! ==================== ! 337 298 z2dt_e = 2. * ( rdt / nn_baro ) 338 IF ( jit == 1 ) z2dt_e = rdt / nn_baro 339 340 ! Time interpolation of open boundary condition data 341 IF( lk_obc ) CALL obc_dta_bt( kt, jit ) 342 IF( lk_bdy .OR. ln_bdy_tides) CALL bdy_dta_bt( kt, jit+1 ) 343 344 ! Horizontal divergence of barotropic transports 345 !-------------------------------------------------- 346 zhdiv(:,:) = 0.e0 347 DO jj = 2, jpjm1 348 DO ji = fs_2, fs_jpim1 ! vector opt. 349 zhdiv(ji,jj) = ( e2u(ji ,jj ) * un_e(ji ,jj) & 350 & -e2u(ji-1,jj ) * un_e(ji-1,jj) & 351 & +e1v(ji ,jj ) * vn_e(ji ,jj) & 352 & -e1v(ji ,jj-1) * vn_e(ji ,jj-1) ) & 353 & / (e1t(ji,jj)*e2t(ji,jj)) 354 END DO 355 END DO 356 299 IF( jn == 1 ) z2dt_e = rdt / nn_baro 300 301 ! !* Update the forcing (OBC, BDY and tides) 302 ! ! ------------------ 303 IF( lk_obc ) CALL obc_dta_bt( kt, jn ) 304 IF( lk_bdy .OR. ln_bdy_tides ) CALL bdy_dta_bt( kt, jn+1 ) 305 306 ! !* after ssh_e 307 ! ! ----------- 308 DO jj = 2, jpjm1 ! Horizontal divergence of barotropic transports 309 DO ji = fs_2, fs_jpim1 ! vector opt. 310 zhdiv(ji,jj) = ( e2u(ji ,jj) * zun_e(ji ,jj) * hu_e(ji ,jj) & 311 & - e2u(ji-1,jj) * zun_e(ji-1,jj) * hu_e(ji-1,jj) & 312 & + e1v(ji,jj ) * zvn_e(ji,jj ) * hv_e(ji,jj ) & 313 & - e1v(ji,jj-1) * zvn_e(ji,jj-1) * hv_e(ji,jj-1) ) / ( e1t(ji,jj) * e2t(ji,jj) ) 314 END DO 315 END DO 316 ! 357 317 #if defined key_obc 358 ! open boundaries (div must be zero behind the open boundary)359 360 IF( lp_obc_east ) zhdiv(nie0p1:nie1p1,nje0 :nje1 )= 0.e0 ! east361 IF( lp_obc_west ) zhdiv(niw0 :niw1 ,njw0 :njw1 )= 0.e0 ! west318 ! ! OBC : zhdiv must be zero behind the open boundary 319 !! mpp remark: The zeroing of hdiv can probably be extended to 1->jpi/jpj for the correct row/column 320 IF( lp_obc_east ) zhdiv(nie0p1:nie1p1,nje0 :nje1 ) = 0.e0 ! east 321 IF( lp_obc_west ) zhdiv(niw0 :niw1 ,njw0 :njw1 ) = 0.e0 ! west 362 322 IF( lp_obc_north ) zhdiv(nin0 :nin1 ,njn0p1:njn1p1) = 0.e0 ! north 363 IF( lp_obc_south ) zhdiv(nis0 :nis1 ,njs0 :njs1 )= 0.e0 ! south323 IF( lp_obc_south ) zhdiv(nis0 :nis1 ,njs0 :njs1 ) = 0.e0 ! south 364 324 #endif 365 366 325 #if defined key_bdy 367 DO jj = 1, jpj 368 DO ji = 1, jpi 369 zhdiv(ji,jj) = zhdiv(ji,jj)*bdytmask(ji,jj) 370 END DO 371 END DO 326 zhdiv(:,:) = zhdiv(:,:) * bdytmask(:,:) ! BDY mask 372 327 #endif 373 374 ! Sea surface height from the barotropic system 375 !---------------------------------------------- 376 DO jj = 2, jpjm1 377 DO ji = fs_2, fs_jpim1 ! vector opt. 378 zssha_e(ji,jj) = ( zsshb_e(ji,jj) - z2dt_e * ( zraur * emp(ji,jj) & 379 & + zhdiv(ji,jj) ) ) * tmask(ji,jj,1) 380 END DO 381 END DO 382 383 ! evolution of the barotropic transport ( following the vorticity scheme used) 384 ! ---------------------------------------------------------------------------- 385 zwx(:,:) = e2u(:,:) * un_e(:,:) 386 zwy(:,:) = e1v(:,:) * vn_e(:,:) 387 388 IF( ln_dynvor_ene .OR. ln_dynvor_mix ) THEN ! energy conserving or mixed scheme 328 ! 329 DO jj = 2, jpjm1 ! leap-frog on ssh_e 330 DO ji = fs_2, fs_jpim1 ! vector opt. 331 ssha_e(ji,jj) = ( zsshb_e(ji,jj) - z2dt_e * ( zraur * emp(ji,jj) + zhdiv(ji,jj) ) ) * tmask(ji,jj,1) 332 END DO 333 END DO 334 335 ! !* after barotropic velocities (vorticity scheme dependent) 336 ! ! --------------------------- 337 zwx(:,:) = e2u(:,:) * zun_e(:,:) * hu_e(:,:) ! now_e transport 338 zwy(:,:) = e1v(:,:) * zvn_e(:,:) * hv_e(:,:) 339 ! 340 IF( ln_dynvor_ene .OR. ln_dynvor_mix ) THEN !== energy conserving or mixed scheme ==! 389 341 DO jj = 2, jpjm1 390 342 DO ji = fs_2, fs_jpim1 ! vector opt. 391 343 ! surface pressure gradient 392 344 IF( lk_vvl) THEN 393 zspgu = -grav * ( ( rhd(ji+1,jj ,1) + 1 ) * sshn_e(ji+1,jj )&394 & - ( rhd(ji ,jj ,1) + 1 ) * sshn_e(ji ,jj ) ) * hu(ji,jj) / e1u(ji,jj)395 zspgv = -grav * ( ( rhd(ji ,jj+1,1) + 1 ) * sshn_e(ji ,jj+1)&396 & - ( rhd(ji ,jj ,1) + 1 ) * sshn_e(ji ,jj ) ) * hv(ji,jj) / e2v(ji,jj)345 zspgu = -grav * ( ( rhd(ji+1,jj ,1) + 1 ) * sshn_e(ji+1,jj ) & 346 & - ( rhd(ji ,jj ,1) + 1 ) * sshn_e(ji ,jj ) ) / e1u(ji,jj) 347 zspgv = -grav * ( ( rhd(ji ,jj+1,1) + 1 ) * sshn_e(ji ,jj+1) & 348 & - ( rhd(ji ,jj ,1) + 1 ) * sshn_e(ji ,jj ) ) / e2v(ji,jj) 397 349 ELSE 398 zspgu = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) * hu(ji,jj)/ e1u(ji,jj)399 zspgv = -grav * ( sshn_e(ji,jj+1) - sshn_e(ji,jj) ) * hv(ji,jj)/ e2v(ji,jj)350 zspgu = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) / e1u(ji,jj) 351 zspgv = -grav * ( sshn_e(ji,jj+1) - sshn_e(ji,jj) ) / e2v(ji,jj) 400 352 ENDIF 401 353 ! energy conserving formulation for planetary vorticity term … … 404 356 zx1 = ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) ) / e2v(ji,jj) 405 357 zx2 = ( zwx(ji ,jj ) + zwx(ji ,jj+1) ) / e2v(ji,jj) 406 zcubt = z fact2 * ( ff(ji ,jj-1) * zy1 + ff(ji,jj) * zy2)407 zcvbt =-z fact2 * ( ff(ji-1,jj ) * zx1 + ff(ji,jj) * zx2)408 ! after transports409 zua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zcubt + zspgu + zua(ji,jj) ) ) * umask(ji,jj,1)410 zva_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zcvbt + zspgv + zva(ji,jj) ) ) * vmask(ji,jj,1)358 zcubt = z1_4 * ( ff(ji ,jj-1) * zy1 + ff(ji,jj) * zy2 ) * hur_e(ji,jj) 359 zcvbt =-z1_4 * ( ff(ji-1,jj ) * zx1 + ff(ji,jj) * zx2 ) * hvr_e(ji,jj) 360 ! after barotropic velocity 361 ua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zcubt + zspgu + zua(ji,jj) ) ) * umask(ji,jj,1) 362 va_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zcvbt + zspgv + zva(ji,jj) ) ) * vmask(ji,jj,1) 411 363 END DO 412 364 END DO 413 365 ! 414 ELSEIF ( ln_dynvor_ens ) THEN ! enstrophy conserving scheme 366 ELSEIF ( ln_dynvor_ens ) THEN !== enstrophy conserving scheme ==! 367 DO jj = 2, jpjm1 368 DO ji = fs_2, fs_jpim1 ! vector opt. 369 ! surface pressure gradient 370 IF( lk_vvl) THEN 371 zspgu = -grav * ( ( rhd(ji+1,jj ,1) + 1 ) * sshn_e(ji+1,jj ) & 372 & - ( rhd(ji ,jj ,1) + 1 ) * sshn_e(ji ,jj ) ) / e1u(ji,jj) 373 zspgv = -grav * ( ( rhd(ji ,jj+1,1) + 1 ) * sshn_e(ji ,jj+1) & 374 & - ( rhd(ji ,jj ,1) + 1 ) * sshn_e(ji ,jj ) ) / e2v(ji,jj) 375 ELSE 376 zspgu = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) / e1u(ji,jj) 377 zspgv = -grav * ( sshn_e(ji,jj+1) - sshn_e(ji,jj) ) / e2v(ji,jj) 378 ENDIF 379 ! enstrophy conserving formulation for planetary vorticity term 380 zy1 = z1_8 * ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) + zwy(ji,jj) + zwy(ji+1,jj ) ) / e1u(ji,jj) 381 zx1 = - z1_8 * ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) + zwx(ji,jj) + zwx(ji ,jj+1) ) / e2v(ji,jj) 382 zcubt = zy1 * ( ff(ji ,jj-1) + ff(ji,jj) ) * hur_e(ji,jj) 383 zcvbt = zx1 * ( ff(ji-1,jj ) + ff(ji,jj) ) * hvr_e(ji,jj) 384 ! after barotropic velocity 385 ua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zcubt + zspgu + zua(ji,jj) ) ) * umask(ji,jj,1) 386 va_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zcvbt + zspgv + zva(ji,jj) ) ) * vmask(ji,jj,1) 387 END DO 388 END DO 389 ! 390 ELSEIF ( ln_dynvor_een ) THEN !== energy and enstrophy conserving scheme ==! 415 391 DO jj = 2, jpjm1 416 392 DO ji = fs_2, fs_jpim1 ! vector opt. 417 393 ! surface pressure gradient 418 394 IF( lk_vvl) THEN 419 zspgu = -grav * ( ( rhd(ji+1,jj ,1) + 1 ) * sshn_e(ji+1,jj )&420 & - ( rhd(ji ,jj ,1) + 1 ) * sshn_e(ji ,jj ) ) * hu(ji,jj) / e1u(ji,jj)421 zspgv = -grav * ( ( rhd(ji ,jj+1,1) + 1 ) * sshn_e(ji ,jj+1)&422 & - ( rhd(ji ,jj ,1) + 1 ) * sshn_e(ji ,jj ) ) * hv(ji,jj) / e2v(ji,jj)395 zspgu = -grav * ( ( rhd(ji+1,jj ,1) + 1 ) * sshn_e(ji+1,jj ) & 396 & - ( rhd(ji ,jj ,1) + 1 ) * sshn_e(ji ,jj ) ) / e1u(ji,jj) 397 zspgv = -grav * ( ( rhd(ji ,jj+1,1) + 1 ) * sshn_e(ji ,jj+1) & 398 & - ( rhd(ji ,jj ,1) + 1 ) * sshn_e(ji ,jj ) ) / e2v(ji,jj) 423 399 ELSE 424 zspgu = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) * hu(ji,jj) / e1u(ji,jj) 425 zspgv = -grav * ( sshn_e(ji,jj+1) - sshn_e(ji,jj) ) * hv(ji,jj) / e2v(ji,jj) 426 ENDIF 427 ! enstrophy conserving formulation for planetary vorticity term 428 zy1 = zfact1 * ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) & 429 + zwy(ji ,jj ) + zwy(ji+1,jj ) ) / e1u(ji,jj) 430 zx1 =-zfact1 * ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) & 431 + zwx(ji ,jj ) + zwx(ji ,jj+1) ) / e2v(ji,jj) 432 zcubt = zy1 * ( ff(ji ,jj-1) + ff(ji,jj) ) 433 zcvbt = zx1 * ( ff(ji-1,jj ) + ff(ji,jj) ) 434 ! after transports 435 zua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zcubt + zspgu + zua(ji,jj) ) ) * umask(ji,jj,1) 436 zva_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zcvbt + zspgv + zva(ji,jj) ) ) * vmask(ji,jj,1) 437 END DO 438 END DO 439 ! 440 ELSEIF ( ln_dynvor_een ) THEN ! energy and enstrophy conserving scheme 441 zfac25 = 0.25 442 DO jj = 2, jpjm1 443 DO ji = fs_2, fs_jpim1 ! vector opt. 444 ! surface pressure gradient 445 IF( lk_vvl) THEN 446 zspgu = -grav * ( ( rhd(ji+1,jj ,1) + 1 ) * sshn_e(ji+1,jj ) & 447 & - ( rhd(ji ,jj ,1) + 1 ) * sshn_e(ji ,jj ) ) * hu(ji,jj) / e1u(ji,jj) 448 zspgv = -grav * ( ( rhd(ji ,jj+1,1) + 1 ) * sshn_e(ji ,jj+1) & 449 & - ( rhd(ji ,jj ,1) + 1 ) * sshn_e(ji ,jj ) ) * hv(ji,jj) / e2v(ji,jj) 450 ELSE 451 zspgu = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) * hu(ji,jj) / e1u(ji,jj) 452 zspgv = -grav * ( sshn_e(ji,jj+1) - sshn_e(ji,jj) ) * hv(ji,jj) / e2v(ji,jj) 400 zspgu = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) / e1u(ji,jj) 401 zspgv = -grav * ( sshn_e(ji,jj+1) - sshn_e(ji,jj) ) / e2v(ji,jj) 453 402 ENDIF 454 403 ! energy/enstrophy conserving formulation for planetary vorticity term 455 zcubt = + z fac25/ e1u(ji,jj) * ( ftne(ji,jj ) * zwy(ji ,jj ) + ftnw(ji+1,jj) * zwy(ji+1,jj ) &456 & + ftse(ji,jj ) * zwy(ji ,jj-1) + ftsw(ji+1,jj) * zwy(ji+1,jj-1))457 zcvbt = - z fac25/ e2v(ji,jj) * ( ftsw(ji,jj+1) * zwx(ji-1,jj+1) + ftse(ji,jj+1) * zwx(ji ,jj+1) &458 & + ftnw(ji,jj ) * zwx(ji-1,jj ) + ftne(ji,jj ) * zwx(ji ,jj ))459 ! after transports460 zua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zcubt + zspgu + zua(ji,jj) ) ) * umask(ji,jj,1)461 zva_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zcvbt + zspgv + zva(ji,jj) ) ) * vmask(ji,jj,1)404 zcubt = + z1_4 / e1u(ji,jj) * ( ftne(ji,jj ) * zwy(ji ,jj ) + ftnw(ji+1,jj) * zwy(ji+1,jj ) & 405 & + ftse(ji,jj ) * zwy(ji ,jj-1) + ftsw(ji+1,jj) * zwy(ji+1,jj-1) ) * hur_e(ji,jj) 406 zcvbt = - z1_4 / e2v(ji,jj) * ( ftsw(ji,jj+1) * zwx(ji-1,jj+1) + ftse(ji,jj+1) * zwx(ji ,jj+1) & 407 & + ftnw(ji,jj ) * zwx(ji-1,jj ) + ftne(ji,jj ) * zwx(ji ,jj ) ) * hvr_e(ji,jj) 408 ! after barotropic velocity 409 ua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zcubt + zspgu + zua(ji,jj) ) ) * umask(ji,jj,1) 410 va_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zcvbt + zspgv + zva(ji,jj) ) ) * vmask(ji,jj,1) 462 411 END DO 463 412 END DO 464 413 ! 465 414 ENDIF 466 467 ! Flather's boundary condition for the barotropic loop : 468 ! - Update sea surface height on each open boundary 469 ! - Correct the barotropic transports 470 IF( lk_obc ) CALL obc_fla_ts 471 IF( lk_bdy .OR. ln_bdy_tides ) CALL bdy_dyn_fla 472 473 ! ... Boundary conditions on ua_e, va_e, ssha_e 474 CALL lbc_lnk( zua_e , 'U', -1. ) 475 CALL lbc_lnk( zva_e , 'V', -1. ) 476 CALL lbc_lnk( zssha_e, 'T', 1. ) 477 478 ! temporal sum 479 !------------- 480 zu_sum (:,:) = zu_sum (:,:) + zua_e (:,:) 481 zv_sum (:,:) = zv_sum (:,:) + zva_e (:,:) 482 zssh_sum(:,:) = zssh_sum(:,:) + zssha_e(:,:) 483 484 ! Time filter and swap of dynamics arrays 485 ! --------------------------------------- 486 IF( jit == 1 ) THEN ! Euler (forward) time stepping 487 zsshb_e(:,:) = sshn_e (:,:) 488 zub_e (:,:) = un_e (:,:) 489 zvb_e (:,:) = vn_e (:,:) 490 sshn_e (:,:) = zssha_e(:,:) 491 un_e (:,:) = zua_e (:,:) 492 vn_e (:,:) = zva_e (:,:) 493 ELSE ! Asselin filtering 494 zsshb_e(:,:) = atfp * ( zsshb_e(:,:) + zssha_e(:,:) ) + atfp1 * sshn_e(:,:) 495 zub_e (:,:) = atfp * ( zub_e (:,:) + zua_e (:,:) ) + atfp1 * un_e (:,:) 496 zvb_e (:,:) = atfp * ( zvb_e (:,:) + zva_e (:,:) ) + atfp1 * vn_e (:,:) 497 sshn_e (:,:) = zssha_e(:,:) 498 un_e (:,:) = zua_e (:,:) 499 vn_e (:,:) = zva_e (:,:) 415 ! !* domain lateral boundary 416 ! ! ----------------------- 417 ! ! Flather's boundary condition for the barotropic loop : 418 ! ! - Update sea surface height on each open boundary 419 ! ! - Correct the velocity 420 421 IF( lk_obc ) CALL obc_fla_ts 422 IF( lk_bdy .OR. ln_bdy_tides ) CALL bdy_dyn_fla( sshn_e ) 423 ! 424 CALL lbc_lnk( ua_e , 'U', -1. ) ! local domain boundaries 425 CALL lbc_lnk( va_e , 'V', -1. ) 426 CALL lbc_lnk( ssha_e, 'T', 1. ) 427 428 zu_sum (:,:) = zu_sum (:,:) + ua_e (:,:) ! Sum over sub-time-steps 429 zv_sum (:,:) = zv_sum (:,:) + va_e (:,:) 430 zssh_sum(:,:) = zssh_sum(:,:) + ssha_e(:,:) 431 432 ! !* Time filter and swap 433 ! ! -------------------- 434 IF( jn == 1 ) THEN ! Swap only (1st Euler time step) 435 zsshb_e(:,:) = sshn_e(:,:) 436 zub_e (:,:) = zun_e (:,:) 437 zvb_e (:,:) = zvn_e (:,:) 438 sshn_e (:,:) = ssha_e(:,:) 439 zun_e (:,:) = ua_e (:,:) 440 zvn_e (:,:) = va_e (:,:) 441 ELSE ! Swap + Filter 442 zsshb_e(:,:) = atfp * ( zsshb_e(:,:) + ssha_e(:,:) ) + atfp1 * sshn_e(:,:) 443 zub_e (:,:) = atfp * ( zub_e (:,:) + ua_e (:,:) ) + atfp1 * zun_e (:,:) 444 zvb_e (:,:) = atfp * ( zvb_e (:,:) + va_e (:,:) ) + atfp1 * zvn_e (:,:) 445 sshn_e (:,:) = ssha_e(:,:) 446 zun_e (:,:) = ua_e (:,:) 447 zvn_e (:,:) = va_e (:,:) 500 448 ENDIF 501 449 502 IF( lk_vvl ) THEN ! Variable volume ! needed for BDY ??? 503 504 ! Sea surface elevation 505 ! --------------------- 506 DO jj = 1, jpjm1 507 DO ji = 1,jpim1 508 509 ! Sea Surface Height at u-point before 510 zsshun_e(ji,jj) = 0.5 * umask(ji,jj,1) / ( e1u(ji,jj) * e2u(ji,jj) ) & 511 & * ( e1t(ji ,jj) * e2t(ji ,jj) * sshn_e(ji ,jj) & 512 & + e1t(ji+1,jj) * e2t(ji+1,jj) * sshn_e(ji+1,jj) ) 513 514 ! Sea Surface Height at v-point before 515 zsshvn_e(ji,jj) = 0.5 * vmask(ji,jj,1) / ( e1v(ji,jj) * e2v(ji,jj) ) & 516 & * ( e1t(ji,jj ) * e2t(ji,jj ) * sshn_e(ji,jj ) & 517 & + e1t(ji,jj+1) * e2t(ji,jj+1) * sshn_e(ji,jj+1) ) 518 519 END DO 520 END DO 521 522 ! Boundaries conditions 523 CALL lbc_lnk( zsshun_e, 'U', 1. ) ; CALL lbc_lnk( zsshvn_e, 'V', 1. ) 524 525 ! Ocean depth at U- and V-points 526 ! ------------------------------------------- 527 hu_e(:,:) = hu_0(:,:) + zsshun_e(:,:) 528 hv_e(:,:) = hv_0(:,:) + zsshvn_e(:,:) 529 450 IF( lk_vvl ) THEN !* Update ocean depth (variable volume case only) 451 ! ! ------------------ 452 DO jj = 1, jpjm1 ! Sea Surface Height at u- & v-points 453 DO ji = 1, fs_jpim1 ! Vector opt. 454 zsshun_e(ji,jj) = 0.5 * umask(ji,jj,1) / ( e1u(ji,jj) * e2u(ji,jj) ) & 455 & * ( e1t(ji ,jj) * e2t(ji ,jj) * sshn_e(ji ,jj) & 456 & + e1t(ji+1,jj) * e2t(ji+1,jj) * sshn_e(ji+1,jj) ) 457 zsshvn_e(ji,jj) = 0.5 * vmask(ji,jj,1) / ( e1v(ji,jj) * e2v(ji,jj) ) & 458 & * ( e1t(ji,jj ) * e2t(ji,jj ) * sshn_e(ji,jj ) & 459 & + e1t(ji,jj+1) * e2t(ji,jj+1) * sshn_e(ji,jj+1) ) 460 END DO 461 END DO 462 CALL lbc_lnk( zsshun_e, 'U', 1. ) ! lateral boundaries conditions 463 CALL lbc_lnk( zsshvn_e, 'V', 1. ) 464 ! 465 hu_e (:,:) = hu_0(:,:) + zsshun_e(:,:) ! Ocean depth at U- and V-points 466 hv_e (:,:) = hv_0(:,:) + zsshvn_e(:,:) 467 hur_e(:,:) = umask(:,:,1) / ( hu_e(:,:) + 1.e0 - umask(:,:,1) ) 468 hvr_e(:,:) = vmask(:,:,1) / ( hv_e(:,:) + 1.e0 - vmask(:,:,1) ) 530 469 ! 531 470 ENDIF … … 534 473 ! ! ==================== ! 535 474 536 ! Time average of after barotropic variables537 zcoef = 1.e0 / ( 2 * nn_baro + 1 )538 un_e (:,:) = zcoef * zu_sum (:,:)539 vn_e (:,:) = zcoef * zv_sum (:,:)540 sshn_e(:,:) = zcoef * zssh_sum(:,:)541 542 475 #if defined key_obc 543 IF( lp_obc_east ) sshfoe_b(:,:) = zcoef * sshfoe_b(:,:) 476 IF( lp_obc_east ) sshfoe_b(:,:) = zcoef * sshfoe_b(:,:) !!gm totally useless ????? 544 477 IF( lp_obc_west ) sshfow_b(:,:) = zcoef * sshfow_b(:,:) 545 478 IF( lp_obc_north ) sshfon_b(:,:) = zcoef * sshfon_b(:,:) … … 548 481 549 482 ! ----------------------------------------------------------------------------- 550 ! Phase 3. Coupling between general trend and barotropic estimates - (2nd step)483 ! Phase 3. update the general trend with the barotropic trend 551 484 ! ----------------------------------------------------------------------------- 552 553 554 555 ! add time averaged barotropic coriolis and surface pressure gradient 556 ! terms to the general momentum trend 557 ! -------------------------------------------------------------------- 485 ! 486 ! !* Time average ==> after barotropic u, v, ssh 487 zcoef = 1.e0 / ( 2 * nn_baro + 1 ) 488 un_b (:,:) = zcoef * zu_sum (:,:) 489 vn_b (:,:) = zcoef * zv_sum (:,:) 490 sshn_b(:,:) = zcoef * zssh_sum(:,:) 491 ! 492 ! !* update the general momentum trend 558 493 DO jk=1,jpkm1 559 ua(:,:,jk) = ua(:,:,jk) + hur(:,:) * ( un_e(:,:) - zub(:,:) ) / z2dt_b560 va(:,:,jk) = va(:,:,jk) + hvr(:,:) * ( vn_e(:,:) - zvb(:,:) ) / z2dt_b494 ua(:,:,jk) = ua(:,:,jk) + ( un_b(:,:) - zub(:,:) ) / z2dt_b 495 va(:,:,jk) = va(:,:,jk) + ( vn_b(:,:) - zvb(:,:) ) / z2dt_b 561 496 END DO 562 563 ! write filtered free surface arrays in restart file 564 ! -------------------------------------------------- 497 ! 498 ! !* write time-spliting arrays in the restart 565 499 IF( lrst_oce ) CALL ts_rst( kt, 'WRITE' ) 566 567 500 ! 568 501 END SUBROUTINE dyn_spg_ts … … 582 515 ! 583 516 IF( TRIM(cdrw) == 'READ' ) THEN 584 IF( iom_varid( numror, 'sshn_e', ldstop = .FALSE. ) > 0 ) THEN 585 CALL iom_get( numror, jpdom_autoglo, 'sshn_e', sshn_e(:,:) ) ! free surface and 586 CALL iom_get( numror, jpdom_autoglo, 'un_e' , un_e (:,:) ) ! horizontal transports issued 587 CALL iom_get( numror, jpdom_autoglo, 'vn_e' , vn_e (:,:) ) ! from barotropic loop 517 IF( iom_varid( numror, 'un_b', ldstop = .FALSE. ) > 0 ) THEN 518 CALL iom_get( numror, jpdom_autoglo, 'un_b' , un_b (:,:) ) ! external velocity issued 519 CALL iom_get( numror, jpdom_autoglo, 'vn_b' , vn_b (:,:) ) ! from barotropic loop 588 520 ELSE 589 sshn_e(:,:) = sshn(:,:) 590 un_e (:,:) = 0.e0 591 vn_e (:,:) = 0.e0 521 un_b (:,:) = 0.e0 522 vn_b (:,:) = 0.e0 592 523 ! vertical sum 593 524 IF( lk_vopt_loop ) THEN ! vector opt., forced unroll 594 525 DO jk = 1, jpkm1 595 526 DO ji = 1, jpij 596 un_ e(ji,1) = un_e(ji,1) + fse3u(ji,1,jk) * un(ji,1,jk)597 vn_ e(ji,1) = vn_e(ji,1) + fse3v(ji,1,jk) * vn(ji,1,jk)527 un_b(ji,1) = un_b(ji,1) + fse3u(ji,1,jk) * un(ji,1,jk) 528 vn_b(ji,1) = vn_b(ji,1) + fse3v(ji,1,jk) * vn(ji,1,jk) 598 529 END DO 599 530 END DO 600 531 ELSE ! No vector opt. 601 532 DO jk = 1, jpkm1 602 un_ e(:,:) = un_e(:,:) + fse3u(:,:,jk) * un(:,:,jk)603 vn_ e(:,:) = vn_e(:,:) + fse3v(:,:,jk) * vn(:,:,jk)533 un_b(:,:) = un_b(:,:) + fse3u(:,:,jk) * un(:,:,jk) 534 vn_b(:,:) = vn_b(:,:) + fse3v(:,:,jk) * vn(:,:,jk) 604 535 END DO 605 536 ENDIF 537 un_b (:,:) = un_b(:,:) * hur(:,:) 538 vn_b (:,:) = vn_b(:,:) * hvr(:,:) 606 539 ENDIF 607 540 ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN 608 CALL iom_rstput( kt, nitrst, numrow, 'sshn_e', sshn_e(:,:) ) ! free surface and 609 CALL iom_rstput( kt, nitrst, numrow, 'un_e' , un_e (:,:) ) ! horizontal transports issued 610 CALL iom_rstput( kt, nitrst, numrow, 'vn_e' , vn_e (:,:) ) ! from barotropic loop 541 CALL iom_rstput( kt, nitrst, numrow, 'un_b' , un_b (:,:) ) ! external velocity issued 542 CALL iom_rstput( kt, nitrst, numrow, 'vn_b' , vn_b (:,:) ) ! from barotropic loop 611 543 ENDIF 612 544 !
Note: See TracChangeset
for help on using the changeset viewer.