- Timestamp:
- 2015-12-17T11:08:49+01:00 (9 years ago)
- Location:
- branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC
- Files:
-
- 2 deleted
- 10 edited
-
. (modified) (1 prop)
-
DYN/dynhpg.F90 (modified) (6 diffs)
-
DYN/dynnxt.F90 (modified) (9 diffs)
-
DYN/dynspg.F90 (modified) (14 diffs)
-
DYN/dynspg_exp.F90 (modified) (3 diffs)
-
DYN/dynspg_flt.F90 (deleted)
-
DYN/dynspg_oce.F90 (deleted)
-
DYN/dynspg_ts.F90 (modified) (27 diffs)
-
DYN/dynvor.F90 (modified) (3 diffs)
-
DYN/dynzdf_exp.F90 (modified) (3 diffs)
-
DYN/dynzdf_imp.F90 (modified) (11 diffs)
-
DYN/sshwzv.F90 (modified) (2 diffs)
Legend:
- Unmodified
- Added
- Removed
-
branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC
- Property svn:mergeinfo deleted
-
branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90
r5901 r6079 36 36 USE trd_oce ! trends: ocean variables 37 37 USE trddyn ! trend manager: dynamics 38 !jc USE zpshde ! partial step: hor. derivative (zps_hde routine) 38 39 ! 39 40 USE in_out_manager ! I/O manager … … 58 59 LOGICAL , PUBLIC :: ln_hpg_prj !: s-coordinate (Pressure Jacobian scheme) 59 60 LOGICAL , PUBLIC :: ln_hpg_isf !: s-coordinate similar to sco modify for isf 60 LOGICAL , PUBLIC :: ln_dynhpg_imp !: semi-implicite hpg flag61 61 62 62 INTEGER , PUBLIC :: nhpg = 0 ! = 0 to 7, type of pressure gradient scheme used ! (deduced from ln_hpg_... flags) (PUBLIC for TAM) … … 132 132 !! 133 133 NAMELIST/namdyn_hpg/ ln_hpg_zco, ln_hpg_zps, ln_hpg_sco, & 134 & ln_hpg_djc, ln_hpg_prj, ln_hpg_isf , ln_dynhpg_imp134 & ln_hpg_djc, ln_hpg_prj, ln_hpg_isf 135 135 !!---------------------------------------------------------------------- 136 136 ! … … 155 155 WRITE(numout,*) ' s-coord. (Density Jacobian: Cubic polynomial) ln_hpg_djc = ', ln_hpg_djc 156 156 WRITE(numout,*) ' s-coord. (Pressure Jacobian: Cubic polynomial) ln_hpg_prj = ', ln_hpg_prj 157 WRITE(numout,*) ' time stepping: centered (F) or semi-implicit (T) ln_dynhpg_imp = ', ln_dynhpg_imp158 157 ENDIF 159 158 ! … … 293 292 ENDIF 294 293 294 ! Partial steps: bottom before horizontal gradient of t, s, rd at the last ocean level 295 !jc CALL zps_hde ( kt, jpts, tsn, gtsu, gtsv, rhd, gru , grv ) 295 296 296 297 ! Local constant initialization … … 491 492 ! iniitialised to 0. zhpi zhpi 492 493 zhpi(:,:,:)=0._wp ; zhpj(:,:,:)=0._wp 494 495 ! Partial steps: top & bottom before horizontal gradient 496 !jc CALL zps_hde_isf( kt, jpts, tsn, gtsu, gtsv, gtui, gtvi, & 497 !jc & rhd, gru , grv , aru , arv , gzu , gzv , ge3ru , ge3rv , & 498 !jc & grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi ) 493 499 494 500 !================================================================================== -
branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90
r5901 r6079 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. 80 !! Note also that in filtered free surface (lk_dynspg_flt=T), 81 !! the time stepping has already been done in dynspg module 75 !! ** Method : * Ensure after velocities transport matches time splitting 76 !! estimate (ln_dynspg_ts=T) 82 77 !! 83 78 !! * Apply lateral boundary conditions on after velocity … … 101 96 INTEGER :: ji, jj, jk ! dummy loop indices 102 97 INTEGER :: iku, ikv ! local integers 103 #if ! defined key_dynspg_flt 104 REAL(wp) :: z2dt ! temporary scalar 105 #endif 106 REAL(wp) :: zue3a, zue3n, zue3b, zuf, zec ! local scalars 98 REAL(wp) :: zue3a, zue3n, zue3b, zuf ! local scalars 107 99 REAL(wp) :: zve3a, zve3n, zve3b, zvf, z1_2dt ! - - 108 100 REAL(wp), POINTER, DIMENSION(:,:) :: zue, zve … … 112 104 IF( nn_timing == 1 ) CALL timing_start('dyn_nxt') 113 105 ! 114 CALL wrk_alloc( jpi,jpj,jpk, ze3u_f, ze3v_f, zua, zva ) 115 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 ) 116 109 ! 117 110 IF( kt == nit000 ) THEN … … 121 114 ENDIF 122 115 123 #if defined key_dynspg_flt 124 ! 125 ! Next velocity : Leap-frog time stepping already done in dynspg_flt.F routine 126 ! ------------- 127 128 ! Update after velocity on domain lateral boundaries (only local domain required) 129 ! -------------------------------------------------- 130 CALL lbc_lnk( ua, 'U', -1. ) ! local domain boundaries 131 CALL lbc_lnk( va, 'V', -1. ) 132 ! 133 #else 134 135 # if defined key_dynspg_exp 136 ! Next velocity : Leap-frog time stepping 137 ! ------------- 138 z2dt = 2. * rdt ! Euler or leap-frog time step 139 IF( neuler == 0 .AND. kt == nit000 ) z2dt = rdt 140 ! 141 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 142 125 DO jk = 1, jpkm1 143 ua(:,:,jk) = ( ub(:,:,jk) + z2dt * ua(:,:,jk) ) * umask(:,:,jk) 144 va(:,:,jk) = ( vb(:,:,jk) + z2dt * va(:,:,jk) ) * vmask(:,:,jk) 145 END DO 146 ELSE ! applied on thickness weighted velocity 147 DO jk = 1, jpkm1 148 ua(:,:,jk) = ( ub(:,:,jk) * fse3u_b(:,:,jk) & 149 & + z2dt * ua(:,:,jk) * fse3u_n(:,:,jk) ) & 150 & / fse3u_a(:,:,jk) * umask(:,:,jk) 151 va(:,:,jk) = ( vb(:,:,jk) * fse3v_b(:,:,jk) & 152 & + z2dt * va(:,:,jk) * fse3v_n(:,:,jk) ) & 153 & / fse3v_a(:,:,jk) * vmask(:,:,jk) 154 END DO 155 ENDIF 156 # endif 157 158 # if defined key_dynspg_ts 159 !!gm IF ( lk_dynspg_ts ) THEN .... 160 ! Ensure below that barotropic velocities match time splitting estimate 161 ! Compute actual transport and replace it with ts estimate at "after" time step 162 zue(:,:) = fse3u_a(:,:,1) * ua(:,:,1) * umask(:,:,1) 163 zve(:,:) = fse3v_a(:,:,1) * va(:,:,1) * vmask(:,:,1) 164 DO jk = 2, jpkm1 165 zue(:,:) = zue(:,:) + fse3u_a(:,:,jk) * ua(:,:,jk) * umask(:,:,jk) 166 zve(:,:) = zve(:,:) + fse3v_a(:,:,jk) * va(:,:,jk) * vmask(:,:,jk) 167 END DO 168 DO jk = 1, jpkm1 169 ua(:,:,jk) = ( ua(:,:,jk) - zue(:,:) * hur_a(:,:) + ua_b(:,:) ) * umask(:,:,jk) 170 va(:,:,jk) = ( va(:,:,jk) - zve(:,:) * hvr_a(:,:) + va_b(:,:) ) * vmask(:,:,jk) 171 END DO 172 173 IF (lk_dynspg_ts.AND.(.NOT.ln_bt_fw)) THEN 174 ! Remove advective velocity from "now velocities" 175 ! prior to asselin filtering 176 ! In the forward case, this is done below after asselin filtering 177 ! so that asselin contribution is removed at the same time 178 DO jk = 1, jpkm1 179 un(:,:,jk) = ( un(:,:,jk) - un_adv(:,:) + un_b(:,:) )*umask(:,:,jk) 180 vn(:,:,jk) = ( vn(:,:,jk) - vn_adv(:,:) + vn_b(:,:) )*vmask(:,:,jk) 181 END DO 182 ENDIF 183 !!gm ENDIF 184 # 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 185 141 186 142 ! Update after velocity on domain lateral boundaries 187 143 ! -------------------------------------------------- 188 CALL lbc_lnk( ua, 'U', -1. ) !* local domain boundaries189 CALL lbc_lnk( va, 'V', -1. )190 !191 # if defined key_bdy192 ! !* BDY open boundaries193 IF( lk_bdy .AND. lk_dynspg_exp ) CALL bdy_dyn( kt )194 IF( lk_bdy .AND. lk_dynspg_ts ) CALL bdy_dyn( kt, dyn3d_only=.true. )195 196 !!$ Do we need a call to bdy_vol here??197 !198 # endif199 !200 144 # if defined key_agrif 201 145 CALL Agrif_dyn( kt ) !* AGRIF zoom boundaries 202 146 # endif 203 #endif 204 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 ! 205 160 IF( l_trddyn ) THEN ! prepare the atf trend computation + some diagnostics 206 161 z1_2dt = 1._wp / (2. * rdt) ! Euler or leap-frog time step … … 259 214 ! (used as a now filtered scale factor until the swap) 260 215 ! ---------------------------------------------------- 261 IF (l k_dynspg_ts.AND.ln_bt_fw) THEN216 IF (ln_dynspg_ts.AND.ln_bt_fw) THEN 262 217 ! No asselin filtering on thicknesses if forward time splitting 263 218 fse3t_b(:,:,:) = fse3t_n(:,:,:) 264 219 ELSE 265 220 fse3t_b(:,:,:) = fse3t_n(:,:,:) + atfp * ( fse3t_b(:,:,:) - 2._wp * fse3t_n(:,:,:) + fse3t_a(:,:,:) ) … … 336 291 ENDIF 337 292 ! 338 IF (l k_dynspg_ts.AND.ln_bt_fw) THEN293 IF (ln_dynspg_ts.AND.ln_bt_fw) THEN 339 294 ! Revert "before" velocities to time split estimate 340 295 ! Doing it here also means that asselin filter contribution is removed … … 400 355 IF(ln_ctl) CALL prt_ctl( tab3d_1=un, clinfo1=' nxt - Un: ', mask1=umask, & 401 356 & tab3d_2=vn, clinfo2=' Vn: ' , mask2=vmask ) 402 ! 403 CALL wrk_dealloc( jpi,jpj,jpk, ze3u_f, ze3v_f, zua, zva ) 404 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 ) 405 361 ! 406 362 IF( nn_timing == 1 ) CALL timing_stop('dyn_nxt') -
branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg.F90
r5901 r6079 6 6 !! History : 1.0 ! 2005-12 (C. Talandier, G. Madec, V. Garnier) Original code 7 7 !! 3.2 ! 2009-07 (R. Benshila) Suppression of rigid-lid option 8 !!---------------------------------------------------------------------- 9 10 !!---------------------------------------------------------------------- 11 !! dyn_spg : update the dynamics trend with the lateral diffusion 12 !! dyn_spg_ctl : initialization, namelist read, and parameters control 8 !! 3.7 ! 2015-11 (J. Chanut) Suppression of filtered free surface 9 !!---------------------------------------------------------------------- 10 11 !!---------------------------------------------------------------------- 12 !! dyn_spg : update the dynamics trend with surface pressure gradient 13 !! dyn_spg_init: initialization, namelist read, and parameters control 13 14 !!---------------------------------------------------------------------- 14 15 USE oce ! ocean dynamics and tracers variables … … 18 19 USE sbc_oce ! surface boundary condition: ocean 19 20 USE sbcapr ! surface boundary condition: atmospheric pressure 20 USE dynspg_oce ! surface pressure gradient variables21 21 USE dynspg_exp ! surface pressure gradient (dyn_spg_exp routine) 22 22 USE dynspg_ts ! surface pressure gradient (dyn_spg_ts routine) 23 USE dynspg_flt ! surface pressure gradient (dyn_spg_flt routine)24 USE dynadv ! dynamics: vector invariant versus flux form25 USE dynhpg, ONLY: ln_dynhpg_imp26 23 USE sbctide 27 24 USE updtide … … 32 29 USE in_out_manager ! I/O manager 33 30 USE lib_mpp ! MPP library 34 USE solver ! solver initialization35 31 USE wrk_nemo ! Memory Allocation 36 32 USE timing ! Timing … … 43 39 PUBLIC dyn_spg_init ! routine called by opa module 44 40 45 INTEGER :: nspg = 0 ! type of surface pressure gradient scheme defined from l k_dynspg_...41 INTEGER :: nspg = 0 ! type of surface pressure gradient scheme defined from ln_dynspg_... 46 42 47 43 !! * Substitutions … … 55 51 CONTAINS 56 52 57 SUBROUTINE dyn_spg( kt , kindic)53 SUBROUTINE dyn_spg( kt ) 58 54 !!---------------------------------------------------------------------- 59 55 !! *** ROUTINE dyn_spg *** … … 66 62 !!gm the after velocity, in the 2 other (ua,va) are still the trends 67 63 !! 68 !! ** Method : T hreeschemes:64 !! ** Method : Two schemes: 69 65 !! - explicit computation : the spg is evaluated at now 70 !! - filtered computation : the Roulet & madec (2000) technique is used71 66 !! - split-explicit computation: a time splitting technique is used 72 67 !! … … 80 75 ! 81 76 INTEGER, INTENT(in ) :: kt ! ocean time-step index 82 INTEGER, INTENT( out) :: kindic ! solver flag83 77 ! 84 78 INTEGER :: ji, jj, jk ! dummy loop indices … … 103 97 104 98 IF( ln_apr_dyn & ! atmos. pressure 105 .OR. ( .NOT.l k_dynspg_ts .AND. (ln_tide_pot .AND. lk_tide) ) & ! tide potential (no time slitting)99 .OR. ( .NOT.ln_dynspg_ts .AND. (ln_tide_pot .AND. lk_tide) ) & ! tide potential (no time slitting) 106 100 .OR. nn_ice_embd == 2 ) THEN ! embedded sea-ice 107 101 ! … … 113 107 END DO 114 108 ! 115 IF( ln_apr_dyn .AND. (.NOT. l k_dynspg_ts) ) THEN !== Atmospheric pressure gradient (added later in time-split case) ==!109 IF( ln_apr_dyn .AND. (.NOT. ln_dynspg_ts) ) THEN !== Atmospheric pressure gradient (added later in time-split case) ==! 116 110 zg_2 = grav * 0.5 117 111 DO jj = 2, jpjm1 ! gradient of Patm using inverse barometer ssh … … 126 120 ! 127 121 ! !== tide potential forcing term ==! 128 IF( .NOT.l k_dynspg_ts .AND. ( ln_tide_pot .AND. lk_tide ) ) THEN ! N.B. added directly at sub-time-step in ts-case122 IF( .NOT.ln_dynspg_ts .AND. ( ln_tide_pot .AND. lk_tide ) ) THEN ! N.B. added directly at sub-time-step in ts-case 129 123 ! 130 124 CALL upd_tide( kt ) ! update tide potential … … 171 165 CASE ( 0 ) ; CALL dyn_spg_exp( kt ) ! explicit 172 166 CASE ( 1 ) ; CALL dyn_spg_ts ( kt ) ! time-splitting 173 CASE ( 2 ) ; CALL dyn_spg_flt( kt, kindic ) ! filtered174 167 ! 175 168 END SELECT 176 169 ! 177 170 IF( l_trddyn ) THEN ! save the surface pressure gradient trends for further diagnostics 178 SELECT CASE ( nspg ) 179 CASE ( 0, 1 ) 180 ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 181 ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 182 CASE( 2 ) 183 z2dt = 2. * rdt 184 IF( neuler == 0 .AND. kt == nit000 ) z2dt = rdt 185 ztrdu(:,:,:) = ( ua(:,:,:) - ub(:,:,:) ) / z2dt - ztrdu(:,:,:) 186 ztrdv(:,:,:) = ( va(:,:,:) - vb(:,:,:) ) / z2dt - ztrdv(:,:,:) 187 END SELECT 171 ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 172 ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 188 173 CALL trd_dyn( ztrdu, ztrdv, jpdyn_spg, kt ) 189 174 ! … … 203 188 !! *** ROUTINE dyn_spg_init *** 204 189 !! 205 !! ** Purpose : Control the consistency between cppoptions for190 !! ** Purpose : Control the consistency between namelist options for 206 191 !! surface pressure gradient schemes 207 192 !!---------------------------------------------------------------------- 208 INTEGER :: ioptio 193 INTEGER :: ioptio, ios 194 ! 195 NAMELIST/namdyn_spg/ ln_dynspg_exp, ln_dynspg_ts, & 196 & ln_bt_fw, ln_bt_av, ln_bt_auto, & 197 & nn_baro, rn_bt_cmax, nn_bt_flt 209 198 !!---------------------------------------------------------------------- 210 199 ! 211 200 IF( nn_timing == 1 ) CALL timing_start('dyn_spg_init') 212 201 ! 213 IF(lwp) THEN ! Control print 202 REWIND( numnam_ref ) ! Namelist namdyn_spg in reference namelist : Free surface 203 READ ( numnam_ref, namdyn_spg, IOSTAT = ios, ERR = 901) 204 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdyn_spg in reference namelist', lwp ) 205 206 REWIND( numnam_cfg ) ! Namelist namdyn_spg in configuration namelist : Free surface 207 READ ( numnam_cfg, namdyn_spg, IOSTAT = ios, ERR = 902 ) 208 902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdyn_spg in configuration namelist', lwp ) 209 IF(lwm) WRITE ( numond, namdyn_spg ) 210 ! 211 IF(lwp) THEN ! Namelist print 214 212 WRITE(numout,*) 215 213 WRITE(numout,*) 'dyn_spg_init : choice of the surface pressure gradient scheme' 216 214 WRITE(numout,*) '~~~~~~~~~~~' 217 WRITE(numout,*) ' Explicit free surface lk_dynspg_exp = ', lk_dynspg_exp 218 WRITE(numout,*) ' Free surface with time splitting lk_dynspg_ts = ', lk_dynspg_ts 219 WRITE(numout,*) ' Filtered free surface cst volume lk_dynspg_flt = ', lk_dynspg_flt 220 ENDIF 221 222 IF( lk_dynspg_ts ) CALL dyn_spg_ts_init( nit000 ) 223 ! (do it now, to set nn_baro, used to allocate some arrays later on) 224 ! ! allocate dyn_spg arrays 225 IF( lk_dynspg_ts ) THEN 226 IF( dynspg_oce_alloc() /= 0 ) CALL ctl_stop('STOP', 'dyn_spg_init: failed to allocate dynspg_oce arrays') 215 WRITE(numout,*) ' Explicit free surface ln_dynspg_exp = ', ln_dynspg_exp 216 WRITE(numout,*) ' Free surface with time splitting ln_dynspg_ts = ', ln_dynspg_ts 217 ENDIF 218 219 IF( ln_dynspg_ts ) THEN 220 CALL dyn_spg_ts_init( nit000 ) ! do it first, to set nn_baro, used to allocate some arrays later on 221 ! ! allocate dyn_spg arrays 227 222 IF( dyn_spg_ts_alloc() /= 0 ) CALL ctl_stop('STOP', 'dyn_spg_init: failed to allocate dynspg_ts arrays') 228 223 IF ((neuler/=0).AND.(ln_bt_fw)) CALL ts_rst( nit000, 'READ' ) … … 231 226 ! ! Control of surface pressure gradient scheme options 232 227 ioptio = 0 233 IF(lk_dynspg_exp) ioptio = ioptio + 1 234 IF(lk_dynspg_ts ) ioptio = ioptio + 1 235 IF(lk_dynspg_flt) ioptio = ioptio + 1 228 IF(ln_dynspg_exp) ioptio = ioptio + 1 229 IF(ln_dynspg_ts ) ioptio = ioptio + 1 236 230 ! 237 231 IF( ioptio > 1 .OR. ( ioptio == 0 .AND. .NOT. lk_c1d ) ) & 238 & CALL ctl_stop( ' Choose only one surface pressure gradient scheme with a key cpp' ) 239 IF( ( lk_dynspg_ts .OR. lk_dynspg_exp ) .AND. ln_isfcav ) & 240 & CALL ctl_stop( ' dynspg_ts and dynspg_exp not tested with ice shelf cavity ' ) 241 ! 242 IF( lk_dynspg_exp) nspg = 0 243 IF( lk_dynspg_ts ) nspg = 1 244 IF( lk_dynspg_flt) nspg = 2 232 & CALL ctl_stop( ' Choose only one surface pressure gradient scheme' ) 233 IF( ln_dynspg_ts .AND. ln_isfcav ) & 234 & CALL ctl_stop( ' dynspg_ts not tested with ice shelf cavity ' ) 235 ! 236 IF( ln_dynspg_exp) nspg = 0 237 IF( ln_dynspg_ts ) nspg = 1 245 238 ! 246 239 IF(lwp) THEN … … 248 241 IF( nspg == 0 ) WRITE(numout,*) ' explicit free surface' 249 242 IF( nspg == 1 ) WRITE(numout,*) ' free surface with time splitting scheme' 250 IF( nspg == 2 ) WRITE(numout,*) ' filtered free surface' 251 ENDIF 252 253 #if defined key_dynspg_flt 254 CALL solver_init( nit000 ) ! Elliptic solver initialisation 255 #endif 256 ! ! Control of hydrostatic pressure choice 257 IF( lk_dynspg_ts .AND. ln_dynhpg_imp ) CALL ctl_stop( 'Semi-implicit hpg not compatible with time splitting' ) 243 ENDIF 258 244 ! 259 245 IF( nn_timing == 1 ) CALL timing_stop('dyn_spg_init') -
branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_exp.F90
r5901 r6079 7 7 !! 3.2 ! 2009-06 (G. Madec, M. Leclair, R. Benshila) introduce sshwzv module 8 8 !!---------------------------------------------------------------------- 9 #if defined key_dynspg_exp10 9 !!---------------------------------------------------------------------- 11 !! 'key_dynspg_exp'explicit free surface10 !! explicit free surface 12 11 !!---------------------------------------------------------------------- 13 12 !! dyn_spg_exp : update the momentum trend with the surface … … 31 30 PRIVATE 32 31 33 PUBLIC dyn_spg_exp ! routine called by step.F9032 PUBLIC dyn_spg_exp ! routine called by dyn_spg 34 33 35 34 !! * Substitutions … … 101 100 END SUBROUTINE dyn_spg_exp 102 101 103 #else104 !!----------------------------------------------------------------------105 !! Default case : Empty module No standart explicit free surface106 !!----------------------------------------------------------------------107 CONTAINS108 SUBROUTINE dyn_spg_exp( kt ) ! Empty routine109 WRITE(*,*) 'dyn_spg_exp: You should not have seen this print! error?', kt110 END SUBROUTINE dyn_spg_exp111 #endif112 113 102 !!====================================================================== 114 103 END MODULE dynspg_exp -
branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90
r5901 r6079 11 11 !! 3.5 ! 2013-07 (J. Chanut) Switch to Forward-backward time stepping 12 12 !! 3.6 ! 2013-11 (A. Coward) Update for z-tilde compatibility 13 !! 3.7 ! 2015-11 (J. Chanut) free surface simplification 13 14 !!--------------------------------------------------------------------- 14 #if defined key_dynspg_ts15 15 !!---------------------------------------------------------------------- 16 !! 'key_dynspg_ts'split explicit free surface16 !! split explicit free surface 17 17 !!---------------------------------------------------------------------- 18 18 !! dyn_spg_ts : compute surface pressure gradient trend using a time- … … 23 23 USE sbc_oce ! surface boundary condition: ocean 24 24 USE sbcisf ! ice shelf variable (fwfisf) 25 USE dynspg_oce ! surface pressure gradient variables26 25 USE phycst ! physical constants 27 26 USE dynvor ! vorticity term 28 27 USE bdy_par ! for lk_bdy 29 USE bdytides ! open boundary condition data 28 USE bdytides ! open boundary condition data 30 29 USE bdydyn2d ! open boundary conditions on barotropic variables 31 30 USE sbctide ! tides … … 68 67 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: ftsw, ftse ! (only used with een vorticity scheme) 69 68 70 ! Arrays below are saved to allow testing of the "no time averaging" option71 ! If this option is not retained, these could be replaced by temporary arrays72 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: sshbb_e, sshb_e, & ! Instantaneous barotropic arrays73 ubb_e, ub_e, &74 vbb_e, vb_e75 76 69 !! * Substitutions 77 70 # include "domzgr_substitute.h90" … … 88 81 !! *** routine dyn_spg_ts_alloc *** 89 82 !!---------------------------------------------------------------------- 90 INTEGER :: ierr( 3)83 INTEGER :: ierr(4) 91 84 !!---------------------------------------------------------------------- 92 85 ierr(:) = 0 93 86 94 ALLOCATE( sshb_e(jpi,jpj), sshbb_e(jpi,jpj), & 95 & ub_e(jpi,jpj) , vb_e(jpi,jpj) , & 96 & ubb_e(jpi,jpj) , vbb_e(jpi,jpj) , STAT= ierr(1) ) 87 ALLOCATE( ssha_e(jpi,jpj), sshn_e(jpi,jpj), sshb_e(jpi,jpj), sshbb_e(jpi,jpj), & 88 & ua_e(jpi,jpj), un_e(jpi,jpj), ub_e(jpi,jpj), ubb_e(jpi,jpj), & 89 & va_e(jpi,jpj), vn_e(jpi,jpj), vb_e(jpi,jpj), vbb_e(jpi,jpj), & 90 & hu_e(jpi,jpj), hur_e(jpi,jpj), hv_e(jpi,jpj), hvr_e(jpi,jpj), STAT= ierr(1) ) 97 91 98 92 ALLOCATE( wgtbtp1(3*nn_baro), wgtbtp2(3*nn_baro), zwz(jpi,jpj), STAT= ierr(2) ) … … 101 95 & ftsw(jpi,jpj) , ftse(jpi,jpj) , STAT=ierr(3) ) 102 96 97 ALLOCATE( ub2_b(jpi,jpj), vb2_b(jpi,jpj), un_adv(jpi,jpj), vn_adv(jpi,jpj), & 98 #if defined key_agrif 99 & ub2_i_b(jpi,jpj), vb2_i_b(jpi,jpj) , & 100 #endif 101 & STAT= ierr(4)) 102 103 103 dyn_spg_ts_alloc = MAXVAL(ierr(:)) 104 104 105 105 IF( lk_mpp ) CALL mpp_sum( dyn_spg_ts_alloc ) 106 IF( dyn_spg_ts_alloc /= 0 ) CALL ctl_warn('dyn spg_oce_alloc: failed to allocate arrays')106 IF( dyn_spg_ts_alloc /= 0 ) CALL ctl_warn('dyn_spg_ts_alloc: failed to allocate arrays') 107 107 ! 108 108 END FUNCTION dyn_spg_ts_alloc … … 146 146 INTEGER :: ikbu, ikbv, noffset ! local integers 147 147 REAL(wp) :: zraur, z1_2dt_b, z2dt_bf ! local scalars 148 REAL(wp) :: zx1, zy1, zx2, zy2 ! - -149 REAL(wp) :: z1_12, z1_8, z1_4, z1_2 150 REAL(wp) :: zu_spg, zv_spg ! - -151 REAL(wp) :: zhura, zhvra 152 REAL(wp) :: za0, za1, za2, za3 153 ! 154 REAL(wp), POINTER, DIMENSION(:,:) :: z un_e, zvn_e, zsshp2_e148 REAL(wp) :: zx1, zy1, zx2, zy2 ! - - 149 REAL(wp) :: z1_12, z1_8, z1_4, z1_2 ! - - 150 REAL(wp) :: zu_spg, zv_spg ! - - 151 REAL(wp) :: zhura, zhvra ! - - 152 REAL(wp) :: za0, za1, za2, za3 ! - - 153 ! 154 REAL(wp), POINTER, DIMENSION(:,:) :: zsshp2_e 155 155 REAL(wp), POINTER, DIMENSION(:,:) :: zu_trd, zv_trd, zu_frc, zv_frc, zssh_frc 156 REAL(wp), POINTER, DIMENSION(:,:) :: z u_sum, zv_sum, zwx, zwy, zhdiv156 REAL(wp), POINTER, DIMENSION(:,:) :: zwx, zwy, zhdiv 157 157 REAL(wp), POINTER, DIMENSION(:,:) :: zhup2_e, zhvp2_e, zhust_e, zhvst_e 158 158 REAL(wp), POINTER, DIMENSION(:,:) :: zsshu_a, zsshv_a … … 164 164 ! !* Allocate temporary arrays 165 165 CALL wrk_alloc( jpi, jpj, zsshp2_e, zhdiv ) 166 CALL wrk_alloc( jpi, jpj, zu_trd, zv_trd , zun_e, zvn_e)167 CALL wrk_alloc( jpi, jpj, zwx, zwy, z u_sum, zv_sum, zssh_frc, zu_frc, zv_frc)166 CALL wrk_alloc( jpi, jpj, zu_trd, zv_trd) 167 CALL wrk_alloc( jpi, jpj, zwx, zwy, zssh_frc, zu_frc, zv_frc) 168 168 CALL wrk_alloc( jpi, jpj, zhup2_e, zhvp2_e, zhust_e, zhvst_e) 169 169 CALL wrk_alloc( jpi, jpj, zsshu_a, zsshv_a ) … … 188 188 ! 189 189 ! time offset in steps for bdy data update 190 IF (.NOT.ln_bt_fw) THEN ; noffset=- 2*nn_baro ; ELSE ; noffset = 0 ; ENDIF190 IF (.NOT.ln_bt_fw) THEN ; noffset=-nn_baro ; ELSE ; noffset = 0 ; ENDIF 191 191 ! 192 192 IF( kt == nit000 ) THEN !* initialisation … … 485 485 IF (ln_bt_fw) THEN ! FORWARD integration: start from NOW fields 486 486 sshn_e(:,:) = sshn (:,:) 487 zun_e(:,:) = un_b (:,:)488 zvn_e(:,:) = vn_b (:,:)487 un_e (:,:) = un_b (:,:) 488 vn_e (:,:) = vn_b (:,:) 489 489 ! 490 490 hu_e (:,:) = hu (:,:) … … 494 494 ELSE ! CENTRED integration: start from BEFORE fields 495 495 sshn_e(:,:) = sshb (:,:) 496 zun_e(:,:) = ub_b (:,:)497 zvn_e(:,:) = vb_b (:,:)496 un_e (:,:) = ub_b (:,:) 497 vn_e (:,:) = vb_b (:,:) 498 498 ! 499 499 hu_e (:,:) = hu_b (:,:) … … 509 509 va_b (:,:) = 0._wp 510 510 ssha (:,:) = 0._wp ! Sum for after averaged sea level 511 zu_sum(:,:) = 0._wp ! Sum for now transport issued from ts loop512 zv_sum(:,:) = 0._wp511 un_adv(:,:) = 0._wp ! Sum for now transport issued from ts loop 512 vn_adv(:,:) = 0._wp 513 513 ! ! ==================== ! 514 514 DO jn = 1, icycle ! sub-time-step loop ! … … 518 518 ! Update only tidal forcing at open boundaries 519 519 #if defined key_tide 520 IF ( lk_bdy .AND. lk_tide ) 521 IF ( ln_tide_pot .AND. lk_tide ) CALL upd_tide( kt, kit=jn, koffset=noffset )520 IF ( lk_bdy .AND. lk_tide ) CALL bdy_dta_tides( kt, kit=jn, time_offset=(noffset+1) ) 521 IF ( ln_tide_pot .AND. lk_tide ) CALL upd_tide( kt, kit=jn, time_offset=noffset ) 522 522 #endif 523 523 ! … … 534 534 535 535 ! Extrapolate barotropic velocities at step jit+0.5: 536 ua_e(:,:) = za1 * zun_e(:,:) + za2 * ub_e(:,:) + za3 * ubb_e(:,:)537 va_e(:,:) = za1 * zvn_e(:,:) + za2 * vb_e(:,:) + za3 * vbb_e(:,:)536 ua_e(:,:) = za1 * un_e(:,:) + za2 * ub_e(:,:) + za3 * ubb_e(:,:) 537 va_e(:,:) = za1 * vn_e(:,:) + za2 * vb_e(:,:) + za3 * vbb_e(:,:) 538 538 539 539 IF( lk_vvl ) THEN !* Update ocean depth (variable volume case only) … … 597 597 ! Sum over sub-time-steps to compute advective velocities 598 598 za2 = wgtbtp2(jn) 599 zu_sum(:,:) = zu_sum(:,:) + za2 * zwx(:,:) * r1_e2u(:,:)600 zv_sum(:,:) = zv_sum(:,:) + za2 * zwy(:,:) * r1_e1v(:,:)599 un_adv(:,:) = un_adv(:,:) + za2 * zwx(:,:) * r1_e2u(:,:) 600 vn_adv(:,:) = vn_adv(:,:) + za2 * zwy(:,:) * r1_e1v(:,:) 601 601 ! 602 602 ! Set next sea level: … … 733 733 ! 734 734 ! Add bottom stresses: 735 zu_trd(:,:) = zu_trd(:,:) + bfrua(:,:) * zun_e(:,:) * hur_e(:,:)736 zv_trd(:,:) = zv_trd(:,:) + bfrva(:,:) * zvn_e(:,:) * hvr_e(:,:)735 zu_trd(:,:) = zu_trd(:,:) + bfrua(:,:) * un_e(:,:) * hur_e(:,:) 736 zv_trd(:,:) = zv_trd(:,:) + bfrva(:,:) * vn_e(:,:) * hvr_e(:,:) 737 737 ! 738 738 ! Surface pressure trend: … … 751 751 DO jj = 2, jpjm1 752 752 DO ji = fs_2, fs_jpim1 ! vector opt. 753 ua_e(ji,jj) = ( zun_e(ji,jj) &753 ua_e(ji,jj) = ( un_e(ji,jj) & 754 754 & + rdtbt * ( zwx(ji,jj) & 755 755 & + zu_trd(ji,jj) & … … 757 757 & ) * umask(ji,jj,1) 758 758 759 va_e(ji,jj) = ( zvn_e(ji,jj) &759 va_e(ji,jj) = ( vn_e(ji,jj) & 760 760 & + rdtbt * ( zwy(ji,jj) & 761 761 & + zv_trd(ji,jj) & … … 772 772 zhvra = vmask(ji,jj,1)/(hv_0(ji,jj) + zsshv_a(ji,jj) + 1._wp - vmask(ji,jj,1)) 773 773 774 ua_e(ji,jj) = ( hu_e(ji,jj) * zun_e(ji,jj) &774 ua_e(ji,jj) = ( hu_e(ji,jj) * un_e(ji,jj) & 775 775 & + rdtbt * ( zhust_e(ji,jj) * zwx(ji,jj) & 776 776 & + zhup2_e(ji,jj) * zu_trd(ji,jj) & … … 778 778 & ) * zhura 779 779 780 va_e(ji,jj) = ( hv_e(ji,jj) * zvn_e(ji,jj) &780 va_e(ji,jj) = ( hv_e(ji,jj) * vn_e(ji,jj) & 781 781 & + rdtbt * ( zhvst_e(ji,jj) * zwy(ji,jj) & 782 782 & + zhvp2_e(ji,jj) * zv_trd(ji,jj) & … … 802 802 #if defined key_bdy 803 803 ! open boundaries 804 IF( lk_bdy ) CALL bdy_dyn2d( jn, ua_e, va_e, zun_e, zvn_e, hur_e, hvr_e, ssha_e )804 IF( lk_bdy ) CALL bdy_dyn2d( jn, ua_e, va_e, un_e, vn_e, hur_e, hvr_e, ssha_e ) 805 805 #endif 806 806 #if defined key_agrif … … 810 810 ! ! ---- 811 811 ubb_e (:,:) = ub_e (:,:) 812 ub_e (:,:) = zun_e(:,:)813 zun_e(:,:) = ua_e (:,:)812 ub_e (:,:) = un_e (:,:) 813 un_e (:,:) = ua_e (:,:) 814 814 ! 815 815 vbb_e (:,:) = vb_e (:,:) 816 vb_e (:,:) = zvn_e(:,:)817 zvn_e(:,:) = va_e (:,:)816 vb_e (:,:) = vn_e (:,:) 817 vn_e (:,:) = va_e (:,:) 818 818 ! 819 819 sshbb_e(:,:) = sshb_e(:,:) … … 840 840 ! ----------------------------------------------------------------------------- 841 841 ! 842 ! At this stage ssha holds a time averaged value 843 ! ! Sea Surface Height at u-,v- and f-points 844 IF( lk_vvl ) THEN ! (required only in key_vvl case) 842 ! Set advection velocity correction: 843 zwx(:,:) = un_adv(:,:) 844 zwy(:,:) = vn_adv(:,:) 845 IF (((kt==nit000).AND.(neuler==0)).OR.(.NOT.ln_bt_fw)) THEN 846 un_adv(:,:) = zwx(:,:)*hur(:,:) 847 vn_adv(:,:) = zwy(:,:)*hvr(:,:) 848 ELSE 849 un_adv(:,:) = z1_2 * ( ub2_b(:,:) + zwx(:,:)) * hur(:,:) 850 vn_adv(:,:) = z1_2 * ( vb2_b(:,:) + zwy(:,:)) * hvr(:,:) 851 END IF 852 853 IF (ln_bt_fw) THEN ! Save integrated transport for next computation 854 ub2_b(:,:) = zwx(:,:) 855 vb2_b(:,:) = zwy(:,:) 856 ENDIF 857 ! 858 ! Update barotropic trend: 859 IF (( ln_dynadv_vec ).OR. (.NOT. lk_vvl)) THEN 860 DO jk=1,jpkm1 861 ua(:,:,jk) = ua(:,:,jk) + ( ua_b(:,:) - ub_b(:,:) ) * z1_2dt_b 862 va(:,:,jk) = va(:,:,jk) + ( va_b(:,:) - vb_b(:,:) ) * z1_2dt_b 863 END DO 864 ELSE 865 ! At this stage, ssha has been corrected: compute new depths at velocity points 845 866 DO jj = 1, jpjm1 846 867 DO ji = 1, jpim1 ! NO Vector Opt. … … 854 875 END DO 855 876 CALL lbc_lnk_multi( zsshu_a, 'U', 1._wp, zsshv_a, 'V', 1._wp ) ! Boundary conditions 856 ENDIF 857 ! 858 ! Set advection velocity correction: 859 IF (((kt==nit000).AND.(neuler==0)).OR.(.NOT.ln_bt_fw)) THEN 860 un_adv(:,:) = zu_sum(:,:)*hur(:,:) 861 vn_adv(:,:) = zv_sum(:,:)*hvr(:,:) 862 ELSE 863 un_adv(:,:) = z1_2 * ( ub2_b(:,:) + zu_sum(:,:)) * hur(:,:) 864 vn_adv(:,:) = z1_2 * ( vb2_b(:,:) + zv_sum(:,:)) * hvr(:,:) 865 END IF 866 867 IF (ln_bt_fw) THEN ! Save integrated transport for next computation 868 ub2_b(:,:) = zu_sum(:,:) 869 vb2_b(:,:) = zv_sum(:,:) 870 ENDIF 871 ! 872 ! Update barotropic trend: 873 IF (( ln_dynadv_vec ).OR. (.NOT. lk_vvl)) THEN 874 DO jk=1,jpkm1 875 ua(:,:,jk) = ua(:,:,jk) + ( ua_b(:,:) - ub_b(:,:) ) * z1_2dt_b 876 va(:,:,jk) = va(:,:,jk) + ( va_b(:,:) - vb_b(:,:) ) * z1_2dt_b 877 END DO 878 ELSE 877 ! 879 878 DO jk=1,jpkm1 880 879 ua(:,:,jk) = ua(:,:,jk) + hur(:,:) * ( ua_b(:,:) - ub_b(:,:) * hu_b(:,:) ) * z1_2dt_b … … 915 914 ! 916 915 CALL wrk_dealloc( jpi, jpj, zsshp2_e, zhdiv ) 917 CALL wrk_dealloc( jpi, jpj, zu_trd, zv_trd , zun_e, zvn_e)918 CALL wrk_dealloc( jpi, jpj, zwx, zwy, z u_sum, zv_sum, zssh_frc, zu_frc, zv_frc )916 CALL wrk_dealloc( jpi, jpj, zu_trd, zv_trd ) 917 CALL wrk_dealloc( jpi, jpj, zwx, zwy, zssh_frc, zu_frc, zv_frc ) 919 918 CALL wrk_dealloc( jpi, jpj, zhup2_e, zhvp2_e, zhust_e, zhvst_e ) 920 919 CALL wrk_dealloc( jpi, jpj, zsshu_a, zsshv_a ) … … 1064 1063 ! 1065 1064 INTEGER :: ji ,jj 1066 INTEGER :: ios ! Local integer output status for namelist read1067 1065 REAL(wp) :: zxr2, zyr2, zcmax 1068 1066 REAL(wp), POINTER, DIMENSION(:,:) :: zcu 1069 1067 !! 1070 NAMELIST/namsplit/ ln_bt_fw, ln_bt_av, ln_bt_nn_auto, &1071 & nn_baro, rn_bt_cmax, nn_bt_flt1072 1068 !!---------------------------------------------------------------------- 1073 1069 ! 1074 REWIND( numnam_ref ) ! Namelist namsplit in reference namelist : time splitting parameters 1075 READ ( numnam_ref, namsplit, IOSTAT = ios, ERR = 901) 1076 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsplit in reference namelist', lwp ) 1077 1078 REWIND( numnam_cfg ) ! Namelist namsplit in configuration namelist : time splitting parameters 1079 READ ( numnam_cfg, namsplit, IOSTAT = ios, ERR = 902 ) 1080 902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsplit in configuration namelist', lwp ) 1081 IF(lwm) WRITE ( numond, namsplit ) 1082 ! 1083 ! ! Max courant number for ext. grav. waves 1070 ! Max courant number for ext. grav. waves 1084 1071 ! 1085 1072 CALL wrk_alloc( jpi, jpj, zcu ) 1086 1073 ! 1087 IF (lk_vvl) THEN 1088 DO jj = 1, jpj 1089 DO ji =1, jpi 1090 zxr2 = r1_e1t(ji,jj) * r1_e1t(ji,jj) 1091 zyr2 = r1_e2t(ji,jj) * r1_e2t(ji,jj) 1092 zcu(ji,jj) = SQRT( grav * ht_0(ji,jj) * (zxr2 + zyr2) ) 1093 END DO 1094 END DO 1095 ELSE 1096 DO jj = 1, jpj 1097 DO ji =1, jpi 1098 zxr2 = r1_e1t(ji,jj) * r1_e1t(ji,jj) 1099 zyr2 = r1_e1t(ji,jj) * r1_e1t(ji,jj) 1100 zcu(ji,jj) = SQRT( grav * ht(ji,jj) * (zxr2 + zyr2) ) 1101 END DO 1102 END DO 1103 ENDIF 1104 1074 DO jj = 1, jpj 1075 DO ji =1, jpi 1076 zxr2 = r1_e1t(ji,jj) * r1_e1t(ji,jj) 1077 zyr2 = r1_e2t(ji,jj) * r1_e2t(ji,jj) 1078 zcu(ji,jj) = SQRT( grav * ht_0(ji,jj) * (zxr2 + zyr2) ) 1079 END DO 1080 END DO 1081 ! 1105 1082 zcmax = MAXVAL( zcu(:,:) ) 1106 1083 IF( lk_mpp ) CALL mpp_max( zcmax ) 1107 1084 1108 1085 ! Estimate number of iterations to satisfy a max courant number= rn_bt_cmax 1109 IF (ln_bt_ nn_auto) nn_baro = CEILING( rdt / rn_bt_cmax * zcmax)1086 IF (ln_bt_auto) nn_baro = CEILING( rdt / rn_bt_cmax * zcmax) 1110 1087 1111 1088 rdtbt = rdt / REAL( nn_baro , wp ) … … 1115 1092 IF(lwp) WRITE(numout,*) 'dyn_spg_ts : split-explicit free surface' 1116 1093 IF(lwp) WRITE(numout,*) '~~~~~~~~~~' 1117 IF( ln_bt_ nn_auto ) THEN1118 IF(lwp) WRITE(numout,*) ' ln_ts_ nn_auto=.true. Automatically set nn_baro '1094 IF( ln_bt_auto ) THEN 1095 IF(lwp) WRITE(numout,*) ' ln_ts_auto=.true. Automatically set nn_baro ' 1119 1096 IF(lwp) WRITE(numout,*) ' Max. courant number allowed: ', rn_bt_cmax 1120 1097 ELSE 1121 IF(lwp) WRITE(numout,*) ' ln_ts_ nn_auto=.false.: Use nn_baro in namelist '1098 IF(lwp) WRITE(numout,*) ' ln_ts_auto=.false.: Use nn_baro in namelist ' 1122 1099 ENDIF 1123 1100 … … 1164 1141 END SUBROUTINE dyn_spg_ts_init 1165 1142 1166 #else1167 !!---------------------------------------------------------------------------1168 !! Default case : Empty module No split explicit free surface1169 !!---------------------------------------------------------------------------1170 CONTAINS1171 INTEGER FUNCTION dyn_spg_ts_alloc() ! Dummy function1172 dyn_spg_ts_alloc = 01173 END FUNCTION dyn_spg_ts_alloc1174 SUBROUTINE dyn_spg_ts( kt ) ! Empty routine1175 INTEGER, INTENT(in) :: kt1176 WRITE(*,*) 'dyn_spg_ts: You should not have seen this print! error?', kt1177 END SUBROUTINE dyn_spg_ts1178 SUBROUTINE ts_rst( kt, cdrw ) ! Empty routine1179 INTEGER , INTENT(in) :: kt ! ocean time-step1180 CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag1181 WRITE(*,*) 'ts_rst : You should not have seen this print! error?', kt, cdrw1182 END SUBROUTINE ts_rst1183 SUBROUTINE dyn_spg_ts_init( kt ) ! Empty routine1184 INTEGER , INTENT(in) :: kt ! ocean time-step1185 WRITE(*,*) 'dyn_spg_ts_init : You should not have seen this print! error?', kt1186 END SUBROUTINE dyn_spg_ts_init1187 #endif1188 1189 1143 !!====================================================================== 1190 1144 END MODULE dynspg_ts -
branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DYN/dynvor.F90
r5901 r6079 32 32 USE trd_oce ! trends: ocean variables 33 33 USE trddyn ! trend manager: dynamics 34 USE c1d ! 1D vertical configuration 34 35 ! 35 36 USE lbclnk ! ocean lateral boundary conditions (or mpp link) … … 547 548 END SELECT 548 549 ! 550 IF( ln_dynvor_msk ) THEN !== mask/unmask vorticity ==! 551 DO jj = 1, jpjm1 552 DO ji = 1, fs_jpim1 ! vector opt. 553 zwz(ji,jj) = zwz(ji,jj) * fmask(ji,jj,jk) 554 END DO 555 END DO 556 ENDIF 557 ! 549 558 CALL lbc_lnk( zwz, 'F', 1. ) 550 !551 IF( ln_dynvor_msk ) THEN !== mask/unmask vorticity ==!552 DO jj = 1, jpjm1553 DO ji = 1, fs_jpim1 ! vector opt.554 zwz(ji,jj) = zwz(ji,jj) * fmask(ji,jj,jk)555 END DO556 END DO557 ENDIF558 559 ! 559 560 ! !== horizontal fluxes ==! … … 662 663 IF( ln_dynvor_een ) THEN ; ioptio = ioptio + 1 ; nvor_scheme = np_EEN ; ENDIF 663 664 ! 664 IF( ioptio /= 1) CALL ctl_stop( ' use ONE and ONLY one vorticity scheme' )665 IF( ( ioptio /= 1).AND.( .NOT.lk_c1d ) ) CALL ctl_stop( ' use ONE and ONLY one vorticity scheme' ) 665 666 ! 666 667 IF(lwp) WRITE(numout,*) ! type of calculated vorticity (set ncor, nrvm, ntot) -
branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf_exp.F90
r3625 r6079 8 8 !! NEMO 0.5 ! 2002-08 (G. Madec) F90: Free form and module 9 9 !! 3.3 ! 2010-04 (M. Leclair, G. Madec) Forcing averaged over 2 time steps 10 !! 3.7 ! 2015-11 (J. Chanut) output velocities instead of trends 10 11 !!---------------------------------------------------------------------- 11 12 … … 18 19 USE phycst ! physical constants 19 20 USE zdf_oce ! ocean vertical physics 21 USE dynadv, ONLY: ln_dynadv_vec ! Momentum advection form 20 22 USE sbc_oce ! surface boundary condition: ocean 21 23 USE lib_mpp ! MPP library … … 118 120 ! 119 121 END DO ! End of time splitting 122 123 ! Time step momentum here to be compliant with what is done in the implicit case 124 ! 125 IF( ln_dynadv_vec .OR. .NOT. lk_vvl ) THEN ! applied on velocity 126 DO jk = 1, jpkm1 127 ua(:,:,jk) = ( ub(:,:,jk) + p2dt * ua(:,:,jk) ) * umask(:,:,jk) 128 va(:,:,jk) = ( vb(:,:,jk) + p2dt * 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 & + p2dt * ua(:,:,jk) * fse3u_n(:,:,jk) ) & 134 & / fse3u_a(:,:,jk) * umask(:,:,jk) 135 va(:,:,jk) = ( vb(:,:,jk) * fse3v_b(:,:,jk) & 136 & + p2dt * va(:,:,jk) * fse3v_n(:,:,jk) ) & 137 & / fse3v_a(:,:,jk) * vmask(:,:,jk) 138 END DO 139 ENDIF 120 140 ! 121 141 CALL wrk_dealloc( jpi,jpj,jpk, zwx, zwy, zwz, zww ) -
branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf_imp.F90
r5901 r6079 25 25 USE wrk_nemo ! Memory Allocation 26 26 USE timing ! Timing 27 USE dynadv ! dynamics: vector invariant versus flux form 28 USE dynspg_oce, ONLY: lk_dynspg_ts 27 USE dynadv, ONLY: ln_dynadv_vec ! Momentum advection form 29 28 30 29 IMPLICIT NONE … … 87 86 ENDIF 88 87 89 ! 0. Local constant initialization 90 ! -------------------------------- 91 z1_p2dt = 1._wp / p2dt ! inverse of the timestep 92 93 ! 1. Apply semi-implicit bottom friction 88 ! 1. Time step dynamics 89 ! --------------------- 90 91 IF( ln_dynadv_vec .OR. .NOT. lk_vvl ) THEN ! applied on velocity 92 DO jk = 1, jpkm1 93 ua(:,:,jk) = ( ub(:,:,jk) + p2dt * ua(:,:,jk) ) * umask(:,:,jk) 94 va(:,:,jk) = ( vb(:,:,jk) + p2dt * va(:,:,jk) ) * vmask(:,:,jk) 95 END DO 96 ELSE ! applied on thickness weighted velocity 97 DO jk = 1, jpkm1 98 ua(:,:,jk) = ( ub(:,:,jk) * fse3u_b(:,:,jk) & 99 & + p2dt * ua(:,:,jk) * fse3u_n(:,:,jk) ) & 100 & / fse3u_a(:,:,jk) * umask(:,:,jk) 101 va(:,:,jk) = ( vb(:,:,jk) * fse3v_b(:,:,jk) & 102 & + p2dt * va(:,:,jk) * fse3v_n(:,:,jk) ) & 103 & / fse3v_a(:,:,jk) * vmask(:,:,jk) 104 END DO 105 ENDIF 106 107 ! 2. Apply semi-implicit bottom friction 94 108 ! -------------------------------------- 95 109 ! Only needed for semi-implicit bottom friction setup. The explicit … … 97 111 ! column vector of the tri-diagonal matrix equation 98 112 ! 99 100 113 IF( ln_bfrimp ) THEN 101 114 DO jj = 2, jpjm1 … … 119 132 ENDIF 120 133 121 #if defined key_dynspg_ts 122 IF( ln_dynadv_vec .OR. .NOT. lk_vvl ) THEN ! applied on velocity 123 DO jk = 1, jpkm1 124 ua(:,:,jk) = ( ub(:,:,jk) + p2dt * ua(:,:,jk) ) * umask(:,:,jk) 125 va(:,:,jk) = ( vb(:,:,jk) + p2dt * va(:,:,jk) ) * vmask(:,:,jk) 126 END DO 127 ELSE ! applied on thickness weighted velocity 128 DO jk = 1, jpkm1 129 ua(:,:,jk) = ( ub(:,:,jk) * fse3u_b(:,:,jk) & 130 & + p2dt * ua(:,:,jk) * fse3u_n(:,:,jk) ) & 131 & / fse3u_a(:,:,jk) * umask(:,:,jk) 132 va(:,:,jk) = ( vb(:,:,jk) * fse3v_b(:,:,jk) & 133 & + p2dt * va(:,:,jk) * fse3v_n(:,:,jk) ) & 134 & / fse3v_a(:,:,jk) * vmask(:,:,jk) 135 END DO 136 ENDIF 137 138 IF ( ln_bfrimp .AND.lk_dynspg_ts ) THEN 134 ! With split-explicit free surface, barotropic stress is treated explicitly 135 ! Update velocities at the bottom. 136 ! J. Chanut: The bottom stress is computed considering after barotropic velocities, which does 137 ! not lead to the effective stress seen over the whole barotropic loop. 138 IF ( ln_bfrimp .AND.ln_dynspg_ts ) THEN 139 139 ! remove barotropic velocities: 140 140 DO jk = 1, jpkm1 … … 166 166 END IF 167 167 ENDIF 168 #endif 169 170 ! 2. Vertical diffusion on u 168 169 ! 3. Vertical diffusion on u 171 170 ! --------------------------- 172 171 ! Matrix and second member construction … … 219 218 DO jj = 2, jpjm1 !== second recurrence: SOLk = RHSk - Lk / Dk-1 Lk-1 == 220 219 DO ji = fs_2, fs_jpim1 ! vector opt. 221 #if defined key_dynspg_ts222 220 ze3ua = ( 1._wp - r_vvl ) * fse3u_n(ji,jj,1) + r_vvl * fse3u_a(ji,jj,1) 223 221 ua(ji,jj,1) = ua(ji,jj,1) + p2dt * 0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) ) & 224 222 & / ( ze3ua * rau0 ) * umask(ji,jj,1) 225 #else226 ua(ji,jj,1) = ub(ji,jj,1) &227 & + p2dt *(ua(ji,jj,1) + 0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) ) &228 & / ( fse3u(ji,jj,1) * rau0 ) * umask(ji,jj,1) )229 #endif230 223 END DO 231 224 END DO … … 233 226 DO jj = 2, jpjm1 234 227 DO ji = fs_2, fs_jpim1 235 #if defined key_dynspg_ts236 228 zrhs = ua(ji,jj,jk) ! zrhs=right hand side 237 #else238 zrhs = ub(ji,jj,jk) + p2dt * ua(ji,jj,jk)239 #endif240 229 ua(ji,jj,jk) = zrhs - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * ua(ji,jj,jk-1) 241 230 END DO … … 256 245 END DO 257 246 258 #if ! defined key_dynspg_ts 259 ! Normalization to obtain the general momentum trend ua 260 DO jk = 1, jpkm1 261 DO jj = 2, jpjm1 262 DO ji = fs_2, fs_jpim1 ! vector opt. 263 ua(ji,jj,jk) = ( ua(ji,jj,jk) - ub(ji,jj,jk) ) * z1_p2dt 264 END DO 265 END DO 266 END DO 267 #endif 268 269 ! 3. Vertical diffusion on v 247 ! 4. Vertical diffusion on v 270 248 ! --------------------------- 271 249 ! Matrix and second member construction … … 317 295 ! 318 296 DO jj = 2, jpjm1 !== second recurrence: SOLk = RHSk - Lk / Dk-1 Lk-1 == 319 DO ji = fs_2, fs_jpim1 ! vector opt. 320 #if defined key_dynspg_ts 297 DO ji = fs_2, fs_jpim1 ! vector opt. 321 298 ze3va = ( 1._wp - r_vvl ) * fse3v_n(ji,jj,1) + r_vvl * fse3v_a(ji,jj,1) 322 299 va(ji,jj,1) = va(ji,jj,1) + p2dt * 0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) ) & 323 300 & / ( ze3va * rau0 ) 324 #else325 va(ji,jj,1) = vb(ji,jj,1) &326 & + p2dt *(va(ji,jj,1) + 0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) ) &327 & / ( fse3v(ji,jj,1) * rau0 ) )328 #endif329 301 END DO 330 302 END DO … … 332 304 DO jj = 2, jpjm1 333 305 DO ji = fs_2, fs_jpim1 ! vector opt. 334 #if defined key_dynspg_ts335 306 zrhs = va(ji,jj,jk) ! zrhs=right hand side 336 #else337 zrhs = vb(ji,jj,jk) + p2dt * va(ji,jj,jk)338 #endif339 307 va(ji,jj,jk) = zrhs - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * va(ji,jj,jk-1) 340 308 END DO … … 355 323 END DO 356 324 357 ! Normalization to obtain the general momentum trend va358 #if ! defined key_dynspg_ts359 DO jk = 1, jpkm1360 DO jj = 2, jpjm1361 DO ji = fs_2, fs_jpim1 ! vector opt.362 va(ji,jj,jk) = ( va(ji,jj,jk) - vb(ji,jj,jk) ) * z1_p2dt363 END DO364 END DO365 END DO366 #endif367 368 325 ! J. Chanut: Lines below are useless ? 369 !! restore bottom layer avmu(v)326 !! restore avmu(v)=0. at bottom (and top if ln_isfcav=T interfaces) 370 327 IF( ln_bfrimp ) THEN 371 328 DO jj = 2, jpjm1 -
branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90
r5901 r6079 106 106 ssha(:,:) = ( sshb(:,:) - z2dt * ( zcoef * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) ) ) * ssmask(:,:) 107 107 108 #if ! defined key_dynspg_ts 109 ! These lines are not necessary with time splitting since110 ! boundary condition on sea level is set during ts loop108 IF ( .NOT.ln_dynspg_ts ) THEN 109 ! These lines are not necessary with time splitting since 110 ! boundary condition on sea level is set during ts loop 111 111 # if defined key_agrif 112 CALL agrif_ssh( kt )112 CALL agrif_ssh( kt ) 113 113 # endif 114 114 # if defined key_bdy 115 IF( lk_bdy ) THEN116 CALL lbc_lnk( ssha, 'T', 1. ) ! Not sure that's necessary117 CALL bdy_ssh( ssha ) ! Duplicate sea level across open boundaries118 ENDIF115 IF( lk_bdy ) THEN 116 CALL lbc_lnk( ssha, 'T', 1. ) ! Not sure that's necessary 117 CALL bdy_ssh( ssha ) ! Duplicate sea level across open boundaries 118 ENDIF 119 119 # endif 120 #endif 120 ENDIF 121 121 122 122 #if defined key_asminc … … 250 250 ENDIF 251 251 252 # if defined key_dynspg_ts 253 IF( ( neuler == 0 .AND. kt == nit000 ) .OR. ln_bt_fw ) THEN !** Euler time-stepping: no filter 254 # else 255 IF ( neuler == 0 .AND. kt == nit000 ) THEN !** Euler time-stepping at first time-step : no filter 256 #endif 252 IF( ( neuler == 0 .AND. kt == nit000 ) .OR. ( ln_bt_fw .AND. ln_dynspg_ts ) ) THEN 253 !** Euler time-stepping: no filter 257 254 sshb(:,:) = sshn(:,:) ! before <-- now 258 255 sshn(:,:) = ssha(:,:) ! now <-- after (before already = now)
Note: See TracChangeset
for help on using the changeset viewer.
