Changeset 1438 for trunk/NEMO/OPA_SRC/DYN/dynspg_exp.F90
- Timestamp:
- 2009-05-11T16:34:47+02:00 (15 years ago)
- File:
-
- 1 edited
-
trunk/NEMO/OPA_SRC/DYN/dynspg_exp.F90 (modified) (7 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMO/OPA_SRC/DYN/dynspg_exp.F90
r1146 r1438 4 4 !! Ocean dynamics: surface pressure gradient trend 5 5 !!====================================================================== 6 !! History : 9.0 ! 2005-11 (V. Garnier, G. Madec, L. Bessieres) Original code 7 !!---------------------------------------------------------------------- 6 8 #if defined key_dynspg_exp || defined key_esopa 7 9 !!---------------------------------------------------------------------- … … 11 13 !! pressure gradient in the free surface constant 12 14 !! volume case with vector optimization 13 !! exp_rst : read/write the explicit restart fields in the ocean restart file14 15 !!---------------------------------------------------------------------- 15 !! * Modules used16 16 USE oce ! ocean dynamics and tracers 17 17 USE dom_oce ! ocean space and time domain … … 31 31 PRIVATE 32 32 33 !! * Accessibility 34 PUBLIC dyn_spg_exp ! routine called by step.F90 35 PUBLIC exp_rst ! routine called by istate.F90 33 PUBLIC dyn_spg_exp ! routine called by step.F90 36 34 37 35 !! * Substitutions … … 39 37 # include "vectopt_loop_substitute.h90" 40 38 !!---------------------------------------------------------------------- 41 !! OPA 9.0 , LOCEAN-IPSL (2005)39 !! NEMO/OPA 3.2 , LOCEAN-IPSL (2009) 42 40 !! $Id$ 43 !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt41 !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 44 42 !!---------------------------------------------------------------------- 45 46 43 47 44 CONTAINS … … 62 59 !! -1- Compute the now surface pressure gradient 63 60 !! -2- Add it to the general trend 64 !! -3- Compute the horizontal divergence of velocities65 !! - the now divergence is given by :66 !! zhdivn = 1/(e1t*e2t*e3t) ( di[e2u*e3u un] + dj[e1v*e3v vn] )67 !! - integrate the horizontal divergence from the bottom68 !! to the surface69 !! - apply lateral boundary conditions on zhdivn70 !! -4- Estimate the after sea surface elevation from the kinematic71 !! surface boundary condition:72 !! zssha = sshb - 2 rdt ( zhdiv + emp )73 !! - Time filter applied on now sea surface elevation to avoid74 !! the divergence of two consecutive time-steps and swap of free75 !! surface arrays to start the next time step:76 !! sshb = sshn + atfp * [ sshb + zssha - 2 sshn ]77 !! sshn = zssha78 !! - apply lateral boundary conditions on sshn79 61 !! 80 62 !! ** Action : - Update (ua,va) with the surf. pressure gradient trend 63 !!--------------------------------------------------------------------- 64 INTEGER, INTENT( in ) :: kt ! ocean time-step index 81 65 !! 82 !! References :83 !!84 !! History :85 !! 9.0 ! 05-11 (V. Garnier, G. Madec, L. Bessieres) Original code86 !!87 !!---------------------------------------------------------------------88 !! * Arguments89 INTEGER, INTENT( in ) :: kt ! ocean time-step index90 91 !! * Local declarations92 66 INTEGER :: ji, jj, jk ! dummy loop indices 93 REAL(wp) :: z2dt, zraur, zssha ! temporary scalars94 REAL(wp), DIMENSION(jpi,jpj) :: & ! temporary arrays95 & zhdiv96 67 !!---------------------------------------------------------------------- 97 68 … … 104 75 spgu(:,:) = 0.e0 ! surface pressure gradient (i-direction) 105 76 spgv(:,:) = 0.e0 ! surface pressure gradient (j-direction) 106 107 CALL exp_rst( nit000, 'READ' ) ! read or initialize the following fields:108 ! ! sshb, sshn109 110 77 ENDIF 111 78 112 IF( .NOT. lk_vvl ) THEN 79 ! read or estimate sea surface height and vertically integrated velocities 80 IF( lk_obc ) CALL obc_dta_bt( kt, 0 ) 113 81 114 ! 0. Initialization115 ! -----------------116 ! read or estimate sea surface height and vertically integrated velocities117 IF( lk_obc ) CALL obc_dta_bt( kt, 0)118 z2dt = 2. * rdt ! time step: leap-frog119 IF( neuler == 0 .AND. kt == nit000 ) z2dt = rdt ! time step: Euler if restart from rest120 zraur = 1.e0 / rauw82 ! Surface pressure gradient (now) 83 DO jj = 2, jpjm1 84 DO ji = fs_2, fs_jpim1 ! vector opt. 85 spgu(ji,jj) = - grav * ( sshn(ji+1,jj) - sshn(ji,jj) ) / e1u(ji,jj) 86 spgv(ji,jj) = - grav * ( sshn(ji,jj+1) - sshn(ji,jj) ) / e2v(ji,jj) 87 END DO 88 END DO 121 89 122 ! 1. Surface pressure gradient (now)123 ! ----------------------------90 ! Add the surface pressure trend to the general trend 91 DO jk = 1, jpkm1 124 92 DO jj = 2, jpjm1 125 93 DO ji = fs_2, fs_jpim1 ! vector opt. 126 spgu(ji,jj) = - grav * ( sshn(ji+1,jj) - sshn(ji,jj) ) / e1u(ji,jj)127 spgv(ji,jj) = - grav * ( sshn(ji,jj+1) - sshn(ji,jj) ) / e2v(ji,jj)94 ua(ji,jj,jk) = ua(ji,jj,jk) + spgu(ji,jj) 95 va(ji,jj,jk) = va(ji,jj,jk) + spgv(ji,jj) 128 96 END DO 129 97 END DO 98 END DO 130 99 131 ! 2. Add the surface pressure trend to the general trend132 ! ------------------------------------------------------133 DO jk = 1, jpkm1134 DO jj = 2, jpjm1135 DO ji = fs_2, fs_jpim1 ! vector opt.136 ua(ji,jj,jk) = ua(ji,jj,jk) + spgu(ji,jj)137 va(ji,jj,jk) = va(ji,jj,jk) + spgv(ji,jj)138 END DO139 END DO140 END DO141 142 ! 3. Vertical integration of horizontal divergence of velocities143 ! --------------------------------144 zhdiv(:,:) = 0.e0145 DO jk = jpkm1, 1, -1146 DO jj = 2, jpjm1147 DO ji = fs_2, fs_jpim1 ! vector opt.148 zhdiv(ji,jj) = zhdiv(ji,jj) + ( e2u(ji ,jj ) * fse3u(ji ,jj ,jk) * un(ji ,jj ,jk) &149 & - e2u(ji-1,jj ) * fse3u(ji-1,jj ,jk) * un(ji-1,jj ,jk) &150 & + e1v(ji ,jj ) * fse3v(ji ,jj ,jk) * vn(ji ,jj ,jk) &151 & - e1v(ji ,jj-1) * fse3v(ji ,jj-1,jk) * vn(ji ,jj-1,jk) ) &152 & / ( e1t(ji,jj) * e2t(ji,jj) )153 END DO154 END DO155 END DO156 157 #if defined key_obc158 ! open boundaries (div must be zero behind the open boundary)159 ! mpp remark: The zeroing of zhdiv can probably be extended to 1->jpi/jpj for the correct row/column160 IF( lp_obc_east ) zhdiv(nie0p1:nie1p1,nje0 :nje1 ) = 0.e0 ! east161 IF( lp_obc_west ) zhdiv(niw0 :niw1 ,njw0 :njw1 ) = 0.e0 ! west162 IF( lp_obc_north ) zhdiv(nin0 :nin1 ,njn0p1:njn1p1) = 0.e0 ! north163 IF( lp_obc_south ) zhdiv(nis0 :nis1 ,njs0 :njs1 ) = 0.e0 ! south164 #endif165 166 ! 4. Sea surface elevation time stepping167 ! --------------------------------------168 ! Euler (forward) time stepping, no time filter169 IF( neuler == 0 .AND. kt == nit000 ) THEN170 DO jj = 1, jpj171 DO ji = 1, jpi172 ! after free surface elevation173 zssha = sshb(ji,jj) - rdt * ( zraur * emp(ji,jj) + zhdiv(ji,jj) ) * tmask(ji,jj,1)174 ! swap of arrays175 sshb(ji,jj) = sshn(ji,jj)176 sshn(ji,jj) = zssha177 END DO178 END DO179 ELSE180 ! Leap-frog time stepping and time filter181 DO jj = 1, jpj182 DO ji = 1, jpi183 ! after free surface elevation184 zssha = sshb(ji,jj) - z2dt * ( zraur * emp(ji,jj) + zhdiv(ji,jj) ) * tmask(ji,jj,1)185 ! time filter and array swap186 sshb(ji,jj) = atfp * ( sshb(ji,jj) + zssha ) + atfp1 * sshn(ji,jj)187 sshn(ji,jj) = zssha188 END DO189 END DO190 ENDIF191 192 ELSE !! Variable volume, ssh time-stepping already done193 194 ! read or estimate sea surface height and vertically integrated velocities195 IF( lk_obc ) CALL obc_dta_bt( kt, 0 )196 197 ! Sea surface elevation swap198 ! -----------------------------199 !200 sshbb(:,:) = sshb(:,:) ! Needed for the dynamics time-stepping201 !202 IF( neuler == 0 .AND. kt == nit000 ) THEN203 ! No time filter204 sshb(:,:) = sshn(:,:)205 sshn(:,:) = ssha(:,:)206 ELSE207 ! Time filter208 sshb(:,:) = atfp * ( sshb(:,:) + ssha(:,:) ) + atfp1 * sshn(:,:)209 sshn(:,:) = ssha(:,:)210 ENDIF211 212 ENDIF213 214 ! Boundary conditions on sshn215 IF( .NOT. lk_obc ) CALL lbc_lnk( sshn, 'T', 1. )216 217 ! write filtered free surface arrays in restart file218 ! --------------------------------------------------219 IF( lrst_oce ) CALL exp_rst( kt, 'WRITE' )220 221 IF(ln_ctl) THEN ! print sum trends (used for debugging)222 CALL prt_ctl(tab2d_1=sshn, clinfo1=' ssh : ', mask1=tmask)223 ENDIF224 225 100 END SUBROUTINE dyn_spg_exp 226 101 227 SUBROUTINE exp_rst( kt, cdrw )228 !!---------------------------------------------------------------------229 !! *** ROUTINE exp_rst ***230 !!231 !! ** Purpose : Read or write explicit arrays in restart file232 !!----------------------------------------------------------------------233 INTEGER , INTENT(in) :: kt ! ocean time-step234 CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag235 !236 !!----------------------------------------------------------------------237 !238 IF( TRIM(cdrw) == 'READ' ) THEN239 IF( iom_varid( numror, 'sshn', ldstop = .FALSE. ) > 0 ) THEN240 CALL iom_get( numror, jpdom_autoglo, 'sshb' , sshb(:,:) )241 CALL iom_get( numror, jpdom_autoglo, 'sshn' , sshn(:,:) )242 IF( neuler == 0 ) sshb(:,:) = sshn(:,:)243 ELSE244 IF( nn_rstssh == 1 ) THEN245 sshb(:,:) = 0.e0246 sshn(:,:) = 0.e0247 ENDIF248 ENDIF249 ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN250 CALL iom_rstput( kt, nitrst, numrow, 'sshb' , sshb (:,:) )251 CALL iom_rstput( kt, nitrst, numrow, 'sshn' , sshn (:,:) )252 ENDIF253 !254 END SUBROUTINE exp_rst255 102 #else 256 103 !!---------------------------------------------------------------------- … … 261 108 WRITE(*,*) 'dyn_spg_exp: You should not have seen this print! error?', kt 262 109 END SUBROUTINE dyn_spg_exp 263 SUBROUTINE exp_rst( kt, cdrw ) ! Empty routine264 INTEGER , INTENT(in) :: kt ! ocean time-step265 CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag266 WRITE(*,*) 'exp_rst: You should not have seen this print! error?', kt, cdrw267 END SUBROUTINE exp_rst268 110 #endif 269 111 270 112 !!====================================================================== 271 113 END MODULE dynspg_exp
Note: See TracChangeset
for help on using the changeset viewer.
