- Timestamp:
- 2015-11-30T12:48:01+01:00 (9 years ago)
- Location:
- branches/2014/dev_r4650_UKMO14.4_OBS_GENERAL_VINTERP/NEMOGCM/NEMO/OPA_SRC/DYN
- Files:
-
- 7 deleted
- 16 edited
- 2 copied
Legend:
- Unmodified
- Added
- Removed
-
branches/2014/dev_r4650_UKMO14.4_OBS_GENERAL_VINTERP/NEMOGCM/NEMO/OPA_SRC/DYN/dynadv.F90
r5950 r5951 76 76 CASE ( 3 ) 77 77 CALL dyn_adv_ubs ( kt ) ! 3rd order UBS scheme 78 !79 CASE (-1 ) ! esopa: test all possibility with control print80 CALL dyn_keg ( kt, nn_dynkeg )81 CALL dyn_zad ( kt )82 CALL dyn_adv_cen2( kt )83 CALL dyn_adv_ubs ( kt )84 78 END SELECT 85 79 ! … … 126 120 IF( ln_dynadv_cen2 ) ioptio = ioptio + 1 127 121 IF( ln_dynadv_ubs ) ioptio = ioptio + 1 128 IF( lk_esopa ) ioptio = 1129 122 130 123 IF( ioptio /= 1 ) CALL ctl_stop( 'Choose ONE advection scheme in namelist namdyn_adv' ) … … 139 132 IF( ln_dynadv_cen2 ) nadv = 2 140 133 IF( ln_dynadv_ubs ) nadv = 3 141 IF( lk_esopa ) nadv = -1142 134 143 135 IF(lwp) THEN ! Print the choice … … 151 143 IF( nadv == 2 ) WRITE(numout,*) ' flux form : 2nd order scheme is used' 152 144 IF( nadv == 3 ) WRITE(numout,*) ' flux form : UBS scheme is used' 153 IF( nadv == -1 ) WRITE(numout,*) ' esopa test: use all advection formulation'154 145 ENDIF 155 146 ! -
branches/2014/dev_r4650_UKMO14.4_OBS_GENERAL_VINTERP/NEMOGCM/NEMO/OPA_SRC/DYN/dynadv_cen2.F90
r5950 r5951 90 90 DO jj = 2, jpjm1 ! divergence of horizontal momentum fluxes 91 91 DO ji = fs_2, fs_jpim1 ! vector opt. 92 zbu = e1 u(ji,jj) *e2u(ji,jj) * fse3u(ji,jj,jk)93 zbv = e1 v(ji,jj) *e2v(ji,jj) * fse3v(ji,jj,jk)92 zbu = e1e2u(ji,jj) * fse3u(ji,jj,jk) 93 zbv = e1e2v(ji,jj) * fse3v(ji,jj,jk) 94 94 ! 95 95 ua(ji,jj,jk) = ua(ji,jj,jk) - ( zfu_t(ji+1,jj ,jk) - zfu_t(ji ,jj ,jk) & … … 114 114 DO jk = 1, jpkm1 ! ==================== ! 115 115 ! ! Vertical volume fluxesÊ 116 zfw(:,:,jk) = 0.25 * e1 t(:,:) *e2t(:,:) * wn(:,:,jk)116 zfw(:,:,jk) = 0.25 * e1e2t(:,:) * wn(:,:,jk) 117 117 ! 118 118 IF( jk == 1 ) THEN ! surface/bottom advective fluxes … … 144 144 DO ji = fs_2, fs_jpim1 ! vector opt. 145 145 ua(ji,jj,jk) = ua(ji,jj,jk) - ( zfu_uw(ji,jj,jk) - zfu_uw(ji,jj,jk+1) ) & 146 & / ( e1 u(ji,jj) *e2u(ji,jj) * fse3u(ji,jj,jk) )146 & / ( e1e2u(ji,jj) * fse3u(ji,jj,jk) ) 147 147 va(ji,jj,jk) = va(ji,jj,jk) - ( zfv_vw(ji,jj,jk) - zfv_vw(ji,jj,jk+1) ) & 148 & / ( e1 v(ji,jj) *e2v(ji,jj) * fse3v(ji,jj,jk) )148 & / ( e1e2v(ji,jj) * fse3v(ji,jj,jk) ) 149 149 END DO 150 150 END DO -
branches/2014/dev_r4650_UKMO14.4_OBS_GENERAL_VINTERP/NEMOGCM/NEMO/OPA_SRC/DYN/dynadv_ubs.F90
r5950 r5951 181 181 DO jj = 2, jpjm1 ! divergence of horizontal momentum fluxes 182 182 DO ji = fs_2, fs_jpim1 ! vector opt. 183 zbu = e1 u(ji,jj) *e2u(ji,jj) * fse3u(ji,jj,jk)184 zbv = e1 v(ji,jj) *e2v(ji,jj) * fse3v(ji,jj,jk)183 zbu = e1e2u(ji,jj) * fse3u(ji,jj,jk) 184 zbv = e1e2v(ji,jj) * fse3v(ji,jj,jk) 185 185 ! 186 186 ua(ji,jj,jk) = ua(ji,jj,jk) - ( zfu_t(ji+1,jj ,jk) - zfu_t(ji ,jj ,jk) & … … 203 203 DO jk = 1, jpkm1 ! ==================== ! 204 204 ! ! Vertical volume fluxesÊ 205 zfw(:,:,jk) = 0.25 * e1 t(:,:) *e2t(:,:) * wn(:,:,jk)205 zfw(:,:,jk) = 0.25 * e1e2t(:,:) * wn(:,:,jk) 206 206 ! 207 207 IF( jk == 1 ) THEN ! surface/bottom advective fluxes … … 233 233 DO ji = fs_2, fs_jpim1 ! vector opt. 234 234 ua(ji,jj,jk) = ua(ji,jj,jk) - ( zfu_uw(ji,jj,jk) - zfu_uw(ji,jj,jk+1) ) & 235 & / ( e1 u(ji,jj) *e2u(ji,jj) * fse3u(ji,jj,jk) )235 & / ( e1e2u(ji,jj) * fse3u(ji,jj,jk) ) 236 236 va(ji,jj,jk) = va(ji,jj,jk) - ( zfv_vw(ji,jj,jk) - zfv_vw(ji,jj,jk+1) ) & 237 & / ( e1 v(ji,jj) *e2v(ji,jj) * fse3v(ji,jj,jk) )237 & / ( e1e2v(ji,jj) * fse3v(ji,jj,jk) ) 238 238 END DO 239 239 END DO -
branches/2014/dev_r4650_UKMO14.4_OBS_GENERAL_VINTERP/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90
r5950 r5951 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 !================================================================================== … … 1046 1052 DO jj = 2, jpjm1 1047 1053 DO ji = 2, jpim1 1048 zsshu_n(ji,jj) = (e1 2u(ji,jj) * sshn(ji,jj) + e12u(ji+1, jj) * sshn(ji+1,jj)) * &1049 & r1_e1 2u(ji,jj) * umask(ji,jj,1) * 0.5_wp1050 zsshv_n(ji,jj) = (e1 2v(ji,jj) * sshn(ji,jj) + e12v(ji+1, jj) * sshn(ji,jj+1)) * &1051 & r1_e1 2v(ji,jj) * vmask(ji,jj,1) * 0.5_wp1054 zsshu_n(ji,jj) = (e1e2u(ji,jj) * sshn(ji,jj) + e1e2u(ji+1, jj) * sshn(ji+1,jj)) * & 1055 & r1_e1e2u(ji,jj) * umask(ji,jj,1) * 0.5_wp 1056 zsshv_n(ji,jj) = (e1e2v(ji,jj) * sshn(ji,jj) + e1e2v(ji+1, jj) * sshn(ji,jj+1)) * & 1057 & r1_e1e2v(ji,jj) * vmask(ji,jj,1) * 0.5_wp 1052 1058 END DO 1053 1059 END DO -
branches/2014/dev_r4650_UKMO14.4_OBS_GENERAL_VINTERP/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf.F90
r5950 r5951 4 4 !! Ocean physics: lateral diffusivity trends 5 5 !!===================================================================== 6 !! History : 9.0 ! 05-11 (G. Madec) Original code (new step architecture) 6 !! History : 2.0 ! 2005-11 (G. Madec) Original code (new step architecture) 7 !! 3.7 ! 2014-01 (F. Lemarie, G. Madec) restructuration/simplification of ahm specification, 8 !! ! add velocity dependent coefficient and optional read in file 7 9 !!---------------------------------------------------------------------- 8 10 … … 14 16 USE dom_oce ! ocean space and time domain 15 17 USE phycst ! physical constants 16 USE ldfdyn_oce ! ocean dynamics lateral physics 17 USE ldftra_oce ! ocean tracers lateral physics 18 USE ldfslp ! lateral mixing: slopes of mixing orientation 19 USE dynldf_bilapg ! lateral mixing (dyn_ldf_bilapg routine) 20 USE dynldf_bilap ! lateral mixing (dyn_ldf_bilap routine) 21 USE dynldf_iso ! lateral mixing (dyn_ldf_iso routine) 22 USE dynldf_lap ! lateral mixing (dyn_ldf_lap routine) 18 USE ldfdyn ! lateral diffusion: eddy viscosity coef. 19 USE ldfslp ! lateral diffusion: slopes of mixing orientation 20 USE dynldf_lap_blp ! lateral mixing (dyn_ldf_lap & dyn_ldf_blp routines) 21 USE dynldf_iso ! lateral mixing (dyn_ldf_iso routine ) 23 22 USE trd_oce ! trends: ocean variables 24 USE trddyn ! trend manager: dynamics (trd_dyn 23 USE trddyn ! trend manager: dynamics (trd_dyn routine) 25 24 ! 26 25 USE prtctl ! Print control … … 28 27 USE lib_mpp ! distribued memory computing library 29 28 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 30 USE wrk_nemo 31 USE timing 29 USE wrk_nemo ! Memory Allocation 30 USE timing ! Timing 32 31 33 32 IMPLICIT NONE … … 37 36 PUBLIC dyn_ldf_init ! called by opa module 38 37 39 INTEGER :: nldf = -2 ! type of lateral diffusion used defined from ln_dynldf_... namlist logicals) 38 ! ! Flag to control the type of lateral viscous operator 39 INTEGER, PARAMETER, PUBLIC :: np_ERROR =-10 ! error in setting the operator 40 INTEGER, PARAMETER, PUBLIC :: np_no_ldf = 00 ! without operator (i.e. no lateral viscous trend) 41 ! !! laplacian ! bilaplacian ! 42 INTEGER, PARAMETER, PUBLIC :: np_lap = 10 , np_blp = 20 ! iso-level operator 43 INTEGER, PARAMETER, PUBLIC :: np_lap_i = 11 ! iso-neutral or geopotential operator 44 45 INTEGER :: nldf ! type of lateral diffusion used defined from ln_dynldf_... (namlist logicals) 40 46 41 47 !! * Substitutions … … 43 49 # include "vectopt_loop_substitute.h90" 44 50 !!---------------------------------------------------------------------- 45 !! NEMO/OPA 3. 3 , NEMO Consortium (2010)51 !! NEMO/OPA 3.7 , NEMO Consortium (2015) 46 52 !! $Id$ 47 53 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) … … 62 68 IF( nn_timing == 1 ) CALL timing_start('dyn_ldf') 63 69 ! 64 IF( l_trddyn ) THEN ! temporary save of ta and satrends65 CALL wrk_alloc( jpi, jpj, jpk,ztrdu, ztrdv )70 IF( l_trddyn ) THEN ! temporary save of momentum trends 71 CALL wrk_alloc( jpi,jpj,jpk, ztrdu, ztrdv ) 66 72 ztrdu(:,:,:) = ua(:,:,:) 67 73 ztrdv(:,:,:) = va(:,:,:) … … 70 76 SELECT CASE ( nldf ) ! compute lateral mixing trend and add it to the general trend 71 77 ! 72 CASE ( 0 ) ; CALL dyn_ldf_lap ( kt ) ! iso-level laplacian 73 CASE ( 1 ) ; CALL dyn_ldf_iso ( kt ) ! rotated laplacian (except dk[ dk[.] ] part) 74 CASE ( 2 ) ; CALL dyn_ldf_bilap ( kt ) ! iso-level bilaplacian 75 CASE ( 3 ) ; CALL dyn_ldf_bilapg ( kt ) ! s-coord. horizontal bilaplacian 76 CASE ( 4 ) ! iso-level laplacian + bilaplacian 77 CALL dyn_ldf_lap ( kt ) 78 CALL dyn_ldf_bilap ( kt ) 79 CASE ( 5 ) ! rotated laplacian + bilaplacian (s-coord) 80 CALL dyn_ldf_iso ( kt ) 81 CALL dyn_ldf_bilapg ( kt ) 78 CASE ( np_lap ) ; CALL dyn_ldf_lap ( kt, ub, vb, ua, va, 1 ) ! iso-level laplacian 79 CASE ( np_lap_i ) ; CALL dyn_ldf_iso ( kt ) ! rotated laplacian 80 CASE ( np_blp ) ; CALL dyn_ldf_blp ( kt, ub, vb, ua, va ) ! iso-level bi-laplacian 82 81 ! 83 CASE ( -1 ) ! esopa: test all possibility with control print84 CALL dyn_ldf_lap ( kt )85 CALL prt_ctl( tab3d_1=ua, clinfo1=' ldf0 - Ua: ', mask1=umask, &86 & tab3d_2=va, clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' )87 CALL dyn_ldf_iso ( kt )88 CALL prt_ctl( tab3d_1=ua, clinfo1=' ldf1 - Ua: ', mask1=umask, &89 & tab3d_2=va, clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' )90 CALL dyn_ldf_bilap ( kt )91 CALL prt_ctl( tab3d_1=ua, clinfo1=' ldf2 - Ua: ', mask1=umask, &92 & tab3d_2=va, clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' )93 CALL dyn_ldf_bilapg ( kt )94 CALL prt_ctl( tab3d_1=ua, clinfo1=' ldf3 - Ua: ', mask1=umask, &95 & tab3d_2=va, clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' )96 !97 CASE ( -2 ) ! neither laplacian nor bilaplacian schemes used98 IF( kt == nit000 ) THEN99 IF(lwp) WRITE(numout,*)100 IF(lwp) WRITE(numout,*) 'dyn_ldf : no lateral diffusion on momentum setup'101 IF(lwp) WRITE(numout,*) '~~~~~~~ '102 ENDIF103 82 END SELECT 104 83 … … 107 86 ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 108 87 CALL trd_dyn( ztrdu, ztrdv, jpdyn_ldf, kt ) 109 CALL wrk_dealloc( jpi, jpj, jpk,ztrdu, ztrdv )88 CALL wrk_dealloc( jpi,jpj,jpk, ztrdu, ztrdv ) 110 89 ENDIF 111 90 ! ! print sum trends (used for debugging) … … 126 105 INTEGER :: ioptio, ierr ! temporary integers 127 106 !!---------------------------------------------------------------------- 128 107 ! 129 108 ! ! Namelist nam_dynldf: already read in ldfdyn module 130 109 ! 131 110 IF(lwp) THEN ! Namelist print 132 111 WRITE(numout,*) … … 134 113 WRITE(numout,*) '~~~~~~~~~~~' 135 114 WRITE(numout,*) ' Namelist nam_dynldf : set lateral mixing parameters (type, direction, coefficients)' 136 WRITE(numout,*) ' laplacian operator ln_dynldf_lap 137 WRITE(numout,*) ' bilaplacian operator ln_dynldf_b ilap = ', ln_dynldf_bilap138 WRITE(numout,*) ' iso-level ln_dynldf_lev el = ', ln_dynldf_level139 WRITE(numout,*) ' horizontal (geopotential) ln_dynldf_hor 140 WRITE(numout,*) ' iso-neutral ln_dynldf_iso 115 WRITE(numout,*) ' laplacian operator ln_dynldf_lap = ', ln_dynldf_lap 116 WRITE(numout,*) ' bilaplacian operator ln_dynldf_blp = ', ln_dynldf_blp 117 WRITE(numout,*) ' iso-level ln_dynldf_lev = ', ln_dynldf_lev 118 WRITE(numout,*) ' horizontal (geopotential) ln_dynldf_hor = ', ln_dynldf_hor 119 WRITE(numout,*) ' iso-neutral ln_dynldf_iso = ', ln_dynldf_iso 141 120 ENDIF 142 143 ! ! control the consistency121 ! ! use of lateral operator or not 122 nldf = np_ERROR 144 123 ioptio = 0 145 IF( ln_dynldf_lap ) ioptio = ioptio + 1 146 IF( ln_dynldf_bilap ) ioptio = ioptio + 1 147 IF( ioptio < 1 ) CALL ctl_warn( ' neither laplacian nor bilaplacian operator set for dynamics' ) 148 ioptio = 0 149 IF( ln_dynldf_level ) ioptio = ioptio + 1 150 IF( ln_dynldf_hor ) ioptio = ioptio + 1 151 IF( ln_dynldf_iso ) ioptio = ioptio + 1 152 IF( ioptio > 1 ) CALL ctl_stop( ' use only ONE direction (level/hor/iso)' ) 153 154 IF( ln_dynldf_iso .AND. ln_traldf_hor ) CALL ctl_stop & 155 & ( 'Not sensible to use geopotential diffusion for tracers with isoneutral diffusion for dynamics' ) 156 157 ! ! Set nldf, the type of lateral diffusion, from ln_dynldf_... logicals 158 ierr = 0 159 IF ( ln_dynldf_lap ) THEN ! laplacian operator 160 IF ( ln_zco ) THEN ! z-coordinate 161 IF ( ln_dynldf_level ) nldf = 0 ! iso-level (no rotation) 162 IF ( ln_dynldf_hor ) nldf = 0 ! horizontal (no rotation) 163 IF ( ln_dynldf_iso ) nldf = 1 ! isoneutral ( rotation) 124 IF( ln_dynldf_lap ) ioptio = ioptio + 1 125 IF( ln_dynldf_blp ) ioptio = ioptio + 1 126 IF( ioptio > 1 ) CALL ctl_stop( 'dyn_ldf_init: use ONE or NONE of the 2 lap/bilap operator type on momentum' ) 127 IF( ioptio == 0 ) nldf = np_no_ldf ! No lateral mixing operator 128 ! 129 IF( nldf /= np_no_ldf ) THEN ! direction ==>> type of operator 130 ioptio = 0 131 IF( ln_dynldf_lev ) ioptio = ioptio + 1 132 IF( ln_dynldf_hor ) ioptio = ioptio + 1 133 IF( ln_dynldf_iso ) ioptio = ioptio + 1 134 IF( ioptio > 1 ) CALL ctl_stop( ' use only ONE direction (level/hor/iso)' ) 135 IF( ioptio == 0 ) CALL ctl_stop( ' use at least ONE direction (level/hor/iso)' ) 136 ! 137 ! ! Set nldf, the type of lateral diffusion, from ln_dynldf_... logicals 138 ierr = 0 139 IF ( ln_dynldf_lap ) THEN ! laplacian operator 140 IF ( ln_zco ) THEN ! z-coordinate 141 IF ( ln_dynldf_lev ) nldf = np_lap ! iso-level = horizontal (no rotation) 142 IF ( ln_dynldf_hor ) nldf = np_lap ! iso-level = horizontal (no rotation) 143 IF ( ln_dynldf_iso ) nldf = np_lap_i ! iso-neutral ( rotation) 144 ENDIF 145 IF ( ln_zps ) THEN ! z-coordinate with partial step 146 IF ( ln_dynldf_lev ) nldf = np_lap ! iso-level (no rotation) 147 IF ( ln_dynldf_hor ) nldf = np_lap ! iso-level (no rotation) 148 IF ( ln_dynldf_iso ) nldf = np_lap_i ! iso-neutral ( rotation) 149 ENDIF 150 IF ( ln_sco ) THEN ! s-coordinate 151 IF ( ln_dynldf_lev ) nldf = np_lap ! iso-level = horizontal (no rotation) 152 IF ( ln_dynldf_hor ) nldf = np_lap_i ! horizontal ( rotation) 153 IF ( ln_dynldf_iso ) nldf = np_lap_i ! iso-neutral ( rotation) 154 ENDIF 164 155 ENDIF 165 IF ( ln_zps ) THEN ! z-coordinate 166 IF ( ln_dynldf_level ) ierr = 1 ! iso-level not allowed 167 IF ( ln_dynldf_hor ) nldf = 0 ! horizontal (no rotation) 168 IF ( ln_dynldf_iso ) nldf = 1 ! isoneutral ( rotation) 156 ! 157 IF( ln_dynldf_blp ) THEN ! bilaplacian operator 158 IF ( ln_zco ) THEN ! z-coordinate 159 IF ( ln_dynldf_lev ) nldf = np_blp ! iso-level = horizontal (no rotation) 160 IF ( ln_dynldf_hor ) nldf = np_blp ! iso-level = horizontal (no rotation) 161 IF ( ln_dynldf_iso ) ierr = 2 ! iso-neutral ( rotation) 162 ENDIF 163 IF ( ln_zps ) THEN ! z-coordinate with partial step 164 IF ( ln_dynldf_lev ) nldf = np_blp ! iso-level (no rotation) 165 IF ( ln_dynldf_hor ) nldf = np_blp ! iso-level (no rotation) 166 IF ( ln_dynldf_iso ) ierr = 2 ! iso-neutral ( rotation) 167 ENDIF 168 IF ( ln_sco ) THEN ! s-coordinate 169 IF ( ln_dynldf_lev ) nldf = np_blp ! iso-level (no rotation) 170 IF ( ln_dynldf_hor ) ierr = 2 ! horizontal ( rotation) 171 IF ( ln_dynldf_iso ) ierr = 2 ! iso-neutral ( rotation) 172 ENDIF 169 173 ENDIF 170 IF ( ln_sco ) THEN ! s-coordinate 171 IF ( ln_dynldf_level ) nldf = 0 ! iso-level (no rotation) 172 IF ( ln_dynldf_hor ) nldf = 1 ! horizontal ( rotation) 173 IF ( ln_dynldf_iso ) nldf = 1 ! isoneutral ( rotation) 174 ENDIF 175 ENDIF 176 177 IF( ln_dynldf_bilap ) THEN ! bilaplacian operator 178 IF ( ln_zco ) THEN ! z-coordinate 179 IF ( ln_dynldf_level ) nldf = 2 ! iso-level (no rotation) 180 IF ( ln_dynldf_hor ) nldf = 2 ! horizontal (no rotation) 181 IF ( ln_dynldf_iso ) ierr = 2 ! isoneutral ( rotation) 182 ENDIF 183 IF ( ln_zps ) THEN ! z-coordinate 184 IF ( ln_dynldf_level ) ierr = 1 ! iso-level not allowed 185 IF ( ln_dynldf_hor ) nldf = 2 ! horizontal (no rotation) 186 IF ( ln_dynldf_iso ) ierr = 2 ! isoneutral ( rotation) 187 ENDIF 188 IF ( ln_sco ) THEN ! s-coordinate 189 IF ( ln_dynldf_level ) nldf = 2 ! iso-level (no rotation) 190 IF ( ln_dynldf_hor ) nldf = 3 ! horizontal ( rotation) 191 IF ( ln_dynldf_iso ) ierr = 2 ! isoneutral ( rotation) 192 ENDIF 193 ENDIF 194 195 IF( ln_dynldf_lap .AND. ln_dynldf_bilap ) THEN ! mixed laplacian and bilaplacian operators 196 IF ( ln_zco ) THEN ! z-coordinate 197 IF ( ln_dynldf_level ) nldf = 4 ! iso-level (no rotation) 198 IF ( ln_dynldf_hor ) nldf = 4 ! horizontal (no rotation) 199 IF ( ln_dynldf_iso ) ierr = 2 ! isoneutral ( rotation) 200 ENDIF 201 IF ( ln_zps ) THEN ! z-coordinate 202 IF ( ln_dynldf_level ) ierr = 1 ! iso-level not allowed 203 IF ( ln_dynldf_hor ) nldf = 4 ! horizontal (no rotation) 204 IF ( ln_dynldf_iso ) ierr = 2 ! isoneutral ( rotation) 205 ENDIF 206 IF ( ln_sco ) THEN ! s-coordinate 207 IF ( ln_dynldf_level ) nldf = 4 ! iso-level (no rotation) 208 IF ( ln_dynldf_hor ) nldf = 5 ! horizontal ( rotation) 209 IF ( ln_dynldf_iso ) ierr = 2 ! isoneutral ( rotation) 210 ENDIF 211 ENDIF 212 213 IF( lk_esopa ) nldf = -1 ! esopa test 214 215 IF( ierr == 1 ) CALL ctl_stop( 'iso-level in z-coordinate - partial step, not allowed' ) 216 IF( ierr == 2 ) CALL ctl_stop( 'isoneutral bilaplacian operator does not exist' ) 217 IF( nldf == 1 .OR. nldf == 3 ) THEN ! rotation 218 IF( .NOT.lk_ldfslp ) CALL ctl_stop( 'the rotation of the diffusive tensor require key_ldfslp' ) 174 ! 175 IF( ierr == 2 ) CALL ctl_stop( 'rotated bi-laplacian operator does not exist' ) 176 ! 177 IF( nldf == np_lap_i ) l_ldfslp = .TRUE. ! rotation require the computation of the slopes 178 ! 219 179 ENDIF 220 180 221 181 IF(lwp) THEN 222 182 WRITE(numout,*) 223 IF( nldf == -2 ) WRITE(numout,*) ' neither laplacian nor bilaplacian schemes used' 224 IF( nldf == -1 ) WRITE(numout,*) ' ESOPA test All scheme used' 225 IF( nldf == 0 ) WRITE(numout,*) ' laplacian operator' 226 IF( nldf == 1 ) WRITE(numout,*) ' rotated laplacian operator' 227 IF( nldf == 2 ) WRITE(numout,*) ' bilaplacian operator' 228 IF( nldf == 3 ) WRITE(numout,*) ' rotated bilaplacian' 229 IF( nldf == 4 ) WRITE(numout,*) ' laplacian and bilaplacian operators' 230 IF( nldf == 5 ) WRITE(numout,*) ' rotated laplacian and bilaplacian operators' 183 IF( nldf == np_no_ldf ) WRITE(numout,*) ' NO lateral viscosity' 184 IF( nldf == np_lap ) WRITE(numout,*) ' iso-level laplacian operator' 185 IF( nldf == np_lap_i ) WRITE(numout,*) ' rotated laplacian operator with iso-level background' 186 IF( nldf == np_blp ) WRITE(numout,*) ' iso-level bi-laplacian operator' 231 187 ENDIF 232 188 ! -
branches/2014/dev_r4650_UKMO14.4_OBS_GENERAL_VINTERP/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf_iso.F90
r5950 r5951 2 2 !!====================================================================== 3 3 !! *** MODULE dynldf_iso *** 4 !! Ocean dynamics: lateral viscosity trend4 !! Ocean dynamics: lateral viscosity trend (rotated laplacian operator) 5 5 !!====================================================================== 6 6 !! History : OPA ! 97-07 (G. Madec) Original code … … 8 8 !! - ! 2004-08 (C. Talandier) New trends organization 9 9 !! 2.0 ! 2005-11 (G. Madec) s-coordinate: horizontal diffusion 10 !! 3.7 ! 2014-01 (F. Lemarie, G. Madec) restructuration/simplification of ahm specification, 11 !! ! add velocity dependent coefficient and optional read in file 10 12 !!---------------------------------------------------------------------- 11 #if defined key_ldfslp || defined key_esopa 12 !!---------------------------------------------------------------------- 13 !! 'key_ldfslp' slopes of the direction of mixing 13 14 14 !!---------------------------------------------------------------------- 15 15 !! dyn_ldf_iso : update the momentum trend with the horizontal part … … 19 19 USE oce ! ocean dynamics and tracers 20 20 USE dom_oce ! ocean space and time domain 21 USE ldfdyn _oce ! ocean dynamics lateral physics22 USE ldftra _oce ! ocean tracer lateral physics21 USE ldfdyn ! lateral diffusion: eddy viscosity coef. 22 USE ldftra ! lateral physics: eddy diffusivity 23 23 USE zdf_oce ! ocean vertical physics 24 24 USE ldfslp ! iso-neutral slopes 25 25 ! 26 USE lbclnk ! ocean lateral boundary conditions (or mpp link)27 26 USE in_out_manager ! I/O manager 28 27 USE lib_mpp ! MPP library 28 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 29 29 USE prtctl ! Print control 30 30 USE wrk_nemo ! Memory Allocation … … 42 42 !! * Substitutions 43 43 # include "domzgr_substitute.h90" 44 # include "ldfdyn_substitute.h90"45 44 # include "vectopt_loop_substitute.h90" 46 45 !!---------------------------------------------------------------------- … … 83 82 !! horizontal fluxes associated with the rotated lateral mixing: 84 83 !! u-component: 85 !! ziut = ( ahmt + ahmb0) e2t * e3t / e1t di[ ub ]86 !! - ahmte2t * mi-1(uslp) dk[ mi(mk(ub)) ]87 !! zjuf = ( ahmf + ahmb0) e1f * e3f / e2f dj[ ub ]88 !! - ahmfe1f * mi(vslp) dk[ mj(mk(ub)) ]84 !! ziut = ( ahmt + rn_ahm_b ) e2t * e3t / e1t di[ ub ] 85 !! - ahmt e2t * mi-1(uslp) dk[ mi(mk(ub)) ] 86 !! zjuf = ( ahmf + rn_ahm_b ) e1f * e3f / e2f dj[ ub ] 87 !! - ahmf e1f * mi(vslp) dk[ mj(mk(ub)) ] 89 88 !! v-component: 90 !! zivf = ( ahmf + ahmb0) e2t * e3t / e1t di[ vb ]91 !! - ahmfe2t * mj(uslp) dk[ mi(mk(vb)) ]92 !! zjvt = ( ahmt + ahmb0) e1f * e3f / e2f dj[ ub ]93 !! - ahmte1f * mj-1(vslp) dk[ mj(mk(vb)) ]89 !! zivf = ( ahmf + rn_ahm_b ) e2t * e3t / e1t di[ vb ] 90 !! - ahmf e2t * mj(uslp) dk[ mi(mk(vb)) ] 91 !! zjvt = ( ahmt + rn_ahm_b ) e1f * e3f / e2f dj[ ub ] 92 !! - ahmt e1f * mj-1(vslp) dk[ mj(mk(vb)) ] 94 93 !! take the horizontal divergence of the fluxes: 95 94 !! diffu = 1/(e1u*e2u*e3u) { di [ ziut ] + dj-1[ zjuf ] } … … 106 105 !! of the rotated operator in dynzdf module 107 106 !!---------------------------------------------------------------------- 108 !109 107 INTEGER, INTENT( in ) :: kt ! ocean time-step index 110 108 ! 111 109 INTEGER :: ji, jj, jk ! dummy loop indices 112 110 REAL(wp) :: zabe1, zabe2, zcof1, zcof2 ! local scalars 113 REAL(wp) :: zmskt, zmskf , zbu, zbv, zuah, zvah! - -111 REAL(wp) :: zmskt, zmskf ! - - 114 112 REAL(wp) :: zcoef0, zcoef3, zcoef4, zmkt, zmkf ! - - 115 113 REAL(wp) :: zuav, zvav, zuwslpi, zuwslpj, zvwslpi, zvwslpj ! - - … … 130 128 ENDIF 131 129 132 ! s-coordinate: Iso-level diffusion on tracer, but geopotential level diffusion on momentum 130 !!gm bug is dyn_ldf_iso called before tra_ldf_iso .... <<<<<===== TO BE CHECKED 131 ! s-coordinate: Iso-level diffusion on momentum but not on tracer 133 132 IF( ln_dynldf_hor .AND. ln_traldf_iso ) THEN 134 133 ! … … 136 135 DO jj = 2, jpjm1 137 136 DO ji = 2, jpim1 138 uslp (ji,jj,jk) = - 1./e1u(ji,jj) * ( fsdept_b(ji+1,jj,jk) - fsdept_b(ji ,jj ,jk)) * umask(ji,jj,jk)139 vslp (ji,jj,jk) = - 1./e2v(ji,jj) * ( fsdept_b(ji,jj+1,jk) - fsdept_b(ji ,jj ,jk)) * vmask(ji,jj,jk)140 wslpi(ji,jj,jk) = - 1./e1t(ji,jj) * ( fsdepw_b(ji+1,jj,jk) - fsdepw_b(ji-1,jj,jk)) * tmask(ji,jj,jk) * 0.5141 wslpj(ji,jj,jk) = - 1./e2t(ji,jj) * ( fsdepw_b(ji,jj+1,jk) - fsdepw_b(ji,jj-1,jk)) * tmask(ji,jj,jk) * 0.5137 uslp (ji,jj,jk) = - ( fsdept_b(ji+1,jj,jk) - fsdept_b(ji ,jj ,jk) ) * r1_e1u(ji,jj) * umask(ji,jj,jk) 138 vslp (ji,jj,jk) = - ( fsdept_b(ji,jj+1,jk) - fsdept_b(ji ,jj ,jk) ) * r1_e2v(ji,jj) * vmask(ji,jj,jk) 139 wslpi(ji,jj,jk) = - ( fsdepw_b(ji+1,jj,jk) - fsdepw_b(ji-1,jj,jk) ) * r1_e1t(ji,jj) * tmask(ji,jj,jk) * 0.5 140 wslpj(ji,jj,jk) = - ( fsdepw_b(ji,jj+1,jk) - fsdepw_b(ji,jj-1,jk) ) * r1_e2t(ji,jj) * tmask(ji,jj,jk) * 0.5 142 141 END DO 143 142 END DO … … 184 183 DO jj = 2, jpjm1 185 184 DO ji = fs_2, jpi ! vector opt. 186 zabe1 = ( fsahmt(ji,jj,jk)+ahmb0) * e2t(ji,jj) * MIN( fse3u(ji,jj,jk), fse3u(ji-1,jj,jk) ) /e1t(ji,jj)187 188 zmskt = 1. /MAX( umask(ji-1,jj,jk )+umask(ji,jj,jk+1)&189 & + umask(ji-1,jj,jk+1)+umask(ji,jj,jk ), 1.)190 191 zcof1 = - aht0 * e2t(ji,jj) * zmskt * 0.5 * ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) )185 zabe1 = ( ahmt(ji,jj,jk)+rn_ahm_b ) * e2t(ji,jj) * MIN( fse3u(ji,jj,jk), fse3u(ji-1,jj,jk) ) * r1_e1t(ji,jj) 186 187 zmskt = 1._wp / MAX( umask(ji-1,jj,jk )+umask(ji,jj,jk+1) & 188 & + umask(ji-1,jj,jk+1)+umask(ji,jj,jk ) , 1._wp ) 189 190 zcof1 = - rn_aht_0 * e2t(ji,jj) * zmskt * 0.5 * ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) ) 192 191 192 ziut(ji,jj) = ( zabe1 * ( ub(ji,jj,jk) - ub(ji-1,jj,jk) ) & 193 & + zcof1 * ( zdku (ji,jj) + zdk1u(ji-1,jj) & 194 & +zdk1u(ji,jj) + zdku (ji-1,jj) ) ) * tmask(ji,jj,jk) 195 END DO 196 END DO 197 ELSE ! other coordinate system (zco or sco) : e3t 198 DO jj = 2, jpjm1 199 DO ji = fs_2, jpi ! vector opt. 200 zabe1 = ( ahmt(ji,jj,jk)+rn_ahm_b ) * e2t(ji,jj) * fse3t(ji,jj,jk) * r1_e1t(ji,jj) 201 202 zmskt = 1._wp / MAX( umask(ji-1,jj,jk ) + umask(ji,jj,jk+1) & 203 & + umask(ji-1,jj,jk+1) + umask(ji,jj,jk ) , 1._wp ) 204 205 zcof1 = - rn_aht_0 * e2t(ji,jj) * zmskt * 0.5 * ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) ) 206 193 207 ziut(ji,jj) = ( zabe1 * ( ub(ji,jj,jk) - ub(ji-1,jj,jk) ) & 194 208 & + zcof1 * ( zdku (ji,jj) + zdk1u(ji-1,jj) & … … 196 210 END DO 197 211 END DO 198 ELSE ! other coordinate system (zco or sco) : e3t199 DO jj = 2, jpjm1200 DO ji = fs_2, jpi ! vector opt.201 zabe1 = (fsahmt(ji,jj,jk)+ahmb0) * e2t(ji,jj) * fse3t(ji,jj,jk) / e1t(ji,jj)202 203 zmskt = 1./MAX( umask(ji-1,jj,jk )+umask(ji,jj,jk+1) &204 & + umask(ji-1,jj,jk+1)+umask(ji,jj,jk ), 1. )205 206 zcof1 = - aht0 * e2t(ji,jj) * zmskt * 0.5 * ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) )207 208 ziut(ji,jj) = ( zabe1 * ( ub(ji,jj,jk) - ub(ji-1,jj,jk) ) &209 & + zcof1 * ( zdku (ji,jj) + zdk1u(ji-1,jj) &210 & +zdk1u(ji,jj) + zdku (ji-1,jj) ) ) * tmask(ji,jj,jk)211 END DO212 END DO213 212 ENDIF 214 213 … … 216 215 DO jj = 1, jpjm1 217 216 DO ji = 1, fs_jpim1 ! vector opt. 218 zabe2 = ( fsahmf(ji,jj,jk) + ahmb0 ) * e1f(ji,jj) * fse3f(ji,jj,jk) /e2f(ji,jj)219 220 zmskf = 1. /MAX( umask(ji,jj+1,jk )+umask(ji,jj,jk+1)&221 & + umask(ji,jj+1,jk+1)+umask(ji,jj,jk ), 1.)222 223 zcof2 = - aht0 * e1f(ji,jj) * zmskf * 0.5 * ( vslp(ji+1,jj,jk) + vslp(ji,jj,jk) )217 zabe2 = ( ahmf(ji,jj,jk) + rn_ahm_b ) * e1f(ji,jj) * fse3f(ji,jj,jk) * r1_e2f(ji,jj) 218 219 zmskf = 1._wp / MAX( umask(ji,jj+1,jk )+umask(ji,jj,jk+1) & 220 & + umask(ji,jj+1,jk+1)+umask(ji,jj,jk ) , 1._wp ) 221 222 zcof2 = - rn_aht_0 * e1f(ji,jj) * zmskf * 0.5 * ( vslp(ji+1,jj,jk) + vslp(ji,jj,jk) ) 224 223 225 224 zjuf(ji,jj) = ( zabe2 * ( ub(ji,jj+1,jk) - ub(ji,jj,jk) ) & … … 237 236 DO jj = 2, jpjm1 238 237 DO ji = 1, fs_jpim1 ! vector opt. 239 zabe1 = ( fsahmf(ji,jj,jk) + ahmb0 ) * e2f(ji,jj) * fse3f(ji,jj,jk) /e1f(ji,jj)240 241 zmskf = 1. /MAX( vmask(ji+1,jj,jk )+vmask(ji,jj,jk+1)&242 & + vmask(ji+1,jj,jk+1)+vmask(ji,jj,jk ), 1.)243 244 zcof1 = - aht0 * e2f(ji,jj) * zmskf * 0.5 * ( uslp(ji,jj+1,jk) + uslp(ji,jj,jk) )245 246 zivf(ji,jj) = ( zabe1 * ( vb(ji+1,jj,jk) - vb(ji,jj,jk) ) &247 & + zcof1 * ( zdkv (ji,jj) + zdk1v(ji+1,jj)&248 & + zdk1v(ji,jj) + zdkv (ji+1,jj) ) ) * fmask(ji,jj,jk)238 zabe1 = ( ahmf(ji,jj,jk) + rn_ahm_b ) * e2f(ji,jj) * fse3f(ji,jj,jk) * r1_e1f(ji,jj) 239 240 zmskf = 1._wp / MAX( vmask(ji+1,jj,jk )+vmask(ji,jj,jk+1) & 241 & + vmask(ji+1,jj,jk+1)+vmask(ji,jj,jk ) , 1._wp ) 242 243 zcof1 = - rn_aht_0 * e2f(ji,jj) * zmskf * 0.5 * ( uslp(ji,jj+1,jk) + uslp(ji,jj,jk) ) 244 245 zivf(ji,jj) = ( zabe1 * ( vb(ji+1,jj,jk) - vb(ji,jj,jk) ) & 246 & + zcof1 * ( zdkv (ji,jj) + zdk1v(ji+1,jj) & 247 & + zdk1v(ji,jj) + zdkv (ji+1,jj) ) ) * fmask(ji,jj,jk) 249 248 END DO 250 249 END DO … … 254 253 DO jj = 2, jpj 255 254 DO ji = 1, fs_jpim1 ! vector opt. 256 zabe2 = (fsahmt(ji,jj,jk)+ahmb0) * e1t(ji,jj) * MIN( fse3v(ji,jj,jk), fse3v(ji,jj-1,jk) ) / e2t(ji,jj) 255 zabe2 = ( ahmt(ji,jj,jk)+rn_ahm_b ) * e1t(ji,jj) * MIN( fse3v(ji,jj,jk), fse3v(ji,jj-1,jk) ) * r1_e2t(ji,jj) 256 257 zmskt = 1._wp / MAX( vmask(ji,jj-1,jk )+vmask(ji,jj,jk+1) & 258 & + vmask(ji,jj-1,jk+1)+vmask(ji,jj,jk ) , 1._wp ) 259 260 zcof2 = - rn_aht_0 * e1t(ji,jj) * zmskt * 0.5 * ( vslp(ji,jj-1,jk) + vslp(ji,jj,jk) ) 261 262 zjvt(ji,jj) = ( zabe2 * ( vb(ji,jj,jk) - vb(ji,jj-1,jk) ) & 263 & + zcof2 * ( zdkv (ji,jj-1) + zdk1v(ji,jj) & 264 & +zdk1v(ji,jj-1) + zdkv (ji,jj) ) ) * tmask(ji,jj,jk) 265 END DO 266 END DO 267 ELSE ! other coordinate system (zco or sco) : e3t 268 DO jj = 2, jpj 269 DO ji = 1, fs_jpim1 ! vector opt. 270 zabe2 = ( ahmt(ji,jj,jk)+rn_ahm_b ) * e1t(ji,jj) * fse3t(ji,jj,jk) * r1_e2t(ji,jj) 257 271 258 272 zmskt = 1./MAX( vmask(ji,jj-1,jk )+vmask(ji,jj,jk+1) & 259 273 & + vmask(ji,jj-1,jk+1)+vmask(ji,jj,jk ), 1. ) 260 274 261 zcof2 = - aht0 * e1t(ji,jj) * zmskt * 0.5 * ( vslp(ji,jj-1,jk) + vslp(ji,jj,jk) )275 zcof2 = - rn_aht_0 * e1t(ji,jj) * zmskt * 0.5 * ( vslp(ji,jj-1,jk) + vslp(ji,jj,jk) ) 262 276 263 277 zjvt(ji,jj) = ( zabe2 * ( vb(ji,jj,jk) - vb(ji,jj-1,jk) ) & … … 266 280 END DO 267 281 END DO 268 ELSE ! other coordinate system (zco or sco) : e3t269 DO jj = 2, jpj270 DO ji = 1, fs_jpim1 ! vector opt.271 zabe2 = (fsahmt(ji,jj,jk)+ahmb0) * e1t(ji,jj) * fse3t(ji,jj,jk) / e2t(ji,jj)272 273 zmskt = 1./MAX( vmask(ji,jj-1,jk )+vmask(ji,jj,jk+1) &274 & + vmask(ji,jj-1,jk+1)+vmask(ji,jj,jk ), 1. )275 276 zcof2 = - aht0 * e1t(ji,jj) * zmskt * 0.5 * ( vslp(ji,jj-1,jk) + vslp(ji,jj,jk) )277 278 zjvt(ji,jj) = ( zabe2 * ( vb(ji,jj,jk) - vb(ji,jj-1,jk) ) &279 & + zcof2 * ( zdkv (ji,jj-1) + zdk1v(ji,jj) &280 & +zdk1v(ji,jj-1) + zdkv (ji,jj) ) ) * tmask(ji,jj,jk)281 END DO282 END DO283 282 ENDIF 284 283 … … 286 285 ! Second derivative (divergence) and add to the general trend 287 286 ! ----------------------------------------------------------- 288 289 287 DO jj = 2, jpjm1 290 DO ji = 2, jpim1 !! Question vectop possible??? !!bug 291 ! volume elements 292 zbu = e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) 293 zbv = e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) 294 ! horizontal component of isopycnal momentum diffusive trends 295 zuah =( ziut (ji+1,jj) - ziut (ji,jj ) + & 296 & zjuf (ji ,jj) - zjuf (ji,jj-1) ) / zbu 297 zvah =( zivf (ji,jj ) - zivf (ji-1,jj) + & 298 & zjvt (ji,jj+1) - zjvt (ji,jj ) ) / zbv 299 ! add the trends to the general trends 300 ua (ji,jj,jk) = ua (ji,jj,jk) + zuah 301 va (ji,jj,jk) = va (ji,jj,jk) + zvah 288 DO ji = 2, jpim1 !!gm Question vectop possible??? !!bug 289 ua(ji,jj,jk) = ua(ji,jj,jk) + ( ziut(ji+1,jj) - ziut(ji,jj ) & 290 & + zjuf(ji ,jj) - zjuf(ji,jj-1) ) / ( e1e2u(ji,jj) * fse3u(ji,jj,jk) ) 291 va(ji,jj,jk) = va(ji,jj,jk) + ( zivf(ji,jj ) - zivf(ji-1,jj) & 292 & + zjvt(ji,jj+1) - zjvt(ji,jj ) ) / ( e1e2v(ji,jj) * fse3v(ji,jj,jk) ) 302 293 END DO 303 294 END DO … … 358 349 DO jk = 2, jpkm1 359 350 DO ji = 2, jpim1 360 zcoef0= 0.5 * aht0 * umask(ji,jj,jk)361 351 zcoef0= 0.5 * rn_aht_0 * umask(ji,jj,jk) 352 ! 362 353 zuwslpi = zcoef0 * ( wslpi(ji+1,jj,jk) + wslpi(ji,jj,jk) ) 363 354 zuwslpj = zcoef0 * ( wslpj(ji+1,jj,jk) + wslpj(ji,jj,jk) ) 364 355 ! 365 356 zmkt = 1./MAX( tmask(ji,jj,jk-1)+tmask(ji+1,jj,jk-1) & 366 357 + tmask(ji,jj,jk )+tmask(ji+1,jj,jk ), 1. ) 367 zmkf = 1./MAX( fmask(ji,jj-1,jk-1) +fmask(ji,jj,jk-1) &368 + fmask(ji,jj-1,jk ) +fmask(ji,jj,jk ), 1. )358 zmkf = 1./MAX( fmask(ji,jj-1,jk-1) + fmask(ji,jj,jk-1) & 359 + fmask(ji,jj-1,jk ) + fmask(ji,jj,jk ), 1. ) 369 360 370 361 zcoef3 = - e2u(ji,jj) * zmkt * zuwslpi … … 376 367 +zdj1u(ji,jk ) + zdju (ji ,jk ) ) 377 368 ! update avmu (add isopycnal vertical coefficient to avmu) 378 ! Caution: zcoef0 include aht0, so divided by aht0 to obtain slp^2 * aht0379 avmu(ji,jj,jk) = avmu(ji,jj,jk) + ( zuwslpi * zuwslpi + zuwslpj * zuwslpj ) / aht0369 ! Caution: zcoef0 include rn_aht_0, so divided by rn_aht_0 to obtain slp^2 * rn_aht_0 370 avmu(ji,jj,jk) = avmu(ji,jj,jk) + ( zuwslpi * zuwslpi + zuwslpj * zuwslpj ) / rn_aht_0 380 371 END DO 381 372 END DO … … 384 375 DO jk = 2, jpkm1 385 376 DO ji = 2, jpim1 386 zcoef0 = 0.5 * aht0 * vmask(ji,jj,jk)377 zcoef0 = 0.5 * rn_aht_0 * vmask(ji,jj,jk) 387 378 388 379 zvwslpi = zcoef0 * ( wslpi(ji,jj+1,jk) + wslpi(ji,jj,jk) ) … … 398 389 ! vertical flux on v field 399 390 zfvw(ji,jk) = zcoef3 * ( zdiv (ji,jk-1) + zdiv (ji-1,jk-1) & 400 401 402 391 & +zdiv (ji,jk ) + zdiv (ji-1,jk ) ) & 392 & + zcoef4 * ( zdjv (ji,jk-1) + zdj1v(ji ,jk-1) & 393 & +zdjv (ji,jk ) + zdj1v(ji ,jk ) ) 403 394 ! update avmv (add isopycnal vertical coefficient to avmv) 404 ! Caution: zcoef0 include aht0, so divided by aht0 to obtain slp^2 * aht0405 avmv(ji,jj,jk) = avmv(ji,jj,jk) + ( zvwslpi * zvwslpi + zvwslpj * zvwslpj ) / aht0395 ! Caution: zcoef0 include rn_aht_0, so divided by rn_aht_0 to obtain slp^2 * rn_aht_0 396 avmv(ji,jj,jk) = avmv(ji,jj,jk) + ( zvwslpi * zvwslpi + zvwslpj * zvwslpj ) / rn_aht_0 406 397 END DO 407 398 END DO … … 412 403 DO jk = 1, jpkm1 413 404 DO ji = 2, jpim1 414 ! volume elements 415 zbu = e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) 416 zbv = e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) 417 ! part of the k-component of isopycnal momentum diffusive trends 418 zuav = ( zfuw(ji,jk) - zfuw(ji,jk+1) ) / zbu 419 zvav = ( zfvw(ji,jk) - zfvw(ji,jk+1) ) / zbv 420 ! add the trends to the general trends 421 ua(ji,jj,jk) = ua(ji,jj,jk) + zuav 422 va(ji,jj,jk) = va(ji,jj,jk) + zvav 405 ua(ji,jj,jk) = ua(ji,jj,jk) + ( zfuw(ji,jk) - zfuw(ji,jk+1) ) / ( e1e2u(ji,jj) * fse3u(ji,jj,jk) ) 406 va(ji,jj,jk) = va(ji,jj,jk) + ( zfvw(ji,jk) - zfvw(ji,jk+1) ) / ( e1e2v(ji,jj) * fse3v(ji,jj,jk) ) 423 407 END DO 424 408 END DO … … 432 416 END SUBROUTINE dyn_ldf_iso 433 417 434 # else435 !!----------------------------------------------------------------------436 !! Dummy module NO rotation of mixing tensor437 !!----------------------------------------------------------------------438 CONTAINS439 SUBROUTINE dyn_ldf_iso( kt ) ! Empty routine440 INTEGER, INTENT(in) :: kt441 WRITE(*,*) 'dyn_ldf_iso: You should not have seen this print! error?', kt442 END SUBROUTINE dyn_ldf_iso443 #endif444 445 418 !!====================================================================== 446 419 END MODULE dynldf_iso -
branches/2014/dev_r4650_UKMO14.4_OBS_GENERAL_VINTERP/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90
r5950 r5951 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(:,:,:) ) 266 221 ! Add volume filter correction: compatibility with tracer advection scheme 267 222 ! => time filter + conservation correction (only at the first level) 268 fse3t_b(:,:,1) = fse3t_b(:,:,1) - atfp * rdt * r1_rau0 * ( emp_b(:,:) - emp(:,:) & 269 & -rnf_b(:,:) + rnf(:,:) ) * tmask(:,:,1) 223 IF ( nn_isf == 0) THEN ! if no ice shelf melting 224 fse3t_b(:,:,1) = fse3t_b(:,:,1) - atfp * rdt * r1_rau0 * ( emp_b(:,:) - emp(:,:) & 225 & -rnf_b(:,:) + rnf(:,:) ) * tmask(:,:,1) 226 ELSE ! if ice shelf melting 227 DO jj = 1,jpj 228 DO ji = 1,jpi 229 jk = mikt(ji,jj) 230 fse3t_b(ji,jj,jk) = fse3t_b(ji,jj,jk) - atfp * rdt * r1_rau0 & 231 & * ( (emp_b(ji,jj) - emp(ji,jj) ) & 232 & - (rnf_b(ji,jj) - rnf(ji,jj) ) & 233 & + (fwfisf_b(ji,jj) - fwfisf(ji,jj)) ) * tmask(ji,jj,jk) 234 END DO 235 END DO 236 END IF 270 237 ENDIF 271 238 ! … … 324 291 ENDIF 325 292 ! 326 IF (l k_dynspg_ts.AND.ln_bt_fw) THEN293 IF (ln_dynspg_ts.AND.ln_bt_fw) THEN 327 294 ! Revert "before" velocities to time split estimate 328 295 ! Doing it here also means that asselin filter contribution is removed … … 388 355 IF(ln_ctl) CALL prt_ctl( tab3d_1=un, clinfo1=' nxt - Un: ', mask1=umask, & 389 356 & tab3d_2=vn, clinfo2=' Vn: ' , mask2=vmask ) 390 ! 391 CALL wrk_dealloc( jpi,jpj,jpk, ze3u_f, ze3v_f, zua, zva ) 392 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 ) 393 361 ! 394 362 IF( nn_timing == 1 ) CALL timing_stop('dyn_nxt') -
branches/2014/dev_r4650_UKMO14.4_OBS_GENERAL_VINTERP/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg.F90
r5950 r5951 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 !! … … 77 72 !! Note that as all external forcing a time averaging over a two rdt 78 73 !! period is used to prevent the divergence of odd and even time step. 79 !!80 !! N.B. : When key_esopa is used all the scheme are tested, regardless81 !! of the physical meaning of the results.82 74 !!---------------------------------------------------------------------- 83 75 ! 84 76 INTEGER, INTENT(in ) :: kt ! ocean time-step index 85 INTEGER, INTENT( out) :: kindic ! solver flag86 77 ! 87 78 INTEGER :: ji, jj, jk ! dummy loop indices … … 106 97 107 98 IF( ln_apr_dyn & ! atmos. pressure 108 .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) 109 100 .OR. nn_ice_embd == 2 ) THEN ! embedded sea-ice 110 101 ! … … 116 107 END DO 117 108 ! 118 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) ==! 119 110 zg_2 = grav * 0.5 120 111 DO jj = 2, jpjm1 ! gradient of Patm using inverse barometer ssh 121 112 DO ji = fs_2, fs_jpim1 ! vector opt. 122 113 spgu(ji,jj) = spgu(ji,jj) + zg_2 * ( ssh_ib (ji+1,jj) - ssh_ib (ji,jj) & 123 & + ssh_ibb(ji+1,jj) - ssh_ibb(ji,jj) ) /e1u(ji,jj)114 & + ssh_ibb(ji+1,jj) - ssh_ibb(ji,jj) ) * r1_e1u(ji,jj) 124 115 spgv(ji,jj) = spgv(ji,jj) + zg_2 * ( ssh_ib (ji,jj+1) - ssh_ib (ji,jj) & 125 & + ssh_ibb(ji,jj+1) - ssh_ibb(ji,jj) ) /e2v(ji,jj)116 & + ssh_ibb(ji,jj+1) - ssh_ibb(ji,jj) ) * r1_e2v(ji,jj) 126 117 END DO 127 118 END DO … … 129 120 ! 130 121 ! !== tide potential forcing term ==! 131 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 132 123 ! 133 124 CALL upd_tide( kt ) ! update tide potential … … 135 126 DO jj = 2, jpjm1 ! add tide potential forcing 136 127 DO ji = fs_2, fs_jpim1 ! vector opt. 137 spgu(ji,jj) = spgu(ji,jj) + grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) /e1u(ji,jj)138 spgv(ji,jj) = spgv(ji,jj) + grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) /e2v(ji,jj)128 spgu(ji,jj) = spgu(ji,jj) + grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) * r1_e1u(ji,jj) 129 spgv(ji,jj) = spgv(ji,jj) + grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) * r1_e2v(ji,jj) 139 130 END DO 140 131 END DO … … 149 140 DO jj = 2, jpjm1 150 141 DO ji = fs_2, fs_jpim1 ! vector opt. 151 spgu(ji,jj) = spgu(ji,jj) + ( zpice(ji+1,jj) - zpice(ji,jj) ) /e1u(ji,jj)152 spgv(ji,jj) = spgv(ji,jj) + ( zpice(ji,jj+1) - zpice(ji,jj) ) /e2v(ji,jj)142 spgu(ji,jj) = spgu(ji,jj) + ( zpice(ji+1,jj) - zpice(ji,jj) ) * r1_e1u(ji,jj) 143 spgv(ji,jj) = spgv(ji,jj) + ( zpice(ji,jj+1) - zpice(ji,jj) ) * r1_e2v(ji,jj) 153 144 END DO 154 145 END DO … … 174 165 CASE ( 0 ) ; CALL dyn_spg_exp( kt ) ! explicit 175 166 CASE ( 1 ) ; CALL dyn_spg_ts ( kt ) ! time-splitting 176 CASE ( 2 ) ; CALL dyn_spg_flt( kt, kindic ) ! filtered177 167 ! 178 CASE ( -1 ) ! esopa: test all possibility with control print179 CALL dyn_spg_exp( kt )180 CALL prt_ctl( tab3d_1=ua, clinfo1=' spg0 - Ua: ', mask1=umask, &181 & tab3d_2=va, clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' )182 CALL dyn_spg_ts ( kt )183 CALL prt_ctl( tab3d_1=ua, clinfo1=' spg1 - Ua: ', mask1=umask, &184 & tab3d_2=va, clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' )185 CALL dyn_spg_flt( kt, kindic )186 CALL prt_ctl( tab3d_1=ua, clinfo1=' spg2 - Ua: ', mask1=umask, &187 & tab3d_2=va, clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' )188 168 END SELECT 189 169 ! 190 170 IF( l_trddyn ) THEN ! save the surface pressure gradient trends for further diagnostics 191 SELECT CASE ( nspg ) 192 CASE ( 0, 1 ) 193 ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 194 ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 195 CASE( 2 ) 196 z2dt = 2. * rdt 197 IF( neuler == 0 .AND. kt == nit000 ) z2dt = rdt 198 ztrdu(:,:,:) = ( ua(:,:,:) - ub(:,:,:) ) / z2dt - ztrdu(:,:,:) 199 ztrdv(:,:,:) = ( va(:,:,:) - vb(:,:,:) ) / z2dt - ztrdv(:,:,:) 200 END SELECT 171 ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 172 ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 201 173 CALL trd_dyn( ztrdu, ztrdv, jpdyn_spg, kt ) 202 174 ! … … 216 188 !! *** ROUTINE dyn_spg_init *** 217 189 !! 218 !! ** Purpose : Control the consistency between cppoptions for190 !! ** Purpose : Control the consistency between namelist options for 219 191 !! surface pressure gradient schemes 220 192 !!---------------------------------------------------------------------- 221 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 222 198 !!---------------------------------------------------------------------- 223 199 ! 224 200 IF( nn_timing == 1 ) CALL timing_start('dyn_spg_init') 225 201 ! 226 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 227 212 WRITE(numout,*) 228 213 WRITE(numout,*) 'dyn_spg_init : choice of the surface pressure gradient scheme' 229 214 WRITE(numout,*) '~~~~~~~~~~~' 230 WRITE(numout,*) ' Explicit free surface lk_dynspg_exp = ', lk_dynspg_exp 231 WRITE(numout,*) ' Free surface with time splitting lk_dynspg_ts = ', lk_dynspg_ts 232 WRITE(numout,*) ' Filtered free surface cst volume lk_dynspg_flt = ', lk_dynspg_flt 233 ENDIF 234 235 IF( lk_dynspg_ts ) CALL dyn_spg_ts_init( nit000 ) 236 ! (do it now, to set nn_baro, used to allocate some arrays later on) 237 ! ! allocate dyn_spg arrays 238 IF( lk_dynspg_ts ) THEN 239 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 240 222 IF( dyn_spg_ts_alloc() /= 0 ) CALL ctl_stop('STOP', 'dyn_spg_init: failed to allocate dynspg_ts arrays') 241 223 IF ((neuler/=0).AND.(ln_bt_fw)) CALL ts_rst( nit000, 'READ' ) … … 244 226 ! ! Control of surface pressure gradient scheme options 245 227 ioptio = 0 246 IF(lk_dynspg_exp) ioptio = ioptio + 1 247 IF(lk_dynspg_ts ) ioptio = ioptio + 1 248 IF(lk_dynspg_flt) ioptio = ioptio + 1 249 ! 250 IF( ( ioptio > 1 .AND. .NOT. lk_esopa ) .OR. ( ioptio == 0 .AND. .NOT. lk_c1d ) ) & 251 & CALL ctl_stop( ' Choose only one surface pressure gradient scheme with a key cpp' ) 252 IF( ( lk_dynspg_ts .OR. lk_dynspg_exp ) .AND. ln_isfcav ) & 253 & CALL ctl_stop( ' dynspg_ts and dynspg_exp not tested with ice shelf cavity ' ) 254 ! 255 IF( lk_esopa ) nspg = -1 256 IF( lk_dynspg_exp) nspg = 0 257 IF( lk_dynspg_ts ) nspg = 1 258 IF( lk_dynspg_flt) nspg = 2 259 ! 260 IF( lk_esopa ) nspg = -1 228 IF(ln_dynspg_exp) ioptio = ioptio + 1 229 IF(ln_dynspg_ts ) ioptio = ioptio + 1 230 ! 231 IF( ioptio > 1 .OR. ( ioptio == 0 .AND. .NOT. lk_c1d ) ) & 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 261 238 ! 262 239 IF(lwp) THEN 263 240 WRITE(numout,*) 264 IF( nspg == -1 ) WRITE(numout,*) ' ESOPA test All scheme used'265 241 IF( nspg == 0 ) WRITE(numout,*) ' explicit free surface' 266 242 IF( nspg == 1 ) WRITE(numout,*) ' free surface with time splitting scheme' 267 IF( nspg == 2 ) WRITE(numout,*) ' filtered free surface'268 ENDIF269 270 #if defined key_dynspg_flt || defined key_esopa271 CALL solver_init( nit000 ) ! Elliptic solver initialisation272 #endif273 274 ! ! Control of timestep choice275 IF( lk_dynspg_ts .OR. lk_dynspg_exp ) THEN276 IF( nn_cla == 1 ) CALL ctl_stop( 'Crossland advection not implemented for this free surface formulation' )277 ENDIF278 279 ! ! Control of hydrostatic pressure choice280 IF( lk_dynspg_ts .AND. ln_dynhpg_imp ) THEN281 CALL ctl_stop( 'Semi-implicit hpg not compatible with time splitting' )282 243 ENDIF 283 244 ! -
branches/2014/dev_r4650_UKMO14.4_OBS_GENERAL_VINTERP/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_exp.F90
r5950 r5951 7 7 !! 3.2 ! 2009-06 (G. Madec, M. Leclair, R. Benshila) introduce sshwzv module 8 8 !!---------------------------------------------------------------------- 9 #if defined key_dynspg_exp || defined key_esopa10 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_r4650_UKMO14.4_OBS_GENERAL_VINTERP/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90
r5950 r5951 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_ts || defined key_esopa15 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) ) 99 93 100 IF( ln_dynvor_een .or. ln_dynvor_een_old ) ALLOCATE( ftnw(jpi,jpj) , ftne(jpi,jpj) , & 101 & ftsw(jpi,jpj) , ftse(jpi,jpj) , STAT=ierr(3) ) 94 IF( ln_dynvor_een ) ALLOCATE( ftnw(jpi,jpj) , ftne(jpi,jpj) , & 95 & ftsw(jpi,jpj) , ftse(jpi,jpj) , STAT=ierr(3) ) 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 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 109 109 110 110 111 SUBROUTINE dyn_spg_ts( kt ) … … 145 146 INTEGER :: ikbu, ikbv, noffset ! local integers 146 147 REAL(wp) :: zraur, z1_2dt_b, z2dt_bf ! local scalars 147 REAL(wp) :: zx1, zy1, zx2, zy2 ! - -148 REAL(wp) :: z1_12, z1_8, z1_4, z1_2 149 REAL(wp) :: zu_spg, zv_spg ! - -150 REAL(wp) :: zhura, zhvra 151 REAL(wp) :: za0, za1, za2, za3 152 ! 153 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 154 155 REAL(wp), POINTER, DIMENSION(:,:) :: zu_trd, zv_trd, zu_frc, zv_frc, zssh_frc 155 REAL(wp), POINTER, DIMENSION(:,:) :: z u_sum, zv_sum, zwx, zwy, zhdiv156 REAL(wp), POINTER, DIMENSION(:,:) :: zwx, zwy, zhdiv 156 157 REAL(wp), POINTER, DIMENSION(:,:) :: zhup2_e, zhvp2_e, zhust_e, zhvst_e 157 158 REAL(wp), POINTER, DIMENSION(:,:) :: zsshu_a, zsshv_a … … 163 164 ! !* Allocate temporary arrays 164 165 CALL wrk_alloc( jpi, jpj, zsshp2_e, zhdiv ) 165 CALL wrk_alloc( jpi, jpj, zu_trd, zv_trd , zun_e, zvn_e)166 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) 167 168 CALL wrk_alloc( jpi, jpj, zhup2_e, zhvp2_e, zhust_e, zhvst_e) 168 169 CALL wrk_alloc( jpi, jpj, zsshu_a, zsshv_a ) … … 187 188 ! 188 189 ! time offset in steps for bdy data update 189 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 190 191 ! 191 192 IF( kt == nit000 ) THEN !* initialisation … … 219 220 ! 220 221 IF ( kt == nit000 .OR. lk_vvl ) THEN 221 IF ( ln_dynvor_een_old ) THEN 222 DO jj = 1, jpjm1 223 DO ji = 1, jpim1 224 zwz(ji,jj) = ( ht(ji ,jj+1) + ht(ji+1,jj+1) + & 225 & ht(ji ,jj ) + ht(ji+1,jj ) ) / 4._wp 226 IF( zwz(ji,jj) /= 0._wp ) zwz(ji,jj) = 1._wp / zwz(ji,jj) 227 END DO 228 END DO 222 IF ( ln_dynvor_een ) THEN !== EEN scheme ==! 223 SELECT CASE( nn_een_e3f ) !* ff/e3 at F-point 224 CASE ( 0 ) ! original formulation (masked averaging of e3t divided by 4) 225 DO jj = 1, jpjm1 226 DO ji = 1, jpim1 227 zwz(ji,jj) = ( ht(ji ,jj+1) + ht(ji+1,jj+1) + & 228 & ht(ji ,jj ) + ht(ji+1,jj ) ) / 4._wp 229 IF( zwz(ji,jj) /= 0._wp ) zwz(ji,jj) = ff(ji,jj) / zwz(ji,jj) 230 END DO 231 END DO 232 CASE ( 1 ) ! new formulation (masked averaging of e3t divided by the sum of mask) 233 DO jj = 1, jpjm1 234 DO ji = 1, jpim1 235 zwz(ji,jj) = ( ht(ji ,jj+1) + ht(ji+1,jj+1) + & 236 & ht(ji ,jj ) + ht(ji+1,jj ) ) & 237 & / ( MAX( 1._wp, tmask(ji ,jj+1, 1) + tmask(ji+1,jj+1, 1) + & 238 & tmask(ji ,jj , 1) + tmask(ji+1,jj , 1) ) ) 239 IF( zwz(ji,jj) /= 0._wp ) zwz(ji,jj) = ff(ji,jj) / zwz(ji,jj) 240 END DO 241 END DO 242 END SELECT 229 243 CALL lbc_lnk( zwz, 'F', 1._wp ) 230 zwz(:,:) = ff(:,:) * zwz(:,:) 231 244 ! 232 245 ftne(1,:) = 0._wp ; ftnw(1,:) = 0._wp ; ftse(1,:) = 0._wp ; ftsw(1,:) = 0._wp 233 246 DO jj = 2, jpj 234 DO ji = fs_2, jpi ! vector opt.247 DO ji = 2, jpi 235 248 ftne(ji,jj) = zwz(ji-1,jj ) + zwz(ji ,jj ) + zwz(ji ,jj-1) 236 249 ftnw(ji,jj) = zwz(ji-1,jj-1) + zwz(ji-1,jj ) + zwz(ji ,jj ) … … 239 252 END DO 240 253 END DO 241 ELSE IF ( ln_dynvor_een ) THEN 242 DO jj = 1, jpjm1 243 DO ji = 1, jpim1 244 zwz(ji,jj) = ( ht(ji ,jj+1) + ht(ji+1,jj+1) + & 245 & ht(ji ,jj ) + ht(ji+1,jj ) ) & 246 & / ( MAX( 1.0_wp, tmask(ji ,jj+1, 1) + tmask(ji+1,jj+1, 1) + & 247 & tmask(ji ,jj , 1) + tmask(ji+1,jj , 1) ) ) 248 IF( zwz(ji,jj) /= 0._wp ) zwz(ji,jj) = 1._wp / zwz(ji,jj) 249 END DO 250 END DO 251 CALL lbc_lnk( zwz, 'F', 1._wp ) 252 zwz(:,:) = ff(:,:) * zwz(:,:) 253 254 ftne(1,:) = 0._wp ; ftnw(1,:) = 0._wp ; ftse(1,:) = 0._wp ; ftsw(1,:) = 0._wp 255 DO jj = 2, jpj 256 DO ji = fs_2, jpi ! vector opt. 257 ftne(ji,jj) = zwz(ji-1,jj ) + zwz(ji ,jj ) + zwz(ji ,jj-1) 258 ftnw(ji,jj) = zwz(ji-1,jj-1) + zwz(ji-1,jj ) + zwz(ji ,jj ) 259 ftse(ji,jj) = zwz(ji ,jj ) + zwz(ji ,jj-1) + zwz(ji-1,jj-1) 260 ftsw(ji,jj) = zwz(ji ,jj-1) + zwz(ji-1,jj-1) + zwz(ji-1,jj ) 261 END DO 262 END DO 263 ELSE 254 ! 255 ELSE !== all other schemes (ENE, ENS, MIX) 264 256 zwz(:,:) = 0._wp 265 zhf(:,:) = 0. 257 zhf(:,:) = 0._wp 266 258 IF ( .not. ln_sco ) THEN 259 260 !!gm agree the JC comment : this should be done in a much clear way 261 267 262 ! JC: It not clear yet what should be the depth at f-points over land in z-coordinate case 268 263 ! Set it to zero for the time being … … 276 271 277 272 DO jj = 1, jpjm1 278 zhf(:,jj) = zhf(:,jj) *(1._wp- umask(:,jj,1) * umask(:,jj+1,1))273 zhf(:,jj) = zhf(:,jj) * (1._wp- umask(:,jj,1) * umask(:,jj+1,1)) 279 274 END DO 280 275 … … 297 292 ! If forward start at previous time step, and centered integration, 298 293 ! then update averaging weights: 299 IF ( (.NOT.ln_bt_fw).AND.((neuler==0).AND.(kt==nit000+1))) THEN294 IF (.NOT.ln_bt_fw .AND.( neuler==0 .AND. kt==nit000+1 ) ) THEN 300 295 ll_fw_start=.FALSE. 301 296 CALL ts_wgt(ln_bt_av, ll_fw_start, icycle, wgtbtp1, wgtbtp2) … … 338 333 DO jj = 2, jpjm1 339 334 DO ji = fs_2, fs_jpim1 ! vector opt. 340 zy1 = ( zwy(ji,jj-1) + zwy(ji+1,jj-1) ) /e1u(ji,jj)341 zy2 = ( zwy(ji,jj ) + zwy(ji+1,jj ) ) /e1u(ji,jj)342 zx1 = ( zwx(ji-1,jj) + zwx(ji-1,jj+1) ) /e2v(ji,jj)343 zx2 = ( zwx(ji ,jj) + zwx(ji ,jj+1) ) /e2v(ji,jj)335 zy1 = ( zwy(ji,jj-1) + zwy(ji+1,jj-1) ) * r1_e1u(ji,jj) 336 zy2 = ( zwy(ji,jj ) + zwy(ji+1,jj ) ) * r1_e1u(ji,jj) 337 zx1 = ( zwx(ji-1,jj) + zwx(ji-1,jj+1) ) * r1_e2v(ji,jj) 338 zx2 = ( zwx(ji ,jj) + zwx(ji ,jj+1) ) * r1_e2v(ji,jj) 344 339 ! energy conserving formulation for planetary vorticity term 345 340 zu_trd(ji,jj) = z1_4 * ( zwz(ji ,jj-1) * zy1 + zwz(ji,jj) * zy2 ) … … 352 347 DO ji = fs_2, fs_jpim1 ! vector opt. 353 348 zy1 = z1_8 * ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) & 354 & + zwy(ji ,jj ) + zwy(ji+1,jj ) ) /e1u(ji,jj)349 & + zwy(ji ,jj ) + zwy(ji+1,jj ) ) * r1_e1u(ji,jj) 355 350 zx1 = - z1_8 * ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) & 356 & + zwx(ji ,jj ) + zwx(ji ,jj+1) ) /e2v(ji,jj)351 & + zwx(ji ,jj ) + zwx(ji ,jj+1) ) * r1_e2v(ji,jj) 357 352 zu_trd(ji,jj) = zy1 * ( zwz(ji ,jj-1) + zwz(ji,jj) ) 358 353 zv_trd(ji,jj) = zx1 * ( zwz(ji-1,jj ) + zwz(ji,jj) ) … … 360 355 END DO 361 356 ! 362 ELSEIF ( ln_dynvor_een .or. ln_dynvor_een_old) THEN ! enstrophy and energy conserving scheme357 ELSEIF ( ln_dynvor_een ) THEN ! enstrophy and energy conserving scheme 363 358 DO jj = 2, jpjm1 364 359 DO ji = fs_2, fs_jpim1 ! vector opt. 365 zu_trd(ji,jj) = + z1_12 /e1u(ji,jj) * ( ftne(ji,jj ) * zwy(ji ,jj ) &366 & + ftnw(ji+1,jj) * zwy(ji+1,jj ) &367 & + ftse(ji,jj ) * zwy(ji ,jj-1) &368 & + ftsw(ji+1,jj) * zwy(ji+1,jj-1) )369 zv_trd(ji,jj) = - z1_12 /e2v(ji,jj) * ( ftsw(ji,jj+1) * zwx(ji-1,jj+1) &370 & + ftse(ji,jj+1) * zwx(ji ,jj+1) &371 & + ftnw(ji,jj ) * zwx(ji-1,jj ) &372 & + ftne(ji,jj ) * zwx(ji ,jj ) )360 zu_trd(ji,jj) = + z1_12 * r1_e1u(ji,jj) * ( ftne(ji,jj ) * zwy(ji ,jj ) & 361 & + ftnw(ji+1,jj) * zwy(ji+1,jj ) & 362 & + ftse(ji,jj ) * zwy(ji ,jj-1) & 363 & + ftsw(ji+1,jj) * zwy(ji+1,jj-1) ) 364 zv_trd(ji,jj) = - z1_12 * r1_e2v(ji,jj) * ( ftsw(ji,jj+1) * zwx(ji-1,jj+1) & 365 & + ftse(ji,jj+1) * zwx(ji ,jj+1) & 366 & + ftnw(ji,jj ) * zwx(ji-1,jj ) & 367 & + ftne(ji,jj ) * zwx(ji ,jj ) ) 373 368 END DO 374 369 END DO … … 381 376 DO jj = 2, jpjm1 382 377 DO ji = fs_2, fs_jpim1 ! vector opt. 383 zu_trd(ji,jj) = zu_trd(ji,jj) - grav * ( sshn(ji+1,jj ) - sshn(ji ,jj ) ) /e1u(ji,jj)384 zv_trd(ji,jj) = zv_trd(ji,jj) - grav * ( sshn(ji ,jj+1) - sshn(ji ,jj ) ) /e2v(ji,jj)378 zu_trd(ji,jj) = zu_trd(ji,jj) - grav * ( sshn(ji+1,jj ) - sshn(ji ,jj ) ) * r1_e1u(ji,jj) 379 zv_trd(ji,jj) = zv_trd(ji,jj) - grav * ( sshn(ji ,jj+1) - sshn(ji ,jj ) ) * r1_e2v(ji,jj) 385 380 END DO 386 381 END DO … … 431 426 DO jj = 2, jpjm1 432 427 DO ji = fs_2, fs_jpim1 ! vector opt. 433 zu_spg = grav * ( ssh_ib (ji+1,jj ) - ssh_ib (ji,jj) ) /e1u(ji,jj)434 zv_spg = grav * ( ssh_ib (ji ,jj+1) - ssh_ib (ji,jj) ) /e2v(ji,jj)428 zu_spg = grav * ( ssh_ib (ji+1,jj ) - ssh_ib (ji,jj) ) * r1_e1u(ji,jj) 429 zv_spg = grav * ( ssh_ib (ji ,jj+1) - ssh_ib (ji,jj) ) * r1_e2v(ji,jj) 435 430 zu_frc(ji,jj) = zu_frc(ji,jj) + zu_spg 436 431 zv_frc(ji,jj) = zv_frc(ji,jj) + zv_spg … … 441 436 DO ji = fs_2, fs_jpim1 ! vector opt. 442 437 zu_spg = grav * z1_2 * ( ssh_ib (ji+1,jj ) - ssh_ib (ji,jj) & 443 & + ssh_ibb(ji+1,jj ) - ssh_ibb(ji,jj) ) /e1u(ji,jj)438 & + ssh_ibb(ji+1,jj ) - ssh_ibb(ji,jj) ) * r1_e1u(ji,jj) 444 439 zv_spg = grav * z1_2 * ( ssh_ib (ji ,jj+1) - ssh_ib (ji,jj) & 445 & + ssh_ibb(ji ,jj+1) - ssh_ibb(ji,jj) ) /e2v(ji,jj)440 & + ssh_ibb(ji ,jj+1) - ssh_ibb(ji,jj) ) * r1_e2v(ji,jj) 446 441 zu_frc(ji,jj) = zu_frc(ji,jj) + zu_spg 447 442 zv_frc(ji,jj) = zv_frc(ji,jj) + zv_spg … … 454 449 ! ! Surface net water flux and rivers 455 450 IF (ln_bt_fw) THEN 456 zssh_frc(:,:) = zraur * ( emp(:,:) - rnf(:,:) + rdivisf *fwfisf(:,:) )451 zssh_frc(:,:) = zraur * ( emp(:,:) - rnf(:,:) + fwfisf(:,:) ) 457 452 ELSE 458 453 zssh_frc(:,:) = zraur * z1_2 * ( emp(:,:) + emp_b(:,:) - rnf(:,:) - rnf_b(:,:) & 459 & + rdivisf * ( fwfisf(:,:) + fwfisf_b(:,:) ))454 & + fwfisf(:,:) + fwfisf_b(:,:) ) 460 455 ENDIF 461 456 #if defined key_asminc … … 465 460 ENDIF 466 461 #endif 467 ! !* Fill boundary data arrays withAGRIF468 ! ! ------------------------------------ -462 ! !* Fill boundary data arrays for AGRIF 463 ! ! ------------------------------------ 469 464 #if defined key_agrif 470 465 IF( .NOT.Agrif_Root() ) CALL agrif_dta_ts( kt ) … … 490 485 IF (ln_bt_fw) THEN ! FORWARD integration: start from NOW fields 491 486 sshn_e(:,:) = sshn (:,:) 492 zun_e(:,:) = un_b (:,:)493 zvn_e(:,:) = vn_b (:,:)487 un_e (:,:) = un_b (:,:) 488 vn_e (:,:) = vn_b (:,:) 494 489 ! 495 490 hu_e (:,:) = hu (:,:) … … 499 494 ELSE ! CENTRED integration: start from BEFORE fields 500 495 sshn_e(:,:) = sshb (:,:) 501 zun_e(:,:) = ub_b (:,:)502 zvn_e(:,:) = vb_b (:,:)496 un_e (:,:) = ub_b (:,:) 497 vn_e (:,:) = vb_b (:,:) 503 498 ! 504 499 hu_e (:,:) = hu_b (:,:) … … 514 509 va_b (:,:) = 0._wp 515 510 ssha (:,:) = 0._wp ! Sum for after averaged sea level 516 zu_sum(:,:) = 0._wp ! Sum for now transport issued from ts loop517 zv_sum(:,:) = 0._wp511 un_adv(:,:) = 0._wp ! Sum for now transport issued from ts loop 512 vn_adv(:,:) = 0._wp 518 513 ! ! ==================== ! 519 514 DO jn = 1, icycle ! sub-time-step loop ! … … 523 518 ! Update only tidal forcing at open boundaries 524 519 #if defined key_tide 525 IF ( lk_bdy .AND. lk_tide ) 526 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 ) 527 522 #endif 528 523 ! … … 539 534 540 535 ! Extrapolate barotropic velocities at step jit+0.5: 541 ua_e(:,:) = za1 * zun_e(:,:) + za2 * ub_e(:,:) + za3 * ubb_e(:,:)542 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(:,:) 543 538 544 539 IF( lk_vvl ) THEN !* Update ocean depth (variable volume case only) … … 549 544 DO jj = 2, jpjm1 ! Sea Surface Height at u- & v-points 550 545 DO ji = 2, fs_jpim1 ! Vector opt. 551 zwx(ji,jj) = z1_2 * umask(ji,jj,1) * r1_e1 2u(ji,jj) &552 & * ( e1 2t(ji ,jj) * zsshp2_e(ji ,jj) &553 & + e1 2t(ji+1,jj) * zsshp2_e(ji+1,jj) )554 zwy(ji,jj) = z1_2 * vmask(ji,jj,1) * r1_e1 2v(ji,jj) &555 & * ( e1 2t(ji,jj ) * zsshp2_e(ji,jj ) &556 & + e1 2t(ji,jj+1) * zsshp2_e(ji,jj+1) )546 zwx(ji,jj) = z1_2 * umask(ji,jj,1) * r1_e1e2u(ji,jj) & 547 & * ( e1e2t(ji ,jj) * zsshp2_e(ji ,jj) & 548 & + e1e2t(ji+1,jj) * zsshp2_e(ji+1,jj) ) 549 zwy(ji,jj) = z1_2 * vmask(ji,jj,1) * r1_e1e2v(ji,jj) & 550 & * ( e1e2t(ji,jj ) * zsshp2_e(ji,jj ) & 551 & + e1e2t(ji,jj+1) * zsshp2_e(ji,jj+1) ) 557 552 END DO 558 553 END DO … … 602 597 ! Sum over sub-time-steps to compute advective velocities 603 598 za2 = wgtbtp2(jn) 604 zu_sum (:,:) = zu_sum (:,:) + za2 * zwx (:,:) / e2u(:,:)605 zv_sum (:,:) = zv_sum (:,:) + za2 * zwy (:,:) / e1v(:,:)599 un_adv(:,:) = un_adv(:,:) + za2 * zwx(:,:) * r1_e2u(:,:) 600 vn_adv(:,:) = vn_adv(:,:) + za2 * zwy(:,:) * r1_e1v(:,:) 606 601 ! 607 602 ! Set next sea level: … … 609 604 DO ji = fs_2, fs_jpim1 ! vector opt. 610 605 zhdiv(ji,jj) = ( zwx(ji,jj) - zwx(ji-1,jj) & 611 & + zwy(ji,jj) - zwy(ji,jj-1) ) * r1_e1 2t(ji,jj)606 & + zwy(ji,jj) - zwy(ji,jj-1) ) * r1_e1e2t(ji,jj) 612 607 END DO 613 608 END DO … … 627 622 DO jj = 2, jpjm1 628 623 DO ji = 2, jpim1 ! NO Vector Opt. 629 zsshu_a(ji,jj) = z1_2 * umask(ji,jj,1) * r1_e1 2u(ji,jj) &630 & * ( e1 2t(ji ,jj ) * ssha_e(ji ,jj ) &631 & + e1 2t(ji+1,jj ) * ssha_e(ji+1,jj ) )632 zsshv_a(ji,jj) = z1_2 * vmask(ji,jj,1) * r1_e1 2v(ji,jj) &633 & * ( e1 2t(ji ,jj ) * ssha_e(ji ,jj ) &634 & + e1 2t(ji ,jj+1) * ssha_e(ji ,jj+1) )624 zsshu_a(ji,jj) = z1_2 * umask(ji,jj,1) * r1_e1e2u(ji,jj) & 625 & * ( e1e2t(ji ,jj ) * ssha_e(ji ,jj ) & 626 & + e1e2t(ji+1,jj ) * ssha_e(ji+1,jj ) ) 627 zsshv_a(ji,jj) = z1_2 * vmask(ji,jj,1) * r1_e1e2v(ji,jj) & 628 & * ( e1e2t(ji ,jj ) * ssha_e(ji ,jj ) & 629 & + e1e2t(ji ,jj+1) * ssha_e(ji ,jj+1) ) 635 630 END DO 636 631 END DO … … 666 661 DO jj = 2, jpjm1 667 662 DO ji = 2, jpim1 668 zx1 = z1_2 * umask(ji ,jj,1) * r1_e1 2u(ji ,jj) &669 & * ( e1 2t(ji ,jj ) * zsshp2_e(ji ,jj) &670 & + e1 2t(ji+1,jj ) * zsshp2_e(ji+1,jj ) )671 zy1 = z1_2 * vmask(ji ,jj,1) * r1_e1 2v(ji ,jj ) &672 & * ( e1 2t(ji ,jj ) * zsshp2_e(ji ,jj ) &673 & + e1 2t(ji ,jj+1) * zsshp2_e(ji ,jj+1) )663 zx1 = z1_2 * umask(ji ,jj,1) * r1_e1e2u(ji ,jj) & 664 & * ( e1e2t(ji ,jj ) * zsshp2_e(ji ,jj) & 665 & + e1e2t(ji+1,jj ) * zsshp2_e(ji+1,jj ) ) 666 zy1 = z1_2 * vmask(ji ,jj,1) * r1_e1e2v(ji ,jj ) & 667 & * ( e1e2t(ji ,jj ) * zsshp2_e(ji ,jj ) & 668 & + e1e2t(ji ,jj+1) * zsshp2_e(ji ,jj+1) ) 674 669 zhust_e(ji,jj) = hu_0(ji,jj) + zx1 675 670 zhvst_e(ji,jj) = hv_0(ji,jj) + zy1 … … 688 683 DO jj = 2, jpjm1 689 684 DO ji = fs_2, fs_jpim1 ! vector opt. 690 zy1 = ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) ) /e1u(ji,jj)691 zy2 = ( zwy(ji ,jj ) + zwy(ji+1,jj ) ) /e1u(ji,jj)692 zx1 = ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) ) /e2v(ji,jj)693 zx2 = ( zwx(ji ,jj ) + zwx(ji ,jj+1) ) /e2v(ji,jj)685 zy1 = ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) ) * r1_e1u(ji,jj) 686 zy2 = ( zwy(ji ,jj ) + zwy(ji+1,jj ) ) * r1_e1u(ji,jj) 687 zx1 = ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) ) * r1_e2v(ji,jj) 688 zx2 = ( zwx(ji ,jj ) + zwx(ji ,jj+1) ) * r1_e2v(ji,jj) 694 689 zu_trd(ji,jj) = z1_4 * ( zwz(ji ,jj-1) * zy1 + zwz(ji,jj) * zy2 ) 695 690 zv_trd(ji,jj) =-z1_4 * ( zwz(ji-1,jj ) * zx1 + zwz(ji,jj) * zx2 ) … … 701 696 DO ji = fs_2, fs_jpim1 ! vector opt. 702 697 zy1 = z1_8 * ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) & 703 & + zwy(ji ,jj ) + zwy(ji+1,jj ) ) /e1u(ji,jj)698 & + zwy(ji ,jj ) + zwy(ji+1,jj ) ) * r1_e1u(ji,jj) 704 699 zx1 = - z1_8 * ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) & 705 & + zwx(ji ,jj ) + zwx(ji ,jj+1) ) /e2v(ji,jj)700 & + zwx(ji ,jj ) + zwx(ji ,jj+1) ) * r1_e2v(ji,jj) 706 701 zu_trd(ji,jj) = zy1 * ( zwz(ji ,jj-1) + zwz(ji,jj) ) 707 702 zv_trd(ji,jj) = zx1 * ( zwz(ji-1,jj ) + zwz(ji,jj) ) … … 709 704 END DO 710 705 ! 711 ELSEIF ( ln_dynvor_een .or. ln_dynvor_een_old) THEN !== energy and enstrophy conserving scheme ==!706 ELSEIF ( ln_dynvor_een ) THEN !== energy and enstrophy conserving scheme ==! 712 707 DO jj = 2, jpjm1 713 708 DO ji = fs_2, fs_jpim1 ! vector opt. 714 zu_trd(ji,jj) = + z1_12 /e1u(ji,jj) * ( ftne(ji,jj ) * zwy(ji ,jj ) &715 & + ftnw(ji+1,jj) * zwy(ji+1,jj ) &716 & + ftse(ji,jj ) * zwy(ji ,jj-1) &717 & + ftsw(ji+1,jj) * zwy(ji+1,jj-1) )718 zv_trd(ji,jj) = - z1_12 /e2v(ji,jj) * ( ftsw(ji,jj+1) * zwx(ji-1,jj+1) &719 & + ftse(ji,jj+1) * zwx(ji ,jj+1) &720 & + ftnw(ji,jj ) * zwx(ji-1,jj ) &721 & + ftne(ji,jj ) * zwx(ji ,jj ) )709 zu_trd(ji,jj) = + z1_12 * r1_e1u(ji,jj) * ( ftne(ji,jj ) * zwy(ji ,jj ) & 710 & + ftnw(ji+1,jj) * zwy(ji+1,jj ) & 711 & + ftse(ji,jj ) * zwy(ji ,jj-1) & 712 & + ftsw(ji+1,jj) * zwy(ji+1,jj-1) ) 713 zv_trd(ji,jj) = - z1_12 * r1_e2v(ji,jj) * ( ftsw(ji,jj+1) * zwx(ji-1,jj+1) & 714 & + ftse(ji,jj+1) * zwx(ji ,jj+1) & 715 & + ftnw(ji,jj ) * zwx(ji-1,jj ) & 716 & + ftne(ji,jj ) * zwx(ji ,jj ) ) 722 717 END DO 723 718 END DO … … 729 724 DO jj = 2, jpjm1 730 725 DO ji = fs_2, fs_jpim1 ! vector opt. 731 zu_spg = grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) /e1u(ji,jj)732 zv_spg = grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) /e2v(ji,jj)726 zu_spg = grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) * r1_e1u(ji,jj) 727 zv_spg = grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) * r1_e2v(ji,jj) 733 728 zu_trd(ji,jj) = zu_trd(ji,jj) + zu_spg 734 729 zv_trd(ji,jj) = zv_trd(ji,jj) + zv_spg … … 738 733 ! 739 734 ! Add bottom stresses: 740 zu_trd(:,:) = zu_trd(:,:) + bfrua(:,:) * zun_e(:,:) * hur_e(:,:)741 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(:,:) 742 737 ! 743 738 ! Surface pressure trend: … … 745 740 DO ji = fs_2, fs_jpim1 ! vector opt. 746 741 ! Add surface pressure gradient 747 zu_spg = - grav * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) /e1u(ji,jj)748 zv_spg = - grav * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) /e2v(ji,jj)742 zu_spg = - grav * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) * r1_e1u(ji,jj) 743 zv_spg = - grav * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) * r1_e2v(ji,jj) 749 744 zwx(ji,jj) = zu_spg 750 745 zwy(ji,jj) = zv_spg … … 756 751 DO jj = 2, jpjm1 757 752 DO ji = fs_2, fs_jpim1 ! vector opt. 758 ua_e(ji,jj) = ( zun_e(ji,jj) &753 ua_e(ji,jj) = ( un_e(ji,jj) & 759 754 & + rdtbt * ( zwx(ji,jj) & 760 755 & + zu_trd(ji,jj) & … … 762 757 & ) * umask(ji,jj,1) 763 758 764 va_e(ji,jj) = ( zvn_e(ji,jj) &759 va_e(ji,jj) = ( vn_e(ji,jj) & 765 760 & + rdtbt * ( zwy(ji,jj) & 766 761 & + zv_trd(ji,jj) & … … 777 772 zhvra = vmask(ji,jj,1)/(hv_0(ji,jj) + zsshv_a(ji,jj) + 1._wp - vmask(ji,jj,1)) 778 773 779 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) & 780 775 & + rdtbt * ( zhust_e(ji,jj) * zwx(ji,jj) & 781 776 & + zhup2_e(ji,jj) * zu_trd(ji,jj) & … … 783 778 & ) * zhura 784 779 785 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) & 786 781 & + rdtbt * ( zhvst_e(ji,jj) * zwy(ji,jj) & 787 782 & + zhvp2_e(ji,jj) * zv_trd(ji,jj) & … … 807 802 #if defined key_bdy 808 803 ! open boundaries 809 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 ) 810 805 #endif 811 806 #if defined key_agrif … … 815 810 ! ! ---- 816 811 ubb_e (:,:) = ub_e (:,:) 817 ub_e (:,:) = zun_e(:,:)818 zun_e(:,:) = ua_e (:,:)812 ub_e (:,:) = un_e (:,:) 813 un_e (:,:) = ua_e (:,:) 819 814 ! 820 815 vbb_e (:,:) = vb_e (:,:) 821 vb_e (:,:) = zvn_e(:,:)822 zvn_e(:,:) = va_e (:,:)816 vb_e (:,:) = vn_e (:,:) 817 vn_e (:,:) = va_e (:,:) 823 818 ! 824 819 sshbb_e(:,:) = sshb_e(:,:) … … 845 840 ! ----------------------------------------------------------------------------- 846 841 ! 847 ! At this stage ssha holds a time averaged value848 ! ! Sea Surface Height at u-,v- and f-points849 IF( lk_vvl ) THEN ! (required only in key_vvl case)850 DO jj = 1, jpjm1851 DO ji = 1, jpim1 ! NO Vector Opt.852 zsshu_a(ji,jj) = z1_2 * umask(ji,jj,1) * r1_e12u(ji,jj) &853 & * ( e12t(ji ,jj) * ssha(ji ,jj) &854 & + e12t(ji+1,jj) * ssha(ji+1,jj) )855 zsshv_a(ji,jj) = z1_2 * vmask(ji,jj,1) * r1_e12v(ji,jj) &856 & * ( e12t(ji,jj ) * ssha(ji,jj ) &857 & + e12t(ji,jj+1) * ssha(ji,jj+1) )858 END DO859 END DO860 CALL lbc_lnk_multi( zsshu_a, 'U', 1._wp, zsshv_a, 'V', 1._wp ) ! Boundary conditions861 ENDIF862 !863 842 ! Set advection velocity correction: 843 zwx(:,:) = un_adv(:,:) 844 zwy(:,:) = vn_adv(:,:) 864 845 IF (((kt==nit000).AND.(neuler==0)).OR.(.NOT.ln_bt_fw)) THEN 865 un_adv(:,:) = z u_sum(:,:)*hur(:,:)866 vn_adv(:,:) = z v_sum(:,:)*hvr(:,:)846 un_adv(:,:) = zwx(:,:)*hur(:,:) 847 vn_adv(:,:) = zwy(:,:)*hvr(:,:) 867 848 ELSE 868 un_adv(:,:) = z1_2 * ( ub2_b(:,:) + z u_sum(:,:)) * hur(:,:)869 vn_adv(:,:) = z1_2 * ( vb2_b(:,:) + z v_sum(:,:)) * hvr(:,:)849 un_adv(:,:) = z1_2 * ( ub2_b(:,:) + zwx(:,:)) * hur(:,:) 850 vn_adv(:,:) = z1_2 * ( vb2_b(:,:) + zwy(:,:)) * hvr(:,:) 870 851 END IF 871 852 872 853 IF (ln_bt_fw) THEN ! Save integrated transport for next computation 873 ub2_b(:,:) = z u_sum(:,:)874 vb2_b(:,:) = z v_sum(:,:)854 ub2_b(:,:) = zwx(:,:) 855 vb2_b(:,:) = zwy(:,:) 875 856 ENDIF 876 857 ! … … 882 863 END DO 883 864 ELSE 865 ! At this stage, ssha has been corrected: compute new depths at velocity points 866 DO jj = 1, jpjm1 867 DO ji = 1, jpim1 ! NO Vector Opt. 868 zsshu_a(ji,jj) = z1_2 * umask(ji,jj,1) * r1_e1e2u(ji,jj) & 869 & * ( e1e2t(ji ,jj) * ssha(ji ,jj) & 870 & + e1e2t(ji+1,jj) * ssha(ji+1,jj) ) 871 zsshv_a(ji,jj) = z1_2 * vmask(ji,jj,1) * r1_e1e2v(ji,jj) & 872 & * ( e1e2t(ji,jj ) * ssha(ji,jj ) & 873 & + e1e2t(ji,jj+1) * ssha(ji,jj+1) ) 874 END DO 875 END DO 876 CALL lbc_lnk_multi( zsshu_a, 'U', 1._wp, zsshv_a, 'V', 1._wp ) ! Boundary conditions 877 ! 884 878 DO jk=1,jpkm1 885 879 ua(:,:,jk) = ua(:,:,jk) + hur(:,:) * ( ua_b(:,:) - ub_b(:,:) * hu_b(:,:) ) * z1_2dt_b … … 900 894 #if defined key_agrif 901 895 ! Save time integrated fluxes during child grid integration 902 ! (used to update coarse grid transports) 903 ! Useless with 2nd order momentum schemes 896 ! (used to update coarse grid transports at next time step) 904 897 ! 905 898 IF ( (.NOT.Agrif_Root()).AND.(ln_bt_fw) ) THEN … … 921 914 ! 922 915 CALL wrk_dealloc( jpi, jpj, zsshp2_e, zhdiv ) 923 CALL wrk_dealloc( jpi, jpj, zu_trd, zv_trd , zun_e, zvn_e)924 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 ) 925 918 CALL wrk_dealloc( jpi, jpj, zhup2_e, zhvp2_e, zhust_e, zhvst_e ) 926 919 CALL wrk_dealloc( jpi, jpj, zsshu_a, zsshv_a ) … … 1070 1063 ! 1071 1064 INTEGER :: ji ,jj 1072 INTEGER :: ios ! Local integer output status for namelist read1073 1065 REAL(wp) :: zxr2, zyr2, zcmax 1074 1066 REAL(wp), POINTER, DIMENSION(:,:) :: zcu 1075 1067 !! 1076 NAMELIST/namsplit/ ln_bt_fw, ln_bt_av, ln_bt_nn_auto, &1077 & nn_baro, rn_bt_cmax, nn_bt_flt1078 1068 !!---------------------------------------------------------------------- 1079 1069 ! 1080 REWIND( numnam_ref ) ! Namelist namsplit in reference namelist : time splitting parameters 1081 READ ( numnam_ref, namsplit, IOSTAT = ios, ERR = 901) 1082 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsplit in reference namelist', lwp ) 1083 1084 REWIND( numnam_cfg ) ! Namelist namsplit in configuration namelist : time splitting parameters 1085 READ ( numnam_cfg, namsplit, IOSTAT = ios, ERR = 902 ) 1086 902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsplit in configuration namelist', lwp ) 1087 IF(lwm) WRITE ( numond, namsplit ) 1088 ! 1089 ! ! Max courant number for ext. grav. waves 1070 ! Max courant number for ext. grav. waves 1090 1071 ! 1091 1072 CALL wrk_alloc( jpi, jpj, zcu ) 1092 1073 ! 1093 IF (lk_vvl) THEN 1094 DO jj = 1, jpj 1095 DO ji =1, jpi 1096 zxr2 = 1./(e1t(ji,jj)*e1t(ji,jj)) 1097 zyr2 = 1./(e2t(ji,jj)*e2t(ji,jj)) 1098 zcu(ji,jj) = sqrt(grav*ht_0(ji,jj)*(zxr2 + zyr2) ) 1099 END DO 1100 END DO 1101 ELSE 1102 DO jj = 1, jpj 1103 DO ji =1, jpi 1104 zxr2 = 1./(e1t(ji,jj)*e1t(ji,jj)) 1105 zyr2 = 1./(e2t(ji,jj)*e2t(ji,jj)) 1106 zcu(ji,jj) = sqrt(grav*ht(ji,jj)*(zxr2 + zyr2) ) 1107 END DO 1108 END DO 1109 ENDIF 1110 1111 zcmax = MAXVAL(zcu(:,:)) 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 ! 1082 zcmax = MAXVAL( zcu(:,:) ) 1112 1083 IF( lk_mpp ) CALL mpp_max( zcmax ) 1113 1084 1114 1085 ! Estimate number of iterations to satisfy a max courant number= rn_bt_cmax 1115 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) 1116 1087 1117 rdtbt = rdt / FLOAT(nn_baro)1088 rdtbt = rdt / REAL( nn_baro , wp ) 1118 1089 zcmax = zcmax * rdtbt 1119 1090 ! Print results … … 1121 1092 IF(lwp) WRITE(numout,*) 'dyn_spg_ts : split-explicit free surface' 1122 1093 IF(lwp) WRITE(numout,*) '~~~~~~~~~~' 1123 IF( ln_bt_ nn_auto ) THEN1124 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 ' 1125 1096 IF(lwp) WRITE(numout,*) ' Max. courant number allowed: ', rn_bt_cmax 1126 1097 ELSE 1127 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 ' 1128 1099 ENDIF 1129 1100 … … 1170 1141 END SUBROUTINE dyn_spg_ts_init 1171 1142 1172 #else1173 !!---------------------------------------------------------------------------1174 !! Default case : Empty module No split explicit free surface1175 !!---------------------------------------------------------------------------1176 CONTAINS1177 INTEGER FUNCTION dyn_spg_ts_alloc() ! Dummy function1178 dyn_spg_ts_alloc = 01179 END FUNCTION dyn_spg_ts_alloc1180 SUBROUTINE dyn_spg_ts( kt ) ! Empty routine1181 INTEGER, INTENT(in) :: kt1182 WRITE(*,*) 'dyn_spg_ts: You should not have seen this print! error?', kt1183 END SUBROUTINE dyn_spg_ts1184 SUBROUTINE ts_rst( kt, cdrw ) ! Empty routine1185 INTEGER , INTENT(in) :: kt ! ocean time-step1186 CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag1187 WRITE(*,*) 'ts_rst : You should not have seen this print! error?', kt, cdrw1188 END SUBROUTINE ts_rst1189 SUBROUTINE dyn_spg_ts_init( kt ) ! Empty routine1190 INTEGER , INTENT(in) :: kt ! ocean time-step1191 WRITE(*,*) 'dyn_spg_ts_init : You should not have seen this print! error?', kt1192 END SUBROUTINE dyn_spg_ts_init1193 #endif1194 1195 1143 !!====================================================================== 1196 1144 END MODULE dynspg_ts 1197 1198 1199 -
branches/2014/dev_r4650_UKMO14.4_OBS_GENERAL_VINTERP/NEMOGCM/NEMO/OPA_SRC/DYN/dynvor.F90
r5950 r5951 16 16 !! 3.3 ! 2010-10 (C. Ethe, G. Madec) reorganisation of initialisation phase 17 17 !! 3.7 ! 2014-04 (G. Madec) trend simplification: suppress jpdyn_trd_dat vorticity 18 !! - ! 2014-06 (G. Madec) suppression of velocity curl from in-core memory 18 19 !!---------------------------------------------------------------------- 19 20 … … 22 23 !! vor_ens : enstrophy conserving scheme (ln_dynvor_ens=T) 23 24 !! vor_ene : energy conserving scheme (ln_dynvor_ene=T) 24 !! vor_mix : mixed enstrophy/energy conserving (ln_dynvor_mix=T)25 25 !! vor_een : energy and enstrophy conserving (ln_dynvor_een=T) 26 26 !! dyn_vor_init : set and control of the different vorticity option … … 32 32 USE trd_oce ! trends: ocean variables 33 33 USE trddyn ! trend manager: dynamics 34 USE c1d ! 1D vertical configuration 35 ! 34 36 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 35 37 USE prtctl ! Print control … … 44 46 45 47 PUBLIC dyn_vor ! routine called by step.F90 46 PUBLIC dyn_vor_init ! routine called by opa.F9048 PUBLIC dyn_vor_init ! routine called by nemogcm.F90 47 49 48 50 ! !!* Namelist namdyn_vor: vorticity term 49 LOGICAL, PUBLIC :: ln_dynvor_ene !: energy conserving scheme 50 LOGICAL, PUBLIC :: ln_dynvor_ens !: enstrophy conserving scheme 51 LOGICAL, PUBLIC :: ln_dynvor_mix !: mixed scheme 52 LOGICAL, PUBLIC :: ln_dynvor_een !: energy and enstrophy conserving scheme 53 LOGICAL, PUBLIC :: ln_dynvor_een_old !: energy and enstrophy conserving scheme (original formulation) 54 55 INTEGER :: nvor = 0 ! type of vorticity trend used 56 INTEGER :: ncor = 1 ! coriolis 57 INTEGER :: nrvm = 2 ! =2 relative vorticity ; =3 metric term 58 INTEGER :: ntot = 4 ! =4 total vorticity (relative + planetary) ; =5 coriolis + metric term 59 51 LOGICAL, PUBLIC :: ln_dynvor_ene !: energy conserving scheme (ENE) 52 LOGICAL, PUBLIC :: ln_dynvor_ens !: enstrophy conserving scheme (ENS) 53 LOGICAL, PUBLIC :: ln_dynvor_mix !: mixed scheme (MIX) 54 LOGICAL, PUBLIC :: ln_dynvor_een !: energy and enstrophy conserving scheme (EEN) 55 INTEGER, PUBLIC :: nn_een_e3f !: e3f=masked averaging of e3t divided by 4 (=0) or by the sum of mask (=1) 56 LOGICAL, PUBLIC :: ln_dynvor_msk !: vorticity multiplied by fmask (=T) or not (=F) (all vorticity schemes) 57 58 INTEGER :: nvor_scheme ! choice of the type of advection scheme 59 ! ! associated indices: 60 INTEGER, PUBLIC, PARAMETER :: np_ENE = 1 ! ENE scheme 61 INTEGER, PUBLIC, PARAMETER :: np_ENS = 2 ! ENS scheme 62 INTEGER, PUBLIC, PARAMETER :: np_MIX = 3 ! MIX scheme 63 INTEGER, PUBLIC, PARAMETER :: np_EEN = 4 ! EEN scheme 64 65 INTEGER :: ncor, nrvm, ntot ! choice of calculated vorticity 66 ! ! associated indices: 67 INTEGER, PARAMETER :: np_COR = 1 ! Coriolis (planetary) 68 INTEGER, PARAMETER :: np_RVO = 2 ! relative vorticity 69 INTEGER, PARAMETER :: np_MET = 3 ! metric term 70 INTEGER, PARAMETER :: np_CRV = 4 ! relative + planetary (total vorticity) 71 INTEGER, PARAMETER :: np_CME = 5 ! Coriolis + metric term 72 73 REAL(wp) :: r1_4 = 0.250_wp ! =1/4 74 REAL(wp) :: r1_8 = 0.125_wp ! =1/8 75 REAL(wp) :: r1_12 = 1._wp / 12._wp ! 1/12 76 60 77 !! * Substitutions 61 78 # include "domzgr_substitute.h90" 62 79 # include "vectopt_loop_substitute.h90" 63 80 !!---------------------------------------------------------------------- 64 !! NEMO/OPA 3. 3 , NEMO Consortium (2010)81 !! NEMO/OPA 3.7 , NEMO Consortium (2014) 65 82 !! $Id$ 66 83 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) … … 87 104 IF( l_trddyn ) CALL wrk_alloc( jpi,jpj,jpk, ztrdu, ztrdv ) 88 105 ! 89 ! ! vorticity term 90 SELECT CASE ( nvor ) ! compute the vorticity trend and add it to the general trend 91 ! 92 CASE ( -1 ) ! esopa: test all possibility with control print 93 CALL vor_ene( kt, ntot, ua, va ) 94 CALL prt_ctl( tab3d_1=ua, clinfo1=' vor0 - Ua: ', mask1=umask, & 95 & tab3d_2=va, clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' ) 96 CALL vor_ens( kt, ntot, ua, va ) 97 CALL prt_ctl( tab3d_1=ua, clinfo1=' vor1 - Ua: ', mask1=umask, & 98 & tab3d_2=va, clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' ) 99 CALL vor_mix( kt ) 100 CALL prt_ctl( tab3d_1=ua, clinfo1=' vor2 - Ua: ', mask1=umask, & 101 & tab3d_2=va, clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' ) 102 CALL vor_een( kt, ntot, ua, va ) 103 CALL prt_ctl( tab3d_1=ua, clinfo1=' vor3 - Ua: ', mask1=umask, & 104 & tab3d_2=va, clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' ) 105 ! 106 CASE ( 0 ) ! energy conserving scheme 107 IF( l_trddyn ) THEN 108 ztrdu(:,:,:) = ua(:,:,:) 109 ztrdv(:,:,:) = va(:,:,:) 110 CALL vor_ene( kt, nrvm, ua, va ) ! relative vorticity or metric trend 106 SELECT CASE ( nvor_scheme ) !== vorticity trend added to the general trend ==! 107 ! 108 CASE ( np_ENE ) !* energy conserving scheme 109 IF( l_trddyn ) THEN ! trend diagnostics: split the trend in two 110 ztrdu(:,:,:) = ua(:,:,:) 111 ztrdv(:,:,:) = va(:,:,:) 112 CALL vor_ene( kt, nrvm, ua, va ) ! relative vorticity or metric trend 111 113 ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 112 114 ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) … … 114 116 ztrdu(:,:,:) = ua(:,:,:) 115 117 ztrdv(:,:,:) = va(:,:,:) 116 CALL vor_ene( kt, ncor, ua, va ) ! planetary vorticity trend118 CALL vor_ene( kt, ncor, ua, va ) ! planetary vorticity trend 117 119 ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 118 120 ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 119 121 CALL trd_dyn( ztrdu, ztrdv, jpdyn_pvo, kt ) 120 122 ELSE 121 CALL vor_ene( kt, ntot, ua, va ) ! total vorticity 122 ENDIF 123 ! 124 CASE ( 1 ) !enstrophy conserving scheme125 IF( l_trddyn ) THEN126 ztrdu(:,:,:) = ua(:,:,:) 127 ztrdv(:,:,:) = va(:,:,:) 128 CALL vor_ens( kt, nrvm, ua, va ) ! relative vorticity or metric trend123 CALL vor_ene( kt, ntot, ua, va ) ! total vorticity trend 124 ENDIF 125 ! 126 CASE ( np_ENS ) !* enstrophy conserving scheme 127 IF( l_trddyn ) THEN ! trend diagnostics: splitthe trend in two 128 ztrdu(:,:,:) = ua(:,:,:) 129 ztrdv(:,:,:) = va(:,:,:) 130 CALL vor_ens( kt, nrvm, ua, va ) ! relative vorticity or metric trend 129 131 ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 130 132 ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) … … 132 134 ztrdu(:,:,:) = ua(:,:,:) 133 135 ztrdv(:,:,:) = va(:,:,:) 134 CALL vor_ens( kt, ncor, ua, va ) ! planetary vorticity trend136 CALL vor_ens( kt, ncor, ua, va ) ! planetary vorticity trend 135 137 ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 136 138 ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 137 139 CALL trd_dyn( ztrdu, ztrdv, jpdyn_pvo, kt ) 138 140 ELSE 139 CALL vor_ens( kt, ntot, ua, va ) ! total vorticity 140 ENDIF 141 ! 142 CASE ( 2 ) !mixed ene-ens scheme143 IF( l_trddyn ) THEN144 ztrdu(:,:,:) = ua(:,:,:) 145 ztrdv(:,:,:) = va(:,:,:) 146 CALL vor_ens( kt, nrvm, ua, va ) ! relative vorticity or metric trend (ens)141 CALL vor_ens( kt, ntot, ua, va ) ! total vorticity trend 142 ENDIF 143 ! 144 CASE ( np_MIX ) !* mixed ene-ens scheme 145 IF( l_trddyn ) THEN ! trend diagnostics: split the trend in two 146 ztrdu(:,:,:) = ua(:,:,:) 147 ztrdv(:,:,:) = va(:,:,:) 148 CALL vor_ens( kt, nrvm, ua, va ) ! relative vorticity or metric trend (ens) 147 149 ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 148 150 ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) … … 150 152 ztrdu(:,:,:) = ua(:,:,:) 151 153 ztrdv(:,:,:) = va(:,:,:) 152 CALL vor_ene( kt, ncor, ua, va ) ! planetary vorticity trend (ene)154 CALL vor_ene( kt, ncor, ua, va ) ! planetary vorticity trend (ene) 153 155 ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 154 156 ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 155 157 CALL trd_dyn( ztrdu, ztrdv, jpdyn_pvo, kt ) 156 158 ELSE 157 CALL vor_mix( kt ) ! total vorticity (mix=ens-ene) 158 ENDIF 159 ! 160 CASE ( 3 ) ! energy and enstrophy conserving scheme 161 IF( l_trddyn ) THEN 162 ztrdu(:,:,:) = ua(:,:,:) 163 ztrdv(:,:,:) = va(:,:,:) 164 CALL vor_een( kt, nrvm, ua, va ) ! relative vorticity or metric trend 159 CALL vor_ens( kt, nrvm, ua, va ) ! relative vorticity or metric trend (ens) 160 CALL vor_ene( kt, ncor, ua, va ) ! planetary vorticity trend (ene) 161 ENDIF 162 ! 163 CASE ( np_EEN ) !* energy and enstrophy conserving scheme 164 IF( l_trddyn ) THEN ! trend diagnostics: split the trend in two 165 ztrdu(:,:,:) = ua(:,:,:) 166 ztrdv(:,:,:) = va(:,:,:) 167 CALL vor_een( kt, nrvm, ua, va ) ! relative vorticity or metric trend 165 168 ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 166 169 ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) … … 168 171 ztrdu(:,:,:) = ua(:,:,:) 169 172 ztrdv(:,:,:) = va(:,:,:) 170 CALL vor_een( kt, ncor, ua, va ) ! planetary vorticity trend173 CALL vor_een( kt, ncor, ua, va ) ! planetary vorticity trend 171 174 ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 172 175 ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 173 176 CALL trd_dyn( ztrdu, ztrdv, jpdyn_pvo, kt ) 174 177 ELSE 175 CALL vor_een( kt, ntot, ua, va ) ! total vorticity 178 CALL vor_een( kt, ntot, ua, va ) ! total vorticity trend 176 179 ENDIF 177 180 ! … … 197 200 !! 198 201 !! ** Method : Trend evaluated using now fields (centered in time) 199 !! and the Sadourny (1975) flux form formulation : conserves the 200 !! horizontal kinetic energy. 201 !! The trend of the vorticity term is given by: 202 !! * s-coordinate (ln_sco=T), the e3. are inside the derivatives: 203 !! voru = 1/e1u mj-1[ (rotn+f)/e3f mi(e1v*e3v vn) ] 204 !! vorv = 1/e2v mi-1[ (rotn+f)/e3f mj(e2u*e3u un) ] 205 !! * z-coordinate (default key), e3t=e3u=e3v, the trend becomes: 206 !! voru = 1/e1u mj-1[ (rotn+f) mi(e1v vn) ] 207 !! vorv = 1/e2v mi-1[ (rotn+f) mj(e2u un) ] 208 !! Add this trend to the general momentum trend (ua,va): 209 !! (ua,va) = (ua,va) + ( voru , vorv ) 202 !! and the Sadourny (1975) flux form formulation : conserves the 203 !! horizontal kinetic energy. 204 !! The general trend of momentum is increased due to the vorticity 205 !! term which is given by: 206 !! voru = 1/e1u mj-1[ (rvor+f)/e3f mi(e1v*e3v vn) ] 207 !! vorv = 1/e2v mi-1[ (rvor+f)/e3f mj(e2u*e3u un) ] 208 !! where rvor is the relative vorticity 210 209 !! 211 210 !! ** Action : - Update (ua,va) with the now vorticity term trend … … 213 212 !! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689. 214 213 !!---------------------------------------------------------------------- 215 !216 214 INTEGER , INTENT(in ) :: kt ! ocean time-step index 217 215 INTEGER , INTENT(in ) :: kvor ! =ncor (planetary) ; =ntot (total) ; … … 220 218 REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: pva ! total v-trend 221 219 ! 222 INTEGER :: ji, jj, jk ! dummy loop indices223 REAL(wp) :: zx1, zy1, z fact2, zx2, zy2 ! local scalars224 REAL(wp), POINTER, DIMENSION(:,:) :: zwx, zwy, zwz220 INTEGER :: ji, jj, jk ! dummy loop indices 221 REAL(wp) :: zx1, zy1, zx2, zy2 ! local scalars 222 REAL(wp), POINTER, DIMENSION(:,:) :: zwx, zwy, zwz ! 2D workspace 225 223 !!---------------------------------------------------------------------- 226 224 ! … … 234 232 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 235 233 ENDIF 236 237 zfact2 = 0.5 * 0.5 ! Local constant initialization 238 239 !CDIR PARALLEL DO PRIVATE( zwx, zwy, zwz ) 234 ! 240 235 ! ! =============== 241 236 DO jk = 1, jpkm1 ! Horizontal slab 242 237 ! ! =============== 243 238 ! 244 ! Potential vorticity and horizontal fluxes 245 ! ----------------------------------------- 246 SELECT CASE( kvor ) ! vorticity considered 247 CASE ( 1 ) ; zwz(:,:) = ff(:,:) ! planetary vorticity (Coriolis) 248 CASE ( 2 ) ; zwz(:,:) = rotn(:,:,jk) ! relative vorticity 249 CASE ( 3 ) ! metric term 239 SELECT CASE( kvor ) !== vorticity considered ==! 240 CASE ( np_COR ) !* Coriolis (planetary vorticity) 241 zwz(:,:) = ff(:,:) 242 CASE ( np_RVO ) !* relative vorticity 243 DO jj = 1, jpjm1 244 DO ji = 1, fs_jpim1 ! vector opt. 245 zwz(ji,jj) = ( e2v(ji+1,jj ) * vn(ji+1,jj ,jk) - e2v(ji,jj) * vn(ji,jj,jk) & 246 & - e1u(ji ,jj+1) * un(ji ,jj+1,jk) + e1u(ji,jj) * un(ji,jj,jk) ) * r1_e1e2f(ji,jj) 247 END DO 248 END DO 249 CASE ( np_MET ) !* metric term 250 250 DO jj = 1, jpjm1 251 251 DO ji = 1, fs_jpim1 ! vector opt. 252 252 zwz(ji,jj) = ( ( vn(ji+1,jj ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) & 253 253 & - ( un(ji ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) ) & 254 & * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) ) 255 END DO 256 END DO 257 CASE ( 4 ) ; zwz(:,:) = ( rotn(:,:,jk) + ff(:,:) ) ! total (relative + planetary vorticity) 258 CASE ( 5 ) ! total (coriolis + metric) 259 DO jj = 1, jpjm1 260 DO ji = 1, fs_jpim1 ! vector opt. 261 zwz(ji,jj) = ( ff (ji,jj) & 262 & + ( ( vn(ji+1,jj ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) & 263 & - ( un(ji ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) ) & 264 & * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) ) & 265 & ) 266 END DO 267 END DO 254 & * 0.5 * r1_e1e2f(ji,jj) 255 END DO 256 END DO 257 CASE ( np_CRV ) !* Coriolis + relative vorticity 258 DO jj = 1, jpjm1 259 DO ji = 1, fs_jpim1 ! vector opt. 260 zwz(ji,jj) = ff(ji,jj) + ( e2v(ji+1,jj ) * vn(ji+1,jj ,jk) - e2v(ji,jj) * vn(ji,jj,jk) & 261 & - e1u(ji ,jj+1) * un(ji ,jj+1,jk) + e1u(ji,jj) * un(ji,jj,jk) ) & 262 & * r1_e1e2f(ji,jj) 263 END DO 264 END DO 265 CASE ( np_CME ) !* Coriolis + metric 266 DO jj = 1, jpjm1 267 DO ji = 1, fs_jpim1 ! vector opt. 268 zwz(ji,jj) = ff(ji,jj) & 269 & + ( ( vn(ji+1,jj ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) & 270 & - ( un(ji ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) ) & 271 & * 0.5 * r1_e1e2f(ji,jj) 272 END DO 273 END DO 274 CASE DEFAULT ! error 275 CALL ctl_stop('STOP','dyn_vor: wrong value for kvor' ) 268 276 END SELECT 277 ! 278 IF( ln_dynvor_msk ) THEN !== mask/unmask vorticity ==! 279 DO jj = 1, jpjm1 280 DO ji = 1, fs_jpim1 ! vector opt. 281 zwz(ji,jj) = zwz(ji,jj) * fmask(ji,jj,jk) 282 END DO 283 END DO 284 ENDIF 269 285 270 286 IF( ln_sco ) THEN … … 276 292 zwy(:,:) = e1v(:,:) * vn(:,:,jk) 277 293 ENDIF 278 279 ! Compute and add the vorticity term trend 280 ! ---------------------------------------- 294 ! !== compute and add the vorticity term trend =! 281 295 DO jj = 2, jpjm1 282 296 DO ji = fs_2, fs_jpim1 ! vector opt. … … 285 299 zx1 = zwx(ji-1,jj) + zwx(ji-1,jj+1) 286 300 zx2 = zwx(ji ,jj) + zwx(ji ,jj+1) 287 pua(ji,jj,jk) = pua(ji,jj,jk) + zfact2 /e1u(ji,jj) * ( zwz(ji ,jj-1) * zy1 + zwz(ji,jj) * zy2 )288 pva(ji,jj,jk) = pva(ji,jj,jk) - zfact2 /e2v(ji,jj) * ( zwz(ji-1,jj ) * zx1 + zwz(ji,jj) * zx2 )301 pua(ji,jj,jk) = pua(ji,jj,jk) + r1_4 * r1_e1u(ji,jj) * ( zwz(ji ,jj-1) * zy1 + zwz(ji,jj) * zy2 ) 302 pva(ji,jj,jk) = pva(ji,jj,jk) - r1_4 * r1_e2v(ji,jj) * ( zwz(ji-1,jj ) * zx1 + zwz(ji,jj) * zx2 ) 289 303 END DO 290 304 END DO … … 299 313 300 314 301 SUBROUTINE vor_mix( kt )302 !!----------------------------------------------------------------------303 !! *** ROUTINE vor_mix ***304 !!305 !! ** Purpose : Compute the now total vorticity trend and add it to306 !! the general trend of the momentum equation.307 !!308 !! ** Method : Trend evaluated using now fields (centered in time)309 !! Mixte formulation : conserves the potential enstrophy of a hori-310 !! zontally non-divergent flow for (rotzu x uh), the relative vor-311 !! ticity term and the horizontal kinetic energy for (f x uh), the312 !! coriolis term. the now trend of the vorticity term is given by:313 !! * s-coordinate (ln_sco=T), the e3. are inside the derivatives:314 !! voru = 1/e1u mj-1(rotn/e3f) mj-1[ mi(e1v*e3v vn) ]315 !! +1/e1u mj-1[ f/e3f mi(e1v*e3v vn) ]316 !! vorv = 1/e2v mi-1(rotn/e3f) mi-1[ mj(e2u*e3u un) ]317 !! +1/e2v mi-1[ f/e3f mj(e2u*e3u un) ]318 !! * z-coordinate (default key), e3t=e3u=e3v, the trend becomes:319 !! voru = 1/e1u mj-1(rotn) mj-1[ mi(e1v vn) ]320 !! +1/e1u mj-1[ f mi(e1v vn) ]321 !! vorv = 1/e2v mi-1(rotn) mi-1[ mj(e2u un) ]322 !! +1/e2v mi-1[ f mj(e2u un) ]323 !! Add this now trend to the general momentum trend (ua,va):324 !! (ua,va) = (ua,va) + ( voru , vorv )325 !!326 !! ** Action : - Update (ua,va) arrays with the now vorticity term trend327 !!328 !! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689.329 !!----------------------------------------------------------------------330 !331 INTEGER, INTENT(in) :: kt ! ocean timestep index332 !333 INTEGER :: ji, jj, jk ! dummy loop indices334 REAL(wp) :: zfact1, zua, zcua, zx1, zy1 ! local scalars335 REAL(wp) :: zfact2, zva, zcva, zx2, zy2 ! - -336 REAL(wp), POINTER, DIMENSION(:,:) :: zwx, zwy, zwz, zww337 !!----------------------------------------------------------------------338 !339 IF( nn_timing == 1 ) CALL timing_start('vor_mix')340 !341 CALL wrk_alloc( jpi, jpj, zwx, zwy, zwz, zww )342 !343 IF( kt == nit000 ) THEN344 IF(lwp) WRITE(numout,*)345 IF(lwp) WRITE(numout,*) 'dyn:vor_mix : vorticity term: mixed energy/enstrophy conserving scheme'346 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'347 ENDIF348 349 zfact1 = 0.5 * 0.25 ! Local constant initialization350 zfact2 = 0.5 * 0.5351 352 !CDIR PARALLEL DO PRIVATE( zwx, zwy, zwz, zww )353 ! ! ===============354 DO jk = 1, jpkm1 ! Horizontal slab355 ! ! ===============356 !357 ! Relative and planetary potential vorticity and horizontal fluxes358 ! ----------------------------------------------------------------359 IF( ln_sco ) THEN360 IF( ln_dynadv_vec ) THEN361 zww(:,:) = rotn(:,:,jk) / fse3f(:,:,jk)362 ELSE363 DO jj = 1, jpjm1364 DO ji = 1, fs_jpim1 ! vector opt.365 zww(ji,jj) = ( ( vn(ji+1,jj ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) &366 & - ( un(ji ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) ) &367 & * 0.5 / ( e1f(ji,jj) * e2f (ji,jj) * fse3f(ji,jj,jk) )368 END DO369 END DO370 ENDIF371 zwz(:,:) = ff (:,:) / fse3f(:,:,jk)372 zwx(:,:) = e2u(:,:) * fse3u(:,:,jk) * un(:,:,jk)373 zwy(:,:) = e1v(:,:) * fse3v(:,:,jk) * vn(:,:,jk)374 ELSE375 IF( ln_dynadv_vec ) THEN376 zww(:,:) = rotn(:,:,jk)377 ELSE378 DO jj = 1, jpjm1379 DO ji = 1, fs_jpim1 ! vector opt.380 zww(ji,jj) = ( ( vn(ji+1,jj ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) &381 & - ( un(ji ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) ) &382 & * 0.5 / ( e1f(ji,jj) * e2f (ji,jj) )383 END DO384 END DO385 ENDIF386 zwz(:,:) = ff (:,:)387 zwx(:,:) = e2u(:,:) * un(:,:,jk)388 zwy(:,:) = e1v(:,:) * vn(:,:,jk)389 ENDIF390 391 ! Compute and add the vorticity term trend392 ! ----------------------------------------393 DO jj = 2, jpjm1394 DO ji = fs_2, fs_jpim1 ! vector opt.395 zy1 = ( zwy(ji,jj-1) + zwy(ji+1,jj-1) ) / e1u(ji,jj)396 zy2 = ( zwy(ji,jj ) + zwy(ji+1,jj ) ) / e1u(ji,jj)397 zx1 = ( zwx(ji-1,jj) + zwx(ji-1,jj+1) ) / e2v(ji,jj)398 zx2 = ( zwx(ji ,jj) + zwx(ji ,jj+1) ) / e2v(ji,jj)399 ! enstrophy conserving formulation for relative vorticity term400 zua = zfact1 * ( zww(ji ,jj-1) + zww(ji,jj) ) * ( zy1 + zy2 )401 zva =-zfact1 * ( zww(ji-1,jj ) + zww(ji,jj) ) * ( zx1 + zx2 )402 ! energy conserving formulation for planetary vorticity term403 zcua = zfact2 * ( zwz(ji ,jj-1) * zy1 + zwz(ji,jj) * zy2 )404 zcva =-zfact2 * ( zwz(ji-1,jj ) * zx1 + zwz(ji,jj) * zx2 )405 ! mixed vorticity trend added to the momentum trends406 ua(ji,jj,jk) = ua(ji,jj,jk) + zcua + zua407 va(ji,jj,jk) = va(ji,jj,jk) + zcva + zva408 END DO409 END DO410 ! ! ===============411 END DO ! End of slab412 ! ! ===============413 CALL wrk_dealloc( jpi, jpj, zwx, zwy, zwz, zww )414 !415 IF( nn_timing == 1 ) CALL timing_stop('vor_mix')416 !417 END SUBROUTINE vor_mix418 419 420 315 SUBROUTINE vor_ens( kt, kvor, pua, pva ) 421 316 !!---------------------------------------------------------------------- … … 429 324 !! potential enstrophy of a horizontally non-divergent flow. the 430 325 !! trend of the vorticity term is given by: 431 !! * s-coordinate (ln_sco=T), the e3. are inside the derivative: 432 !! voru = 1/e1u mj-1[ (rotn+f)/e3f ] mj-1[ mi(e1v*e3v vn) ] 433 !! vorv = 1/e2v mi-1[ (rotn+f)/e3f ] mi-1[ mj(e2u*e3u un) ] 434 !! * z-coordinate (default key), e3t=e3u=e3v, the trend becomes: 435 !! voru = 1/e1u mj-1[ rotn+f ] mj-1[ mi(e1v vn) ] 436 !! vorv = 1/e2v mi-1[ rotn+f ] mi-1[ mj(e2u un) ] 326 !! voru = 1/e1u mj-1[ (rvor+f)/e3f ] mj-1[ mi(e1v*e3v vn) ] 327 !! vorv = 1/e2v mi-1[ (rvor+f)/e3f ] mi-1[ mj(e2u*e3u un) ] 437 328 !! Add this trend to the general momentum trend (ua,va): 438 329 !! (ua,va) = (ua,va) + ( voru , vorv ) … … 442 333 !! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689. 443 334 !!---------------------------------------------------------------------- 444 !445 335 INTEGER , INTENT(in ) :: kt ! ocean time-step index 446 336 INTEGER , INTENT(in ) :: kvor ! =ncor (planetary) ; =ntot (total) ; … … 449 339 REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: pva ! total v-trend 450 340 ! 451 INTEGER :: ji, jj, jk 452 REAL(wp) :: z fact1, zuav, zvau ! temporaryscalars453 REAL(wp), POINTER, DIMENSION(:,:) :: zwx, zwy, zwz, zww341 INTEGER :: ji, jj, jk ! dummy loop indices 342 REAL(wp) :: zuav, zvau ! local scalars 343 REAL(wp), POINTER, DIMENSION(:,:) :: zwx, zwy, zwz, zww ! 2D workspace 454 344 !!---------------------------------------------------------------------- 455 345 ! … … 463 353 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 464 354 ENDIF 465 466 zfact1 = 0.5 * 0.25 ! Local constant initialization467 468 !CDIR PARALLEL DO PRIVATE( zwx, zwy, zwz )469 355 ! ! =============== 470 356 DO jk = 1, jpkm1 ! Horizontal slab 471 357 ! ! =============== 472 358 ! 473 ! Potential vorticity and horizontal fluxes 474 ! ----------------------------------------- 475 SELECT CASE( kvor ) ! vorticity considered 476 CASE ( 1 ) ; zwz(:,:) = ff(:,:) ! planetary vorticity (Coriolis) 477 CASE ( 2 ) ; zwz(:,:) = rotn(:,:,jk) ! relative vorticity 478 CASE ( 3 ) ! metric term 359 SELECT CASE( kvor ) !== vorticity considered ==! 360 CASE ( np_COR ) !* Coriolis (planetary vorticity) 361 zwz(:,:) = ff(:,:) 362 CASE ( np_RVO ) !* relative vorticity 363 DO jj = 1, jpjm1 364 DO ji = 1, fs_jpim1 ! vector opt. 365 zwz(ji,jj) = ( e2v(ji+1,jj ) * vn(ji+1,jj ,jk) - e2v(ji,jj) * vn(ji,jj,jk) & 366 & - e1u(ji ,jj+1) * un(ji ,jj+1,jk) + e1u(ji,jj) * un(ji,jj,jk) ) * r1_e1e2f(ji,jj) 367 END DO 368 END DO 369 CASE ( np_MET ) !* metric term 479 370 DO jj = 1, jpjm1 480 371 DO ji = 1, fs_jpim1 ! vector opt. 481 372 zwz(ji,jj) = ( ( vn(ji+1,jj ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) & 482 373 & - ( un(ji ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) ) & 483 & * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) ) 484 END DO 485 END DO 486 CASE ( 4 ) ; zwz(:,:) = ( rotn(:,:,jk) + ff(:,:) ) ! total (relative + planetary vorticity) 487 CASE ( 5 ) ! total (coriolis + metric) 488 DO jj = 1, jpjm1 489 DO ji = 1, fs_jpim1 ! vector opt. 490 zwz(ji,jj) = ( ff (ji,jj) & 491 & + ( ( vn(ji+1,jj ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) & 492 & - ( un(ji ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) ) & 493 & * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) ) & 494 & ) 495 END DO 496 END DO 374 & * 0.5 * r1_e1e2f(ji,jj) 375 END DO 376 END DO 377 CASE ( np_CRV ) !* Coriolis + relative vorticity 378 DO jj = 1, jpjm1 379 DO ji = 1, fs_jpim1 ! vector opt. 380 zwz(ji,jj) = ff(ji,jj) + ( e2v(ji+1,jj ) * vn(ji+1,jj ,jk) - e2v(ji,jj) * vn(ji,jj,jk) & 381 & - e1u(ji ,jj+1) * un(ji ,jj+1,jk) + e1u(ji,jj) * un(ji,jj,jk) ) & 382 & * r1_e1e2f(ji,jj) 383 END DO 384 END DO 385 CASE ( np_CME ) !* Coriolis + metric 386 DO jj = 1, jpjm1 387 DO ji = 1, fs_jpim1 ! vector opt. 388 zwz(ji,jj) = ff(ji,jj) & 389 & + ( ( vn(ji+1,jj ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) & 390 & - ( un(ji ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) ) & 391 & * 0.5 * r1_e1e2f(ji,jj) 392 END DO 393 END DO 394 CASE DEFAULT ! error 395 CALL ctl_stop('STOP','dyn_vor: wrong value for kvor' ) 497 396 END SELECT 498 397 ! 499 IF( ln_sco ) THEN 500 DO jj = 1, jpj ! caution: don't use (:,:) for this loop 501 DO ji = 1, jpi ! it causes optimization problems on NEC in auto-tasking 502 zwz(ji,jj) = zwz(ji,jj) / fse3f(ji,jj,jk) 503 zwx(ji,jj) = e2u(ji,jj) * fse3u(ji,jj,jk) * un(ji,jj,jk) 504 zwy(ji,jj) = e1v(ji,jj) * fse3v(ji,jj,jk) * vn(ji,jj,jk) 505 END DO 506 END DO 398 IF( ln_dynvor_msk ) THEN !== mask/unmask vorticity ==! 399 DO jj = 1, jpjm1 400 DO ji = 1, fs_jpim1 ! vector opt. 401 zwz(ji,jj) = zwz(ji,jj) * fmask(ji,jj,jk) 402 END DO 403 END DO 404 ENDIF 405 ! 406 IF( ln_sco ) THEN !== horizontal fluxes ==! 407 zwz(:,:) = zwz(:,:) / fse3f(:,:,jk) 408 zwx(:,:) = e2u(:,:) * fse3u(:,:,jk) * un(:,:,jk) 409 zwy(:,:) = e1v(:,:) * fse3v(:,:,jk) * vn(:,:,jk) 507 410 ELSE 508 DO jj = 1, jpj ! caution: don't use (:,:) for this loop 509 DO ji = 1, jpi ! it causes optimization problems on NEC in auto-tasking 510 zwx(ji,jj) = e2u(ji,jj) * un(ji,jj,jk) 511 zwy(ji,jj) = e1v(ji,jj) * vn(ji,jj,jk) 512 END DO 513 END DO 514 ENDIF 515 ! 516 ! Compute and add the vorticity term trend 517 ! ---------------------------------------- 411 zwx(:,:) = e2u(:,:) * un(:,:,jk) 412 zwy(:,:) = e1v(:,:) * vn(:,:,jk) 413 ENDIF 414 ! !== compute and add the vorticity term trend =! 518 415 DO jj = 2, jpjm1 519 416 DO ji = fs_2, fs_jpim1 ! vector opt. 520 zuav = zfact1 /e1u(ji,jj) * ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) &521 & 522 zvau =- zfact1 /e2v(ji,jj) * ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) &523 & 417 zuav = r1_8 * r1_e1u(ji,jj) * ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) & 418 & + zwy(ji ,jj ) + zwy(ji+1,jj ) ) 419 zvau =-r1_8 * r1_e2v(ji,jj) * ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) & 420 & + zwx(ji ,jj ) + zwx(ji ,jj+1) ) 524 421 pua(ji,jj,jk) = pua(ji,jj,jk) + zuav * ( zwz(ji ,jj-1) + zwz(ji,jj) ) 525 422 pva(ji,jj,jk) = pva(ji,jj,jk) + zvau * ( zwz(ji-1,jj ) + zwz(ji,jj) ) … … 553 450 !! References : Arakawa and Lamb 1980, Mon. Wea. Rev., 109, 18-36 554 451 !!---------------------------------------------------------------------- 555 !556 452 INTEGER , INTENT(in ) :: kt ! ocean time-step index 557 INTEGER , INTENT(in ) :: kvor ! =ncor (planetary) ; =ntot (total) ; 558 ! ! =nrvm (relative vorticity or metric) 453 INTEGER , INTENT(in ) :: kvor ! =ncor (planetary) ; =ntot (total) ; =nrvm (relative or metric) 559 454 REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: pua ! total u-trend 560 455 REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: pva ! total v-trend 561 !! 562 INTEGER :: ji, jj, jk ! dummy loop indices 563 INTEGER :: ierr ! local integer 564 REAL(wp) :: zfac12, zua, zva ! local scalars 565 REAL(wp) :: zmsk, ze3 ! local scalars 566 ! ! 3D workspace 567 REAL(wp), POINTER , DIMENSION(:,: ) :: zwx, zwy, zwz 568 REAL(wp), POINTER , DIMENSION(:,: ) :: ztnw, ztne, ztsw, ztse 569 #if defined key_vvl 570 REAL(wp), POINTER , DIMENSION(:,:,:) :: ze3f ! 3D workspace (lk_vvl=T) 571 #else 572 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), SAVE :: ze3f ! lk_vvl=F, ze3f=1/e3f saved one for all 573 #endif 456 ! 457 INTEGER :: ji, jj, jk ! dummy loop indices 458 INTEGER :: ierr ! local integer 459 REAL(wp) :: zua, zva ! local scalars 460 REAL(wp) :: zmsk, ze3 ! local scalars 461 ! 462 REAL(wp), POINTER, DIMENSION(:,:) :: zwx, zwy, zwz, z1_e3f 463 REAL(wp), POINTER, DIMENSION(:,:) :: ztnw, ztne, ztsw, ztse 574 464 !!---------------------------------------------------------------------- 575 465 ! 576 466 IF( nn_timing == 1 ) CALL timing_start('vor_een') 577 467 ! 578 CALL wrk_alloc( jpi, jpj, zwx , zwy , zwz ) 579 CALL wrk_alloc( jpi, jpj, ztnw, ztne, ztsw, ztse ) 580 #if defined key_vvl 581 CALL wrk_alloc( jpi, jpj, jpk, ze3f ) 582 #endif 468 CALL wrk_alloc( jpi,jpj, zwx , zwy , zwz , z1_e3f ) 469 CALL wrk_alloc( jpi,jpj, ztnw, ztne, ztsw, ztse ) 583 470 ! 584 471 IF( kt == nit000 ) THEN … … 586 473 IF(lwp) WRITE(numout,*) 'dyn:vor_een : vorticity term: energy and enstrophy conserving scheme' 587 474 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 588 #if ! defined key_vvl589 IF( .NOT.ALLOCATED(ze3f) ) THEN590 ALLOCATE( ze3f(jpi,jpj,jpk) , STAT=ierr )591 IF( lk_mpp ) CALL mpp_sum ( ierr )592 IF( ierr /= 0 ) CALL ctl_stop( 'STOP', 'dyn:vor_een : unable to allocate arrays' )593 ENDIF594 ze3f(:,:,:) = 0._wp595 #endif596 475 ENDIF 597 598 IF( kt == nit000 .OR. lk_vvl ) THEN ! reciprocal of e3 at F-point (masked averaging of e3t over ocean points) 599 600 IF( ln_dynvor_een_old ) THEN ! original formulation 601 DO jk = 1, jpk 602 DO jj = 1, jpjm1 603 DO ji = 1, jpim1 604 ze3 = ( fse3t(ji,jj+1,jk)*tmask(ji,jj+1,jk) + fse3t(ji+1,jj+1,jk)*tmask(ji+1,jj+1,jk) & 605 & + fse3t(ji,jj ,jk)*tmask(ji,jj ,jk) + fse3t(ji+1,jj ,jk)*tmask(ji+1,jj ,jk) ) 606 IF( ze3 /= 0._wp ) ze3f(ji,jj,jk) = 4.0_wp / ze3 607 END DO 608 END DO 609 END DO 610 ELSE ! new formulation from NEMO 3.6 611 DO jk = 1, jpk 612 DO jj = 1, jpjm1 613 DO ji = 1, jpim1 614 ze3 = ( fse3t(ji,jj+1,jk)*tmask(ji,jj+1,jk) + fse3t(ji+1,jj+1,jk)*tmask(ji+1,jj+1,jk) & 615 & + fse3t(ji,jj ,jk)*tmask(ji,jj ,jk) + fse3t(ji+1,jj ,jk)*tmask(ji+1,jj ,jk) ) 616 zmsk = ( tmask(ji,jj+1,jk) + tmask(ji+1,jj+1,jk) & 617 & + tmask(ji,jj ,jk) + tmask(ji+1,jj ,jk) ) 618 IF( ze3 /= 0._wp ) ze3f(ji,jj,jk) = zmsk / ze3 619 END DO 620 END DO 621 END DO 622 ENDIF 623 624 CALL lbc_lnk( ze3f, 'F', 1. ) 625 ENDIF 626 627 zfac12 = 1._wp / 12._wp ! Local constant initialization 628 629 630 !CDIR PARALLEL DO PRIVATE( zwx, zwy, zwz, ztnw, ztne, ztsw, ztse ) 476 ! 631 477 ! ! =============== 632 478 DO jk = 1, jpkm1 ! Horizontal slab 633 479 ! ! =============== 634 635 ! Potential vorticity and horizontal fluxes 636 ! ----------------------------------------- 637 SELECT CASE( kvor ) ! vorticity considered 638 CASE ( 1 ) ! planetary vorticity (Coriolis) 639 zwz(:,:) = ff(:,:) * ze3f(:,:,jk) 640 CASE ( 2 ) ! relative vorticity 641 zwz(:,:) = rotn(:,:,jk) * ze3f(:,:,jk) 642 CASE ( 3 ) ! metric term 480 ! 481 SELECT CASE( nn_een_e3f ) ! == reciprocal of e3 at F-point 482 CASE ( 0 ) ! original formulation (masked averaging of e3t divided by 4) 483 DO jj = 1, jpjm1 484 DO ji = 1, fs_jpim1 ! vector opt. 485 ze3 = ( fse3t(ji,jj+1,jk)*tmask(ji,jj+1,jk) + fse3t(ji+1,jj+1,jk)*tmask(ji+1,jj+1,jk) & 486 & + fse3t(ji,jj ,jk)*tmask(ji,jj ,jk) + fse3t(ji+1,jj ,jk)*tmask(ji+1,jj ,jk) ) 487 IF( ze3 /= 0._wp ) THEN ; z1_e3f(ji,jj) = 4.0_wp / ze3 488 ELSE ; z1_e3f(ji,jj) = 0.0_wp 489 ENDIF 490 END DO 491 END DO 492 CASE ( 1 ) ! new formulation (masked averaging of e3t divided by the sum of mask) 493 DO jj = 1, jpjm1 494 DO ji = 1, fs_jpim1 ! vector opt. 495 ze3 = ( fse3t(ji,jj+1,jk)*tmask(ji,jj+1,jk) + fse3t(ji+1,jj+1,jk)*tmask(ji+1,jj+1,jk) & 496 & + fse3t(ji,jj ,jk)*tmask(ji,jj ,jk) + fse3t(ji+1,jj ,jk)*tmask(ji+1,jj ,jk) ) 497 zmsk = ( tmask(ji,jj+1,jk) + tmask(ji+1,jj+1,jk) & 498 & + tmask(ji,jj ,jk) + tmask(ji+1,jj ,jk) ) 499 IF( ze3 /= 0._wp ) THEN ; z1_e3f(ji,jj) = zmsk / ze3 500 ELSE ; z1_e3f(ji,jj) = 0.0_wp 501 ENDIF 502 END DO 503 END DO 504 END SELECT 505 ! 506 SELECT CASE( kvor ) !== vorticity considered ==! 507 CASE ( np_COR ) !* Coriolis (planetary vorticity) 508 DO jj = 1, jpjm1 509 DO ji = 1, fs_jpim1 ! vector opt. 510 zwz(ji,jj) = ff(ji,jj) * z1_e3f(ji,jj) 511 END DO 512 END DO 513 CASE ( np_RVO ) !* relative vorticity 514 DO jj = 1, jpjm1 515 DO ji = 1, fs_jpim1 ! vector opt. 516 zwz(ji,jj) = ( e2v(ji+1,jj ) * vn(ji+1,jj ,jk) - e2v(ji,jj) * vn(ji,jj,jk) & 517 & - e1u(ji ,jj+1) * un(ji ,jj+1,jk) + e1u(ji,jj) * un(ji,jj,jk) ) & 518 & * r1_e1e2f(ji,jj) * z1_e3f(ji,jj) 519 END DO 520 END DO 521 CASE ( np_MET ) !* metric term 643 522 DO jj = 1, jpjm1 644 523 DO ji = 1, fs_jpim1 ! vector opt. 645 524 zwz(ji,jj) = ( ( vn(ji+1,jj ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) & 646 525 & - ( un(ji ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) ) & 647 & * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) ) * ze3f(ji,jj,jk) 648 END DO 649 END DO 650 CALL lbc_lnk( zwz, 'F', 1. ) 651 CASE ( 4 ) ! total (relative + planetary vorticity) 652 zwz(:,:) = ( rotn(:,:,jk) + ff(:,:) ) * ze3f(:,:,jk) 653 CASE ( 5 ) ! total (coriolis + metric) 654 DO jj = 1, jpjm1 655 DO ji = 1, fs_jpim1 ! vector opt. 656 zwz(ji,jj) = ( ff (ji,jj) & 657 & + ( ( vn(ji+1,jj ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) & 658 & - ( un(ji ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) ) & 659 & * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) ) & 660 & ) * ze3f(ji,jj,jk) 661 END DO 662 END DO 663 CALL lbc_lnk( zwz, 'F', 1. ) 526 & * 0.5 * r1_e1e2f(ji,jj) * z1_e3f(ji,jj) 527 END DO 528 END DO 529 CASE ( np_CRV ) !* Coriolis + relative vorticity 530 DO jj = 1, jpjm1 531 DO ji = 1, fs_jpim1 ! vector opt. 532 zwz(ji,jj) = ( ff(ji,jj) + ( e2v(ji+1,jj ) * vn(ji+1,jj ,jk) - e2v(ji,jj) * vn(ji,jj,jk) & 533 & - e1u(ji ,jj+1) * un(ji ,jj+1,jk) + e1u(ji,jj) * un(ji,jj,jk) ) & 534 & * r1_e1e2f(ji,jj) ) * z1_e3f(ji,jj) 535 END DO 536 END DO 537 CASE ( np_CME ) !* Coriolis + metric 538 DO jj = 1, jpjm1 539 DO ji = 1, fs_jpim1 ! vector opt. 540 zwz(ji,jj) = ( ff(ji,jj) & 541 & + ( ( vn(ji+1,jj ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) & 542 & - ( un(ji ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) ) & 543 & * 0.5 * r1_e1e2f(ji,jj) ) * z1_e3f(ji,jj) 544 END DO 545 END DO 546 CASE DEFAULT ! error 547 CALL ctl_stop('STOP','dyn_vor: wrong value for kvor' ) 664 548 END SELECT 665 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 ! 558 CALL lbc_lnk( zwz, 'F', 1. ) 559 ! 560 ! !== horizontal fluxes ==! 666 561 zwx(:,:) = e2u(:,:) * fse3u(:,:,jk) * un(:,:,jk) 667 562 zwy(:,:) = e1v(:,:) * fse3v(:,:,jk) * vn(:,:,jk) 668 563 669 ! Compute and add the vorticity term trend 670 ! ---------------------------------------- 564 ! !== compute and add the vorticity term trend =! 671 565 jj = 2 672 566 ztne(1,:) = 0 ; ztnw(1,:) = 0 ; ztse(1,:) = 0 ; ztsw(1,:) = 0 673 DO ji = 2, jpi 567 DO ji = 2, jpi ! split in 2 parts due to vector opt. 674 568 ztne(ji,jj) = zwz(ji-1,jj ) + zwz(ji ,jj ) + zwz(ji ,jj-1) 675 569 ztnw(ji,jj) = zwz(ji-1,jj-1) + zwz(ji-1,jj ) + zwz(ji ,jj ) … … 687 581 DO jj = 2, jpjm1 688 582 DO ji = fs_2, fs_jpim1 ! vector opt. 689 zua = + zfac12 /e1u(ji,jj) * ( ztne(ji,jj ) * zwy(ji ,jj ) + ztnw(ji+1,jj) * zwy(ji+1,jj ) &690 & + ztse(ji,jj ) * zwy(ji ,jj-1) + ztsw(ji+1,jj) * zwy(ji+1,jj-1) )691 zva = - zfac12 /e2v(ji,jj) * ( ztsw(ji,jj+1) * zwx(ji-1,jj+1) + ztse(ji,jj+1) * zwx(ji ,jj+1) &692 & + ztnw(ji,jj ) * zwx(ji-1,jj ) + ztne(ji,jj ) * zwx(ji ,jj ) )583 zua = + r1_12 * r1_e1u(ji,jj) * ( ztne(ji,jj ) * zwy(ji ,jj ) + ztnw(ji+1,jj) * zwy(ji+1,jj ) & 584 & + ztse(ji,jj ) * zwy(ji ,jj-1) + ztsw(ji+1,jj) * zwy(ji+1,jj-1) ) 585 zva = - r1_12 * r1_e2v(ji,jj) * ( ztsw(ji,jj+1) * zwx(ji-1,jj+1) + ztse(ji,jj+1) * zwx(ji ,jj+1) & 586 & + ztnw(ji,jj ) * zwx(ji-1,jj ) + ztne(ji,jj ) * zwx(ji ,jj ) ) 693 587 pua(ji,jj,jk) = pua(ji,jj,jk) + zua 694 588 pva(ji,jj,jk) = pva(ji,jj,jk) + zva … … 698 592 END DO ! End of slab 699 593 ! ! =============== 700 CALL wrk_dealloc( jpi, jpj, zwx , zwy , zwz ) 701 CALL wrk_dealloc( jpi, jpj, ztnw, ztne, ztsw, ztse ) 702 #if defined key_vvl 703 CALL wrk_dealloc( jpi, jpj, jpk, ze3f ) 704 #endif 594 ! 595 CALL wrk_dealloc( jpi,jpj, zwx , zwy , zwz , z1_e3f ) 596 CALL wrk_dealloc( jpi,jpj, ztnw, ztne, ztsw, ztse ) 705 597 ! 706 598 IF( nn_timing == 1 ) CALL timing_stop('vor_een') … … 720 612 INTEGER :: ios ! Local integer output status for namelist read 721 613 !! 722 NAMELIST/namdyn_vor/ ln_dynvor_ens, ln_dynvor_ene, ln_dynvor_mix, ln_dynvor_een, ln_dynvor_een_old614 NAMELIST/namdyn_vor/ ln_dynvor_ens, ln_dynvor_ene, ln_dynvor_mix, ln_dynvor_een, nn_een_e3f, ln_dynvor_msk 723 615 !!---------------------------------------------------------------------- 724 616 … … 737 629 WRITE(numout,*) '~~~~~~~~~~~~' 738 630 WRITE(numout,*) ' Namelist namdyn_vor : choice of the vorticity term scheme' 739 WRITE(numout,*) ' energy conserving scheme ln_dynvor_ene = ', ln_dynvor_ene 740 WRITE(numout,*) ' enstrophy conserving scheme ln_dynvor_ens = ', ln_dynvor_ens 741 WRITE(numout,*) ' mixed enstrophy/energy conserving scheme ln_dynvor_mix = ', ln_dynvor_mix 742 WRITE(numout,*) ' enstrophy and energy conserving scheme ln_dynvor_een = ', ln_dynvor_een 743 WRITE(numout,*) ' enstrophy and energy conserving scheme (old) ln_dynvor_een_old= ', ln_dynvor_een_old 631 WRITE(numout,*) ' energy conserving scheme ln_dynvor_ene = ', ln_dynvor_ene 632 WRITE(numout,*) ' enstrophy conserving scheme ln_dynvor_ens = ', ln_dynvor_ens 633 WRITE(numout,*) ' mixed enstrophy/energy conserving scheme ln_dynvor_mix = ', ln_dynvor_mix 634 WRITE(numout,*) ' enstrophy and energy conserving scheme ln_dynvor_een = ', ln_dynvor_een 635 WRITE(numout,*) ' e3f = averaging /4 (=0) or /sum(tmask) (=1) nn_een_e3f = ', nn_een_e3f 636 WRITE(numout,*) ' masked (=1) or unmasked(=0) vorticity ln_dynvor_msk = ', ln_dynvor_msk 744 637 ENDIF 745 638 639 !!gm this should be removed when choosing a unique strategy for fmask at the coast 746 640 ! If energy, enstrophy or mixed advection of momentum in vector form change the value for masks 747 641 ! at angles with three ocean points and one land point 642 IF(lwp) WRITE(numout,*) 643 IF(lwp) WRITE(numout,*) ' namlbc: change fmask value in the angles (T) ln_vorlat = ', ln_vorlat 748 644 IF( ln_vorlat .AND. ( ln_dynvor_ene .OR. ln_dynvor_ens .OR. ln_dynvor_mix ) ) THEN 749 645 DO jk = 1, jpk … … 759 655 ! 760 656 ENDIF 761 762 ioptio = 0 ! Control of vorticity scheme options 763 IF( ln_dynvor_ene ) ioptio = ioptio + 1 764 IF( ln_dynvor_ens ) ioptio = ioptio + 1 765 IF( ln_dynvor_mix ) ioptio = ioptio + 1 766 IF( ln_dynvor_een ) ioptio = ioptio + 1 767 IF( ln_dynvor_een_old ) ioptio = ioptio + 1 768 IF( lk_esopa ) ioptio = 1 769 770 IF( ioptio /= 1 ) CALL ctl_stop( ' use ONE and ONLY one vorticity scheme' ) 771 772 ! ! Set nvor (type of scheme for vorticity) 773 IF( ln_dynvor_ene ) nvor = 0 774 IF( ln_dynvor_ens ) nvor = 1 775 IF( ln_dynvor_mix ) nvor = 2 776 IF( ln_dynvor_een .or. ln_dynvor_een_old ) nvor = 3 777 IF( lk_esopa ) nvor = -1 778 779 ! ! Set ncor, nrvm, ntot (type of vorticity) 780 IF(lwp) WRITE(numout,*) 781 ncor = 1 657 !!gm end 658 659 ioptio = 0 ! type of scheme for vorticity (set nvor_scheme) 660 IF( ln_dynvor_ene ) THEN ; ioptio = ioptio + 1 ; nvor_scheme = np_ENE ; ENDIF 661 IF( ln_dynvor_ens ) THEN ; ioptio = ioptio + 1 ; nvor_scheme = np_ENS ; ENDIF 662 IF( ln_dynvor_mix ) THEN ; ioptio = ioptio + 1 ; nvor_scheme = np_MIX ; ENDIF 663 IF( ln_dynvor_een ) THEN ; ioptio = ioptio + 1 ; nvor_scheme = np_EEN ; ENDIF 664 ! 665 IF( ( ioptio /= 1).AND.( .NOT.lk_c1d ) ) CALL ctl_stop( ' use ONE and ONLY one vorticity scheme' ) 666 ! 667 IF(lwp) WRITE(numout,*) ! type of calculated vorticity (set ncor, nrvm, ntot) 668 ncor = np_COR 782 669 IF( ln_dynadv_vec ) THEN 783 670 IF(lwp) WRITE(numout,*) ' Vector form advection : vorticity = Coriolis + relative vorticity' 784 nrvm = 2785 ntot = 4671 nrvm = np_RVO ! relative vorticity 672 ntot = np_CRV ! relative + planetary vorticity 786 673 ELSE 787 674 IF(lwp) WRITE(numout,*) ' Flux form advection : vorticity = Coriolis + metric term' 788 nrvm = 3789 ntot = 5675 nrvm = np_MET ! metric term 676 ntot = np_CME ! Coriolis + metric term 790 677 ENDIF 791 678 792 679 IF(lwp) THEN ! Print the choice 793 680 WRITE(numout,*) 794 IF( nvor == 0 ) WRITE(numout,*) ' vorticity scheme : energy conserving scheme' 795 IF( nvor == 1 ) WRITE(numout,*) ' vorticity scheme : enstrophy conserving scheme' 796 IF( nvor == 2 ) WRITE(numout,*) ' vorticity scheme : mixed enstrophy/energy conserving scheme' 797 IF( nvor == 3 ) WRITE(numout,*) ' vorticity scheme : energy and enstrophy conserving scheme' 798 IF( nvor == -1 ) WRITE(numout,*) ' esopa test: use all lateral physics options' 681 IF( nvor_scheme == np_ENE ) WRITE(numout,*) ' vorticity scheme ==>> energy conserving scheme' 682 IF( nvor_scheme == np_ENS ) WRITE(numout,*) ' vorticity scheme ==>> enstrophy conserving scheme' 683 IF( nvor_scheme == np_MIX ) WRITE(numout,*) ' vorticity scheme ==>> mixed enstrophy/energy conserving scheme' 684 IF( nvor_scheme == np_EEN ) WRITE(numout,*) ' vorticity scheme ==>> energy and enstrophy conserving scheme' 799 685 ENDIF 800 686 ! -
branches/2014/dev_r4650_UKMO14.4_OBS_GENERAL_VINTERP/NEMOGCM/NEMO/OPA_SRC/DYN/dynzad.F90
r5950 r5951 49 49 !! 50 50 !! ** Method : The now vertical advection of momentum is given by: 51 !! w dz(u) = ua + 1/(e1 u*e2u*e3u) mk+1[ mi(e1t*e2t*wn) dk(un) ]52 !! w dz(v) = va + 1/(e1 v*e2v*e3v) mk+1[ mj(e1t*e2t*wn) dk(vn) ]51 !! w dz(u) = ua + 1/(e1e2u*e3u) mk+1[ mi(e1e2t*wn) dk(un) ] 52 !! w dz(v) = va + 1/(e1e2v*e3v) mk+1[ mj(e1e2t*wn) dk(vn) ] 53 53 !! Add this trend to the general trend (ua,va): 54 54 !! (ua,va) = (ua,va) + w dz(u,v) … … 85 85 DO jj = 2, jpj ! vertical fluxes 86 86 DO ji = fs_2, jpi ! vector opt. 87 zww(ji,jj) = 0.25_wp * e1 t(ji,jj) *e2t(ji,jj) * wn(ji,jj,jk)87 zww(ji,jj) = 0.25_wp * e1e2t(ji,jj) * wn(ji,jj,jk) 88 88 END DO 89 89 END DO … … 121 121 DO ji = fs_2, fs_jpim1 ! vector opt. 122 122 ! ! vertical momentum advective trends 123 zua = - ( zwuw(ji,jj,jk) + zwuw(ji,jj,jk+1) ) / ( e1 u(ji,jj) *e2u(ji,jj) * fse3u(ji,jj,jk) )124 zva = - ( zwvw(ji,jj,jk) + zwvw(ji,jj,jk+1) ) / ( e1 v(ji,jj) *e2v(ji,jj) * fse3v(ji,jj,jk) )123 zua = - ( zwuw(ji,jj,jk) + zwuw(ji,jj,jk+1) ) / ( e1e2u(ji,jj) * fse3u(ji,jj,jk) ) 124 zva = - ( zwvw(ji,jj,jk) + zwvw(ji,jj,jk+1) ) / ( e1e2v(ji,jj) * fse3v(ji,jj,jk) ) 125 125 ! ! add the trends to the general momentum trends 126 126 ua(ji,jj,jk) = ua(ji,jj,jk) + zua … … 146 146 ! 147 147 END SUBROUTINE dyn_zad 148 148 149 149 150 SUBROUTINE dyn_zad_zts ( kt ) … … 182 183 IF( nn_timing == 1 ) CALL timing_start('dyn_zad_zts') 183 184 ! 184 CALL wrk_alloc( jpi,jpj,jpk, zwuw, zwvw, zww )185 CALL wrk_alloc( jpi,jpj,jpk,3, zus, zvs )185 CALL wrk_alloc( jpi,jpj,jpk, zwuw, zwvw, zww ) 186 CALL wrk_alloc( jpi,jpj,jpk,3, zus , zvs ) 186 187 ! 187 188 IF( kt == nit000 ) THEN … … 205 206 DO jj = 2, jpj 206 207 DO ji = fs_2, jpi ! vector opt. 207 zww(ji,jj,jk) = 0.25_wp * e1 t(ji,jj) *e2t(ji,jj) * wn(ji,jj,jk)208 zww(ji,jj,jk) = 0.25_wp * e1e2t(ji,jj) * wn(ji,jj,jk) 208 209 END DO 209 210 END DO 210 211 END DO 211 ! 212 ! Surface and bottom advective fluxes set to zero 213 DO jj = 2, jpjm1 212 213 DO jj = 2, jpjm1 ! Surface and bottom advective fluxes set to zero 214 214 DO ji = fs_2, fs_jpim1 ! vector opt. 215 !!gm missing ISF boundary condition 215 216 zwuw(ji,jj, 1 ) = 0._wp 216 217 zwvw(ji,jj, 1 ) = 0._wp … … 251 252 DO ji = fs_2, fs_jpim1 ! vector opt. 252 253 ! ! vertical momentum advective trends 253 zua = - ( zwuw(ji,jj,jk) + zwuw(ji,jj,jk+1) ) / ( e1 u(ji,jj) *e2u(ji,jj) * fse3u(ji,jj,jk) )254 zva = - ( zwvw(ji,jj,jk) + zwvw(ji,jj,jk+1) ) / ( e1 v(ji,jj) *e2v(ji,jj) * fse3v(ji,jj,jk) )254 zua = - ( zwuw(ji,jj,jk) + zwuw(ji,jj,jk+1) ) / ( e1e2u(ji,jj) * fse3u(ji,jj,jk) ) 255 zva = - ( zwvw(ji,jj,jk) + zwvw(ji,jj,jk+1) ) / ( e1e2v(ji,jj) * fse3v(ji,jj,jk) ) 255 256 zus(ji,jj,jk,jta) = zus(ji,jj,jk,jtb) + zua * zts 256 257 zvs(ji,jj,jk,jta) = zvs(ji,jj,jk,jtb) + zva * zts … … 283 284 & tab3d_2=va, clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' ) 284 285 ! 285 CALL wrk_dealloc( jpi,jpj,jpk, zwuw, zwvw, zww )286 CALL wrk_dealloc( jpi,jpj,jpk,3, zus, zvs )286 CALL wrk_dealloc( jpi,jpj,jpk, zwuw, zwvw, zww ) 287 CALL wrk_dealloc( jpi,jpj,jpk,3, zus , zvs ) 287 288 ! 288 289 IF( nn_timing == 1 ) CALL timing_stop('dyn_zad_zts') -
branches/2014/dev_r4650_UKMO14.4_OBS_GENERAL_VINTERP/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf.F90
r5950 r5951 9 9 10 10 !!---------------------------------------------------------------------- 11 !! dyn_zdf : Update the momentum trend with the vertical diffusion12 !! dyn_zdf_init : initializations of the vertical diffusion scheme11 !! dyn_zdf : Update the momentum trend with the vertical diffusion 12 !! dyn_zdf_init : initializations of the vertical diffusion scheme 13 13 !!---------------------------------------------------------------------- 14 USE oce ! ocean dynamics and tracers variables 15 USE dom_oce ! ocean space and time domain variables 16 USE zdf_oce ! ocean vertical physics variables 17 18 USE dynzdf_exp ! vertical diffusion: explicit (dyn_zdf_exp routine) 19 USE dynzdf_imp ! vertical diffusion: implicit (dyn_zdf_imp routine) 20 21 USE ldfdyn_oce ! ocean dynamics: lateral physics 22 USE trd_oce ! trends: ocean variables 23 USE trddyn ! trend manager: dynamics 24 USE in_out_manager ! I/O manager 25 USE lib_mpp ! MPP library 26 USE prtctl ! Print control 27 USE wrk_nemo ! Memory Allocation 28 USE timing ! Timing 14 USE oce ! ocean dynamics and tracers variables 15 USE dom_oce ! ocean space and time domain variables 16 USE zdf_oce ! ocean vertical physics variables 17 USE dynzdf_exp ! vertical diffusion: explicit (dyn_zdf_exp routine) 18 USE dynzdf_imp ! vertical diffusion: implicit (dyn_zdf_imp routine) 19 USE ldfdyn ! lateral diffusion: eddy viscosity coef. 20 USE trd_oce ! trends: ocean variables 21 USE trddyn ! trend manager: dynamics 22 ! 23 USE in_out_manager ! I/O manager 24 USE lib_mpp ! MPP library 25 USE prtctl ! Print control 26 USE wrk_nemo ! Memory Allocation 27 USE timing ! Timing 29 28 30 29 IMPLICIT NONE … … 61 60 !!--------------------------------------------------------------------- 62 61 ! 63 IF( nn_timing == 1 ) CALL timing_start('dyn_zdf')62 IF( nn_timing == 1 ) CALL timing_start('dyn_zdf') 64 63 ! 65 64 ! ! set time step … … 79 78 CASE ( 1 ) ; CALL dyn_zdf_imp( kt, r2dt ) ! implicit scheme 80 79 ! 81 CASE ( -1 ) ! esopa: test all possibility with control print82 CALL dyn_zdf_exp( kt, r2dt )83 CALL prt_ctl( tab3d_1=ua, clinfo1=' zdf0 - Ua: ', mask1=umask, &84 & tab3d_2=va, clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' )85 CALL dyn_zdf_imp( kt, r2dt )86 CALL prt_ctl( tab3d_1=ua, clinfo1=' zdf1 - Ua: ', mask1=umask, &87 & tab3d_2=va, clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' )88 80 END SELECT 89 81 … … 96 88 ! ! print mean trends (used for debugging) 97 89 IF(ln_ctl) CALL prt_ctl( tab3d_1=ua, clinfo1=' zdf - Ua: ', mask1=umask, & 98 &tab3d_2=va, clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' )99 !100 IF( nn_timing == 1 ) CALL timing_stop('dyn_zdf')90 & tab3d_2=va, clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' ) 91 ! 92 IF( nn_timing == 1 ) CALL timing_stop('dyn_zdf') 101 93 ! 102 94 END SUBROUTINE dyn_zdf … … 114 106 USE zdftke 115 107 USE zdfgls 116 USE zdfkpp117 108 !!---------------------------------------------------------------------- 118 109 ! … … 123 114 ! 124 115 ! Force implicit schemes 125 IF( lk_zdftke .OR. lk_zdfgls .OR. lk_zdfkpp ) nzdf = 1 ! TKE, GLS or KPP physics 126 IF( ln_dynldf_iso ) nzdf = 1 ! iso-neutral lateral physics 127 IF( ln_dynldf_hor .AND. ln_sco ) nzdf = 1 ! horizontal lateral physics in s-coordinate 128 ! 129 IF( lk_esopa ) nzdf = -1 ! Esopa key: All schemes used 116 IF( lk_zdftke .OR. lk_zdfgls ) nzdf = 1 ! TKE or GLS physics 117 IF( ln_dynldf_iso ) nzdf = 1 ! iso-neutral lateral physics 118 IF( ln_dynldf_hor .AND. ln_sco ) nzdf = 1 ! horizontal lateral physics in s-coordinate 130 119 ! 131 120 IF(lwp) THEN ! Print the choice … … 133 122 WRITE(numout,*) 'dyn_zdf_init : vertical dynamics physics scheme' 134 123 WRITE(numout,*) '~~~~~~~~~~~' 135 IF( nzdf == -1 ) WRITE(numout,*) ' ESOPA test All scheme used'136 124 IF( nzdf == 0 ) WRITE(numout,*) ' Explicit time-splitting scheme' 137 125 IF( nzdf == 1 ) WRITE(numout,*) ' Implicit (euler backward) scheme' -
branches/2014/dev_r4650_UKMO14.4_OBS_GENERAL_VINTERP/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf_exp.F90
r5950 r5951 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_r4650_UKMO14.4_OBS_GENERAL_VINTERP/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf_imp.F90
r5950 r5951 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 … … 209 208 !----------------------------------------------------------------------- 210 209 ! 211 !== First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 (increasing k) == 212 DO jk = 2, jpkm1 210 DO jk = 2, jpkm1 !== First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 (increasing k) == 213 211 DO jj = 2, jpjm1 214 212 DO ji = fs_2, fs_jpim1 ! vector opt. … … 220 218 DO jj = 2, jpjm1 !== second recurrence: SOLk = RHSk - Lk / Dk-1 Lk-1 == 221 219 DO ji = fs_2, fs_jpim1 ! vector opt. 222 #if defined key_dynspg_ts223 220 ze3ua = ( 1._wp - r_vvl ) * fse3u_n(ji,jj,1) + r_vvl * fse3u_a(ji,jj,1) 224 221 ua(ji,jj,1) = ua(ji,jj,1) + p2dt * 0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) ) & 225 222 & / ( ze3ua * rau0 ) * umask(ji,jj,1) 226 #else227 ua(ji,jj,1) = ub(ji,jj,1) &228 & + p2dt *(ua(ji,jj,1) + 0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) ) &229 & / ( fse3u(ji,jj,1) * rau0 ) * umask(ji,jj,1) )230 #endif231 223 END DO 232 224 END DO … … 234 226 DO jj = 2, jpjm1 235 227 DO ji = fs_2, fs_jpim1 236 #if defined key_dynspg_ts237 228 zrhs = ua(ji,jj,jk) ! zrhs=right hand side 238 #else239 zrhs = ub(ji,jj,jk) + p2dt * ua(ji,jj,jk)240 #endif241 229 ua(ji,jj,jk) = zrhs - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * ua(ji,jj,jk-1) 242 230 END DO … … 257 245 END DO 258 246 259 #if ! defined key_dynspg_ts 260 ! Normalization to obtain the general momentum trend ua 261 DO jk = 1, jpkm1 262 DO jj = 2, jpjm1 263 DO ji = fs_2, fs_jpim1 ! vector opt. 264 ua(ji,jj,jk) = ( ua(ji,jj,jk) - ub(ji,jj,jk) ) * z1_p2dt 265 END DO 266 END DO 267 END DO 268 #endif 269 270 ! 3. Vertical diffusion on v 247 ! 4. Vertical diffusion on v 271 248 ! --------------------------- 272 249 ! Matrix and second member construction … … 309 286 !----------------------------------------------------------------------- 310 287 ! 311 !== First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 (increasing k) == 312 DO jk = 2, jpkm1 288 DO jk = 2, jpkm1 !== First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 (increasing k) == 313 289 DO jj = 2, jpjm1 314 290 DO ji = fs_2, fs_jpim1 ! vector opt. … … 319 295 ! 320 296 DO jj = 2, jpjm1 !== second recurrence: SOLk = RHSk - Lk / Dk-1 Lk-1 == 321 DO ji = fs_2, fs_jpim1 ! vector opt. 322 #if defined key_dynspg_ts 297 DO ji = fs_2, fs_jpim1 ! vector opt. 323 298 ze3va = ( 1._wp - r_vvl ) * fse3v_n(ji,jj,1) + r_vvl * fse3v_a(ji,jj,1) 324 299 va(ji,jj,1) = va(ji,jj,1) + p2dt * 0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) ) & 325 300 & / ( ze3va * rau0 ) 326 #else327 va(ji,jj,1) = vb(ji,jj,1) &328 & + p2dt *(va(ji,jj,1) + 0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) ) &329 & / ( fse3v(ji,jj,1) * rau0 ) )330 #endif331 301 END DO 332 302 END DO … … 334 304 DO jj = 2, jpjm1 335 305 DO ji = fs_2, fs_jpim1 ! vector opt. 336 #if defined key_dynspg_ts337 306 zrhs = va(ji,jj,jk) ! zrhs=right hand side 338 #else339 zrhs = vb(ji,jj,jk) + p2dt * va(ji,jj,jk)340 #endif341 307 va(ji,jj,jk) = zrhs - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * va(ji,jj,jk-1) 342 308 END DO … … 357 323 END DO 358 324 359 ! Normalization to obtain the general momentum trend va360 #if ! defined key_dynspg_ts361 DO jk = 1, jpkm1362 DO jj = 2, jpjm1363 DO ji = fs_2, fs_jpim1 ! vector opt.364 va(ji,jj,jk) = ( va(ji,jj,jk) - vb(ji,jj,jk) ) * z1_p2dt365 END DO366 END DO367 END DO368 #endif369 370 325 ! J. Chanut: Lines below are useless ? 371 !! restore bottom layer avmu(v)326 !! restore avmu(v)=0. at bottom (and top if ln_isfcav=T interfaces) 372 327 IF( ln_bfrimp ) THEN 373 328 DO jj = 2, jpjm1 -
branches/2014/dev_r4650_UKMO14.4_OBS_GENERAL_VINTERP/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90
r5950 r5951 20 20 USE sbc_oce ! surface boundary condition: ocean 21 21 USE domvvl ! Variable volume 22 USE divcur ! hor. divergence and curl (div & cur routines) 23 USE restart ! only for lrst_oce 24 USE in_out_manager ! I/O manager 25 USE prtctl ! Print control 26 USE phycst 27 USE lbclnk ! ocean lateral boundary condition (or mpp link) 28 USE lib_mpp ! MPP library 22 USE divhor ! horizontal divergence 23 USE phycst ! physical constants 29 24 USE bdy_oce 30 25 USE bdy_par 31 26 USE bdydyn2d ! bdy_ssh routine 32 27 #if defined key_agrif 33 USE agrif_opa_update34 28 USE agrif_opa_interp 35 29 #endif … … 37 31 USE asminc ! Assimilation increment 38 32 #endif 33 USE in_out_manager ! I/O manager 34 USE restart ! only for lrst_oce 35 USE prtctl ! Print control 36 USE lbclnk ! ocean lateral boundary condition (or mpp link) 37 USE lib_mpp ! MPP library 39 38 USE wrk_nemo ! Memory Allocation 40 39 USE timing ! Timing … … 67 66 !! by the time step. 68 67 !! 69 !! ** action : ssha :after sea surface height68 !! ** action : ssha, after sea surface height 70 69 !! 71 70 !! Reference : Leclair, M., and G. Madec, 2009, Ocean Modelling. 72 71 !!---------------------------------------------------------------------- 73 ! 74 REAL(wp), POINTER, DIMENSION(:,: ) :: zhdiv 75 INTEGER, INTENT(in) :: kt ! time step 72 INTEGER, INTENT(in) :: kt ! time step 76 73 ! 77 INTEGER :: jk ! dummy loop indice 78 REAL(wp) :: z2dt, z1_rau0 ! local scalars 79 !!---------------------------------------------------------------------- 80 ! 81 IF( nn_timing == 1 ) CALL timing_start('ssh_nxt') 82 ! 83 CALL wrk_alloc( jpi, jpj, zhdiv ) 74 INTEGER :: jk ! dummy loop indice 75 REAL(wp) :: z2dt, zcoef ! local scalars 76 REAL(wp), POINTER, DIMENSION(:,: ) :: zhdiv ! 2D workspace 77 !!---------------------------------------------------------------------- 78 ! 79 IF( nn_timing == 1 ) CALL timing_start('ssh_nxt') 80 ! 81 CALL wrk_alloc( jpi,jpj, zhdiv ) 84 82 ! 85 83 IF( kt == nit000 ) THEN 86 !87 84 IF(lwp) WRITE(numout,*) 88 85 IF(lwp) WRITE(numout,*) 'ssh_nxt : after sea surface height' 89 86 IF(lwp) WRITE(numout,*) '~~~~~~~ ' 90 ! 91 ENDIF 92 ! 93 CALL div_cur( kt ) ! Horizontal divergence & Relative vorticity 87 ENDIF 88 ! 89 CALL div_hor( kt ) ! Horizontal divergence 94 90 ! 95 91 z2dt = 2._wp * rdt ! set time step size (Euler/Leapfrog) … … 107 103 ! compute the vertical velocity which can be used to compute the non-linear terms of the momentum equations. 108 104 ! 109 z 1_rau0= 0.5_wp * r1_rau0110 ssha(:,:) = ( sshb(:,:) - z2dt * ( z 1_rau0* ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) ) ) * ssmask(:,:)111 112 #if ! defined key_dynspg_ts 113 ! These lines are not necessary with time splitting since114 ! boundary condition on sea level is set during ts loop115 # if defined key_agrif116 CALL agrif_ssh( kt )117 # endif118 # if defined key_bdy119 IF (lk_bdy) THEN120 CALL lbc_lnk( ssha, 'T', 1. )! Not sure that's necessary121 CALL bdy_ssh( ssha )! Duplicate sea level across open boundaries122 ENDIF123 # endif124 #endif 105 zcoef = 0.5_wp * r1_rau0 106 ssha(:,:) = ( sshb(:,:) - z2dt * ( zcoef * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) ) ) * ssmask(:,:) 107 108 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 # if defined key_agrif 112 CALL agrif_ssh( kt ) 113 # endif 114 # if defined key_bdy 115 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 # endif 120 ENDIF 125 121 126 122 #if defined key_asminc 127 ! ! Include the IAU weighted SSH increment 128 IF( lk_asminc .AND. ln_sshinc .AND. ln_asmiau ) THEN 123 IF( lk_asminc .AND. ln_sshinc .AND. ln_asmiau ) THEN ! Include the IAU weighted SSH increment 129 124 CALL ssh_asm_inc( kt ) 130 125 ssha(:,:) = ssha(:,:) + z2dt * ssh_iau(:,:) 131 126 ENDIF 132 127 #endif 133 134 128 ! !------------------------------! 135 129 ! ! outputs ! … … 160 154 !! Reference : Leclair, M., and G. Madec, 2009, Ocean Modelling. 161 155 !!---------------------------------------------------------------------- 162 ! 163 INTEGER, INTENT(in) :: kt ! time step 156 INTEGER, INTENT(in) :: kt ! time step 157 ! 158 INTEGER :: ji, jj, jk ! dummy loop indices 159 REAL(wp) :: z1_2dt ! local scalars 164 160 REAL(wp), POINTER, DIMENSION(:,: ) :: z2d 165 161 REAL(wp), POINTER, DIMENSION(:,:,:) :: z3d, zhdiv 166 ! 167 INTEGER :: ji, jj, jk ! dummy loop indices 168 REAL(wp) :: z1_2dt ! local scalars 169 !!---------------------------------------------------------------------- 170 171 IF( nn_timing == 1 ) CALL timing_start('wzv') 162 !!---------------------------------------------------------------------- 163 ! 164 IF( nn_timing == 1 ) CALL timing_start('wzv') 172 165 ! 173 166 IF( kt == nit000 ) THEN 174 !175 167 IF(lwp) WRITE(numout,*) 176 168 IF(lwp) WRITE(numout,*) 'wzv : now vertical velocity ' … … 178 170 ! 179 171 wn(:,:,jpk) = 0._wp ! bottom boundary condition: w=0 (set once for all) 180 !181 172 ENDIF 182 173 ! !------------------------------! … … 194 185 DO jj = 2, jpjm1 195 186 DO ji = fs_2, fs_jpim1 ! vector opt. 196 zhdiv(ji,jj,jk) = r1_e1 2t(ji,jj) * ( un_td(ji,jj,jk) - un_td(ji-1,jj,jk) + vn_td(ji,jj,jk) - vn_td(ji,jj-1,jk) )187 zhdiv(ji,jj,jk) = r1_e1e2t(ji,jj) * ( un_td(ji,jj,jk) - un_td(ji-1,jj,jk) + vn_td(ji,jj,jk) - vn_td(ji,jj-1,jk) ) 197 188 END DO 198 189 END DO … … 217 208 218 209 #if defined key_bdy 219 IF (lk_bdy) THEN210 IF( lk_bdy ) THEN 220 211 DO jk = 1, jpkm1 221 212 wn(:,:,jk) = wn(:,:,jk) * bdytmask(:,:) … … 225 216 ! 226 217 IF( nn_timing == 1 ) CALL timing_stop('wzv') 227 228 218 ! 229 219 END SUBROUTINE wzv 220 230 221 231 222 SUBROUTINE ssh_swp( kt ) … … 259 250 ENDIF 260 251 261 # if defined key_dynspg_ts 262 IF( ( neuler == 0 .AND. kt == nit000 ) .OR. ln_bt_fw ) THEN !** Euler time-stepping: no filter 263 # else 264 IF ( neuler == 0 .AND. kt == nit000 ) THEN !** Euler time-stepping at first time-step : no filter 265 #endif 252 IF( ( neuler == 0 .AND. kt == nit000 ) .OR. ( ln_bt_fw .AND. ln_dynspg_ts ) ) THEN 253 !** Euler time-stepping: no filter 266 254 sshb(:,:) = sshn(:,:) ! before <-- now 267 255 sshn(:,:) = ssha(:,:) ! now <-- after (before already = now) 256 ! 268 257 ELSE !** Leap-Frog time-stepping: Asselin filter + swap 269 258 sshb(:,:) = sshn(:,:) + atfp * ( sshb(:,:) - 2 * sshn(:,:) + ssha(:,:) ) ! before <-- now filtered 270 IF( lk_vvl ) sshb(:,:) = sshb(:,:) - atfp * rdt / rau0 * ( emp_b(:,:) - emp(:,:) - rnf_b(:,:) + rnf(:,:) ) * ssmask(:,:) 259 IF( lk_vvl ) sshb(:,:) = sshb(:,:) - atfp * rdt / rau0 * ( emp_b(:,:) - emp(:,:) & 260 & - rnf_b(:,:) + rnf(:,:) & 261 & + fwfisf_b(:,:) - fwfisf(:,:) ) * ssmask(:,:) 271 262 sshn(:,:) = ssha(:,:) ! now <-- after 272 263 ENDIF 273 !274 ! Update velocity at AGRIF zoom boundaries275 #if defined key_agrif276 IF ( .NOT.Agrif_Root() ) CALL Agrif_Update_Dyn( kt )277 #endif278 264 ! 279 265 IF(ln_ctl) CALL prt_ctl( tab2d_1=sshb, clinfo1=' sshb - : ', mask1=tmask, ovlap=1 )
Note: See TracChangeset
for help on using the changeset viewer.