- 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/dynzad.F90
r7753 r9019 5 5 !!====================================================================== 6 6 !! History : OPA ! 1991-01 (G. Madec) Original code 7 !! 7.0 ! 1991-11 (G. Madec)8 !! 7.5 ! 1996-01 (G. Madec) statement function for e39 7 !! NEMO 0.5 ! 2002-07 (G. Madec) Free form, F90 10 8 !!---------------------------------------------------------------------- … … 22 20 USE lib_mpp ! MPP library 23 21 USE prtctl ! Print control 24 USE wrk_nemo ! Memory Allocation25 22 USE timing ! Timing 26 23 … … 29 26 30 27 PUBLIC dyn_zad ! routine called by dynadv.F90 31 PUBLIC dyn_zad_zts ! routine called by dynadv.F9032 28 33 29 !! * Substitutions 34 30 # include "vectopt_loop_substitute.h90" 35 31 !!---------------------------------------------------------------------- 36 !! NEMO/OPA 3.3 , NEMO Consortium (2010)32 !! NEMO/OPA 4.0 , NEMO Consortium (2017) 37 33 !! $Id$ 38 34 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) … … 58 54 INTEGER, INTENT(in) :: kt ! ocean time-step inedx 59 55 ! 60 INTEGER :: ji, jj, jk 61 REAL(wp) :: zua, zva ! temporaryscalars62 REAL(wp), POINTER, DIMENSION(:,:,:) :: zwuw , zwvw63 REAL(wp), POINTER, DIMENSION(:,: ) :: zww64 REAL(wp), POINTER, DIMENSION(:,:,:) ::ztrdu, ztrdv56 INTEGER :: ji, jj, jk ! dummy loop indices 57 REAL(wp) :: zua, zva ! local scalars 58 REAL(wp), DIMENSION(jpi,jpj) :: zww 59 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zwuw, zwvw 60 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ztrdu, ztrdv 65 61 !!---------------------------------------------------------------------- 66 62 ! 67 IF( nn_timing == 1 ) CALL timing_start('dyn_zad') 68 ! 69 CALL wrk_alloc( jpi,jpj, zww ) 70 CALL wrk_alloc( jpi,jpj,jpk, zwuw , zwvw ) 63 IF( ln_timing ) CALL timing_start('dyn_zad') 71 64 ! 72 65 IF( kt == nit000 ) THEN 73 IF(lwp) WRITE(numout,*)74 IF(lwp) WRITE(numout,*) 'dyn_zad : arakawaadvection scheme'66 IF(lwp) WRITE(numout,*) 67 IF(lwp) WRITE(numout,*) 'dyn_zad : 2nd order vertical advection scheme' 75 68 ENDIF 76 69 77 70 IF( l_trddyn ) THEN ! Save ua and va trends 78 CALL wrk_alloc( jpi, jpj, jpk, ztrdu, ztrdv)71 ALLOCATE( ztrdu(jpi,jpj,jpk) , ztrdv(jpi,jpj,jpk) ) 79 72 ztrdu(:,:,:) = ua(:,:,:) 80 73 ztrdv(:,:,:) = va(:,:,:) … … 96 89 ! 97 90 ! Surface and bottom advective fluxes set to zero 98 IF 91 IF( ln_isfcav ) THEN 99 92 DO jj = 2, jpjm1 100 93 DO ji = fs_2, fs_jpim1 ! vector opt. … … 119 112 DO jj = 2, jpjm1 120 113 DO ji = fs_2, fs_jpim1 ! vector opt. 121 ! ! vertical momentum advective trends 122 zua = - ( zwuw(ji,jj,jk) + zwuw(ji,jj,jk+1) ) * r1_e1e2u(ji,jj) / e3u_n(ji,jj,jk) 123 zva = - ( zwvw(ji,jj,jk) + zwvw(ji,jj,jk+1) ) * r1_e1e2v(ji,jj) / e3v_n(ji,jj,jk) 124 ! ! add the trends to the general momentum trends 125 ua(ji,jj,jk) = ua(ji,jj,jk) + zua 126 va(ji,jj,jk) = va(ji,jj,jk) + zva 114 ua(ji,jj,jk) = ua(ji,jj,jk) - ( zwuw(ji,jj,jk) + zwuw(ji,jj,jk+1) ) * r1_e1e2u(ji,jj) / e3u_n(ji,jj,jk) 115 va(ji,jj,jk) = va(ji,jj,jk) - ( zwvw(ji,jj,jk) + zwvw(ji,jj,jk+1) ) * r1_e1e2v(ji,jj) / e3v_n(ji,jj,jk) 127 116 END DO 128 117 END DO … … 133 122 ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 134 123 CALL trd_dyn( ztrdu, ztrdv, jpdyn_zad, kt ) 135 CALL wrk_dealloc( jpi, jpj, jpk,ztrdu, ztrdv )124 DEALLOCATE( ztrdu, ztrdv ) 136 125 ENDIF 137 126 ! ! Control print … … 139 128 & tab3d_2=va, clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' ) 140 129 ! 141 CALL wrk_dealloc( jpi,jpj, zww ) 142 CALL wrk_dealloc( jpi,jpj,jpk, zwuw , zwvw ) 143 ! 144 IF( nn_timing == 1 ) CALL timing_stop('dyn_zad') 130 IF( ln_timing ) CALL timing_stop('dyn_zad') 145 131 ! 146 132 END SUBROUTINE dyn_zad 147 133 148 149 SUBROUTINE dyn_zad_zts ( kt )150 !!----------------------------------------------------------------------151 !! *** ROUTINE dynzad_zts ***152 !!153 !! ** Purpose : Compute the now vertical momentum advection trend and154 !! add it to the general trend of momentum equation. This version155 !! uses sub-timesteps for improved numerical stability with small156 !! vertical grid sizes. This is especially relevant when using157 !! embedded ice with thin surface boxes.158 !!159 !! ** Method : The now vertical advection of momentum is given by:160 !! w dz(u) = ua + 1/(e1u*e2u*e3u) mk+1[ mi(e1t*e2t*wn) dk(un) ]161 !! w dz(v) = va + 1/(e1v*e2v*e3v) mk+1[ mj(e1t*e2t*wn) dk(vn) ]162 !! Add this trend to the general trend (ua,va):163 !! (ua,va) = (ua,va) + w dz(u,v)164 !!165 !! ** Action : - Update (ua,va) with the vert. momentum adv. trends166 !! - Save the trends in (ztrdu,ztrdv) ('key_trddyn')167 !!----------------------------------------------------------------------168 INTEGER, INTENT(in) :: kt ! ocean time-step inedx169 !170 INTEGER :: ji, jj, jk, jl ! dummy loop indices171 INTEGER :: jnzts = 5 ! number of sub-timesteps for vertical advection172 INTEGER :: jtb, jtn, jta ! sub timestep pointers for leap-frog/euler forward steps173 REAL(wp) :: zua, zva ! temporary scalars174 REAL(wp) :: zr_rdt ! temporary scalar175 REAL(wp) :: z2dtzts ! length of Euler forward sub-timestep for vertical advection176 REAL(wp) :: zts ! length of sub-timestep for vertical advection177 REAL(wp), POINTER, DIMENSION(:,:,:) :: zwuw , zwvw, zww178 REAL(wp), POINTER, DIMENSION(:,:,:) :: ztrdu, ztrdv179 REAL(wp), POINTER, DIMENSION(:,:,:,:) :: zus , zvs180 !!----------------------------------------------------------------------181 !182 IF( nn_timing == 1 ) CALL timing_start('dyn_zad_zts')183 !184 CALL wrk_alloc( jpi,jpj,jpk, zwuw, zwvw, zww )185 CALL wrk_alloc( jpi,jpj,jpk,3, zus , zvs )186 !187 IF( kt == nit000 ) THEN188 IF(lwp)WRITE(numout,*)189 IF(lwp)WRITE(numout,*) 'dyn_zad_zts : arakawa advection scheme with sub-timesteps'190 ENDIF191 192 IF( l_trddyn ) THEN ! Save ua and va trends193 CALL wrk_alloc( jpi, jpj, jpk, ztrdu, ztrdv )194 ztrdu(:,:,:) = ua(:,:,:)195 ztrdv(:,:,:) = va(:,:,:)196 ENDIF197 198 IF( neuler == 0 .AND. kt == nit000 ) THEN199 z2dtzts = rdt / REAL( jnzts, wp ) ! = rdt (restart with Euler time stepping)200 ELSE201 z2dtzts = 2._wp * rdt / REAL( jnzts, wp ) ! = 2 rdt (leapfrog)202 ENDIF203 204 DO jk = 2, jpkm1 ! Calculate and store vertical fluxes205 DO jj = 2, jpj206 DO ji = fs_2, jpi ! vector opt.207 zww(ji,jj,jk) = 0.25_wp * e1e2t(ji,jj) * wn(ji,jj,jk)208 END DO209 END DO210 END DO211 212 DO jj = 2, jpjm1 ! Surface and bottom advective fluxes set to zero213 DO ji = fs_2, fs_jpim1 ! vector opt.214 !!gm missing ISF boundary condition215 zwuw(ji,jj, 1 ) = 0._wp216 zwvw(ji,jj, 1 ) = 0._wp217 zwuw(ji,jj,jpk) = 0._wp218 zwvw(ji,jj,jpk) = 0._wp219 END DO220 END DO221 222 ! Start with before values and use sub timestepping to reach after values223 224 zus(:,:,:,1) = ub(:,:,:)225 zvs(:,:,:,1) = vb(:,:,:)226 227 DO jl = 1, jnzts ! Start of sub timestepping loop228 229 IF( jl == 1 ) THEN ! Euler forward to kick things off230 jtb = 1 ; jtn = 1 ; jta = 2231 zts = z2dtzts232 ELSEIF( jl == 2 ) THEN ! First leapfrog step233 jtb = 1 ; jtn = 2 ; jta = 3234 zts = 2._wp * z2dtzts235 ELSE ! Shuffle pointers for subsequent leapfrog steps236 jtb = MOD(jtb,3) + 1237 jtn = MOD(jtn,3) + 1238 jta = MOD(jta,3) + 1239 ENDIF240 241 DO jk = 2, jpkm1 ! Vertical momentum advection at level w and u- and v- vertical242 DO jj = 2, jpjm1 ! vertical momentum advection at w-point243 DO ji = fs_2, fs_jpim1 ! vector opt.244 zwuw(ji,jj,jk) = ( zww(ji+1,jj ,jk) + zww(ji,jj,jk) ) * ( zus(ji,jj,jk-1,jtn)-zus(ji,jj,jk,jtn) ) !* wumask(ji,jj,jk)245 zwvw(ji,jj,jk) = ( zww(ji ,jj+1,jk) + zww(ji,jj,jk) ) * ( zvs(ji,jj,jk-1,jtn)-zvs(ji,jj,jk,jtn) ) !* wvmask(ji,jj,jk)246 END DO247 END DO248 END DO249 DO jk = 1, jpkm1 ! Vertical momentum advection at u- and v-points250 DO jj = 2, jpjm1251 DO ji = fs_2, fs_jpim1 ! vector opt.252 ! ! vertical momentum advective trends253 zua = - ( zwuw(ji,jj,jk) + zwuw(ji,jj,jk+1) ) * r1_e1e2u(ji,jj) / e3u_n(ji,jj,jk)254 zva = - ( zwvw(ji,jj,jk) + zwvw(ji,jj,jk+1) ) * r1_e1e2v(ji,jj) / e3v_n(ji,jj,jk)255 zus(ji,jj,jk,jta) = zus(ji,jj,jk,jtb) + zua * zts256 zvs(ji,jj,jk,jta) = zvs(ji,jj,jk,jtb) + zva * zts257 END DO258 END DO259 END DO260 261 END DO ! End of sub timestepping loop262 263 zr_rdt = 1._wp / ( REAL( jnzts, wp ) * z2dtzts )264 DO jk = 1, jpkm1 ! Recover trends over the outer timestep265 DO jj = 2, jpjm1266 DO ji = fs_2, fs_jpim1 ! vector opt.267 ! ! vertical momentum advective trends268 ! ! add the trends to the general momentum trends269 ua(ji,jj,jk) = ua(ji,jj,jk) + ( zus(ji,jj,jk,jta) - ub(ji,jj,jk)) * zr_rdt270 va(ji,jj,jk) = va(ji,jj,jk) + ( zvs(ji,jj,jk,jta) - vb(ji,jj,jk)) * zr_rdt271 END DO272 END DO273 END DO274 275 IF( l_trddyn ) THEN ! save the vertical advection trends for diagnostic276 ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)277 ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)278 CALL trd_dyn( ztrdu, ztrdv, jpdyn_zad, kt )279 CALL wrk_dealloc( jpi, jpj, jpk, ztrdu, ztrdv )280 ENDIF281 ! ! Control print282 IF(ln_ctl) CALL prt_ctl( tab3d_1=ua, clinfo1=' zad - Ua: ', mask1=umask, &283 & tab3d_2=va, clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' )284 !285 CALL wrk_dealloc( jpi,jpj,jpk, zwuw, zwvw, zww )286 CALL wrk_dealloc( jpi,jpj,jpk,3, zus , zvs )287 !288 IF( nn_timing == 1 ) CALL timing_stop('dyn_zad_zts')289 !290 END SUBROUTINE dyn_zad_zts291 292 134 !!====================================================================== 293 135 END MODULE dynzad
Note: See TracChangeset
for help on using the changeset viewer.
