Changeset 5902 for branches/2015/dev_r5847_MERCATOR9_solveur_simplification/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90
- Timestamp:
- 2015-11-20T10:58:48+01:00 (8 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2015/dev_r5847_MERCATOR9_solveur_simplification/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90
r5868 r5902 19 19 !! 3.5 ! 2013-07 (J. Chanut) Compliant with time splitting changes 20 20 !! 3.7 ! 2014-04 (G. Madec) add the diagnostic of the time filter trends 21 !! 3.7 ! 2015-11 (J. Chanut) Free surface simplification 21 22 !!------------------------------------------------------------------------- 22 23 … … 28 29 USE sbc_oce ! Surface boundary condition: ocean fields 29 30 USE phycst ! physical constants 30 USE dynspg_oce ! type of surface pressure gradient31 31 USE dynadv ! dynamics: vector invariant versus flux form 32 32 USE domvvl ! variable volume … … 68 68 !! *** ROUTINE dyn_nxt *** 69 69 !! 70 !! ** Purpose : Compute the after horizontal velocity. Apply the boundary71 !! condition on the after velocity, achieve dthe time stepping70 !! ** Purpose : Finalize after horizontal velocity. Apply the boundary 71 !! condition on the after velocity, achieve the time stepping 72 72 !! by applying the Asselin filter on now fields and swapping 73 73 !! the fields. 74 74 !! 75 !! ** Method : * After velocity is compute using a leap-frog scheme: 76 !! (ua,va) = (ub,vb) + 2 rdt (ua,va) 77 !! Note that with flux form advection and variable volume layer 78 !! (lk_vvl=T), the leap-frog is applied on thickness weighted 79 !! velocity. 75 !! ** Method : * Ensure after velocities transport matches time splitting 76 !! estimate (ln_dynspg_ts=T) 80 77 !! 81 78 !! * Apply lateral boundary conditions on after velocity … … 99 96 INTEGER :: ji, jj, jk ! dummy loop indices 100 97 INTEGER :: iku, ikv ! local integers 101 REAL(wp) :: z2dt ! temporary scalar 102 REAL(wp) :: zue3a, zue3n, zue3b, zuf, zec ! local scalars 98 REAL(wp) :: zue3a, zue3n, zue3b, zuf ! local scalars 103 99 REAL(wp) :: zve3a, zve3n, zve3b, zvf, z1_2dt ! - - 104 100 REAL(wp), POINTER, DIMENSION(:,:) :: zue, zve … … 108 104 IF( nn_timing == 1 ) CALL timing_start('dyn_nxt') 109 105 ! 110 CALL wrk_alloc( jpi,jpj,jpk, ze3u_f, ze3v_f, zua, zva ) 111 IF( lk_dynspg_ts ) CALL wrk_alloc( jpi,jpj, zue, zve ) 106 IF( ln_dynspg_ts ) CALL wrk_alloc( jpi,jpj, zue, zve ) 107 IF( lk_vvl.AND.(.NOT.ln_dynadv_vec ) ) CALL wrk_alloc( jpi,jpj,jpk, ze3u_f, ze3v_f ) 108 IF( l_trddyn ) CALL wrk_alloc( jpi,jpj,jpk, zua, zva ) 112 109 ! 113 110 IF( kt == nit000 ) THEN … … 117 114 ENDIF 118 115 119 # if defined key_dynspg_exp 120 ! Next velocity : Leap-frog time stepping 121 ! ------------- 122 z2dt = 2. * rdt ! Euler or leap-frog time step 123 IF( neuler == 0 .AND. kt == nit000 ) z2dt = rdt 124 ! 125 IF( ln_dynadv_vec .OR. .NOT. lk_vvl ) THEN ! applied on velocity 116 IF ( ln_dynspg_ts ) THEN 117 ! Ensure below that barotropic velocities match time splitting estimate 118 ! Compute actual transport and replace it with ts estimate at "after" time step 119 zue(:,:) = fse3u_a(:,:,1) * ua(:,:,1) * umask(:,:,1) 120 zve(:,:) = fse3v_a(:,:,1) * va(:,:,1) * vmask(:,:,1) 121 DO jk = 2, jpkm1 122 zue(:,:) = zue(:,:) + fse3u_a(:,:,jk) * ua(:,:,jk) * umask(:,:,jk) 123 zve(:,:) = zve(:,:) + fse3v_a(:,:,jk) * va(:,:,jk) * vmask(:,:,jk) 124 END DO 126 125 DO jk = 1, jpkm1 127 ua(:,:,jk) = ( ub(:,:,jk) + z2dt * ua(:,:,jk) ) * umask(:,:,jk) 128 va(:,:,jk) = ( vb(:,:,jk) + z2dt * va(:,:,jk) ) * vmask(:,:,jk) 129 END DO 130 ELSE ! applied on thickness weighted velocity 131 DO jk = 1, jpkm1 132 ua(:,:,jk) = ( ub(:,:,jk) * fse3u_b(:,:,jk) & 133 & + z2dt * ua(:,:,jk) * fse3u_n(:,:,jk) ) & 134 & / fse3u_a(:,:,jk) * umask(:,:,jk) 135 va(:,:,jk) = ( vb(:,:,jk) * fse3v_b(:,:,jk) & 136 & + z2dt * va(:,:,jk) * fse3v_n(:,:,jk) ) & 137 & / fse3v_a(:,:,jk) * vmask(:,:,jk) 138 END DO 139 ENDIF 140 # endif 141 142 # if defined key_dynspg_ts 143 !!gm IF ( lk_dynspg_ts ) THEN .... 144 ! Ensure below that barotropic velocities match time splitting estimate 145 ! Compute actual transport and replace it with ts estimate at "after" time step 146 zue(:,:) = fse3u_a(:,:,1) * ua(:,:,1) * umask(:,:,1) 147 zve(:,:) = fse3v_a(:,:,1) * va(:,:,1) * vmask(:,:,1) 148 DO jk = 2, jpkm1 149 zue(:,:) = zue(:,:) + fse3u_a(:,:,jk) * ua(:,:,jk) * umask(:,:,jk) 150 zve(:,:) = zve(:,:) + fse3v_a(:,:,jk) * va(:,:,jk) * vmask(:,:,jk) 151 END DO 152 DO jk = 1, jpkm1 153 ua(:,:,jk) = ( ua(:,:,jk) - zue(:,:) * hur_a(:,:) + ua_b(:,:) ) * umask(:,:,jk) 154 va(:,:,jk) = ( va(:,:,jk) - zve(:,:) * hvr_a(:,:) + va_b(:,:) ) * vmask(:,:,jk) 155 END DO 156 157 IF (lk_dynspg_ts.AND.(.NOT.ln_bt_fw)) THEN 158 ! Remove advective velocity from "now velocities" 159 ! prior to asselin filtering 160 ! In the forward case, this is done below after asselin filtering 161 ! so that asselin contribution is removed at the same time 162 DO jk = 1, jpkm1 163 un(:,:,jk) = ( un(:,:,jk) - un_adv(:,:) + un_b(:,:) )*umask(:,:,jk) 164 vn(:,:,jk) = ( vn(:,:,jk) - vn_adv(:,:) + vn_b(:,:) )*vmask(:,:,jk) 165 END DO 166 ENDIF 167 !!gm ENDIF 168 # endif 126 ua(:,:,jk) = ( ua(:,:,jk) - zue(:,:) * hur_a(:,:) + ua_b(:,:) ) * umask(:,:,jk) 127 va(:,:,jk) = ( va(:,:,jk) - zve(:,:) * hvr_a(:,:) + va_b(:,:) ) * vmask(:,:,jk) 128 END DO 129 130 IF ( .NOT.ln_bt_fw ) THEN 131 ! Remove advective velocity from "now velocities" 132 ! prior to asselin filtering 133 ! In the forward case, this is done below after asselin filtering 134 ! so that asselin contribution is removed at the same time 135 DO jk = 1, jpkm1 136 un(:,:,jk) = ( un(:,:,jk) - un_adv(:,:) + un_b(:,:) )*umask(:,:,jk) 137 vn(:,:,jk) = ( vn(:,:,jk) - vn_adv(:,:) + vn_b(:,:) )*vmask(:,:,jk) 138 END DO 139 ENDIF 140 ENDIF 169 141 170 142 ! Update after velocity on domain lateral boundaries 171 143 ! -------------------------------------------------- 172 CALL lbc_lnk( ua, 'U', -1. ) !* local domain boundaries173 CALL lbc_lnk( va, 'V', -1. )174 !175 # if defined key_bdy176 ! !* BDY open boundaries177 IF( lk_bdy .AND. lk_dynspg_exp ) CALL bdy_dyn( kt )178 IF( lk_bdy .AND. lk_dynspg_ts ) CALL bdy_dyn( kt, dyn3d_only=.true. )179 180 !!$ Do we need a call to bdy_vol here??181 !182 # endif183 !184 144 # if defined key_agrif 185 145 CALL Agrif_dyn( kt ) !* AGRIF zoom boundaries 186 146 # endif 187 147 ! 148 CALL lbc_lnk( ua, 'U', -1. ) !* local domain boundaries 149 CALL lbc_lnk( va, 'V', -1. ) 150 ! 151 # if defined key_bdy 152 ! !* BDY open boundaries 153 IF( lk_bdy .AND. ln_dynspg_exp ) CALL bdy_dyn( kt ) 154 IF( lk_bdy .AND. ln_dynspg_ts ) CALL bdy_dyn( kt, dyn3d_only=.true. ) 155 156 !!$ Do we need a call to bdy_vol here?? 157 ! 158 # endif 159 ! 188 160 IF( l_trddyn ) THEN ! prepare the atf trend computation + some diagnostics 189 161 z1_2dt = 1._wp / (2. * rdt) ! Euler or leap-frog time step … … 242 214 ! (used as a now filtered scale factor until the swap) 243 215 ! ---------------------------------------------------- 244 IF (l k_dynspg_ts.AND.ln_bt_fw) THEN216 IF (ln_dynspg_ts.AND.ln_bt_fw) THEN 245 217 ! No asselin filtering on thicknesses if forward time splitting 246 218 fse3t_b(:,:,:) = fse3t_n(:,:,:) 247 219 ELSE 248 220 fse3t_b(:,:,:) = fse3t_n(:,:,:) + atfp * ( fse3t_b(:,:,:) - 2._wp * fse3t_n(:,:,:) + fse3t_a(:,:,:) ) … … 319 291 ENDIF 320 292 ! 321 IF (l k_dynspg_ts.AND.ln_bt_fw) THEN293 IF (ln_dynspg_ts.AND.ln_bt_fw) THEN 322 294 ! Revert "before" velocities to time split estimate 323 295 ! Doing it here also means that asselin filter contribution is removed … … 383 355 IF(ln_ctl) CALL prt_ctl( tab3d_1=un, clinfo1=' nxt - Un: ', mask1=umask, & 384 356 & tab3d_2=vn, clinfo2=' Vn: ' , mask2=vmask ) 385 ! 386 CALL wrk_dealloc( jpi,jpj,jpk, ze3u_f, ze3v_f, zua, zva ) 387 IF( lk_dynspg_ts ) CALL wrk_dealloc( jpi,jpj, zue, zve ) 357 ! 358 IF( ln_dynspg_ts ) CALL wrk_dealloc( jpi,jpj, zue, zve ) 359 IF( lk_vvl.AND.(.NOT.ln_dynadv_vec ) ) CALL wrk_dealloc( jpi,jpj,jpk, ze3u_f, ze3v_f ) 360 IF( l_trddyn ) CALL wrk_dealloc( jpi,jpj,jpk, zua, zva ) 388 361 ! 389 362 IF( nn_timing == 1 ) CALL timing_stop('dyn_nxt')
Note: See TracChangeset
for help on using the changeset viewer.
