- Timestamp:
- 2020-11-02T10:56:42+01:00 (4 years ago)
- Location:
- NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves
- Files:
-
- 17 edited
- 1 copied
-
. (modified) (1 prop)
-
src/OCE/DYN/divhor.F90 (modified) (4 diffs)
-
src/OCE/DYN/dynadv_cen2.F90 (modified) (5 diffs)
-
src/OCE/DYN/dynadv_ubs.F90 (modified) (8 diffs)
-
src/OCE/DYN/dynatf.F90 (modified) (11 diffs)
-
src/OCE/DYN/dynatf_qco.F90 (copied) (copied from NEMO/trunk/src/OCE/DYN/dynatf_qco.F90)
-
src/OCE/DYN/dynhpg.F90 (modified) (29 diffs)
-
src/OCE/DYN/dynkeg.F90 (modified) (3 diffs)
-
src/OCE/DYN/dynldf_iso.F90 (modified) (11 diffs)
-
src/OCE/DYN/dynldf_lap_blp.F90 (modified) (3 diffs)
-
src/OCE/DYN/dynspg.F90 (modified) (7 diffs)
-
src/OCE/DYN/dynspg_exp.F90 (modified) (1 diff)
-
src/OCE/DYN/dynspg_ts.F90 (modified) (52 diffs)
-
src/OCE/DYN/dynvor.F90 (modified) (31 diffs)
-
src/OCE/DYN/dynzad.F90 (modified) (7 diffs)
-
src/OCE/DYN/dynzdf.F90 (modified) (20 diffs)
-
src/OCE/DYN/sshwzv.F90 (modified) (15 diffs)
-
src/OCE/DYN/wet_dry.F90 (modified) (12 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves
- Property svn:externals
-
old new 3 3 ^/utils/build/mk@HEAD mk 4 4 ^/utils/tools@HEAD tools 5 ^/vendors/AGRIF/dev @HEADext/AGRIF5 ^/vendors/AGRIF/dev_r12970_AGRIF_CMEMS ext/AGRIF 6 6 ^/vendors/FCM@HEAD ext/FCM 7 7 ^/vendors/IOIPSL@HEAD ext/IOIPSL 8 8 9 9 # SETTE 10 ^/utils/CI/sette@ HEADsette10 ^/utils/CI/sette@13559 sette
-
- Property svn:externals
-
NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/OCE/DYN/divhor.F90
r12377 r13710 40 40 !! * Substitutions 41 41 # include "do_loop_substitute.h90" 42 # include "domzgr_substitute.h90" 42 43 !!---------------------------------------------------------------------- 43 44 !! NEMO/OCE 4.0 , NEMO Consortium (2018) … … 76 77 ENDIF 77 78 ! 78 DO_3D _00_00( 1, jpkm1 )79 hdiv(ji,jj,jk) = ( e2u(ji ,jj) * e3u(ji ,jj,jk,Kmm) * uu(ji ,jj,jk,Kmm) &79 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) !== Horizontal divergence ==! 80 hdiv(ji,jj,jk) = ( e2u(ji ,jj) * e3u(ji ,jj,jk,Kmm) * uu(ji ,jj,jk,Kmm) & 80 81 & - e2u(ji-1,jj) * e3u(ji-1,jj,jk,Kmm) * uu(ji-1,jj,jk,Kmm) & 81 82 & + e1v(ji,jj ) * e3v(ji,jj ,jk,Kmm) * vv(ji,jj ,jk,Kmm) & … … 83 84 & * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm) 84 85 END_3D 85 !86 #if defined key_agrif87 IF( .NOT. Agrif_Root() ) THEN88 IF( nbondi == -1 .OR. nbondi == 2 ) hdiv( 2 , : ,:) = 0._wp ! west89 IF( nbondi == 1 .OR. nbondi == 2 ) hdiv( nlci-1, : ,:) = 0._wp ! east90 IF( nbondj == -1 .OR. nbondj == 2 ) hdiv( : , 2 ,:) = 0._wp ! south91 IF( nbondj == 1 .OR. nbondj == 2 ) hdiv( : ,nlcj-1,:) = 0._wp ! north92 ENDIF93 #endif94 86 ! 95 87 IF( ln_rnf ) CALL sbc_rnf_div( hdiv, Kmm ) !== runoffs ==! (update hdiv field) … … 102 94 IF( ln_isf ) CALL isf_hdiv( kt, Kmm, hdiv ) !== ice shelf ==! (update hdiv field) 103 95 ! 104 CALL lbc_lnk( 'divhor', hdiv, 'T', 1. ) ! (no sign change)96 CALL lbc_lnk( 'divhor', hdiv, 'T', 1.0_wp ) ! (no sign change) 105 97 ! 106 98 IF( ln_timing ) CALL timing_stop('div_hor') -
NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/OCE/DYN/dynadv_cen2.F90
r12377 r13710 28 28 !! * Substitutions 29 29 # include "do_loop_substitute.h90" 30 # include "domzgr_substitute.h90" 30 31 !!---------------------------------------------------------------------- 31 32 !! NEMO/OCE 4.0 , NEMO Consortium (2018) … … 71 72 zfu(:,:,jk) = 0.25_wp * e2u(:,:) * e3u(:,:,jk,Kmm) * puu(:,:,jk,Kmm) 72 73 zfv(:,:,jk) = 0.25_wp * e1v(:,:) * e3v(:,:,jk,Kmm) * pvv(:,:,jk,Kmm) 73 DO_2D _10_1074 DO_2D( 1, 0, 1, 0 ) ! horizontal momentum fluxes (at T- and F-point) 74 75 zfu_t(ji+1,jj ,jk) = ( zfu(ji,jj,jk) + zfu(ji+1,jj,jk) ) * ( puu(ji,jj,jk,Kmm) + puu(ji+1,jj ,jk,Kmm) ) 75 76 zfv_f(ji ,jj ,jk) = ( zfv(ji,jj,jk) + zfv(ji+1,jj,jk) ) * ( puu(ji,jj,jk,Kmm) + puu(ji ,jj+1,jk,Kmm) ) … … 77 78 zfv_t(ji ,jj+1,jk) = ( zfv(ji,jj,jk) + zfv(ji,jj+1,jk) ) * ( pvv(ji,jj,jk,Kmm) + pvv(ji ,jj+1,jk,Kmm) ) 78 79 END_2D 79 DO_2D _00_0080 DO_2D( 0, 0, 0, 0 ) ! divergence of horizontal momentum fluxes 80 81 puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) - ( zfu_t(ji+1,jj,jk) - zfu_t(ji,jj ,jk) & 81 & + zfv_f(ji ,jj,jk) - zfv_f(ji,jj-1,jk) ) * r1_e1e2u(ji,jj) / e3u(ji,jj,jk,Kmm) 82 & + zfv_f(ji ,jj,jk) - zfv_f(ji,jj-1,jk) ) * r1_e1e2u(ji,jj) & 83 & / e3u(ji,jj,jk,Kmm) 82 84 pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) - ( zfu_f(ji,jj ,jk) - zfu_f(ji-1,jj,jk) & 83 & + zfv_t(ji,jj+1,jk) - zfv_t(ji ,jj,jk) ) * r1_e1e2v(ji,jj) / e3v(ji,jj,jk,Kmm) 85 & + zfv_t(ji,jj+1,jk) - zfv_t(ji ,jj,jk) ) * r1_e1e2v(ji,jj) & 86 & / e3v(ji,jj,jk,Kmm) 84 87 END_2D 85 88 END DO … … 95 98 ! !== Vertical advection ==! 96 99 ! 97 DO_2D _00_00100 DO_2D( 0, 0, 0, 0 ) ! surface/bottom advective fluxes set to zero 98 101 zfu_uw(ji,jj,jpk) = 0._wp ; zfv_vw(ji,jj,jpk) = 0._wp 99 102 zfu_uw(ji,jj, 1 ) = 0._wp ; zfv_vw(ji,jj, 1 ) = 0._wp 100 103 END_2D 101 104 IF( ln_linssh ) THEN ! linear free surface: advection through the surface 102 DO_2D _00_00105 DO_2D( 0, 0, 0, 0 ) 103 106 zfu_uw(ji,jj,1) = 0.5_wp * ( e1e2t(ji,jj) * ww(ji,jj,1) + e1e2t(ji+1,jj) * ww(ji+1,jj,1) ) * puu(ji,jj,1,Kmm) 104 107 zfv_vw(ji,jj,1) = 0.5_wp * ( e1e2t(ji,jj) * ww(ji,jj,1) + e1e2t(ji,jj+1) * ww(ji,jj+1,1) ) * pvv(ji,jj,1,Kmm) … … 106 109 ENDIF 107 110 DO jk = 2, jpkm1 ! interior advective fluxes 108 DO_2D _01_01111 DO_2D( 0, 1, 0, 1 ) ! 1/4 * Vertical transport 109 112 zfw(ji,jj,jk) = 0.25_wp * e1e2t(ji,jj) * ww(ji,jj,jk) 110 113 END_2D 111 DO_2D _00_00114 DO_2D( 0, 0, 0, 0 ) 112 115 zfu_uw(ji,jj,jk) = ( zfw(ji,jj,jk) + zfw(ji+1,jj ,jk) ) * ( puu(ji,jj,jk,Kmm) + puu(ji,jj,jk-1,Kmm) ) 113 116 zfv_vw(ji,jj,jk) = ( zfw(ji,jj,jk) + zfw(ji ,jj+1,jk) ) * ( pvv(ji,jj,jk,Kmm) + pvv(ji,jj,jk-1,Kmm) ) 114 117 END_2D 115 118 END DO 116 DO_3D_00_00( 1, jpkm1 ) 117 puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) - ( zfu_uw(ji,jj,jk) - zfu_uw(ji,jj,jk+1) ) * r1_e1e2u(ji,jj) / e3u(ji,jj,jk,Kmm) 118 pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) - ( zfv_vw(ji,jj,jk) - zfv_vw(ji,jj,jk+1) ) * r1_e1e2v(ji,jj) / e3v(ji,jj,jk,Kmm) 119 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) ! divergence of vertical momentum flux divergence 120 puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) - ( zfu_uw(ji,jj,jk) - zfu_uw(ji,jj,jk+1) ) * r1_e1e2u(ji,jj) & 121 & / e3u(ji,jj,jk,Kmm) 122 pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) - ( zfv_vw(ji,jj,jk) - zfv_vw(ji,jj,jk+1) ) * r1_e1e2v(ji,jj) & 123 & / e3v(ji,jj,jk,Kmm) 119 124 END_3D 120 125 ! -
NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/OCE/DYN/dynadv_ubs.F90
r12377 r13710 34 34 !! * Substitutions 35 35 # include "do_loop_substitute.h90" 36 # include "domzgr_substitute.h90" 36 37 !!---------------------------------------------------------------------- 37 38 !! NEMO/OCE 4.0 , NEMO Consortium (2018) … … 107 108 zfv(:,:,jk) = e1v(:,:) * e3v(:,:,jk,Kmm) * pvv(:,:,jk,Kmm) 108 109 ! 109 DO_2D _00_00110 DO_2D( 0, 0, 0, 0 ) ! laplacian 110 111 zlu_uu(ji,jj,jk,1) = ( puu (ji+1,jj ,jk,Kbb) - 2.*puu (ji,jj,jk,Kbb) + puu (ji-1,jj ,jk,Kbb) ) * umask(ji,jj,jk) 111 112 zlv_vv(ji,jj,jk,1) = ( pvv (ji ,jj+1,jk,Kbb) - 2.*pvv (ji,jj,jk,Kbb) + pvv (ji ,jj-1,jk,Kbb) ) * vmask(ji,jj,jk) … … 123 124 END_2D 124 125 END DO 125 CALL lbc_lnk_multi( 'dynadv_ubs', zlu_uu(:,:,:,1), 'U', 1. , zlu_uv(:,:,:,1), 'U', 1., &126 & zlu_uu(:,:,:,2), 'U', 1. , zlu_uv(:,:,:,2), 'U', 1., &127 & zlv_vv(:,:,:,1), 'V', 1. , zlv_vu(:,:,:,1), 'V', 1., &128 & zlv_vv(:,:,:,2), 'V', 1. , zlv_vu(:,:,:,2), 'V', 1.)126 CALL lbc_lnk_multi( 'dynadv_ubs', zlu_uu(:,:,:,1), 'U', 1.0_wp , zlu_uv(:,:,:,1), 'U', 1.0_wp, & 127 & zlu_uu(:,:,:,2), 'U', 1.0_wp , zlu_uv(:,:,:,2), 'U', 1.0_wp, & 128 & zlv_vv(:,:,:,1), 'V', 1.0_wp , zlv_vu(:,:,:,1), 'V', 1.0_wp, & 129 & zlv_vv(:,:,:,2), 'V', 1.0_wp , zlv_vu(:,:,:,2), 'V', 1.0_wp ) 129 130 ! 130 131 ! ! ====================== ! … … 135 136 zfv(:,:,jk) = 0.25_wp * e1v(:,:) * e3v(:,:,jk,Kmm) * pvv(:,:,jk,Kmm) 136 137 ! 137 DO_2D _10_10138 DO_2D( 1, 0, 1, 0 ) ! horizontal momentum fluxes at T- and F-point 138 139 zui = ( puu(ji,jj,jk,Kmm) + puu(ji+1,jj ,jk,Kmm) ) 139 140 zvj = ( pvv(ji,jj,jk,Kmm) + pvv(ji ,jj+1,jk,Kmm) ) … … 167 168 & * ( pvv(ji,jj,jk,Kmm) + pvv(ji+1,jj ,jk,Kmm) - gamma1 * zl_v ) 168 169 END_2D 169 DO_2D _00_00170 DO_2D( 0, 0, 0, 0 ) ! divergence of horizontal momentum fluxes 170 171 puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) - ( zfu_t(ji+1,jj,jk) - zfu_t(ji,jj ,jk) & 171 & + zfv_f(ji ,jj,jk) - zfv_f(ji,jj-1,jk) ) * r1_e1e2u(ji,jj) / e3u(ji,jj,jk,Kmm) 172 & + zfv_f(ji ,jj,jk) - zfv_f(ji,jj-1,jk) ) * r1_e1e2u(ji,jj) & 173 & / e3u(ji,jj,jk,Kmm) 172 174 pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) - ( zfu_f(ji,jj ,jk) - zfu_f(ji-1,jj,jk) & 173 & + zfv_t(ji,jj+1,jk) - zfv_t(ji ,jj,jk) ) * r1_e1e2v(ji,jj) / e3v(ji,jj,jk,Kmm) 175 & + zfv_t(ji,jj+1,jk) - zfv_t(ji ,jj,jk) ) * r1_e1e2v(ji,jj) & 176 & / e3v(ji,jj,jk,Kmm) 174 177 END_2D 175 178 END DO … … 184 187 ! ! Vertical advection ! 185 188 ! ! ==================== ! 186 DO_2D _00_00189 DO_2D( 0, 0, 0, 0 ) ! surface/bottom advective fluxes set to zero 187 190 zfu_uw(ji,jj,jpk) = 0._wp 188 191 zfv_vw(ji,jj,jpk) = 0._wp … … 191 194 END_2D 192 195 IF( ln_linssh ) THEN ! constant volume : advection through the surface 193 DO_2D _00_00196 DO_2D( 0, 0, 0, 0 ) 194 197 zfu_uw(ji,jj,1) = 0.5_wp * ( e1e2t(ji,jj) * ww(ji,jj,1) + e1e2t(ji+1,jj) * ww(ji+1,jj,1) ) * puu(ji,jj,1,Kmm) 195 198 zfv_vw(ji,jj,1) = 0.5_wp * ( e1e2t(ji,jj) * ww(ji,jj,1) + e1e2t(ji,jj+1) * ww(ji,jj+1,1) ) * pvv(ji,jj,1,Kmm) … … 197 200 ENDIF 198 201 DO jk = 2, jpkm1 ! interior fluxes 199 DO_2D _01_01202 DO_2D( 0, 1, 0, 1 ) 200 203 zfw(ji,jj,jk) = 0.25_wp * e1e2t(ji,jj) * ww(ji,jj,jk) 201 204 END_2D 202 DO_2D _00_00205 DO_2D( 0, 0, 0, 0 ) 203 206 zfu_uw(ji,jj,jk) = ( zfw(ji,jj,jk)+ zfw(ji+1,jj,jk) ) * ( puu(ji,jj,jk,Kmm) + puu(ji,jj,jk-1,Kmm) ) 204 207 zfv_vw(ji,jj,jk) = ( zfw(ji,jj,jk)+ zfw(ji,jj+1,jk) ) * ( pvv(ji,jj,jk,Kmm) + pvv(ji,jj,jk-1,Kmm) ) 205 208 END_2D 206 209 END DO 207 DO_3D_00_00( 1, jpkm1 ) 208 puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) - ( zfu_uw(ji,jj,jk) - zfu_uw(ji,jj,jk+1) ) * r1_e1e2u(ji,jj) / e3u(ji,jj,jk,Kmm) 209 pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) - ( zfv_vw(ji,jj,jk) - zfv_vw(ji,jj,jk+1) ) * r1_e1e2v(ji,jj) / e3v(ji,jj,jk,Kmm) 210 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) ! divergence of vertical momentum flux divergence 211 puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) - ( zfu_uw(ji,jj,jk) - zfu_uw(ji,jj,jk+1) ) * r1_e1e2u(ji,jj) & 212 & / e3u(ji,jj,jk,Kmm) 213 pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) - ( zfv_vw(ji,jj,jk) - zfv_vw(ji,jj,jk+1) ) * r1_e1e2v(ji,jj) & 214 & / e3v(ji,jj,jk,Kmm) 210 215 END_3D 211 216 ! -
NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/OCE/DYN/dynatf.F90
r12489 r13710 34 34 USE dynspg_ts ! surface pressure gradient: split-explicit scheme 35 35 USE domvvl ! variable volume 36 USE bdy_oce , ONLY: ln_bdy36 USE bdy_oce , ONLY : ln_bdy 37 37 USE bdydta ! ocean open boundary conditions 38 38 USE bdydyn ! ocean open boundary conditions … … 50 50 USE prtctl ! Print control 51 51 USE timing ! Timing 52 USE zdfdrg , ONLY : ln_drgice_imp, rCdU_top 52 53 #if defined key_agrif 53 54 USE agrif_oce_interp … … 58 59 59 60 PUBLIC dyn_atf ! routine called by step.F90 61 62 #if defined key_qco 63 !!---------------------------------------------------------------------- 64 !! 'key_qco' EMPTY ROUTINE Quasi-Eulerian vertical coordonate 65 !!---------------------------------------------------------------------- 66 CONTAINS 67 68 SUBROUTINE dyn_atf ( kt, Kbb, Kmm, Kaa, puu, pvv, pe3t, pe3u, pe3v ) 69 INTEGER , INTENT(in ) :: kt ! ocean time-step index 70 INTEGER , INTENT(in ) :: Kbb, Kmm, Kaa ! before and after time level indices 71 REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: puu, pvv ! velocities to be time filtered 72 REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: pe3t, pe3u, pe3v ! scale factors to be time filtered 73 74 WRITE(*,*) 'dyn_atf: You should not have seen this print! error?', kt 75 END SUBROUTINE dyn_atf 76 77 #else 60 78 61 79 !! * Substitutions … … 103 121 REAL(wp) :: zve3a, zve3n, zve3b, z1_2dt ! - - 104 122 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zue, zve, zwfld 123 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zutau, zvtau 105 124 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: ze3t_f, ze3u_f, ze3v_f, zua, zva 106 125 !!---------------------------------------------------------------------- … … 148 167 # endif 149 168 ! 150 CALL lbc_lnk_multi( 'dynatf', puu(:,:,:,Kaa), 'U', -1. , pvv(:,:,:,Kaa), 'V', -1.) !* local domain boundaries169 CALL lbc_lnk_multi( 'dynatf', puu(:,:,:,Kaa), 'U', -1.0_wp, pvv(:,:,:,Kaa), 'V', -1.0_wp ) !* local domain boundaries 151 170 ! 152 171 ! !* BDY open boundaries … … 180 199 IF( ln_linssh ) THEN ! Fixed volume ! 181 200 ! ! =============! 182 DO_3D _11_11(1, jpkm1 )201 DO_3D( 1, 1, 1, 1, 1, jpkm1 ) 183 202 puu(ji,jj,jk,Kmm) = puu(ji,jj,jk,Kmm) + rn_atfp * ( puu(ji,jj,jk,Kbb) - 2._wp * puu(ji,jj,jk,Kmm) + puu(ji,jj,jk,Kaa) ) 184 203 pvv(ji,jj,jk,Kmm) = pvv(ji,jj,jk,Kmm) + rn_atfp * ( pvv(ji,jj,jk,Kbb) - 2._wp * pvv(ji,jj,jk,Kmm) + pvv(ji,jj,jk,Kaa) ) … … 198 217 zwfld(:,:) = emp_b(:,:) - emp(:,:) 199 218 IF ( ln_rnf ) zwfld(:,:) = zwfld(:,:) - ( rnf_b(:,:) - rnf(:,:) ) 219 200 220 DO jk = 1, jpkm1 201 221 ze3t_f(:,:,jk) = ze3t_f(:,:,jk) - zcoef * zwfld(:,:) * tmask(:,:,jk) & … … 215 235 CALL dom_vvl_interpol( pe3t(:,:,:,Kmm), pe3u(:,:,:,Kmm), 'U' ) 216 236 CALL dom_vvl_interpol( pe3t(:,:,:,Kmm), pe3v(:,:,:,Kmm), 'V' ) 217 DO_3D _11_11(1, jpkm1 )237 DO_3D( 1, 1, 1, 1, 1, jpkm1 ) 218 238 puu(ji,jj,jk,Kmm) = puu(ji,jj,jk,Kmm) + rn_atfp * ( puu(ji,jj,jk,Kbb) - 2._wp * puu(ji,jj,jk,Kmm) + puu(ji,jj,jk,Kaa) ) 219 239 pvv(ji,jj,jk,Kmm) = pvv(ji,jj,jk,Kmm) + rn_atfp * ( pvv(ji,jj,jk,Kbb) - 2._wp * pvv(ji,jj,jk,Kmm) + pvv(ji,jj,jk,Kaa) ) … … 226 246 CALL dom_vvl_interpol( pe3t(:,:,:,Kmm), ze3u_f, 'U' ) 227 247 CALL dom_vvl_interpol( pe3t(:,:,:,Kmm), ze3v_f, 'V' ) 228 DO_3D _11_11(1, jpkm1 )248 DO_3D( 1, 1, 1, 1, 1, jpkm1 ) 229 249 zue3a = pe3u(ji,jj,jk,Kaa) * puu(ji,jj,jk,Kaa) 230 250 zve3a = pe3v(ji,jj,jk,Kaa) * pvv(ji,jj,jk,Kaa) … … 303 323 ENDIF 304 324 ! 325 IF ( iom_use("utau") ) THEN 326 IF ( ln_drgice_imp.OR.ln_isfcav ) THEN 327 ALLOCATE(zutau(jpi,jpj)) 328 DO_2D( 0, 0, 0, 0 ) 329 jk = miku(ji,jj) 330 zutau(ji,jj) = utau(ji,jj) + 0.5_wp * rho0 * ( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * puu(ji,jj,jk,Kaa) 331 END_2D 332 CALL iom_put( "utau", zutau(:,:) ) 333 DEALLOCATE(zutau) 334 ELSE 335 CALL iom_put( "utau", utau(:,:) ) 336 ENDIF 337 ENDIF 338 ! 339 IF ( iom_use("vtau") ) THEN 340 IF ( ln_drgice_imp.OR.ln_isfcav ) THEN 341 ALLOCATE(zvtau(jpi,jpj)) 342 DO_2D( 0, 0, 0, 0 ) 343 jk = mikv(ji,jj) 344 zvtau(ji,jj) = vtau(ji,jj) + 0.5_wp * rho0 * ( rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) * pvv(ji,jj,jk,Kaa) 345 END_2D 346 CALL iom_put( "vtau", zvtau(:,:) ) 347 DEALLOCATE(zvtau) 348 ELSE 349 CALL iom_put( "vtau", vtau(:,:) ) 350 ENDIF 351 ENDIF 352 ! 305 353 IF(sn_cfctl%l_prtctl) CALL prt_ctl( tab3d_1=puu(:,:,:,Kaa), clinfo1=' nxt - puu(:,:,:,Kaa): ', mask1=umask, & 306 354 & tab3d_2=pvv(:,:,:,Kaa), clinfo2=' pvv(:,:,:,Kaa): ' , mask2=vmask ) … … 312 360 END SUBROUTINE dyn_atf 313 361 362 #endif 363 314 364 !!========================================================================= 315 365 END MODULE dynatf -
NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/OCE/DYN/dynhpg.F90
r12377 r13710 76 76 !! * Substitutions 77 77 # include "do_loop_substitute.h90" 78 # include "domzgr_substitute.h90" 79 78 80 !!---------------------------------------------------------------------- 79 81 !! NEMO/OCE 4.0 , NEMO Consortium (2018) … … 255 257 256 258 ! Surface value 257 DO_2D _00_00259 DO_2D( 0, 0, 0, 0 ) 258 260 zcoef1 = zcoef0 * e3w(ji,jj,1,Kmm) 259 261 ! hydrostatic pressure gradient … … 267 269 ! 268 270 ! interior value (2=<jk=<jpkm1) 269 DO_3D _00_00(2, jpkm1 )271 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 270 272 zcoef1 = zcoef0 * e3w(ji,jj,jk,Kmm) 271 273 ! hydrostatic pressure gradient … … 317 319 318 320 ! Surface value (also valid in partial step case) 319 DO_2D _00_00321 DO_2D( 0, 0, 0, 0 ) 320 322 zcoef1 = zcoef0 * e3w(ji,jj,1,Kmm) 321 323 ! hydrostatic pressure gradient … … 328 330 329 331 ! interior value (2=<jk=<jpkm1) 330 DO_3D _00_00(2, jpkm1 )332 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 331 333 zcoef1 = zcoef0 * e3w(ji,jj,jk,Kmm) 332 334 ! hydrostatic pressure gradient … … 344 346 345 347 ! partial steps correction at the last level (use zgru & zgrv computed in zpshde.F90) 346 DO_2D _00_00348 DO_2D( 0, 0, 0, 0 ) 347 349 iku = mbku(ji,jj) 348 350 ikv = mbkv(ji,jj) … … 409 411 ! 410 412 IF( ln_wd_il ) THEN 411 DO_2D _00_00413 DO_2D( 0, 0, 0, 0 ) 412 414 ll_tmp1 = MIN( ssh(ji,jj,Kmm) , ssh(ji+1,jj,Kmm) ) > & 413 415 & MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) .AND. & … … 446 448 END IF 447 449 END_2D 448 CALL lbc_lnk_multi( 'dynhpg', zcpx, 'U', 1. , zcpy, 'V', 1.)450 CALL lbc_lnk_multi( 'dynhpg', zcpx, 'U', 1.0_wp, zcpy, 'V', 1.0_wp ) 449 451 END IF 450 452 451 453 ! Surface value 452 DO_2D _00_00454 DO_2D( 0, 0, 0, 0 ) 453 455 ! hydrostatic pressure gradient along s-surfaces 454 zhpi(ji,jj,1) = zcoef0 * ( e3w(ji+1,jj ,1,Kmm) * ( znad + rhd(ji+1,jj ,1) ) & 455 & - e3w(ji ,jj ,1,Kmm) * ( znad + rhd(ji ,jj ,1) ) ) * r1_e1u(ji,jj) 456 zhpj(ji,jj,1) = zcoef0 * ( e3w(ji ,jj+1,1,Kmm) * ( znad + rhd(ji ,jj+1,1) ) & 457 & - e3w(ji ,jj ,1,Kmm) * ( znad + rhd(ji ,jj ,1) ) ) * r1_e2v(ji,jj) 456 zhpi(ji,jj,1) = & 457 & zcoef0 * ( e3w(ji+1,jj ,1,Kmm) * ( znad + rhd(ji+1,jj ,1) ) & 458 & - e3w(ji ,jj ,1,Kmm) * ( znad + rhd(ji ,jj ,1) ) ) & 459 & * r1_e1u(ji,jj) 460 zhpj(ji,jj,1) = & 461 & zcoef0 * ( e3w(ji ,jj+1,1,Kmm) * ( znad + rhd(ji ,jj+1,1) ) & 462 & - e3w(ji ,jj ,1,Kmm) * ( znad + rhd(ji ,jj ,1) ) ) & 463 & * r1_e2v(ji,jj) 458 464 ! s-coordinate pressure gradient correction 459 465 zuap = -zcoef0 * ( rhd (ji+1,jj,1) + rhd (ji,jj,1) + 2._wp * znad ) & … … 475 481 476 482 ! interior value (2=<jk=<jpkm1) 477 DO_3D _00_00(2, jpkm1 )483 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 478 484 ! hydrostatic pressure gradient along s-surfaces 479 485 zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + zcoef0 * r1_e1u(ji,jj) & … … 557 563 !===== Compute surface value ===================================================== 558 564 !================================================================================== 559 DO_2D _00_00565 DO_2D( 0, 0, 0, 0 ) 560 566 ikt = mikt(ji,jj) 561 567 iktp1i = mikt(ji+1,jj) … … 586 592 !================================================================================== 587 593 ! interior value (2=<jk=<jpkm1) 588 DO_3D _00_00(2, jpkm1 )594 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 589 595 ! hydrostatic pressure gradient along s-surfaces 590 596 zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + zcoef0 / e1u(ji,jj) & 591 & * ( e3w(ji+1,jj,jk,Kmm) * ( rhd(ji+1,jj,jk) + rhd(ji+1,jj,jk-1) + 2*znad ) * wmask(ji+1,jj,jk) & 592 & - e3w(ji ,jj,jk,Kmm) * ( rhd(ji ,jj,jk) + rhd(ji ,jj,jk-1) + 2*znad ) * wmask(ji ,jj,jk) ) 597 & * ( e3w(ji+1,jj,jk,Kmm) & 598 & * ( rhd(ji+1,jj,jk) + rhd(ji+1,jj,jk-1) + 2*znad ) * wmask(ji+1,jj,jk) & 599 & - e3w(ji ,jj,jk,Kmm) & 600 & * ( rhd(ji ,jj,jk) + rhd(ji ,jj,jk-1) + 2*znad ) * wmask(ji ,jj,jk) ) 593 601 zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) + zcoef0 / e2v(ji,jj) & 594 & * ( e3w(ji,jj+1,jk,Kmm) * ( rhd(ji,jj+1,jk) + rhd(ji,jj+1,jk-1) + 2*znad ) * wmask(ji,jj+1,jk) & 595 & - e3w(ji,jj ,jk,Kmm) * ( rhd(ji,jj, jk) + rhd(ji,jj ,jk-1) + 2*znad ) * wmask(ji,jj ,jk) ) 602 & * ( e3w(ji,jj+1,jk,Kmm) & 603 & * ( rhd(ji,jj+1,jk) + rhd(ji,jj+1,jk-1) + 2*znad ) * wmask(ji,jj+1,jk) & 604 & - e3w(ji,jj ,jk,Kmm) & 605 & * ( rhd(ji,jj, jk) + rhd(ji,jj ,jk-1) + 2*znad ) * wmask(ji,jj ,jk) ) 596 606 ! s-coordinate pressure gradient correction 597 607 zuap = -zcoef0 * ( rhd (ji+1,jj ,jk) + rhd (ji,jj,jk) + 2._wp * znad ) & … … 633 643 IF( ln_wd_il ) THEN 634 644 ALLOCATE( zcpx(jpi,jpj) , zcpy(jpi,jpj) ) 635 DO_2D _00_00645 DO_2D( 0, 0, 0, 0 ) 636 646 ll_tmp1 = MIN( ssh(ji,jj,Kmm) , ssh(ji+1,jj,Kmm) ) > & 637 647 & MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) .AND. & … … 669 679 END IF 670 680 END_2D 671 CALL lbc_lnk_multi( 'dynhpg', zcpx, 'U', 1. , zcpy, 'V', 1.)681 CALL lbc_lnk_multi( 'dynhpg', zcpx, 'U', 1.0_wp, zcpy, 'V', 1.0_wp ) 672 682 END IF 673 683 … … 689 699 !!bug gm Not a true bug, but... dzz=e3w for dzx, dzy verify what it is really 690 700 691 DO_3D _00_00(2, jpkm1 )701 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 692 702 drhoz(ji,jj,jk) = rhd (ji ,jj ,jk) - rhd (ji,jj,jk-1) 693 703 dzz (ji,jj,jk) = gde3w(ji ,jj ,jk) - gde3w(ji,jj,jk-1) … … 706 716 !!bug gm idem for drhox, drhoy et ji=jpi and jj=jpj 707 717 708 DO_3D _00_00(2, jpkm1 )718 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 709 719 cffw = 2._wp * drhoz(ji ,jj ,jk) * drhoz(ji,jj,jk-1) 710 720 … … 771 781 !------------------------------------------------------------- 772 782 773 !!bug gm : e3w-gde3w = 0.5*e3w .... and gde3w(2)-gde3w(1)=e3w(2) .... to be verified774 ! true if gde3w is really defined as the sum of the e3w scale factors as, it seems to me, it should be775 776 DO_2D _00_00783 !!bug gm : e3w-gde3w(:,:,:) = 0.5*e3w .... and gde3w(:,:,2)-gde3w(:,:,1)=e3w(:,:,2,Kmm) .... to be verified 784 ! true if gde3w(:,:,:) is really defined as the sum of the e3w scale factors as, it seems to me, it should be 785 786 DO_2D( 0, 0, 0, 0 ) 777 787 rho_k(ji,jj,1) = -grav * ( e3w(ji,jj,1,Kmm) - gde3w(ji,jj,1) ) & 778 788 & * ( rhd(ji,jj,1) & … … 785 795 !!bug gm : optimisation: 1/10 and 1/12 the division should be done before the loop 786 796 787 DO_3D _00_00(2, jpkm1 )797 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 788 798 789 799 rho_k(ji,jj,jk) = zcoef0 * ( rhd (ji,jj,jk) + rhd (ji,jj,jk-1) ) & … … 815 825 816 826 END_3D 817 CALL lbc_lnk_multi( 'dynhpg', rho_k, 'W', 1. , rho_i, 'U', 1., rho_j, 'V', 1.)827 CALL lbc_lnk_multi( 'dynhpg', rho_k, 'W', 1.0_wp, rho_i, 'U', 1.0_wp, rho_j, 'V', 1.0_wp ) 818 828 819 829 ! --------------- 820 830 ! Surface value 821 831 ! --------------- 822 DO_2D _00_00832 DO_2D( 0, 0, 0, 0 ) 823 833 zhpi(ji,jj,1) = ( rho_k(ji+1,jj ,1) - rho_k(ji,jj,1) - rho_i(ji,jj,1) ) * r1_e1u(ji,jj) 824 834 zhpj(ji,jj,1) = ( rho_k(ji ,jj+1,1) - rho_k(ji,jj,1) - rho_j(ji,jj,1) ) * r1_e2v(ji,jj) … … 835 845 ! interior value (2=<jk=<jpkm1) 836 846 ! ---------------- 837 DO_3D _00_00(2, jpkm1 )847 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 838 848 ! hydrostatic pressure gradient along s-surfaces 839 849 zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) & … … 901 911 IF( ln_wd_il ) THEN 902 912 ALLOCATE( zcpx(jpi,jpj) , zcpy(jpi,jpj) ) 903 DO_2D _00_00913 DO_2D( 0, 0, 0, 0 ) 904 914 ll_tmp1 = MIN( ssh(ji,jj,Kmm) , ssh(ji+1,jj,Kmm) ) > & 905 915 & MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) .AND. & … … 942 952 ENDIF 943 953 END_2D 944 CALL lbc_lnk_multi( 'dynhpg', zcpx, 'U', 1. , zcpy, 'V', 1.)954 CALL lbc_lnk_multi( 'dynhpg', zcpx, 'U', 1.0_wp, zcpy, 'V', 1.0_wp ) 945 955 ENDIF 946 956 … … 950 960 951 961 ! Preparing vertical density profile "zrhh(:,:,:)" for hybrid-sco coordinate 952 DO_2D _11_11953 jk = mbkt(ji,jj) +1954 IF( jk <= 0) THEN ; zrhh(ji,jj, : ) = 0._wp955 ELSEIF( jk == 1) THEN ; zrhh(ji,jj,jk+1:jpk) = rhd(ji,jj,jk)962 DO_2D( 1, 1, 1, 1 ) 963 jk = mbkt(ji,jj) 964 IF( jk <= 1 ) THEN ; zrhh(ji,jj, : ) = 0._wp 965 ELSEIF( jk == 2 ) THEN ; zrhh(ji,jj,jk+1:jpk) = rhd(ji,jj,jk) 956 966 ELSEIF( jk < jpkm1 ) THEN 957 967 DO jkk = jk+1, jpk 958 968 zrhh(ji,jj,jkk) = interp1(gde3w(ji,jj,jkk ), gde3w(ji,jj,jkk-1), & 959 & gde3w(ji,jj,jkk-2), rhd (ji,jj,jkk-1), rhd(ji,jj,jkk-2))969 & gde3w(ji,jj,jkk-2), zrhh (ji,jj,jkk-1), zrhh(ji,jj,jkk-2)) 960 970 END DO 961 971 ENDIF … … 963 973 964 974 ! Transfer the depth of "T(:,:,:)" to vertical coordinate "zdept(:,:,:)" 965 DO_2D _11_11975 DO_2D( 1, 1, 1, 1 ) 966 976 zdept(ji,jj,1) = 0.5_wp * e3w(ji,jj,1,Kmm) - ssh(ji,jj,Kmm) * znad 967 977 END_2D 968 978 969 DO_3D _11_11(2, jpk )979 DO_3D( 1, 1, 1, 1, 2, jpk ) 970 980 zdept(ji,jj,jk) = zdept(ji,jj,jk-1) + e3w(ji,jj,jk,Kmm) 971 981 END_3D … … 980 990 981 991 ! Integrate the hydrostatic pressure "zhpi(:,:,:)" at "T(ji,jj,1)" 982 DO_2D _01_01992 DO_2D( 0, 1, 0, 1 ) 983 993 zrhdt1 = zrhh(ji,jj,1) - interp3( zdept(ji,jj,1), asp(ji,jj,1), bsp(ji,jj,1), & 984 994 & csp(ji,jj,1), dsp(ji,jj,1) ) * 0.25_wp * e3w(ji,jj,1,Kmm) … … 989 999 990 1000 ! Calculate the pressure "zhpi(:,:,:)" at "T(ji,jj,2:jpkm1)" 991 DO_3D _01_01(2, jpkm1 )1001 DO_3D( 0, 1, 0, 1, 2, jpkm1 ) 992 1002 zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + & 993 1003 & integ_spline( zdept(ji,jj,jk-1), zdept(ji,jj,jk), & … … 999 1009 1000 1010 ! Prepare zsshu_n and zsshv_n 1001 DO_2D _00_001011 DO_2D( 0, 0, 0, 0 ) 1002 1012 !!gm BUG ? if it is ssh at u- & v-point then it should be: 1003 1013 ! zsshu_n(ji,jj) = (e1e2t(ji,jj) * ssh(ji,jj,Kmm) + e1e2t(ji+1,jj) * ssh(ji+1,jj,Kmm)) * & … … 1012 1022 END_2D 1013 1023 1014 CALL lbc_lnk_multi ('dynhpg', zsshu_n, 'U', 1. , zsshv_n, 'V', 1.)1015 1016 DO_2D _00_001024 CALL lbc_lnk_multi ('dynhpg', zsshu_n, 'U', 1.0_wp, zsshv_n, 'V', 1.0_wp ) 1025 1026 DO_2D( 0, 0, 0, 0 ) 1017 1027 zu(ji,jj,1) = - ( e3u(ji,jj,1,Kmm) - zsshu_n(ji,jj) * znad) 1018 1028 zv(ji,jj,1) = - ( e3v(ji,jj,1,Kmm) - zsshv_n(ji,jj) * znad) 1019 1029 END_2D 1020 1030 1021 DO_3D _00_00(2, jpkm1 )1031 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 1022 1032 zu(ji,jj,jk) = zu(ji,jj,jk-1) - e3u(ji,jj,jk,Kmm) 1023 1033 zv(ji,jj,jk) = zv(ji,jj,jk-1) - e3v(ji,jj,jk,Kmm) 1024 1034 END_3D 1025 1035 1026 DO_3D _00_00(1, jpkm1 )1036 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 1027 1037 zu(ji,jj,jk) = zu(ji,jj,jk) + 0.5_wp * e3u(ji,jj,jk,Kmm) 1028 1038 zv(ji,jj,jk) = zv(ji,jj,jk) + 0.5_wp * e3v(ji,jj,jk,Kmm) 1029 1039 END_3D 1030 1040 1031 DO_3D _00_00(1, jpkm1 )1041 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 1032 1042 zu(ji,jj,jk) = MIN( zu(ji,jj,jk) , MAX( -zdept(ji,jj,jk) , -zdept(ji+1,jj,jk) ) ) 1033 1043 zu(ji,jj,jk) = MAX( zu(ji,jj,jk) , MIN( -zdept(ji,jj,jk) , -zdept(ji+1,jj,jk) ) ) … … 1037 1047 1038 1048 1039 DO_3D _00_00(1, jpkm1 )1049 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 1040 1050 zpwes = 0._wp; zpwed = 0._wp 1041 1051 zpnss = 0._wp; zpnsd = 0._wp … … 1359 1369 !!====================================================================== 1360 1370 END MODULE dynhpg 1361 -
NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/OCE/DYN/dynkeg.F90
r12377 r13710 101 101 ! 102 102 CASE ( nkeg_C2 ) !-- Standard scheme --! 103 DO_3D _01_01(1, jpkm1 )103 DO_3D( 0, 1, 0, 1, 1, jpkm1 ) 104 104 zu = puu(ji-1,jj ,jk,Kmm) * puu(ji-1,jj ,jk,Kmm) & 105 105 & + puu(ji ,jj ,jk,Kmm) * puu(ji ,jj ,jk,Kmm) … … 109 109 END_3D 110 110 CASE ( nkeg_HW ) !-- Hollingsworth scheme --! 111 DO_3D _00_00(1, jpkm1 )111 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 112 112 zu = 8._wp * ( puu(ji-1,jj ,jk,Kmm) * puu(ji-1,jj ,jk,Kmm) & 113 113 & + puu(ji ,jj ,jk,Kmm) * puu(ji ,jj ,jk,Kmm) ) & … … 121 121 zhke(ji,jj,jk) = r1_48 * ( zv + zu ) 122 122 END_3D 123 CALL lbc_lnk( 'dynkeg', zhke, 'T', 1. )123 CALL lbc_lnk( 'dynkeg', zhke, 'T', 1.0_wp ) 124 124 ! 125 125 END SELECT 126 126 ! 127 DO_3D _00_00( 1, jpkm1 )127 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) !== grad( KE ) added to the general momentum trends ==! 128 128 puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) - ( zhke(ji+1,jj ,jk) - zhke(ji,jj,jk) ) / e1u(ji,jj) 129 129 pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) - ( zhke(ji ,jj+1,jk) - zhke(ji,jj,jk) ) / e2v(ji,jj) -
NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/OCE/DYN/dynldf_iso.F90
r12377 r13710 42 42 !! * Substitutions 43 43 # include "do_loop_substitute.h90" 44 # include "domzgr_substitute.h90" 44 45 !!---------------------------------------------------------------------- 45 46 !! NEMO/OCE 4.0 , NEMO Consortium (2018) … … 127 128 IF( ln_dynldf_hor .AND. ln_traldf_iso ) THEN 128 129 ! 129 DO_3D _00_00( 1, jpk )130 DO_3D( 0, 0, 0, 0, 1, jpk ) ! set the slopes of iso-level 130 131 uslp (ji,jj,jk) = - ( gdept(ji+1,jj,jk,Kbb) - gdept(ji ,jj ,jk,Kbb) ) * r1_e1u(ji,jj) * umask(ji,jj,jk) 131 132 vslp (ji,jj,jk) = - ( gdept(ji,jj+1,jk,Kbb) - gdept(ji ,jj ,jk,Kbb) ) * r1_e2v(ji,jj) * vmask(ji,jj,jk) … … 134 135 END_3D 135 136 ! Lateral boundary conditions on the slopes 136 CALL lbc_lnk_multi( 'dynldf_iso', uslp , 'U', -1. , vslp , 'V', -1., wslpi, 'W', -1., wslpj, 'W', -1.)137 CALL lbc_lnk_multi( 'dynldf_iso', uslp , 'U', -1.0_wp, vslp , 'V', -1.0_wp, wslpi, 'W', -1.0_wp, wslpj, 'W', -1.0_wp ) 137 138 ! 138 139 ENDIF … … 167 168 168 169 IF( ln_zps ) THEN ! z-coordinate - partial steps : min(e3u) 169 DO_2D_00_01 170 zabe1 = ( ahmt(ji,jj,jk)+rn_ahm_b ) * e2t(ji,jj) * MIN( e3u(ji,jj,jk,Kmm), e3u(ji-1,jj,jk,Kmm) ) * r1_e1t(ji,jj) 170 DO_2D( 0, 0, 0, 1 ) 171 zabe1 = ( ahmt(ji,jj,jk)+rn_ahm_b ) * e2t(ji,jj) & 172 & * MIN( e3u(ji ,jj,jk,Kmm), & 173 & e3u(ji-1,jj,jk,Kmm) ) * r1_e1t(ji,jj) 171 174 172 175 zmskt = 1._wp / MAX( umask(ji-1,jj,jk )+umask(ji,jj,jk+1) & … … 180 183 END_2D 181 184 ELSE ! other coordinate system (zco or sco) : e3t 182 DO_2D_00_01 183 zabe1 = ( ahmt(ji,jj,jk)+rn_ahm_b ) * e2t(ji,jj) * e3t(ji,jj,jk,Kmm) * r1_e1t(ji,jj) 185 DO_2D( 0, 0, 0, 1 ) 186 zabe1 = ( ahmt(ji,jj,jk)+rn_ahm_b ) & 187 & * e2t(ji,jj) * e3t(ji,jj,jk,Kmm) * r1_e1t(ji,jj) 184 188 185 189 zmskt = 1._wp / MAX( umask(ji-1,jj,jk ) + umask(ji,jj,jk+1) & … … 195 199 196 200 ! j-flux at f-point 197 DO_2D_10_10 198 zabe2 = ( ahmf(ji,jj,jk) + rn_ahm_b ) * e1f(ji,jj) * e3f(ji,jj,jk) * r1_e2f(ji,jj) 201 DO_2D( 1, 0, 1, 0 ) 202 zabe2 = ( ahmf(ji,jj,jk) + rn_ahm_b ) & 203 & * e1f(ji,jj) * e3f(ji,jj,jk) * r1_e2f(ji,jj) 199 204 200 205 zmskf = 1._wp / MAX( umask(ji,jj+1,jk )+umask(ji,jj,jk+1) & … … 214 219 ! i-flux at f-point | t | 215 220 216 DO_2D_00_10 217 zabe1 = ( ahmf(ji,jj,jk) + rn_ahm_b ) * e2f(ji,jj) * e3f(ji,jj,jk) * r1_e1f(ji,jj) 221 DO_2D( 0, 0, 1, 0 ) 222 zabe1 = ( ahmf(ji,jj,jk) + rn_ahm_b ) & 223 & * e2f(ji,jj) * e3f(ji,jj,jk) * r1_e1f(ji,jj) 218 224 219 225 zmskf = 1._wp / MAX( vmask(ji+1,jj,jk )+vmask(ji,jj,jk+1) & … … 229 235 ! j-flux at t-point 230 236 IF( ln_zps ) THEN ! z-coordinate - partial steps : min(e3u) 231 DO_2D_01_10 232 zabe2 = ( ahmt(ji,jj,jk)+rn_ahm_b ) * e1t(ji,jj) * MIN( e3v(ji,jj,jk,Kmm), e3v(ji,jj-1,jk,Kmm) ) * r1_e2t(ji,jj) 237 DO_2D( 0, 1, 1, 0 ) 238 zabe2 = ( ahmt(ji,jj,jk)+rn_ahm_b ) * e1t(ji,jj) & 239 & * MIN( e3v(ji,jj ,jk,Kmm), & 240 & e3v(ji,jj-1,jk,Kmm) ) * r1_e2t(ji,jj) 233 241 234 242 zmskt = 1._wp / MAX( vmask(ji,jj-1,jk )+vmask(ji,jj,jk+1) & … … 242 250 END_2D 243 251 ELSE ! other coordinate system (zco or sco) : e3t 244 DO_2D_01_10 245 zabe2 = ( ahmt(ji,jj,jk)+rn_ahm_b ) * e1t(ji,jj) * e3t(ji,jj,jk,Kmm) * r1_e2t(ji,jj) 252 DO_2D( 0, 1, 1, 0 ) 253 zabe2 = ( ahmt(ji,jj,jk)+rn_ahm_b ) & 254 & * e1t(ji,jj) * e3t(ji,jj,jk,Kmm) * r1_e2t(ji,jj) 246 255 247 256 zmskt = 1./MAX( vmask(ji,jj-1,jk )+vmask(ji,jj,jk+1) & … … 259 268 ! Second derivative (divergence) and add to the general trend 260 269 ! ----------------------------------------------------------- 261 DO_2D _00_00270 DO_2D( 0, 0, 0, 0 ) !!gm Question vectop possible??? !!bug 262 271 puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + ( ziut(ji+1,jj) - ziut(ji,jj ) & 263 & + zjuf(ji ,jj) - zjuf(ji,jj-1) ) * r1_e1e2u(ji,jj) / e3u(ji,jj,jk,Kmm) 272 & + zjuf(ji ,jj) - zjuf(ji,jj-1) ) * r1_e1e2u(ji,jj) & 273 & / e3u(ji,jj,jk,Kmm) 264 274 pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + ( zivf(ji,jj ) - zivf(ji-1,jj) & 265 & + zjvt(ji,jj+1) - zjvt(ji,jj ) ) * r1_e1e2v(ji,jj) / e3v(ji,jj,jk,Kmm) 275 & + zjvt(ji,jj+1) - zjvt(ji,jj ) ) * r1_e1e2v(ji,jj) & 276 & / e3v(ji,jj,jk,Kmm) 266 277 END_2D 267 278 ! ! =============== … … 375 386 DO jk = 1, jpkm1 376 387 DO ji = 2, jpim1 377 puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + ( zfuw(ji,jk) - zfuw(ji,jk+1) ) * r1_e1e2u(ji,jj) / e3u(ji,jj,jk,Kmm) 378 pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + ( zfvw(ji,jk) - zfvw(ji,jk+1) ) * r1_e1e2v(ji,jj) / e3v(ji,jj,jk,Kmm) 388 puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + ( zfuw(ji,jk) - zfuw(ji,jk+1) ) * r1_e1e2u(ji,jj) & 389 & / e3u(ji,jj,jk,Kmm) 390 pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + ( zfvw(ji,jk) - zfvw(ji,jk+1) ) * r1_e1e2v(ji,jj) & 391 & / e3v(ji,jj,jk,Kmm) 379 392 END DO 380 393 END DO -
NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/OCE/DYN/dynldf_lap_blp.F90
r12377 r13710 28 28 !! * Substitutions 29 29 # include "do_loop_substitute.h90" 30 # include "domzgr_substitute.h90" 30 31 !!---------------------------------------------------------------------- 31 32 !! NEMO/OCE 4.0 , NEMO Consortium (2018) … … 72 73 DO jk = 1, jpkm1 ! Horizontal slab 73 74 ! ! =============== 74 DO_2D _01_0175 DO_2D( 0, 1, 0, 1 ) 75 76 ! ! ahm * e3 * curl (computed from 1 to jpim1/jpjm1) 76 !!gm open question here : e3f at before or now ? probably now... 77 !!gm note that ahmf has already been multiplied by fmask 78 zcur(ji-1,jj-1) = ahmf(ji-1,jj-1,jk) * e3f(ji-1,jj-1,jk) * r1_e1e2f(ji-1,jj-1) & 77 zcur(ji-1,jj-1) = ahmf(ji-1,jj-1,jk) * e3f(ji-1,jj-1,jk) * r1_e1e2f(ji-1,jj-1) & ! ahmf already * by fmask 79 78 & * ( e2v(ji ,jj-1) * pv(ji ,jj-1,jk) - e2v(ji-1,jj-1) * pv(ji-1,jj-1,jk) & 80 79 & - e1u(ji-1,jj ) * pu(ji-1,jj ,jk) + e1u(ji-1,jj-1) * pu(ji-1,jj-1,jk) ) 81 80 ! ! ahm * div (computed from 2 to jpi/jpj) 82 !!gm note that ahmt has already been multiplied by tmask 83 zdiv(ji,jj) = ahmt(ji,jj,jk) * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kbb) & 81 zdiv(ji,jj) = ahmt(ji,jj,jk) * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kbb) & ! ahmt already * by tmask 84 82 & * ( e2u(ji,jj)*e3u(ji,jj,jk,Kbb) * pu(ji,jj,jk) - e2u(ji-1,jj)*e3u(ji-1,jj,jk,Kbb) * pu(ji-1,jj,jk) & 85 83 & + e1v(ji,jj)*e3v(ji,jj,jk,Kbb) * pv(ji,jj,jk) - e1v(ji,jj-1)*e3v(ji,jj-1,jk,Kbb) * pv(ji,jj-1,jk) ) 86 84 END_2D 87 85 ! 88 DO_2D _00_0089 pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + zsign * ( &86 DO_2D( 0, 0, 0, 0 ) ! - curl( curl) + grad( div ) 87 pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + zsign * umask(ji,jj,jk) * ( & ! * by umask is mandatory for dyn_ldf_blp use 90 88 & - ( zcur(ji ,jj) - zcur(ji,jj-1) ) * r1_e2u(ji,jj) / e3u(ji,jj,jk,Kmm) & 91 & + ( zdiv(ji+1,jj) - zdiv(ji,jj ) ) * r1_e1u(ji,jj) )89 & + ( zdiv(ji+1,jj) - zdiv(ji,jj ) ) * r1_e1u(ji,jj) ) 92 90 ! 93 pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) + zsign * ( &91 pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) + zsign * vmask(ji,jj,jk) * ( & ! * by vmask is mandatory for dyn_ldf_blp use 94 92 & ( zcur(ji,jj ) - zcur(ji-1,jj) ) * r1_e1v(ji,jj) / e3v(ji,jj,jk,Kmm) & 95 & + ( zdiv(ji,jj+1) - zdiv(ji ,jj) ) * r1_e2v(ji,jj) )93 & + ( zdiv(ji,jj+1) - zdiv(ji ,jj) ) * r1_e2v(ji,jj) ) 96 94 END_2D 97 95 ! ! =============== … … 134 132 CALL dyn_ldf_lap( kt, Kbb, Kmm, pu, pv, zulap, zvlap, 1 ) ! rotated laplacian applied to pt (output in zlap,Kbb) 135 133 ! 136 CALL lbc_lnk_multi( 'dynldf_lap_blp', zulap, 'U', -1. , zvlap, 'V', -1.) ! Lateral boundary conditions134 CALL lbc_lnk_multi( 'dynldf_lap_blp', zulap, 'U', -1.0_wp, zvlap, 'V', -1.0_wp ) ! Lateral boundary conditions 137 135 ! 138 136 CALL dyn_ldf_lap( kt, Kbb, Kmm, zulap, zvlap, pu_rhs, pv_rhs, 2 ) ! rotated laplacian applied to zlap (output in pt(:,:,:,:,Krhs)) -
NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/OCE/DYN/dynspg.F90
r12991 r13710 96 96 .OR. ln_ice_embd ) THEN ! embedded sea-ice 97 97 ! 98 DO_2D _00_0098 DO_2D( 0, 0, 0, 0 ) 99 99 spgu(ji,jj) = 0._wp 100 100 spgv(ji,jj) = 0._wp … … 103 103 IF( ln_apr_dyn .AND. .NOT.ln_dynspg_ts ) THEN !== Atmospheric pressure gradient (added later in time-split case) ==! 104 104 zg_2 = grav * 0.5 105 DO_2D _00_00105 DO_2D( 0, 0, 0, 0 ) ! gradient of Patm using inverse barometer ssh 106 106 spgu(ji,jj) = spgu(ji,jj) + zg_2 * ( ssh_ib (ji+1,jj) - ssh_ib (ji,jj) & 107 107 & + ssh_ibb(ji+1,jj) - ssh_ibb(ji,jj) ) * r1_e1u(ji,jj) … … 118 118 CALL upd_tide(zt0step, Kmm) 119 119 ! 120 DO_2D _00_00120 DO_2D( 0, 0, 0, 0 ) ! add tide potential forcing 121 121 spgu(ji,jj) = spgu(ji,jj) + grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) * r1_e1u(ji,jj) 122 122 spgv(ji,jj) = spgv(ji,jj) + grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) * r1_e2v(ji,jj) … … 125 125 IF (ln_scal_load) THEN 126 126 zld = rn_scal_load * grav 127 DO_2D _00_00127 DO_2D( 0, 0, 0, 0 ) ! add scalar approximation for load potential 128 128 spgu(ji,jj) = spgu(ji,jj) + zld * ( pssh(ji+1,jj,Kmm) - pssh(ji,jj,Kmm) ) * r1_e1u(ji,jj) 129 129 spgv(ji,jj) = spgv(ji,jj) + zld * ( pssh(ji,jj+1,Kmm) - pssh(ji,jj,Kmm) ) * r1_e2v(ji,jj) … … 137 137 zgrho0r = - grav * r1_rho0 138 138 zpice(:,:) = ( zintp * snwice_mass(:,:) + ( 1.- zintp ) * snwice_mass_b(:,:) ) * zgrho0r 139 DO_2D _00_00139 DO_2D( 0, 0, 0, 0 ) 140 140 spgu(ji,jj) = spgu(ji,jj) + ( zpice(ji+1,jj) - zpice(ji,jj) ) * r1_e1u(ji,jj) 141 141 spgv(ji,jj) = spgv(ji,jj) + ( zpice(ji,jj+1) - zpice(ji,jj) ) * r1_e2v(ji,jj) … … 145 145 ! 146 146 IF( ln_wave .and. ln_bern_srfc ) THEN !== Add J terms: depth-independent Bernoulli head 147 DO_2D _00_00147 DO_2D( 0, 0, 0, 0 ) 148 148 spgu(ji,jj) = spgu(ji,jj) + ( bhd_wave(ji+1,jj) - bhd_wave(ji,jj) ) / e1u(ji,jj) !++ bhd_wave from wave model in m2/s2 [BHD parameters in WW3] 149 149 spgv(ji,jj) = spgv(ji,jj) + ( bhd_wave(ji,jj+1) - bhd_wave(ji,jj) ) / e2v(ji,jj) … … 151 151 ENDIF 152 152 ! 153 DO_3D _00_00( 1, jpkm1 )153 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) !== Add all terms to the general trend 154 154 puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + spgu(ji,jj) 155 155 pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + spgv(ji,jj) -
NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/OCE/DYN/dynspg_exp.F90
r12489 r13710 74 74 IF( ln_linssh ) THEN !* linear free surface : add the surface pressure gradient trend 75 75 ! 76 DO_2D _00_0076 DO_2D( 0, 0, 0, 0 ) ! now surface pressure gradient 77 77 spgu(ji,jj) = - grav * ( ssh(ji+1,jj,Kmm) - ssh(ji,jj,Kmm) ) * r1_e1u(ji,jj) 78 78 spgv(ji,jj) = - grav * ( ssh(ji,jj+1,Kmm) - ssh(ji,jj,Kmm) ) * r1_e2v(ji,jj) 79 79 END_2D 80 80 ! 81 DO_3D _00_00( 1, jpkm1 )81 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) ! Add it to the general trend 82 82 puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + spgu(ji,jj) 83 83 pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + spgv(ji,jj) -
NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/OCE/DYN/dynspg_ts.F90
r12489 r13710 87 87 !! * Substitutions 88 88 # include "do_loop_substitute.h90" 89 # include "domzgr_substitute.h90" 89 90 !!---------------------------------------------------------------------- 90 91 !! NEMO/OCE 4.0 , NEMO Consortium (2018) … … 161 162 REAL(wp), DIMENSION(jpi,jpj) :: zCdU_u, zCdU_v ! top/bottom stress at u- & v-points 162 163 REAL(wp), DIMENSION(jpi,jpj) :: zhU, zhV ! fluxes 164 REAL(wp), DIMENSION(jpi, jpj, jpk) :: ze3u, ze3v 163 165 ! 164 166 REAL(wp) :: zwdramp ! local scalar - only used if ln_wd_dl = .True. … … 227 229 ! != zu_frc = 1/H e3*d/dt(Ua) =! (Vertical mean of Ua, the 3D trends) 228 230 ! ! --------------------------- ! 229 zu_frc(:,:) = SUM( e3u(:,:,:,Kmm) * uu(:,:,:,Krhs) * umask(:,:,:) , DIM=3 ) * r1_hu(:,:,Kmm) 230 zv_frc(:,:) = SUM( e3v(:,:,:,Kmm) * vv(:,:,:,Krhs) * vmask(:,:,:) , DIM=3 ) * r1_hv(:,:,Kmm) 231 DO jk = 1 , jpk 232 ze3u(:,:,jk) = e3u(:,:,jk,Kmm) 233 ze3v(:,:,jk) = e3v(:,:,jk,Kmm) 234 END DO 235 ! 236 zu_frc(:,:) = SUM( ze3u(:,:,:) * uu(:,:,:,Krhs) * umask(:,:,:) , DIM=3 ) * r1_hu(:,:,Kmm) 237 zv_frc(:,:) = SUM( ze3v(:,:,:) * vv(:,:,:,Krhs) * vmask(:,:,:) , DIM=3 ) * r1_hv(:,:,Kmm) 231 238 ! 232 239 ! … … 250 257 zhV(:,:) = pvv_b(:,:,Kmm) * hv(:,:,Kmm) * e1v(:,:) ! NB: FULL domain : put a value in last row and column 251 258 ! 252 CALL dyn_cor_2d( h u(:,:,Kmm), hv(:,:,Kmm), puu_b(:,:,Kmm), pvv_b(:,:,Kmm), zhU, zhV, & ! <<== in259 CALL dyn_cor_2d( ht(:,:), hu(:,:,Kmm), hv(:,:,Kmm), puu_b(:,:,Kmm), pvv_b(:,:,Kmm), zhU, zhV, & ! <<== in 253 260 & zu_trd, zv_trd ) ! ==>> out 254 261 ! … … 257 264 IF( ln_wd_il ) THEN ! W/D : limiter applied to spgspg 258 265 CALL wad_spg( pssh(:,:,Kmm), zcpx, zcpy ) ! Calculating W/D gravity filters, zcpx and zcpy 259 DO_2D _00_00266 DO_2D( 0, 0, 0, 0 ) ! SPG with the application of W/D gravity filters 260 267 zu_trd(ji,jj) = zu_trd(ji,jj) - grav * ( pssh(ji+1,jj ,Kmm) - pssh(ji ,jj ,Kmm) ) & 261 268 & * r1_e1u(ji,jj) * zcpx(ji,jj) * wdrampu(ji,jj) !jth … … 264 271 END_2D 265 272 ELSE ! now suface pressure gradient 266 DO_2D _00_00273 DO_2D( 0, 0, 0, 0 ) 267 274 zu_trd(ji,jj) = zu_trd(ji,jj) - grav * ( pssh(ji+1,jj ,Kmm) - pssh(ji ,jj ,Kmm) ) * r1_e1u(ji,jj) 268 275 zv_trd(ji,jj) = zv_trd(ji,jj) - grav * ( pssh(ji ,jj+1,Kmm) - pssh(ji ,jj ,Kmm) ) * r1_e2v(ji,jj) … … 272 279 ENDIF 273 280 ! 274 DO_2D _00_00281 DO_2D( 0, 0, 0, 0 ) ! Remove coriolis term (and possibly spg) from barotropic trend 275 282 zu_frc(ji,jj) = zu_frc(ji,jj) - zu_trd(ji,jj) * ssumask(ji,jj) 276 283 zv_frc(ji,jj) = zv_frc(ji,jj) - zv_trd(ji,jj) * ssvmask(ji,jj) … … 284 291 IF( ln_apr_dyn ) THEN 285 292 IF( ln_bt_fw ) THEN ! FORWARD integration: use kt+1/2 pressure (NOW+1/2) 286 DO_2D _00_00293 DO_2D( 0, 0, 0, 0 ) 287 294 zu_frc(ji,jj) = zu_frc(ji,jj) + grav * ( ssh_ib (ji+1,jj ) - ssh_ib (ji,jj) ) * r1_e1u(ji,jj) 288 295 zv_frc(ji,jj) = zv_frc(ji,jj) + grav * ( ssh_ib (ji ,jj+1) - ssh_ib (ji,jj) ) * r1_e2v(ji,jj) … … 290 297 ELSE ! CENTRED integration: use kt-1/2 + kt+1/2 pressure (NOW) 291 298 zztmp = grav * r1_2 292 DO_2D _00_00299 DO_2D( 0, 0, 0, 0 ) 293 300 zu_frc(ji,jj) = zu_frc(ji,jj) + zztmp * ( ssh_ib (ji+1,jj ) - ssh_ib (ji,jj) & 294 301 & + ssh_ibb(ji+1,jj ) - ssh_ibb(ji,jj) ) * r1_e1u(ji,jj) … … 302 309 ! ! ---------------------------------- ! 303 310 IF( ln_bt_fw ) THEN ! Add wind forcing 304 DO_2D _00_00311 DO_2D( 0, 0, 0, 0 ) 305 312 zu_frc(ji,jj) = zu_frc(ji,jj) + r1_rho0 * utau(ji,jj) * r1_hu(ji,jj,Kmm) 306 313 zv_frc(ji,jj) = zv_frc(ji,jj) + r1_rho0 * vtau(ji,jj) * r1_hv(ji,jj,Kmm) … … 308 315 ELSE 309 316 zztmp = r1_rho0 * r1_2 310 DO_2D _00_00317 DO_2D( 0, 0, 0, 0 ) 311 318 zu_frc(ji,jj) = zu_frc(ji,jj) + zztmp * ( utau_b(ji,jj) + utau(ji,jj) ) * r1_hu(ji,jj,Kmm) 312 319 zv_frc(ji,jj) = zv_frc(ji,jj) + zztmp * ( vtau_b(ji,jj) + vtau(ji,jj) ) * r1_hv(ji,jj,Kmm) … … 468 475 ! 469 476 ! ! ocean u- and v-depth at mid-step (separate DO-loops remove the need of a lbc_lnk) 470 DO_2D _11_10477 DO_2D( 1, 1, 1, 0 ) ! not jpi-column 471 478 zhup2_e(ji,jj) = hu_0(ji,jj) + r1_2 * r1_e1e2u(ji,jj) & 472 479 & * ( e1e2t(ji ,jj) * zsshp2_e(ji ,jj) & 473 480 & + e1e2t(ji+1,jj) * zsshp2_e(ji+1,jj) ) * ssumask(ji,jj) 474 481 END_2D 475 DO_2D _10_11482 DO_2D( 1, 0, 1, 1 ) ! not jpj-row 476 483 zhvp2_e(ji,jj) = hv_0(ji,jj) + r1_2 * r1_e1e2v(ji,jj) & 477 484 & * ( e1e2t(ji,jj ) * zsshp2_e(ji,jj ) & … … 508 515 !-- ssh = ssh - delta_t' * [ frc + div( flux ) ] --! 509 516 !-------------------------------------------------------------------------! 510 DO_2D _00_00517 DO_2D( 0, 0, 0, 0 ) 511 518 zhdiv = ( zhU(ji,jj) - zhU(ji-1,jj) + zhV(ji,jj) - zhV(ji,jj-1) ) * r1_e1e2t(ji,jj) 512 519 ssha_e(ji,jj) = ( sshn_e(ji,jj) - rDt_e * ( zssh_frc(ji,jj) + zhdiv ) ) * ssmask(ji,jj) … … 514 521 ! 515 522 CALL lbc_lnk_multi( 'dynspg_ts', ssha_e, 'T', 1._wp, zhU, 'U', -1._wp, zhV, 'V', -1._wp ) 523 ! 524 ! Duplicate sea level across open boundaries (this is only cosmetic if linssh=T) 525 IF( ln_bdy ) CALL bdy_ssh( ssha_e ) 526 #if defined key_agrif 527 IF( .NOT.Agrif_Root() ) CALL agrif_ssh_ts( jn ) 528 #endif 516 529 ! 517 530 ! ! Sum over sub-time-steps to compute advective velocities … … 525 538 END IF 526 539 ! 527 ! Duplicate sea level across open boundaries (this is only cosmetic if linssh=T)528 IF( ln_bdy ) CALL bdy_ssh( ssha_e )529 #if defined key_agrif530 IF( .NOT.Agrif_Root() ) CALL agrif_ssh_ts( jn )531 #endif532 540 ! 533 541 ! Sea Surface Height at u-,v-points (vvl case only) 534 542 IF( .NOT.ln_linssh ) THEN 535 DO_2D _00_00543 DO_2D( 0, 0, 0, 0 ) 536 544 zsshu_a(ji,jj) = r1_2 * ssumask(ji,jj) * r1_e1e2u(ji,jj) & 537 545 & * ( e1e2t(ji ,jj ) * ssha_e(ji ,jj ) & … … 553 561 ! ! Surface pressure gradient 554 562 zldg = ( 1._wp - rn_scal_load ) * grav ! local factor 555 DO_2D _00_00563 DO_2D( 0, 0, 0, 0 ) 556 564 zu_spg(ji,jj) = - zldg * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) * r1_e1u(ji,jj) 557 565 zv_spg(ji,jj) = - zldg * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) * r1_e2v(ji,jj) … … 567 575 ! at each time step. We however keep them constant here for optimization. 568 576 ! Recall that zhU and zhV hold fluxes at jn+0.5 (extrapolated not backward interpolated) 569 CALL dyn_cor_2d( zh up2_e, zhvp2_e, ua_e, va_e, zhU, zhV, zu_trd, zv_trd )577 CALL dyn_cor_2d( zhtp2_e, zhup2_e, zhvp2_e, ua_e, va_e, zhU, zhV, zu_trd, zv_trd ) 570 578 ! 571 579 ! Add tidal astronomical forcing if defined 572 580 IF ( ln_tide .AND. ln_tide_pot ) THEN 573 DO_2D _00_00581 DO_2D( 0, 0, 0, 0 ) 574 582 zu_trd(ji,jj) = zu_trd(ji,jj) + grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) * r1_e1u(ji,jj) 575 583 zv_trd(ji,jj) = zv_trd(ji,jj) + grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) * r1_e2v(ji,jj) … … 580 588 !jth do implicitly instead 581 589 IF ( .NOT. ll_wd ) THEN ! Revert to explicit for bit comparison tests in non wad runs 582 DO_2D _00_00590 DO_2D( 0, 0, 0, 0 ) 583 591 zu_trd(ji,jj) = zu_trd(ji,jj) + zCdU_u(ji,jj) * un_e(ji,jj) * hur_e(ji,jj) 584 592 zv_trd(ji,jj) = zv_trd(ji,jj) + zCdU_v(ji,jj) * vn_e(ji,jj) * hvr_e(ji,jj) … … 598 606 !------------------------------------------------------------------------------------------------------------------------! 599 607 IF( ln_dynadv_vec .OR. ln_linssh ) THEN !* Vector form 600 DO_2D _00_00608 DO_2D( 0, 0, 0, 0 ) 601 609 ua_e(ji,jj) = ( un_e(ji,jj) & 602 610 & + rDt_e * ( zu_spg(ji,jj) & … … 613 621 ! 614 622 ELSE !* Flux form 615 DO_2D _00_00623 DO_2D( 0, 0, 0, 0 ) 616 624 ! ! hu_e, hv_e hold depth at jn, zhup2_e, zhvp2_e hold extrapolated depth at jn+1/2 617 625 ! ! backward interpolated depth used in spg terms at jn+1/2 … … 637 645 !jth implicit bottom friction: 638 646 IF ( ll_wd ) THEN ! revert to explicit for bit comparison tests in non wad runs 639 DO_2D _00_00647 DO_2D( 0, 0, 0, 0 ) 640 648 ua_e(ji,jj) = ua_e(ji,jj) /(1.0 - rDt_e * zCdU_u(ji,jj) * hur_e(ji,jj)) 641 649 va_e(ji,jj) = va_e(ji,jj) /(1.0 - rDt_e * zCdU_v(ji,jj) * hvr_e(ji,jj)) … … 643 651 ENDIF 644 652 645 IF( .NOT.ln_linssh ) THEN 653 IF( .NOT.ln_linssh ) THEN !* Update ocean depth (variable volume case only) 646 654 hu_e (2:jpim1,2:jpjm1) = hu_0(2:jpim1,2:jpjm1) + zsshu_a(2:jpim1,2:jpjm1) 647 655 hur_e(2:jpim1,2:jpjm1) = ssumask(2:jpim1,2:jpjm1) / ( hu_e(2:jpim1,2:jpjm1) + 1._wp - ssumask(2:jpim1,2:jpjm1) ) 648 656 hv_e (2:jpim1,2:jpjm1) = hv_0(2:jpim1,2:jpjm1) + zsshv_a(2:jpim1,2:jpjm1) 649 657 hvr_e(2:jpim1,2:jpjm1) = ssvmask(2:jpim1,2:jpjm1) / ( hv_e(2:jpim1,2:jpjm1) + 1._wp - ssvmask(2:jpim1,2:jpjm1) ) 658 ENDIF 659 ! 660 IF( .NOT.ln_linssh ) THEN !* Update ocean depth (variable volume case only) 650 661 CALL lbc_lnk_multi( 'dynspg_ts', ua_e , 'U', -1._wp, va_e , 'V', -1._wp & 651 662 & , hu_e , 'U', 1._wp, hv_e , 'V', 1._wp & … … 654 665 CALL lbc_lnk_multi( 'dynspg_ts', ua_e , 'U', -1._wp, va_e , 'V', -1._wp ) 655 666 ENDIF 656 !657 !658 667 ! ! open boundaries 659 668 IF( ln_bdy ) CALL bdy_dyn2d( jn, ua_e, va_e, un_e, vn_e, hur_e, hvr_e, ssha_e ) … … 703 712 IF (ln_bt_fw) THEN 704 713 IF( .NOT.( kt == nit000 .AND. l_1st_euler ) ) THEN 705 DO_2D _11_11714 DO_2D( 1, 1, 1, 1 ) 706 715 zun_save = un_adv(ji,jj) 707 716 zvn_save = vn_adv(ji,jj) … … 734 743 ELSE 735 744 ! At this stage, pssh(:,:,:,Krhs) has been corrected: compute new depths at velocity points 736 DO_2D _10_10745 DO_2D( 1, 0, 1, 0 ) 737 746 zsshu_a(ji,jj) = r1_2 * ssumask(ji,jj) * r1_e1e2u(ji,jj) & 738 747 & * ( e1e2t(ji ,jj) * pssh(ji ,jj,Kaa) & … … 891 900 ! ! --------------- 892 901 IF( ln_rstart .AND. ln_bt_fw .AND. (.NOT.l_1st_euler) ) THEN !* Read the restart file 893 CALL iom_get( numror, jpdom_auto glo, 'ub2_b' , ub2_b (:,:), ldxios = lrxios )894 CALL iom_get( numror, jpdom_auto glo, 'vb2_b' , vb2_b (:,:), ldxios = lrxios )895 CALL iom_get( numror, jpdom_auto glo, 'un_bf' , un_bf (:,:), ldxios = lrxios )896 CALL iom_get( numror, jpdom_auto glo, 'vn_bf' , vn_bf (:,:), ldxios = lrxios )902 CALL iom_get( numror, jpdom_auto, 'ub2_b' , ub2_b (:,:), cd_type = 'U', psgn = -1._wp, ldxios = lrxios ) 903 CALL iom_get( numror, jpdom_auto, 'vb2_b' , vb2_b (:,:), cd_type = 'V', psgn = -1._wp, ldxios = lrxios ) 904 CALL iom_get( numror, jpdom_auto, 'un_bf' , un_bf (:,:), cd_type = 'U', psgn = -1._wp, ldxios = lrxios ) 905 CALL iom_get( numror, jpdom_auto, 'vn_bf' , vn_bf (:,:), cd_type = 'V', psgn = -1._wp, ldxios = lrxios ) 897 906 IF( .NOT.ln_bt_av ) THEN 898 CALL iom_get( numror, jpdom_auto glo, 'sshbb_e' , sshbb_e(:,:), ldxios = lrxios )899 CALL iom_get( numror, jpdom_auto glo, 'ubb_e' , ubb_e(:,:), ldxios = lrxios )900 CALL iom_get( numror, jpdom_auto glo, 'vbb_e' , vbb_e(:,:), ldxios = lrxios )901 CALL iom_get( numror, jpdom_auto glo, 'sshb_e' , sshb_e(:,:), ldxios = lrxios )902 CALL iom_get( numror, jpdom_auto glo, 'ub_e' , ub_e(:,:), ldxios = lrxios )903 CALL iom_get( numror, jpdom_auto glo, 'vb_e' , vb_e(:,:), ldxios = lrxios )907 CALL iom_get( numror, jpdom_auto, 'sshbb_e' , sshbb_e(:,:), cd_type = 'T', psgn = 1._wp, ldxios = lrxios ) 908 CALL iom_get( numror, jpdom_auto, 'ubb_e' , ubb_e(:,:), cd_type = 'U', psgn = -1._wp, ldxios = lrxios ) 909 CALL iom_get( numror, jpdom_auto, 'vbb_e' , vbb_e(:,:), cd_type = 'V', psgn = -1._wp, ldxios = lrxios ) 910 CALL iom_get( numror, jpdom_auto, 'sshb_e' , sshb_e(:,:), cd_type = 'T', psgn = 1._wp, ldxios = lrxios ) 911 CALL iom_get( numror, jpdom_auto, 'ub_e' , ub_e(:,:), cd_type = 'U', psgn = -1._wp, ldxios = lrxios ) 912 CALL iom_get( numror, jpdom_auto, 'vb_e' , vb_e(:,:), cd_type = 'V', psgn = -1._wp, ldxios = lrxios ) 904 913 ENDIF 905 914 #if defined key_agrif 906 915 ! Read time integrated fluxes 907 916 IF ( .NOT.Agrif_Root() ) THEN 908 CALL iom_get( numror, jpdom_autoglo, 'ub2_i_b' , ub2_i_b(:,:), ldxios = lrxios ) 909 CALL iom_get( numror, jpdom_autoglo, 'vb2_i_b' , vb2_i_b(:,:), ldxios = lrxios ) 917 CALL iom_get( numror, jpdom_auto, 'ub2_i_b' , ub2_i_b(:,:), cd_type = 'U', psgn = -1._wp, ldxios = lrxios ) 918 CALL iom_get( numror, jpdom_auto, 'vb2_i_b' , vb2_i_b(:,:), cd_type = 'V', psgn = -1._wp, ldxios = lrxios ) 919 ELSE 920 ub2_i_b(:,:) = 0._wp ; vb2_i_b(:,:) = 0._wp ! used in the 1st update of agrif 910 921 ENDIF 911 922 #endif … … 913 924 IF(lwp) WRITE(numout,*) 914 925 IF(lwp) WRITE(numout,*) ' ==>>> start from rest: set barotropic values to 0' 915 ub2_b (:,:) = 0._wp ; vb2_b(:,:) = 0._wp ! used in the 1st interpol of agrif916 un_adv (:,:) = 0._wp ; vn_adv(:,:) = 0._wp ! used in the 1st interpol of agrif917 un_bf (:,:) = 0._wp ; vn_bf(:,:) = 0._wp ! used in the 1st update of agrif926 ub2_b (:,:) = 0._wp ; vb2_b (:,:) = 0._wp ! used in the 1st interpol of agrif 927 un_adv (:,:) = 0._wp ; vn_adv (:,:) = 0._wp ! used in the 1st interpol of agrif 928 un_bf (:,:) = 0._wp ; vn_bf (:,:) = 0._wp ! used in the 1st update of agrif 918 929 #if defined key_agrif 919 IF ( .NOT.Agrif_Root() ) THEN 920 ub2_i_b(:,:) = 0._wp ; vb2_i_b(:,:) = 0._wp ! used in the 1st update of agrif 921 ENDIF 930 ub2_i_b(:,:) = 0._wp ; vb2_i_b(:,:) = 0._wp ! used in the 1st update of agrif 922 931 #endif 923 932 ENDIF … … 966 975 ! Max courant number for ext. grav. waves 967 976 ! 968 DO_2D _11_11977 DO_2D( 0, 0, 0, 0 ) 969 978 zxr2 = r1_e1t(ji,jj) * r1_e1t(ji,jj) 970 979 zyr2 = r1_e2t(ji,jj) * r1_e2t(ji,jj) … … 972 981 END_2D 973 982 ! 974 zcmax = MAXVAL( zcu( :,:) )983 zcmax = MAXVAL( zcu(Nis0:Nie0,Njs0:Nje0) ) 975 984 CALL mpp_max( 'dynspg_ts', zcmax ) 976 985 … … 1088 1097 ! 1089 1098 SELECT CASE( nvor_scheme ) 1090 CASE( np_EEN ) != EEN scheme using e3f (energy & enstrophy scheme)1099 CASE( np_EEN ) != EEN scheme using e3f energy & enstrophy scheme 1091 1100 SELECT CASE( nn_een_e3f ) !* ff_f/e3 at F-point 1092 1101 CASE ( 0 ) ! original formulation (masked averaging of e3t divided by 4) 1093 DO_2D _10_101102 DO_2D( 1, 0, 1, 0 ) 1094 1103 zwz(ji,jj) = ( ht(ji ,jj+1) + ht(ji+1,jj+1) + & 1095 1104 & ht(ji ,jj ) + ht(ji+1,jj ) ) * 0.25_wp … … 1097 1106 END_2D 1098 1107 CASE ( 1 ) ! new formulation (masked averaging of e3t divided by the sum of mask) 1099 DO_2D _10_101108 DO_2D( 1, 0, 1, 0 ) 1100 1109 zwz(ji,jj) = ( ht (ji ,jj+1) + ht (ji+1,jj+1) & 1101 1110 & + ht (ji ,jj ) + ht (ji+1,jj ) ) & … … 1108 1117 ! 1109 1118 ftne(1,:) = 0._wp ; ftnw(1,:) = 0._wp ; ftse(1,:) = 0._wp ; ftsw(1,:) = 0._wp 1110 DO_2D _01_011119 DO_2D( 0, 1, 0, 1 ) 1111 1120 ftne(ji,jj) = zwz(ji-1,jj ) + zwz(ji ,jj ) + zwz(ji ,jj-1) 1112 1121 ftnw(ji,jj) = zwz(ji-1,jj-1) + zwz(ji-1,jj ) + zwz(ji ,jj ) … … 1115 1124 END_2D 1116 1125 ! 1117 CASE( np_EET ) != EEN scheme using e3t (energy conserving scheme)1126 CASE( np_EET ) != EEN scheme using e3t energy conserving scheme 1118 1127 ftne(1,:) = 0._wp ; ftnw(1,:) = 0._wp ; ftse(1,:) = 0._wp ; ftsw(1,:) = 0._wp 1119 DO_2D _01_011128 DO_2D( 0, 1, 0, 1 ) 1120 1129 z1_ht = ssmask(ji,jj) / ( ht(ji,jj) + 1._wp - ssmask(ji,jj) ) 1121 1130 ftne(ji,jj) = ( ff_f(ji-1,jj ) + ff_f(ji ,jj ) + ff_f(ji ,jj-1) ) * z1_ht … … 1150 1159 ! 1151 1160 !zhf(:,:) = hbatf(:,:) 1152 DO_2D _10_101161 DO_2D( 1, 0, 1, 0 ) 1153 1162 zhf(ji,jj) = ( ht_0 (ji,jj ) + ht_0 (ji+1,jj ) & 1154 1163 & + ht_0 (ji,jj+1) + ht_0 (ji+1,jj+1) ) & … … 1169 1178 CALL lbc_lnk( 'dynspg_ts', zhf, 'F', 1._wp ) 1170 1179 ! JC: TBC. hf should be greater than 0 1171 DO_2D _11_111180 DO_2D( 1, 1, 1, 1 ) 1172 1181 IF( zhf(ji,jj) /= 0._wp ) zwz(ji,jj) = 1._wp / zhf(ji,jj) 1173 1182 END_2D … … 1179 1188 1180 1189 1181 SUBROUTINE dyn_cor_2d( ph u, phv, punb, pvnb, zhU, zhV, zu_trd, zv_trd )1190 SUBROUTINE dyn_cor_2d( pht, phu, phv, punb, pvnb, zhU, zhV, zu_trd, zv_trd ) 1182 1191 !!--------------------------------------------------------------------- 1183 1192 !! *** ROUTINE dyn_cor_2d *** … … 1187 1196 INTEGER :: ji ,jj ! dummy loop indices 1188 1197 REAL(wp) :: zx1, zx2, zy1, zy2, z1_hu, z1_hv ! - - 1189 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: ph u, phv, punb, pvnb, zhU, zhV1198 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pht, phu, phv, punb, pvnb, zhU, zhV 1190 1199 REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: zu_trd, zv_trd 1191 1200 !!---------------------------------------------------------------------- 1192 1201 SELECT CASE( nvor_scheme ) 1193 1202 CASE( np_ENT ) ! enstrophy conserving scheme (f-point) 1194 DO_2D _00_001203 DO_2D( 0, 0, 0, 0 ) 1195 1204 z1_hu = ssumask(ji,jj) / ( phu(ji,jj) + 1._wp - ssumask(ji,jj) ) 1196 1205 z1_hv = ssvmask(ji,jj) / ( phv(ji,jj) + 1._wp - ssvmask(ji,jj) ) 1197 1206 zu_trd(ji,jj) = + r1_4 * r1_e1e2u(ji,jj) * z1_hu & 1198 & * ( e1e2t(ji+1,jj)* ht(ji+1,jj)*ff_t(ji+1,jj) * ( pvnb(ji+1,jj) + pvnb(ji+1,jj-1) ) &1199 & + e1e2t(ji ,jj)* ht(ji ,jj)*ff_t(ji ,jj) * ( pvnb(ji ,jj) + pvnb(ji ,jj-1) ) )1207 & * ( e1e2t(ji+1,jj)*pht(ji+1,jj)*ff_t(ji+1,jj) * ( pvnb(ji+1,jj) + pvnb(ji+1,jj-1) ) & 1208 & + e1e2t(ji ,jj)*pht(ji ,jj)*ff_t(ji ,jj) * ( pvnb(ji ,jj) + pvnb(ji ,jj-1) ) ) 1200 1209 ! 1201 1210 zv_trd(ji,jj) = - r1_4 * r1_e1e2v(ji,jj) * z1_hv & 1202 & * ( e1e2t(ji,jj+1)* ht(ji,jj+1)*ff_t(ji,jj+1) * ( punb(ji,jj+1) + punb(ji-1,jj+1) ) &1203 & + e1e2t(ji,jj )* ht(ji,jj )*ff_t(ji,jj ) * ( punb(ji,jj ) + punb(ji-1,jj ) ) )1211 & * ( e1e2t(ji,jj+1)*pht(ji,jj+1)*ff_t(ji,jj+1) * ( punb(ji,jj+1) + punb(ji-1,jj+1) ) & 1212 & + e1e2t(ji,jj )*pht(ji,jj )*ff_t(ji,jj ) * ( punb(ji,jj ) + punb(ji-1,jj ) ) ) 1204 1213 END_2D 1205 1214 ! 1206 1215 CASE( np_ENE , np_MIX ) ! energy conserving scheme (t-point) ENE or MIX 1207 DO_2D _00_001216 DO_2D( 0, 0, 0, 0 ) 1208 1217 zy1 = ( zhV(ji,jj-1) + zhV(ji+1,jj-1) ) * r1_e1u(ji,jj) 1209 1218 zy2 = ( zhV(ji,jj ) + zhV(ji+1,jj ) ) * r1_e1u(ji,jj) … … 1216 1225 ! 1217 1226 CASE( np_ENS ) ! enstrophy conserving scheme (f-point) 1218 DO_2D _00_001227 DO_2D( 0, 0, 0, 0 ) 1219 1228 zy1 = r1_8 * ( zhV(ji ,jj-1) + zhV(ji+1,jj-1) & 1220 1229 & + zhV(ji ,jj ) + zhV(ji+1,jj ) ) * r1_e1u(ji,jj) … … 1226 1235 ! 1227 1236 CASE( np_EET , np_EEN ) ! energy & enstrophy scheme (using e3t or e3f) 1228 DO_2D _00_001237 DO_2D( 0, 0, 0, 0 ) 1229 1238 zu_trd(ji,jj) = + r1_12 * r1_e1u(ji,jj) * ( ftne(ji,jj ) * zhV(ji ,jj ) & 1230 1239 & + ftnw(ji+1,jj) * zhV(ji+1,jj ) & … … 1260 1269 ! 1261 1270 IF( ln_wd_dl_rmp ) THEN 1262 DO_2D _11_111271 DO_2D( 1, 1, 1, 1 ) 1263 1272 IF ( pssh(ji,jj) + ht_0(ji,jj) > 2._wp * rn_wdmin1 ) THEN 1264 1273 ! IF ( pssh(ji,jj) + ht_0(ji,jj) > rn_wdmin2 ) THEN … … 1271 1280 END_2D 1272 1281 ELSE 1273 DO_2D _11_111282 DO_2D( 1, 1, 1, 1 ) 1274 1283 IF ( pssh(ji,jj) + ht_0(ji,jj) > rn_wdmin1 ) THEN ; ptmsk(ji,jj) = 1._wp 1275 1284 ELSE ; ptmsk(ji,jj) = 0._wp … … 1299 1308 !!---------------------------------------------------------------------- 1300 1309 ! 1301 DO_2D _11_101310 DO_2D( 1, 1, 1, 0 ) ! not jpi-column 1302 1311 IF ( phU(ji,jj) > 0._wp ) THEN ; pUmsk(ji,jj) = pTmsk(ji ,jj) 1303 1312 ELSE ; pUmsk(ji,jj) = pTmsk(ji+1,jj) … … 1307 1316 END_2D 1308 1317 ! 1309 DO_2D _10_111318 DO_2D( 1, 0, 1, 1 ) ! not jpj-row 1310 1319 IF ( phV(ji,jj) > 0._wp ) THEN ; pVmsk(ji,jj) = pTmsk(ji,jj ) 1311 1320 ELSE ; pVmsk(ji,jj) = pTmsk(ji,jj+1) … … 1329 1338 REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: zcpx, zcpy 1330 1339 !!---------------------------------------------------------------------- 1331 DO_2D _00_001340 DO_2D( 0, 0, 0, 0 ) 1332 1341 ll_tmp1 = MIN( pshn(ji,jj) , pshn(ji+1,jj) ) > & 1333 1342 & MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) .AND. & … … 1396 1405 ! !== Set the barotropic drag coef. ==! 1397 1406 ! 1398 IF( ln_isfcav ) THEN ! top+bottom friction (ocean cavities)1407 IF( ln_isfcav.OR.ln_drgice_imp ) THEN ! top+bottom friction (ocean cavities) 1399 1408 1400 DO_2D _00_001409 DO_2D( 0, 0, 0, 0 ) 1401 1410 pCdU_u(ji,jj) = r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) + rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) 1402 1411 pCdU_v(ji,jj) = r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) + rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) 1403 1412 END_2D 1404 1413 ELSE ! bottom friction only 1405 DO_2D _00_001414 DO_2D( 0, 0, 0, 0 ) 1406 1415 pCdU_u(ji,jj) = r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) 1407 1416 pCdU_v(ji,jj) = r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) … … 1413 1422 IF( ln_bt_fw ) THEN ! FORWARD integration: use NOW bottom baroclinic velocities 1414 1423 1415 DO_2D _00_001424 DO_2D( 0, 0, 0, 0 ) 1416 1425 ikbu = mbku(ji,jj) 1417 1426 ikbv = mbkv(ji,jj) … … 1421 1430 ELSE ! CENTRED integration: use BEFORE bottom baroclinic velocities 1422 1431 1423 DO_2D _00_001432 DO_2D( 0, 0, 0, 0 ) 1424 1433 ikbu = mbku(ji,jj) 1425 1434 ikbv = mbkv(ji,jj) … … 1431 1440 IF( ln_wd_il ) THEN ! W/D : use the "clipped" bottom friction !!gm explain WHY, please ! 1432 1441 zztmp = -1._wp / rDt_e 1433 DO_2D _00_001442 DO_2D( 0, 0, 0, 0 ) 1434 1443 pu_RHSi(ji,jj) = pu_RHSi(ji,jj) + zu_i(ji,jj) * wdrampu(ji,jj) * MAX( & 1435 1444 & r1_hu(ji,jj,Kmm) * r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) , zztmp ) … … 1439 1448 ELSE ! use "unclipped" drag (even if explicit friction is used in 3D calculation) 1440 1449 1441 DO_2D _00_001450 DO_2D( 0, 0, 0, 0 ) 1442 1451 pu_RHSi(ji,jj) = pu_RHSi(ji,jj) + r1_hu(ji,jj,Kmm) * r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) * zu_i(ji,jj) 1443 1452 pv_RHSi(ji,jj) = pv_RHSi(ji,jj) + r1_hv(ji,jj,Kmm) * r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) * zv_i(ji,jj) … … 1447 1456 ! !== TOP stress contribution from baroclinic velocities ==! (no W/D case) 1448 1457 ! 1449 IF( ln_isfcav ) THEN1458 IF( ln_isfcav.OR.ln_drgice_imp ) THEN 1450 1459 ! 1451 1460 IF( ln_bt_fw ) THEN ! FORWARD integration: use NOW top baroclinic velocity 1452 1461 1453 DO_2D _00_001462 DO_2D( 0, 0, 0, 0 ) 1454 1463 iktu = miku(ji,jj) 1455 1464 iktv = mikv(ji,jj) … … 1459 1468 ELSE ! CENTRED integration: use BEFORE top baroclinic velocity 1460 1469 1461 DO_2D _00_001470 DO_2D( 0, 0, 0, 0 ) 1462 1471 iktu = miku(ji,jj) 1463 1472 iktv = mikv(ji,jj) … … 1469 1478 ! ! use "unclipped" top drag (even if explicit friction is used in 3D calculation) 1470 1479 1471 DO_2D _00_001480 DO_2D( 0, 0, 0, 0 ) 1472 1481 pu_RHSi(ji,jj) = pu_RHSi(ji,jj) + r1_hu(ji,jj,Kmm) * r1_2*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * zu_i(ji,jj) 1473 1482 pv_RHSi(ji,jj) = pv_RHSi(ji,jj) + r1_hv(ji,jj,Kmm) * r1_2*( rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) * zv_i(ji,jj) -
NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/OCE/DYN/dynvor.F90
r12991 r13710 81 81 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: di_e2u_2 ! = di(e2u)/2 used in T-point metric term calculation 82 82 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: dj_e1v_2 ! = dj(e1v)/2 - - - - 83 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: di_e2v_2e1e2f ! = di(e2 u)/(2*e1e2f) used in F-point metric term calculation84 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: dj_e1u_2e1e2f ! = dj(e1 v)/(2*e1e2f) - - - -83 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: di_e2v_2e1e2f ! = di(e2v)/(2*e1e2f) used in F-point metric term calculation 84 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: dj_e1u_2e1e2f ! = dj(e1u)/(2*e1e2f) - - - - 85 85 86 86 REAL(wp) :: r1_4 = 0.250_wp ! =1/4 … … 90 90 !! * Substitutions 91 91 # include "do_loop_substitute.h90" 92 # include "domzgr_substitute.h90" 93 92 94 !!---------------------------------------------------------------------- 93 95 !! NEMO/OCE 4.0 , NEMO Consortium (2018) … … 233 235 INTEGER :: ji, jj, jk ! dummy loop indices 234 236 REAL(wp) :: zx1, zy1, zx2, zy2 ! local scalars 235 REAL(wp), DIMENSION(jpi,jpj) :: zwx, zwy, zwt ! 2D workspace236 REAL(wp), DIMENSION(jpi,jpj,jpk ) :: zwz ! 3D workspace237 REAL(wp), DIMENSION(jpi,jpj) :: zwx, zwy, zwt ! 2D workspace 238 REAL(wp), DIMENSION(jpi,jpj,jpkm1) :: zwz ! 3D workspace, jpkm1 -> avoid lbc_lnk on jpk that is not defined 237 239 !!---------------------------------------------------------------------- 238 240 ! … … 247 249 CASE ( np_RVO ) !* relative vorticity 248 250 DO jk = 1, jpkm1 ! Horizontal slab 249 DO_2D _10_10251 DO_2D( 1, 0, 1, 0 ) 250 252 zwz(ji,jj,jk) = ( e2v(ji+1,jj) * pv(ji+1,jj,jk) - e2v(ji,jj) * pv(ji,jj,jk) & 251 253 & - e1u(ji,jj+1) * pu(ji,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk) ) * r1_e1e2f(ji,jj) 252 254 END_2D 253 255 IF( ln_dynvor_msk ) THEN ! mask/unmask relative vorticity 254 DO_2D _10_10256 DO_2D( 1, 0, 1, 0 ) 255 257 zwz(ji,jj,jk) = zwz(ji,jj,jk) * fmask(ji,jj,jk) 256 258 END_2D … … 258 260 END DO 259 261 260 CALL lbc_lnk( 'dynvor', zwz, 'F', 1. )262 CALL lbc_lnk( 'dynvor', zwz, 'F', 1.0_wp ) 261 263 262 264 CASE ( np_CRV ) !* Coriolis + relative vorticity 263 265 DO jk = 1, jpkm1 ! Horizontal slab 264 DO_2D _10_10266 DO_2D( 1, 0, 1, 0 ) ! relative vorticity 265 267 zwz(ji,jj,jk) = ( e2v(ji+1,jj) * pv(ji+1,jj,jk) - e2v(ji,jj) * pv(ji,jj,jk) & 266 268 & - e1u(ji,jj+1) * pu(ji,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk) ) * r1_e1e2f(ji,jj) 267 269 END_2D 268 270 IF( ln_dynvor_msk ) THEN ! mask/unmask relative vorticity 269 DO_2D _10_10271 DO_2D( 1, 0, 1, 0 ) 270 272 zwz(ji,jj,jk) = zwz(ji,jj,jk) * fmask(ji,jj,jk) 271 273 END_2D … … 273 275 END DO 274 276 275 CALL lbc_lnk( 'dynvor', zwz, 'F', 1. )277 CALL lbc_lnk( 'dynvor', zwz, 'F', 1.0_wp ) 276 278 277 279 END SELECT … … 285 287 zwt(:,:) = ff_t(:,:) * e1e2t(:,:)*e3t(:,:,jk,Kmm) 286 288 CASE ( np_RVO ) !* relative vorticity 287 DO_2D _01_01289 DO_2D( 0, 1, 0, 1 ) 288 290 zwt(ji,jj) = r1_4 * ( zwz(ji-1,jj ,jk) + zwz(ji,jj ,jk) & 289 & + zwz(ji-1,jj-1,jk) + zwz(ji,jj-1,jk) ) * e1e2t(ji,jj)*e3t(ji,jj,jk,Kmm) 291 & + zwz(ji-1,jj-1,jk) + zwz(ji,jj-1,jk) ) & 292 & * e1e2t(ji,jj)*e3t(ji,jj,jk,Kmm) 290 293 END_2D 291 294 CASE ( np_MET ) !* metric term 292 DO_2D _01_01295 DO_2D( 0, 1, 0, 1 ) 293 296 zwt(ji,jj) = ( ( pv(ji,jj,jk) + pv(ji,jj-1,jk) ) * di_e2u_2(ji,jj) & 294 & - ( pu(ji,jj,jk) + pu(ji-1,jj,jk) ) * dj_e1v_2(ji,jj) ) * e3t(ji,jj,jk,Kmm) 297 & - ( pu(ji,jj,jk) + pu(ji-1,jj,jk) ) * dj_e1v_2(ji,jj) ) & 298 & * e3t(ji,jj,jk,Kmm) 295 299 END_2D 296 300 CASE ( np_CRV ) !* Coriolis + relative vorticity 297 DO_2D _01_01301 DO_2D( 0, 1, 0, 1 ) 298 302 zwt(ji,jj) = ( ff_t(ji,jj) + r1_4 * ( zwz(ji-1,jj ,jk) + zwz(ji,jj ,jk) & 299 & + zwz(ji-1,jj-1,jk) + zwz(ji,jj-1,jk) ) ) * e1e2t(ji,jj)*e3t(ji,jj,jk,Kmm) 303 & + zwz(ji-1,jj-1,jk) + zwz(ji,jj-1,jk) ) ) & 304 & * e1e2t(ji,jj)*e3t(ji,jj,jk,Kmm) 300 305 END_2D 301 306 CASE ( np_CME ) !* Coriolis + metric 302 DO_2D _01_01307 DO_2D( 0, 1, 0, 1 ) 303 308 zwt(ji,jj) = ( ff_t(ji,jj) * e1e2t(ji,jj) & 304 309 & + ( pv(ji,jj,jk) + pv(ji,jj-1,jk) ) * di_e2u_2(ji,jj) & 305 & - ( pu(ji,jj,jk) + pu(ji-1,jj,jk) ) * dj_e1v_2(ji,jj) ) * e3t(ji,jj,jk,Kmm) 310 & - ( pu(ji,jj,jk) + pu(ji-1,jj,jk) ) * dj_e1v_2(ji,jj) ) & 311 & * e3t(ji,jj,jk,Kmm) 306 312 END_2D 307 313 CASE DEFAULT ! error … … 310 316 ! 311 317 ! !== compute and add the vorticity term trend =! 312 DO_2D _00_00318 DO_2D( 0, 0, 0, 0 ) 313 319 pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + r1_4 * r1_e1e2u(ji,jj) / e3u(ji,jj,jk,Kmm) & 314 320 & * ( zwt(ji+1,jj) * ( pv(ji+1,jj,jk) + pv(ji+1,jj-1,jk) ) & … … 370 376 zwz(:,:) = ff_f(:,:) 371 377 CASE ( np_RVO ) !* relative vorticity 372 DO_2D _10_10378 DO_2D( 1, 0, 1, 0 ) 373 379 zwz(ji,jj) = ( e2v(ji+1,jj ) * pv(ji+1,jj ,jk) - e2v(ji,jj) * pv(ji,jj,jk) & 374 380 & - e1u(ji ,jj+1) * pu(ji ,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk) ) * r1_e1e2f(ji,jj) 375 381 END_2D 376 382 CASE ( np_MET ) !* metric term 377 DO_2D _10_10383 DO_2D( 1, 0, 1, 0 ) 378 384 zwz(ji,jj) = ( pv(ji+1,jj ,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj) & 379 385 & - ( pu(ji ,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj) 380 386 END_2D 381 387 CASE ( np_CRV ) !* Coriolis + relative vorticity 382 DO_2D _10_10388 DO_2D( 1, 0, 1, 0 ) 383 389 zwz(ji,jj) = ff_f(ji,jj) + ( e2v(ji+1,jj) * pv(ji+1,jj,jk) - e2v(ji,jj) * pv(ji,jj,jk) & 384 390 & - e1u(ji,jj+1) * pu(ji,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk) ) * r1_e1e2f(ji,jj) 385 391 END_2D 386 392 CASE ( np_CME ) !* Coriolis + metric 387 DO_2D _10_10393 DO_2D( 1, 0, 1, 0 ) 388 394 zwz(ji,jj) = ff_f(ji,jj) + ( pv(ji+1,jj ,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj) & 389 395 & - ( pu(ji ,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj) … … 394 400 ! 395 401 IF( ln_dynvor_msk ) THEN !== mask/unmask vorticity ==! 396 DO_2D _10_10402 DO_2D( 1, 0, 1, 0 ) 397 403 zwz(ji,jj) = zwz(ji,jj) * fmask(ji,jj,jk) 398 404 END_2D … … 408 414 ENDIF 409 415 ! !== compute and add the vorticity term trend =! 410 DO_2D _00_00416 DO_2D( 0, 0, 0, 0 ) 411 417 zy1 = zwy(ji,jj-1) + zwy(ji+1,jj-1) 412 418 zy2 = zwy(ji,jj ) + zwy(ji+1,jj ) … … 466 472 zwz(:,:) = ff_f(:,:) 467 473 CASE ( np_RVO ) !* relative vorticity 468 DO_2D _10_10474 DO_2D( 1, 0, 1, 0 ) 469 475 zwz(ji,jj) = ( e2v(ji+1,jj ) * pv(ji+1,jj ,jk) - e2v(ji,jj) * pv(ji,jj,jk) & 470 476 & - e1u(ji ,jj+1) * pu(ji ,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk) ) * r1_e1e2f(ji,jj) 471 477 END_2D 472 478 CASE ( np_MET ) !* metric term 473 DO_2D _10_10479 DO_2D( 1, 0, 1, 0 ) 474 480 zwz(ji,jj) = ( pv(ji+1,jj ,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj) & 475 481 & - ( pu(ji ,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj) 476 482 END_2D 477 483 CASE ( np_CRV ) !* Coriolis + relative vorticity 478 DO_2D _10_10484 DO_2D( 1, 0, 1, 0 ) 479 485 zwz(ji,jj) = ff_f(ji,jj) + ( e2v(ji+1,jj ) * pv(ji+1,jj ,jk) - e2v(ji,jj) * pv(ji,jj,jk) & 480 486 & - e1u(ji ,jj+1) * pu(ji ,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk) ) * r1_e1e2f(ji,jj) 481 487 END_2D 482 488 CASE ( np_CME ) !* Coriolis + metric 483 DO_2D _10_10489 DO_2D( 1, 0, 1, 0 ) 484 490 zwz(ji,jj) = ff_f(ji,jj) + ( pv(ji+1,jj ,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj) & 485 491 & - ( pu(ji ,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj) … … 490 496 ! 491 497 IF( ln_dynvor_msk ) THEN !== mask/unmask vorticity ==! 492 DO_2D _10_10498 DO_2D( 1, 0, 1, 0 ) 493 499 zwz(ji,jj) = zwz(ji,jj) * fmask(ji,jj,jk) 494 500 END_2D … … 504 510 ENDIF 505 511 ! !== compute and add the vorticity term trend =! 506 DO_2D _00_00512 DO_2D( 0, 0, 0, 0 ) 507 513 zuav = r1_8 * r1_e1u(ji,jj) * ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) & 508 514 & + zwy(ji ,jj ) + zwy(ji+1,jj ) ) … … 545 551 REAL(wp) :: zua, zva ! local scalars 546 552 REAL(wp) :: zmsk, ze3f ! local scalars 547 REAL(wp), DIMENSION(jpi,jpj) :: zwx , zwy , z1_e3f548 REAL(wp), DIMENSION(jpi,jpj) :: ztnw, ztne, ztsw, ztse549 REAL(wp), DIMENSION(jpi,jpj,jpk ) :: zwz553 REAL(wp), DIMENSION(jpi,jpj) :: zwx , zwy , z1_e3f 554 REAL(wp), DIMENSION(jpi,jpj) :: ztnw, ztne, ztsw, ztse 555 REAL(wp), DIMENSION(jpi,jpj,jpkm1) :: zwz ! 3D workspace, jpkm1 -> jpkm1 -> avoid lbc_lnk on jpk that is not defined 550 556 !!---------------------------------------------------------------------- 551 557 ! … … 562 568 SELECT CASE( nn_een_e3f ) ! == reciprocal of e3 at F-point 563 569 CASE ( 0 ) ! original formulation (masked averaging of e3t divided by 4) 564 DO_2D_10_10 565 ze3f = ( e3t(ji,jj+1,jk,Kmm)*tmask(ji,jj+1,jk) + e3t(ji+1,jj+1,jk,Kmm)*tmask(ji+1,jj+1,jk) & 566 & + e3t(ji,jj ,jk,Kmm)*tmask(ji,jj ,jk) + e3t(ji+1,jj ,jk,Kmm)*tmask(ji+1,jj ,jk) ) 570 DO_2D( 1, 0, 1, 0 ) 571 ze3f = ( e3t(ji ,jj+1,jk,Kmm)*tmask(ji ,jj+1,jk) & 572 & + e3t(ji+1,jj+1,jk,Kmm)*tmask(ji+1,jj+1,jk) & 573 & + e3t(ji ,jj ,jk,Kmm)*tmask(ji ,jj ,jk) & 574 & + e3t(ji+1,jj ,jk,Kmm)*tmask(ji+1,jj ,jk) ) 567 575 IF( ze3f /= 0._wp ) THEN ; z1_e3f(ji,jj) = 4._wp / ze3f 568 576 ELSE ; z1_e3f(ji,jj) = 0._wp … … 570 578 END_2D 571 579 CASE ( 1 ) ! new formulation (masked averaging of e3t divided by the sum of mask) 572 DO_2D_10_10 573 ze3f = ( e3t(ji,jj+1,jk,Kmm)*tmask(ji,jj+1,jk) + e3t(ji+1,jj+1,jk,Kmm)*tmask(ji+1,jj+1,jk) & 574 & + e3t(ji,jj ,jk,Kmm)*tmask(ji,jj ,jk) + e3t(ji+1,jj ,jk,Kmm)*tmask(ji+1,jj ,jk) ) 580 DO_2D( 1, 0, 1, 0 ) 581 ze3f = ( e3t(ji ,jj+1,jk,Kmm)*tmask(ji ,jj+1,jk) & 582 & + e3t(ji+1,jj+1,jk,Kmm)*tmask(ji+1,jj+1,jk) & 583 & + e3t(ji ,jj ,jk,Kmm)*tmask(ji ,jj ,jk) & 584 & + e3t(ji+1,jj ,jk,Kmm)*tmask(ji+1,jj ,jk) ) 575 585 zmsk = ( tmask(ji,jj+1,jk) + tmask(ji+1,jj+1,jk) & 576 586 & + tmask(ji,jj ,jk) + tmask(ji+1,jj ,jk) ) … … 583 593 SELECT CASE( kvor ) !== vorticity considered ==! 584 594 CASE ( np_COR ) !* Coriolis (planetary vorticity) 585 DO_2D _10_10595 DO_2D( 1, 0, 1, 0 ) 586 596 zwz(ji,jj,jk) = ff_f(ji,jj) * z1_e3f(ji,jj) 587 597 END_2D 588 598 CASE ( np_RVO ) !* relative vorticity 589 DO_2D _10_10599 DO_2D( 1, 0, 1, 0 ) 590 600 zwz(ji,jj,jk) = ( e2v(ji+1,jj ) * pv(ji+1,jj,jk) - e2v(ji,jj) * pv(ji,jj,jk) & 591 601 & - e1u(ji ,jj+1) * pu(ji,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk) ) * r1_e1e2f(ji,jj)*z1_e3f(ji,jj) 592 602 END_2D 593 603 CASE ( np_MET ) !* metric term 594 DO_2D _10_10604 DO_2D( 1, 0, 1, 0 ) 595 605 zwz(ji,jj,jk) = ( ( pv(ji+1,jj,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj) & 596 606 & - ( pu(ji,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj) ) * z1_e3f(ji,jj) 597 607 END_2D 598 608 CASE ( np_CRV ) !* Coriolis + relative vorticity 599 DO_2D _10_10609 DO_2D( 1, 0, 1, 0 ) 600 610 zwz(ji,jj,jk) = ( ff_f(ji,jj) + ( e2v(ji+1,jj ) * pv(ji+1,jj,jk) - e2v(ji,jj) * pv(ji,jj,jk) & 601 611 & - e1u(ji ,jj+1) * pu(ji,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk) ) & … … 603 613 END_2D 604 614 CASE ( np_CME ) !* Coriolis + metric 605 DO_2D _10_10615 DO_2D( 1, 0, 1, 0 ) 606 616 zwz(ji,jj,jk) = ( ff_f(ji,jj) + ( pv(ji+1,jj ,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj) & 607 617 & - ( pu(ji ,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj) ) * z1_e3f(ji,jj) … … 612 622 ! 613 623 IF( ln_dynvor_msk ) THEN !== mask/unmask vorticity ==! 614 DO_2D _10_10624 DO_2D( 1, 0, 1, 0 ) 615 625 zwz(ji,jj,jk) = zwz(ji,jj,jk) * fmask(ji,jj,jk) 616 626 END_2D … … 618 628 END DO ! End of slab 619 629 ! 620 CALL lbc_lnk( 'dynvor', zwz, 'F', 1. )630 CALL lbc_lnk( 'dynvor', zwz, 'F', 1.0_wp ) 621 631 622 632 DO jk = 1, jpkm1 ! Horizontal slab … … 643 653 END DO 644 654 END DO 645 DO_2D _00_00655 DO_2D( 0, 0, 0, 0 ) 646 656 zua = + r1_12 * r1_e1u(ji,jj) * ( ztne(ji,jj ) * zwy(ji ,jj ) + ztnw(ji+1,jj) * zwy(ji+1,jj ) & 647 657 & + ztse(ji,jj ) * zwy(ji ,jj-1) + ztsw(ji+1,jj) * zwy(ji+1,jj-1) ) … … 685 695 REAL(wp) :: zua, zva ! local scalars 686 696 REAL(wp) :: zmsk, z1_e3t ! local scalars 687 REAL(wp), DIMENSION(jpi,jpj) :: zwx , zwy688 REAL(wp), DIMENSION(jpi,jpj) :: ztnw, ztne, ztsw, ztse689 REAL(wp), DIMENSION(jpi,jpj,jpk ) :: zwz697 REAL(wp), DIMENSION(jpi,jpj) :: zwx , zwy 698 REAL(wp), DIMENSION(jpi,jpj) :: ztnw, ztne, ztsw, ztse 699 REAL(wp), DIMENSION(jpi,jpj,jpkm1) :: zwz ! 3D workspace, avoid lbc_lnk on jpk that is not defined 690 700 !!---------------------------------------------------------------------- 691 701 ! … … 703 713 SELECT CASE( kvor ) !== vorticity considered ==! 704 714 CASE ( np_COR ) !* Coriolis (planetary vorticity) 705 DO_2D _10_10715 DO_2D( 1, 0, 1, 0 ) 706 716 zwz(ji,jj,jk) = ff_f(ji,jj) 707 717 END_2D 708 718 CASE ( np_RVO ) !* relative vorticity 709 DO_2D _10_10719 DO_2D( 1, 0, 1, 0 ) 710 720 zwz(ji,jj,jk) = ( e2v(ji+1,jj ) * pv(ji+1,jj ,jk) - e2v(ji,jj) * pv(ji,jj,jk) & 711 721 & - e1u(ji ,jj+1) * pu(ji ,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk) ) & … … 713 723 END_2D 714 724 CASE ( np_MET ) !* metric term 715 DO_2D _10_10725 DO_2D( 1, 0, 1, 0 ) 716 726 zwz(ji,jj,jk) = ( pv(ji+1,jj ,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj) & 717 727 & - ( pu(ji ,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj) 718 728 END_2D 719 729 CASE ( np_CRV ) !* Coriolis + relative vorticity 720 DO_2D _10_10730 DO_2D( 1, 0, 1, 0 ) 721 731 zwz(ji,jj,jk) = ( ff_f(ji,jj) + ( e2v(ji+1,jj ) * pv(ji+1,jj ,jk) - e2v(ji,jj) * pv(ji,jj,jk) & 722 732 & - e1u(ji ,jj+1) * pu(ji ,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk) ) & … … 724 734 END_2D 725 735 CASE ( np_CME ) !* Coriolis + metric 726 DO_2D _10_10736 DO_2D( 1, 0, 1, 0 ) 727 737 zwz(ji,jj,jk) = ff_f(ji,jj) + ( pv(ji+1,jj ,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj) & 728 738 & - ( pu(ji ,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj) … … 733 743 ! 734 744 IF( ln_dynvor_msk ) THEN !== mask/unmask vorticity ==! 735 DO_2D _10_10745 DO_2D( 1, 0, 1, 0 ) 736 746 zwz(ji,jj,jk) = zwz(ji,jj,jk) * fmask(ji,jj,jk) 737 747 END_2D … … 739 749 END DO 740 750 ! 741 CALL lbc_lnk( 'dynvor', zwz, 'F', 1. )751 CALL lbc_lnk( 'dynvor', zwz, 'F', 1.0_wp ) 742 752 ! 743 753 DO jk = 1, jpkm1 ! Horizontal slab … … 766 776 END DO 767 777 END DO 768 DO_2D _00_00778 DO_2D( 0, 0, 0, 0 ) 769 779 zua = + r1_12 * r1_e1u(ji,jj) * ( ztne(ji,jj ) * zwy(ji ,jj ) + ztnw(ji+1,jj) * zwy(ji+1,jj ) & 770 780 & + ztse(ji,jj ) * zwy(ji ,jj-1) + ztsw(ji+1,jj) * zwy(ji+1,jj-1) ) … … 826 836 IF(lwp) WRITE(numout,*) ' change fmask value in the angles (T) ln_vorlat = ', ln_vorlat 827 837 IF( ln_vorlat .AND. ( ln_dynvor_ene .OR. ln_dynvor_ens .OR. ln_dynvor_mix ) ) THEN 828 DO_3D _10_10(1, jpk )838 DO_3D( 1, 0, 1, 0, 1, jpk ) 829 839 IF( tmask(ji,jj+1,jk) + tmask(ji+1,jj+1,jk) & 830 & + tmask(ji,jj ,jk) + tmask(ji+1,jj +1,jk) == 3._wp ) fmask(ji,jj,jk) = 1._wp840 & + tmask(ji,jj ,jk) + tmask(ji+1,jj ,jk) == 3._wp ) fmask(ji,jj,jk) = 1._wp 831 841 END_3D 832 842 ! … … 865 875 CASE( np_ENT ) !* T-point metric term : pre-compute di(e2u)/2 and dj(e1v)/2 866 876 ALLOCATE( di_e2u_2(jpi,jpj), dj_e1v_2(jpi,jpj) ) 867 DO_2D _00_00877 DO_2D( 0, 0, 0, 0 ) 868 878 di_e2u_2(ji,jj) = ( e2u(ji,jj) - e2u(ji-1,jj ) ) * 0.5_wp 869 879 dj_e1v_2(ji,jj) = ( e1v(ji,jj) - e1v(ji ,jj-1) ) * 0.5_wp 870 880 END_2D 871 CALL lbc_lnk_multi( 'dynvor', di_e2u_2, 'T', -1. , dj_e1v_2, 'T', -1.) ! Lateral boundary conditions881 CALL lbc_lnk_multi( 'dynvor', di_e2u_2, 'T', -1.0_wp , dj_e1v_2, 'T', -1.0_wp ) ! Lateral boundary conditions 872 882 ! 873 883 CASE DEFAULT !* F-point metric term : pre-compute di(e2u)/(2*e1e2f) and dj(e1v)/(2*e1e2f) 874 884 ALLOCATE( di_e2v_2e1e2f(jpi,jpj), dj_e1u_2e1e2f(jpi,jpj) ) 875 DO_2D _10_10885 DO_2D( 1, 0, 1, 0 ) 876 886 di_e2v_2e1e2f(ji,jj) = ( e2v(ji+1,jj ) - e2v(ji,jj) ) * 0.5 * r1_e1e2f(ji,jj) 877 887 dj_e1u_2e1e2f(ji,jj) = ( e1u(ji ,jj+1) - e1u(ji,jj) ) * 0.5 * r1_e1e2f(ji,jj) 878 888 END_2D 879 CALL lbc_lnk_multi( 'dynvor', di_e2v_2e1e2f, 'F', -1. , dj_e1u_2e1e2f, 'F', -1.) ! Lateral boundary conditions889 CALL lbc_lnk_multi( 'dynvor', di_e2v_2e1e2f, 'F', -1.0_wp , dj_e1u_2e1e2f, 'F', -1.0_wp ) ! Lateral boundary conditions 880 890 END SELECT 881 891 ! -
NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/OCE/DYN/dynzad.F90
r12991 r13710 30 30 !! * Substitutions 31 31 # include "do_loop_substitute.h90" 32 # include "domzgr_substitute.h90" 32 33 !!---------------------------------------------------------------------- 33 34 !! NEMO/OCE 4.0 , NEMO Consortium (2018) … … 71 72 ENDIF 72 73 73 IF( l_trddyn ) THEN ! Save puu(:,:,:,Krhs) and pvv(:,:,:,Krhs) trends74 IF( l_trddyn ) THEN ! Save puu(:,:,:,Krhs) and pvv(:,:,:,Krhs) trends 74 75 ALLOCATE( ztrdu(jpi,jpj,jpk) , ztrdv(jpi,jpj,jpk) ) 75 76 ztrdu(:,:,:) = puu(:,:,:,Krhs) … … 77 78 ENDIF 78 79 79 DO jk = 2, jpkm1 ! Vertical momentum advection at level w and u- and v- vertical80 DO_2D _01_0180 DO jk = 2, jpkm1 ! Vertical momentum advection at level w and u- and v- vertical 81 DO_2D( 0, 1, 0, 1 ) ! vertical fluxes 81 82 IF( ln_vortex_force ) THEN 82 83 zww(ji,jj) = 0.25_wp * e1e2t(ji,jj) * ( ww(ji,jj,jk) + wsd(ji,jj,jk) ) … … 85 86 ENDIF 86 87 END_2D 87 DO_2D _00_0088 DO_2D( 0, 0, 0, 0 ) ! vertical momentum advection at w-point 88 89 zwuw(ji,jj,jk) = ( zww(ji+1,jj ) + zww(ji,jj) ) * ( puu(ji,jj,jk-1,Kmm) - puu(ji,jj,jk,Kmm) ) 89 90 zwvw(ji,jj,jk) = ( zww(ji ,jj+1) + zww(ji,jj) ) * ( pvv(ji,jj,jk-1,Kmm) - pvv(ji,jj,jk,Kmm) ) … … 92 93 ! 93 94 ! Surface and bottom advective fluxes set to zero 94 DO_2D _00_0095 DO_2D( 0, 0, 0, 0 ) 95 96 zwuw(ji,jj, 1 ) = 0._wp 96 97 zwvw(ji,jj, 1 ) = 0._wp … … 99 100 END_2D 100 101 ! 101 DO_3D_00_00( 1, jpkm1 ) 102 puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) - ( zwuw(ji,jj,jk) + zwuw(ji,jj,jk+1) ) * r1_e1e2u(ji,jj) / e3u(ji,jj,jk,Kmm) 103 pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) - ( zwvw(ji,jj,jk) + zwvw(ji,jj,jk+1) ) * r1_e1e2v(ji,jj) / e3v(ji,jj,jk,Kmm) 102 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) ! Vertical momentum advection at u- and v-points 103 puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) - ( zwuw(ji,jj,jk) + zwuw(ji,jj,jk+1) ) * r1_e1e2u(ji,jj) & 104 & / e3u(ji,jj,jk,Kmm) 105 pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) - ( zwvw(ji,jj,jk) + zwvw(ji,jj,jk+1) ) * r1_e1e2v(ji,jj) & 106 & / e3v(ji,jj,jk,Kmm) 104 107 END_3D 105 108 106 IF( l_trddyn ) THEN ! save the vertical advection trends for diagnostic109 IF( l_trddyn ) THEN ! save the vertical advection trends for diagnostic 107 110 ztrdu(:,:,:) = puu(:,:,:,Krhs) - ztrdu(:,:,:) 108 111 ztrdv(:,:,:) = pvv(:,:,:,Krhs) - ztrdv(:,:,:) … … 110 113 DEALLOCATE( ztrdu, ztrdv ) 111 114 ENDIF 112 ! ! Control print115 ! ! Control print 113 116 IF(sn_cfctl%l_prtctl) CALL prt_ctl( tab3d_1=puu(:,:,:,Krhs), clinfo1=' zad - Ua: ', mask1=umask, & 114 117 & tab3d_2=pvv(:,:,:,Krhs), clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' ) -
NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/OCE/DYN/dynzdf.F90
r12489 r13710 38 38 !! * Substitutions 39 39 # include "do_loop_substitute.h90" 40 # include "domzgr_substitute.h90" 40 41 !!---------------------------------------------------------------------- 41 42 !! NEMO/OCE 4.0 , NEMO Consortium (2018) … … 55 56 !! ** Method : - Leap-Frog time stepping on all trends but the vertical mixing 56 57 !! u(after) = u(before) + 2*dt * u(rhs) vector form or linear free surf. 57 !! u(after) = ( e3u_b*u(before) + 2*dt * e3u_n*u(rhs) ) / e3u (after)otherwise58 !! u(after) = ( e3u_b*u(before) + 2*dt * e3u_n*u(rhs) ) / e3u_after otherwise 58 59 !! - update the after velocity with the implicit vertical mixing. 59 60 !! This requires to solver the following system: 60 !! u(after) = u(after) + 1/e3u (after) dk+1[ mi(avm) / e3uw(after)dk[ua] ]61 !! u(after) = u(after) + 1/e3u_after dk+1[ mi(avm) / e3uw_after dk[ua] ] 61 62 !! with the following surface/top/bottom boundary condition: 62 63 !! surface: wind stress input (averaged over kt-1/2 & kt+1/2) … … 106 107 ! ! time stepping except vertical diffusion 107 108 IF( ln_dynadv_vec .OR. ln_linssh ) THEN ! applied on velocity 108 DO jk = 1, jpkm1109 puu( :,:,jk,Kaa) = ( puu(:,:,jk,Kbb) + rDt * puu(:,:,jk,Krhs) ) * umask(:,:,jk)110 pvv( :,:,jk,Kaa) = ( pvv(:,:,jk,Kbb) + rDt * pvv(:,:,jk,Krhs) ) * vmask(:,:,jk)111 END DO109 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 110 puu(ji,jj,jk,Kaa) = ( puu(ji,jj,jk,Kbb) + rDt * puu(ji,jj,jk,Krhs) ) * umask(ji,jj,jk) 111 pvv(ji,jj,jk,Kaa) = ( pvv(ji,jj,jk,Kbb) + rDt * pvv(ji,jj,jk,Krhs) ) * vmask(ji,jj,jk) 112 END_3D 112 113 ELSE ! applied on thickness weighted velocity 113 DO jk = 1, jpkm1 114 puu(:,:,jk,Kaa) = ( e3u(:,:,jk,Kbb) * puu(:,:,jk,Kbb) & 115 & + rDt * e3u(:,:,jk,Kmm) * puu(:,:,jk,Krhs) ) / e3u(:,:,jk,Kaa) * umask(:,:,jk) 116 pvv(:,:,jk,Kaa) = ( e3v(:,:,jk,Kbb) * pvv(:,:,jk,Kbb) & 117 & + rDt * e3v(:,:,jk,Kmm) * pvv(:,:,jk,Krhs) ) / e3v(:,:,jk,Kaa) * vmask(:,:,jk) 118 END DO 114 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 115 puu(ji,jj,jk,Kaa) = ( e3u(ji,jj,jk,Kbb) * puu(ji,jj,jk,Kbb ) & 116 & + rDt * e3u(ji,jj,jk,Kmm) * puu(ji,jj,jk,Krhs) ) & 117 & / e3u(ji,jj,jk,Kaa) * umask(ji,jj,jk) 118 pvv(ji,jj,jk,Kaa) = ( e3v(ji,jj,jk,Kbb) * pvv(ji,jj,jk,Kbb ) & 119 & + rDt * e3v(ji,jj,jk,Kmm) * pvv(ji,jj,jk,Krhs) ) & 120 & / e3v(ji,jj,jk,Kaa) * vmask(ji,jj,jk) 121 END_3D 119 122 ENDIF 120 123 ! ! add top/bottom friction … … 124 127 ! G. Madec : in linear free surface, e3u(:,:,:,Kaa) = e3u(:,:,:,Kmm) = e3u_0, so systematic use of e3u(:,:,:,Kaa) 125 128 IF( ln_drgimp .AND. ln_dynspg_ts ) THEN 126 DO jk = 1, jpkm1! remove barotropic velocities127 puu( :,:,jk,Kaa) = ( puu(:,:,jk,Kaa) - uu_b(:,:,Kaa) ) * umask(:,:,jk)128 pvv( :,:,jk,Kaa) = ( pvv(:,:,jk,Kaa) - vv_b(:,:,Kaa) ) * vmask(:,:,jk)129 END DO130 DO_2D _00_00129 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) ! remove barotropic velocities 130 puu(ji,jj,jk,Kaa) = ( puu(ji,jj,jk,Kaa) - uu_b(ji,jj,Kaa) ) * umask(ji,jj,jk) 131 pvv(ji,jj,jk,Kaa) = ( pvv(ji,jj,jk,Kaa) - vv_b(ji,jj,Kaa) ) * vmask(ji,jj,jk) 132 END_3D 133 DO_2D( 0, 0, 0, 0 ) ! Add bottom/top stress due to barotropic component only 131 134 iku = mbku(ji,jj) ! ocean bottom level at u- and v-points 132 135 ikv = mbkv(ji,jj) ! (deepest ocean u- and v-points) 133 ze3ua = ( 1._wp - r_vvl ) * e3u(ji,jj,iku,Kmm) + r_vvl * e3u(ji,jj,iku,Kaa) 134 ze3va = ( 1._wp - r_vvl ) * e3v(ji,jj,ikv,Kmm) + r_vvl * e3v(ji,jj,ikv,Kaa) 136 ze3ua = ( 1._wp - r_vvl ) * e3u(ji,jj,iku,Kmm) & 137 & + r_vvl * e3u(ji,jj,iku,Kaa) 138 ze3va = ( 1._wp - r_vvl ) * e3v(ji,jj,ikv,Kmm) & 139 & + r_vvl * e3v(ji,jj,ikv,Kaa) 135 140 puu(ji,jj,iku,Kaa) = puu(ji,jj,iku,Kaa) + rDt * 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) * uu_b(ji,jj,Kaa) / ze3ua 136 141 pvv(ji,jj,ikv,Kaa) = pvv(ji,jj,ikv,Kaa) + rDt * 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) * vv_b(ji,jj,Kaa) / ze3va 137 142 END_2D 138 IF( ln_isfcav ) THEN ! Ocean cavities (ISF)139 DO_2D _00_00143 IF( ln_isfcav.OR.ln_drgice_imp ) THEN ! Ocean cavities (ISF) 144 DO_2D( 0, 0, 0, 0 ) 140 145 iku = miku(ji,jj) ! top ocean level at u- and v-points 141 146 ikv = mikv(ji,jj) ! (first wet ocean u- and v-points) 142 ze3ua = ( 1._wp - r_vvl ) * e3u(ji,jj,iku,Kmm) + r_vvl * e3u(ji,jj,iku,Kaa) 143 ze3va = ( 1._wp - r_vvl ) * e3v(ji,jj,ikv,Kmm) + r_vvl * e3v(ji,jj,ikv,Kaa) 147 ze3ua = ( 1._wp - r_vvl ) * e3u(ji,jj,iku,Kmm) & 148 & + r_vvl * e3u(ji,jj,iku,Kaa) 149 ze3va = ( 1._wp - r_vvl ) * e3v(ji,jj,ikv,Kmm) & 150 & + r_vvl * e3v(ji,jj,ikv,Kaa) 144 151 puu(ji,jj,iku,Kaa) = puu(ji,jj,iku,Kaa) + rDt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * uu_b(ji,jj,Kaa) / ze3ua 145 152 pvv(ji,jj,ikv,Kaa) = pvv(ji,jj,ikv,Kaa) + rDt * 0.5*( rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) * vv_b(ji,jj,Kaa) / ze3va … … 155 162 SELECT CASE( nldf_dyn ) 156 163 CASE( np_lap_i ) ! rotated lateral mixing: add its vertical mixing (akzu) 157 DO_3D_00_00( 1, jpkm1 ) 158 ze3ua = ( 1._wp - r_vvl ) * e3u(ji,jj,jk,Kmm) + r_vvl * e3u(ji,jj,jk,Kaa) ! after scale factor at U-point 164 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 165 ze3ua = ( 1._wp - r_vvl ) * e3u(ji,jj,jk,Kmm) & 166 & + r_vvl * e3u(ji,jj,jk,Kaa) ! after scale factor at U-point 159 167 zzwi = - zdt * ( avm(ji+1,jj,jk ) + avm(ji,jj,jk ) + akzu(ji,jj,jk ) ) & 160 168 & / ( ze3ua * e3uw(ji,jj,jk ,Kmm) ) * wumask(ji,jj,jk ) … … 168 176 END_3D 169 177 CASE DEFAULT ! iso-level lateral mixing 170 DO_3D_00_00( 1, jpkm1 ) 171 ze3ua = ( 1._wp - r_vvl ) * e3u(ji,jj,jk,Kmm) + r_vvl * e3u(ji,jj,jk,Kaa) ! after scale factor at U-point 172 zzwi = - zdt * ( avm(ji+1,jj,jk ) + avm(ji,jj,jk ) ) / ( ze3ua * e3uw(ji,jj,jk ,Kmm) ) * wumask(ji,jj,jk ) 173 zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) ) / ( ze3ua * e3uw(ji,jj,jk+1,Kmm) ) * wumask(ji,jj,jk+1) 178 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 179 ze3ua = ( 1._wp - r_vvl ) * e3u(ji,jj,jk,Kmm) & ! after scale factor at U-point 180 & + r_vvl * e3u(ji,jj,jk,Kaa) 181 zzwi = - zdt * ( avm(ji+1,jj,jk ) + avm(ji,jj,jk ) ) & 182 & / ( ze3ua * e3uw(ji,jj,jk ,Kmm) ) * wumask(ji,jj,jk ) 183 zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) ) & 184 & / ( ze3ua * e3uw(ji,jj,jk+1,Kmm) ) * wumask(ji,jj,jk+1) 174 185 zWui = ( wi(ji,jj,jk ) + wi(ji+1,jj,jk ) ) / ze3ua 175 186 zWus = ( wi(ji,jj,jk+1) + wi(ji+1,jj,jk+1) ) / ze3ua … … 179 190 END_3D 180 191 END SELECT 181 DO_2D _00_00192 DO_2D( 0, 0, 0, 0 ) !* Surface boundary conditions 182 193 zwi(ji,jj,1) = 0._wp 183 ze3ua = ( 1._wp - r_vvl ) * e3u(ji,jj,1,Kmm) + r_vvl * e3u(ji,jj,1,Kaa) 184 zzws = - zdt * ( avm(ji+1,jj,2) + avm(ji ,jj,2) ) / ( ze3ua * e3uw(ji,jj,2,Kmm) ) * wumask(ji,jj,2) 194 ze3ua = ( 1._wp - r_vvl ) * e3u(ji,jj,1,Kmm) & 195 & + r_vvl * e3u(ji,jj,1,Kaa) 196 zzws = - zdt * ( avm(ji+1,jj,2) + avm(ji ,jj,2) ) & 197 & / ( ze3ua * e3uw(ji,jj,2,Kmm) ) * wumask(ji,jj,2) 185 198 zWus = ( wi(ji ,jj,2) + wi(ji+1,jj,2) ) / ze3ua 186 199 zws(ji,jj,1 ) = zzws - zdt * MAX( zWus, 0._wp ) … … 190 203 SELECT CASE( nldf_dyn ) 191 204 CASE( np_lap_i ) ! rotated lateral mixing: add its vertical mixing (akzu) 192 DO_3D_00_00( 1, jpkm1 ) 193 ze3ua = ( 1._wp - r_vvl ) * e3u(ji,jj,jk,Kmm) + r_vvl * e3u(ji,jj,jk,Kaa) ! after scale factor at U-point 205 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 206 ze3ua = ( 1._wp - r_vvl ) * e3u(ji,jj,jk,Kmm) & 207 & + r_vvl * e3u(ji,jj,jk,Kaa) ! after scale factor at U-point 194 208 zzwi = - zdt * ( avm(ji+1,jj,jk ) + avm(ji,jj,jk ) + akzu(ji,jj,jk ) ) & 195 209 & / ( ze3ua * e3uw(ji,jj,jk ,Kmm) ) * wumask(ji,jj,jk ) … … 201 215 END_3D 202 216 CASE DEFAULT ! iso-level lateral mixing 203 DO_3D_00_00( 1, jpkm1 ) 204 ze3ua = ( 1._wp - r_vvl ) * e3u(ji,jj,jk,Kmm) + r_vvl * e3u(ji,jj,jk,Kaa) ! after scale factor at U-point 205 zzwi = - zdt * ( avm(ji+1,jj,jk ) + avm(ji,jj,jk ) ) / ( ze3ua * e3uw(ji,jj,jk ,Kmm) ) * wumask(ji,jj,jk ) 206 zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) ) / ( ze3ua * e3uw(ji,jj,jk+1,Kmm) ) * wumask(ji,jj,jk+1) 217 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 218 ze3ua = ( 1._wp - r_vvl ) * e3u(ji,jj,jk,Kmm) & 219 & + r_vvl * e3u(ji,jj,jk,Kaa) ! after scale factor at U-point 220 zzwi = - zdt * ( avm(ji+1,jj,jk ) + avm(ji,jj,jk ) ) & 221 & / ( ze3ua * e3uw(ji,jj,jk ,Kmm) ) * wumask(ji,jj,jk ) 222 zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) ) & 223 & / ( ze3ua * e3uw(ji,jj,jk+1,Kmm) ) * wumask(ji,jj,jk+1) 207 224 zwi(ji,jj,jk) = zzwi 208 225 zws(ji,jj,jk) = zzws … … 210 227 END_3D 211 228 END SELECT 212 DO_2D _00_00229 DO_2D( 0, 0, 0, 0 ) !* Surface boundary conditions 213 230 zwi(ji,jj,1) = 0._wp 214 231 zwd(ji,jj,1) = 1._wp - zws(ji,jj,1) … … 224 241 ! 225 242 IF ( ln_drgimp ) THEN ! implicit bottom friction 226 DO_2D _00_00243 DO_2D( 0, 0, 0, 0 ) 227 244 iku = mbku(ji,jj) ! ocean bottom level at u- and v-points 228 ze3ua = ( 1._wp - r_vvl ) * e3u(ji,jj,iku,Kmm) + r_vvl * e3u(ji,jj,iku,Kaa) ! after scale factor at T-point 245 ze3ua = ( 1._wp - r_vvl ) * e3u(ji,jj,iku,Kmm) & 246 & + r_vvl * e3u(ji,jj,iku,Kaa) ! after scale factor at T-point 229 247 zwd(ji,jj,iku) = zwd(ji,jj,iku) - rDt * 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) / ze3ua 230 248 END_2D 231 IF ( ln_isfcav ) THEN ! top friction (always implicit)232 DO_2D _00_00249 IF ( ln_isfcav.OR.ln_drgice_imp ) THEN ! top friction (always implicit) 250 DO_2D( 0, 0, 0, 0 ) 233 251 !!gm top Cd is masked (=0 outside cavities) no need of test on mik>=2 ==>> it has been suppressed 234 252 iku = miku(ji,jj) ! ocean top level at u- and v-points 235 ze3ua = ( 1._wp - r_vvl ) * e3u(ji,jj,iku,Kmm) + r_vvl * e3u(ji,jj,iku,Kaa) ! after scale factor at T-point 253 ze3ua = ( 1._wp - r_vvl ) * e3u(ji,jj,iku,Kmm) & 254 & + r_vvl * e3u(ji,jj,iku,Kaa) ! after scale factor at T-point 236 255 zwd(ji,jj,iku) = zwd(ji,jj,iku) - rDt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) / ze3ua 237 256 END_2D … … 254 273 !----------------------------------------------------------------------- 255 274 ! 256 DO_3D _00_00( 2, jpkm1 )275 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) !== First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 (increasing k) == 257 276 zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1) 258 277 END_3D 259 278 ! 260 DO_2D_00_00 261 ze3ua = ( 1._wp - r_vvl ) * e3u(ji,jj,1,Kmm) + r_vvl * e3u(ji,jj,1,Kaa) 279 DO_2D( 0, 0, 0, 0 ) !== second recurrence: SOLk = RHSk - Lk / Dk-1 Lk-1 ==! 280 ze3ua = ( 1._wp - r_vvl ) * e3u(ji,jj,1,Kmm) & 281 & + r_vvl * e3u(ji,jj,1,Kaa) 262 282 puu(ji,jj,1,Kaa) = puu(ji,jj,1,Kaa) + rDt * 0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) ) & 263 283 & / ( ze3ua * rho0 ) * umask(ji,jj,1) 264 284 END_2D 265 DO_3D _00_00(2, jpkm1 )285 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 266 286 puu(ji,jj,jk,Kaa) = puu(ji,jj,jk,Kaa) - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * puu(ji,jj,jk-1,Kaa) 267 287 END_3D 268 288 ! 269 DO_2D _00_00289 DO_2D( 0, 0, 0, 0 ) !== thrid recurrence : SOLk = ( Lk - Uk * Ek+1 ) / Dk ==! 270 290 puu(ji,jj,jpkm1,Kaa) = puu(ji,jj,jpkm1,Kaa) / zwd(ji,jj,jpkm1) 271 291 END_2D 272 DO_3DS _00_00(jpk-2, 1, -1 )292 DO_3DS( 0, 0, 0, 0, jpk-2, 1, -1 ) 273 293 puu(ji,jj,jk,Kaa) = ( puu(ji,jj,jk,Kaa) - zws(ji,jj,jk) * puu(ji,jj,jk+1,Kaa) ) / zwd(ji,jj,jk) 274 294 END_3D … … 281 301 SELECT CASE( nldf_dyn ) 282 302 CASE( np_lap_i ) ! rotated lateral mixing: add its vertical mixing (akzv) 283 DO_3D_00_00( 1, jpkm1 ) 284 ze3va = ( 1._wp - r_vvl ) * e3v(ji,jj,jk,Kmm) + r_vvl * e3v(ji,jj,jk,Kaa) ! after scale factor at V-point 303 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 304 ze3va = ( 1._wp - r_vvl ) * e3v(ji,jj,jk,Kmm) & 305 & + r_vvl * e3v(ji,jj,jk,Kaa) ! after scale factor at V-point 285 306 zzwi = - zdt * ( avm(ji,jj+1,jk ) + avm(ji,jj,jk ) + akzv(ji,jj,jk ) ) & 286 307 & / ( ze3va * e3vw(ji,jj,jk ,Kmm) ) * wvmask(ji,jj,jk ) … … 294 315 END_3D 295 316 CASE DEFAULT ! iso-level lateral mixing 296 DO_3D_00_00( 1, jpkm1 ) 297 ze3va = ( 1._wp - r_vvl ) * e3v(ji,jj,jk,Kmm) + r_vvl * e3v(ji,jj,jk,Kaa) ! after scale factor at V-point 298 zzwi = - zdt * ( avm(ji,jj+1,jk ) + avm(ji,jj,jk ) ) / ( ze3va * e3vw(ji,jj,jk ,Kmm) ) * wvmask(ji,jj,jk ) 299 zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) ) / ( ze3va * e3vw(ji,jj,jk+1,Kmm) ) * wvmask(ji,jj,jk+1) 317 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 318 ze3va = ( 1._wp - r_vvl ) * e3v(ji,jj,jk,Kmm) & 319 & + r_vvl * e3v(ji,jj,jk,Kaa) ! after scale factor at V-point 320 zzwi = - zdt * ( avm(ji,jj+1,jk ) + avm(ji,jj,jk ) ) & 321 & / ( ze3va * e3vw(ji,jj,jk ,Kmm) ) * wvmask(ji,jj,jk ) 322 zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) ) & 323 & / ( ze3va * e3vw(ji,jj,jk+1,Kmm) ) * wvmask(ji,jj,jk+1) 300 324 zWvi = ( wi(ji,jj,jk ) + wi(ji,jj+1,jk ) ) / ze3va 301 325 zWvs = ( wi(ji,jj,jk+1) + wi(ji,jj+1,jk+1) ) / ze3va … … 305 329 END_3D 306 330 END SELECT 307 DO_2D _00_00331 DO_2D( 0, 0, 0, 0 ) !* Surface boundary conditions 308 332 zwi(ji,jj,1) = 0._wp 309 ze3va = ( 1._wp - r_vvl ) * e3v(ji,jj,1,Kmm) + r_vvl * e3v(ji,jj,1,Kaa) 310 zzws = - zdt * ( avm(ji,jj+1,2) + avm(ji,jj,2) ) / ( ze3va * e3vw(ji,jj,2,Kmm) ) * wvmask(ji,jj,2) 333 ze3va = ( 1._wp - r_vvl ) * e3v(ji,jj,1,Kmm) & 334 & + r_vvl * e3v(ji,jj,1,Kaa) 335 zzws = - zdt * ( avm(ji,jj+1,2) + avm(ji,jj,2) ) & 336 & / ( ze3va * e3vw(ji,jj,2,Kmm) ) * wvmask(ji,jj,2) 311 337 zWvs = ( wi(ji,jj ,2) + wi(ji,jj+1,2) ) / ze3va 312 338 zws(ji,jj,1 ) = zzws - zdt * MAX( zWvs, 0._wp ) … … 316 342 SELECT CASE( nldf_dyn ) 317 343 CASE( np_lap_i ) ! rotated lateral mixing: add its vertical mixing (akzu) 318 DO_3D_00_00( 1, jpkm1 ) 319 ze3va = ( 1._wp - r_vvl ) * e3v(ji,jj,jk,Kmm) + r_vvl * e3v(ji,jj,jk,Kaa) ! after scale factor at V-point 344 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 345 ze3va = ( 1._wp - r_vvl ) * e3v(ji,jj,jk,Kmm) & 346 & + r_vvl * e3v(ji,jj,jk,Kaa) ! after scale factor at V-point 320 347 zzwi = - zdt * ( avm(ji,jj+1,jk ) + avm(ji,jj,jk ) + akzv(ji,jj,jk ) ) & 321 348 & / ( ze3va * e3vw(ji,jj,jk ,Kmm) ) * wvmask(ji,jj,jk ) … … 327 354 END_3D 328 355 CASE DEFAULT ! iso-level lateral mixing 329 DO_3D_00_00( 1, jpkm1 ) 330 ze3va = ( 1._wp - r_vvl ) * e3v(ji,jj,jk,Kmm) + r_vvl * e3v(ji,jj,jk,Kaa) ! after scale factor at V-point 331 zzwi = - zdt * ( avm(ji,jj+1,jk ) + avm(ji,jj,jk ) ) / ( ze3va * e3vw(ji,jj,jk ,Kmm) ) * wvmask(ji,jj,jk ) 332 zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) ) / ( ze3va * e3vw(ji,jj,jk+1,Kmm) ) * wvmask(ji,jj,jk+1) 356 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 357 ze3va = ( 1._wp - r_vvl ) * e3v(ji,jj,jk,Kmm) & 358 & + r_vvl * e3v(ji,jj,jk,Kaa) ! after scale factor at V-point 359 zzwi = - zdt * ( avm(ji,jj+1,jk ) + avm(ji,jj,jk ) ) & 360 & / ( ze3va * e3vw(ji,jj,jk ,Kmm) ) * wvmask(ji,jj,jk ) 361 zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) ) & 362 & / ( ze3va * e3vw(ji,jj,jk+1,Kmm) ) * wvmask(ji,jj,jk+1) 333 363 zwi(ji,jj,jk) = zzwi 334 364 zws(ji,jj,jk) = zzws … … 336 366 END_3D 337 367 END SELECT 338 DO_2D _00_00368 DO_2D( 0, 0, 0, 0 ) !* Surface boundary conditions 339 369 zwi(ji,jj,1) = 0._wp 340 370 zwd(ji,jj,1) = 1._wp - zws(ji,jj,1) … … 349 379 ! 350 380 IF( ln_drgimp ) THEN 351 DO_2D _00_00381 DO_2D( 0, 0, 0, 0 ) 352 382 ikv = mbkv(ji,jj) ! (deepest ocean u- and v-points) 353 ze3va = ( 1._wp - r_vvl ) * e3v(ji,jj,ikv,Kmm) + r_vvl * e3v(ji,jj,ikv,Kaa) ! after scale factor at T-point 383 ze3va = ( 1._wp - r_vvl ) * e3v(ji,jj,ikv,Kmm) & 384 & + r_vvl * e3v(ji,jj,ikv,Kaa) ! after scale factor at T-point 354 385 zwd(ji,jj,ikv) = zwd(ji,jj,ikv) - rDt * 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) / ze3va 355 386 END_2D 356 IF ( ln_isfcav ) THEN357 DO_2D _00_00387 IF ( ln_isfcav.OR.ln_drgice_imp ) THEN 388 DO_2D( 0, 0, 0, 0 ) 358 389 ikv = mikv(ji,jj) ! (first wet ocean u- and v-points) 359 ze3va = ( 1._wp - r_vvl ) * e3v(ji,jj,ikv,Kmm) + r_vvl * e3v(ji,jj,ikv,Kaa) ! after scale factor at T-point 390 ze3va = ( 1._wp - r_vvl ) * e3v(ji,jj,ikv,Kmm) & 391 & + r_vvl * e3v(ji,jj,ikv,Kaa) ! after scale factor at T-point 360 392 zwd(ji,jj,ikv) = zwd(ji,jj,ikv) - rDt * 0.5*( rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) / ze3va 361 393 END_2D … … 378 410 !----------------------------------------------------------------------- 379 411 ! 380 DO_3D _00_00( 2, jpkm1 )412 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) !== First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 (increasing k) == 381 413 zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1) 382 414 END_3D 383 415 ! 384 DO_2D_00_00 385 ze3va = ( 1._wp - r_vvl ) * e3v(ji,jj,1,Kmm) + r_vvl * e3v(ji,jj,1,Kaa) 416 DO_2D( 0, 0, 0, 0 ) !== second recurrence: SOLk = RHSk - Lk / Dk-1 Lk-1 ==! 417 ze3va = ( 1._wp - r_vvl ) * e3v(ji,jj,1,Kmm) & 418 & + r_vvl * e3v(ji,jj,1,Kaa) 386 419 pvv(ji,jj,1,Kaa) = pvv(ji,jj,1,Kaa) + rDt * 0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) ) & 387 420 & / ( ze3va * rho0 ) * vmask(ji,jj,1) 388 421 END_2D 389 DO_3D _00_00(2, jpkm1 )422 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 390 423 pvv(ji,jj,jk,Kaa) = pvv(ji,jj,jk,Kaa) - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * pvv(ji,jj,jk-1,Kaa) 391 424 END_3D 392 425 ! 393 DO_2D _00_00426 DO_2D( 0, 0, 0, 0 ) !== third recurrence : SOLk = ( Lk - Uk * SOLk+1 ) / Dk ==! 394 427 pvv(ji,jj,jpkm1,Kaa) = pvv(ji,jj,jpkm1,Kaa) / zwd(ji,jj,jpkm1) 395 428 END_2D 396 DO_3DS _00_00(jpk-2, 1, -1 )429 DO_3DS( 0, 0, 0, 0, jpk-2, 1, -1 ) 397 430 pvv(ji,jj,jk,Kaa) = ( pvv(ji,jj,jk,Kaa) - zws(ji,jj,jk) * pvv(ji,jj,jk+1,Kaa) ) / zwd(ji,jj,jk) 398 431 END_3D -
NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/OCE/DYN/sshwzv.F90
r12489 r13710 28 28 USE bdydyn2d ! bdy_ssh routine 29 29 #if defined key_agrif 30 USE agrif_oce 30 31 USE agrif_oce_interp 31 32 #endif … … 50 51 !! * Substitutions 51 52 # include "do_loop_substitute.h90" 53 # include "domzgr_substitute.h90" 54 52 55 !!---------------------------------------------------------------------- 53 56 !! NEMO/OCE 4.0 , NEMO Consortium (2018) … … 110 113 ! 111 114 #if defined key_agrif 112 Kbb_a = Kbb; Kmm_a = Kmm; Krhs_a = Kaa; CALL agrif_ssh( kt ) 115 Kbb_a = Kbb ; Kmm_a = Kmm ; Krhs_a = Kaa 116 CALL agrif_ssh( kt ) 113 117 #endif 114 118 ! 115 119 IF ( .NOT.ln_dynspg_ts ) THEN 116 120 IF( ln_bdy ) THEN 117 CALL lbc_lnk( 'sshwzv', pssh(:,:,Kaa), 'T', 1. ) ! Not sure that's necessary121 CALL lbc_lnk( 'sshwzv', pssh(:,:,Kaa), 'T', 1.0_wp ) ! Not sure that's necessary 118 122 CALL bdy_ssh( pssh(:,:,Kaa) ) ! Duplicate sea level across open boundaries 119 123 ENDIF … … 130 134 131 135 132 SUBROUTINE wzv( kt, Kbb, Kmm, pww, Kaa)136 SUBROUTINE wzv( kt, Kbb, Kmm, Kaa, pww ) 133 137 !!---------------------------------------------------------------------- 134 138 !! *** ROUTINE wzv *** … … 147 151 INTEGER , INTENT(in) :: kt ! time step 148 152 INTEGER , INTENT(in) :: Kbb, Kmm, Kaa ! time level indices 149 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pww ! now vertical velocity153 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pww ! vertical velocity at Kmm 150 154 ! 151 155 INTEGER :: ji, jj, jk ! dummy loop indices … … 166 170 ! !------------------------------! 167 171 ! 168 IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN ! z_tilde and layer cases 172 ! !===============================! 173 IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN !== z_tilde and layer cases ==! 174 ! !===============================! 169 175 ALLOCATE( zhdiv(jpi,jpj,jpk) ) 170 176 ! … … 172 178 ! horizontal divergence of thickness diffusion transport ( velocity multiplied by e3t) 173 179 ! - ML - note: computation already done in dom_vvl_sf_nxt. Could be optimized (not critical and clearer this way) 174 DO_2D _00_00180 DO_2D( 0, 0, 0, 0 ) 175 181 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) ) 176 182 END_2D 177 183 END DO 178 CALL lbc_lnk('sshwzv', zhdiv, 'T', 1. ) ! - ML - Perhaps not necessary: not used for horizontal "connexions"184 CALL lbc_lnk('sshwzv', zhdiv, 'T', 1.0_wp) ! - ML - Perhaps not necessary: not used for horizontal "connexions" 179 185 ! ! Is it problematic to have a wrong vertical velocity in boundary cells? 180 186 ! ! Same question holds for hdiv. Perhaps just for security 181 187 DO jk = jpkm1, 1, -1 ! integrate from the bottom the hor. divergence 182 188 ! computation of w 183 pww(:,:,jk) = pww(:,:,jk+1) - ( e3t(:,:,jk,Kmm) * hdiv(:,:,jk) + zhdiv(:,:,jk) & 184 & + r1_Dt * ( e3t(:,:,jk,Kaa) - e3t(:,:,jk,Kbb) ) ) * tmask(:,:,jk) 189 pww(:,:,jk) = pww(:,:,jk+1) - ( e3t(:,:,jk,Kmm) * hdiv(:,:,jk) & 190 & + zhdiv(:,:,jk) & 191 & + r1_Dt * ( e3t(:,:,jk,Kaa) & 192 & - e3t(:,:,jk,Kbb) ) ) * tmask(:,:,jk) 185 193 END DO 186 194 ! IF( ln_vvl_layer ) pww(:,:,:) = 0.e0 187 195 DEALLOCATE( zhdiv ) 188 ELSE ! z_star and linear free surface cases 189 DO jk = jpkm1, 1, -1 ! integrate from the bottom the hor. divergence 190 ! computation of w 196 ! !=================================! 197 ELSEIF( ln_linssh ) THEN !== linear free surface cases ==! 198 ! !=================================! 199 DO jk = jpkm1, 1, -1 ! integrate from the bottom the hor. divergence 200 pww(:,:,jk) = pww(:,:,jk+1) - ( e3t(:,:,jk,Kmm) * hdiv(:,:,jk) ) * tmask(:,:,jk) 201 END DO 202 ! !==========================================! 203 ELSE !== Quasi-Eulerian vertical coordinate ==! ('key_qco') 204 ! !==========================================! 205 DO jk = jpkm1, 1, -1 ! integrate from the bottom the hor. divergence 191 206 pww(:,:,jk) = pww(:,:,jk+1) - ( e3t(:,:,jk,Kmm) * hdiv(:,:,jk) & 192 & + r1_Dt * ( e3t(:,:,jk,Kaa) - e3t(:,:,jk,Kbb) ) ) * tmask(:,:,jk) 207 & + r1_Dt * ( e3t(:,:,jk,Kaa) & 208 & - e3t(:,:,jk,Kbb) ) ) * tmask(:,:,jk) 193 209 END DO 194 210 ENDIF … … 200 216 ENDIF 201 217 ! 202 #if defined key_agrif 203 IF( .NOT. AGRIF_Root() ) THEN 204 IF ((nbondi == 1).OR.(nbondi == 2)) pww(nlci-1 , : ,:) = 0.e0 ! east 205 IF ((nbondi == -1).OR.(nbondi == 2)) pww(2 , : ,:) = 0.e0 ! west 206 IF ((nbondj == 1).OR.(nbondj == 2)) pww(: ,nlcj-1 ,:) = 0.e0 ! north 207 IF ((nbondj == -1).OR.(nbondj == 2)) pww(: ,2 ,:) = 0.e0 ! south 218 #if defined key_agrif 219 IF( .NOT. AGRIF_Root() ) THEN 220 ! 221 ! Mask vertical velocity at first/last columns/row 222 ! inside computational domain (cosmetic) 223 DO jk = 1, jpkm1 224 IF( lk_west ) THEN ! --- West --- ! 225 DO ji = mi0(2+nn_hls), mi1(2+nn_hls) 226 DO jj = 1, jpj 227 pww(ji,jj,jk) = 0._wp 228 END DO 229 END DO 230 ENDIF 231 IF( lk_east ) THEN ! --- East --- ! 232 DO ji = mi0(jpiglo-1-nn_hls), mi1(jpiglo-1-nn_hls) 233 DO jj = 1, jpj 234 pww(ji,jj,jk) = 0._wp 235 END DO 236 END DO 237 ENDIF 238 IF( lk_south ) THEN ! --- South --- ! 239 DO jj = mj0(2+nn_hls), mj1(2+nn_hls) 240 DO ji = 1, jpi 241 pww(ji,jj,jk) = 0._wp 242 END DO 243 END DO 244 ENDIF 245 IF( lk_north ) THEN ! --- North --- ! 246 DO jj = mj0(jpjglo-1-nn_hls), mj1(jpjglo-1-nn_hls) 247 DO ji = 1, jpi 248 pww(ji,jj,jk) = 0._wp 249 END DO 250 END DO 251 ENDIF 252 ! 253 END DO 254 ! 208 255 ENDIF 209 #endif 256 #endif 210 257 ! 211 258 IF( ln_timing ) CALL timing_stop('wzv') … … 214 261 215 262 216 SUBROUTINE ssh_atf( kt, Kbb, Kmm, Kaa, pssh )263 SUBROUTINE ssh_atf( kt, Kbb, Kmm, Kaa, pssh, pssh_f ) 217 264 !!---------------------------------------------------------------------- 218 265 !! *** ROUTINE ssh_atf *** … … 231 278 INTEGER , INTENT(in ) :: kt ! ocean time-step index 232 279 INTEGER , INTENT(in ) :: Kbb, Kmm, Kaa ! ocean time level indices 233 REAL(wp), DIMENSION(jpi,jpj,jpt), INTENT(inout) :: pssh ! SSH field 280 REAL(wp), DIMENSION(jpi,jpj,jpt) , TARGET, INTENT(inout) :: pssh ! SSH field 281 REAL(wp), DIMENSION(jpi,jpj ), OPTIONAL, TARGET, INTENT( out) :: pssh_f ! filtered SSH field 234 282 ! 235 283 REAL(wp) :: zcoef ! local scalar 284 REAL(wp), POINTER, DIMENSION(:,:) :: zssh ! pointer for filtered SSH 236 285 !!---------------------------------------------------------------------- 237 286 ! … … 245 294 ! !== Euler time-stepping: no filter, just swap ==! 246 295 IF ( .NOT.( l_1st_euler ) ) THEN ! Only do time filtering for leapfrog timesteps 296 IF( PRESENT( pssh_f ) ) THEN ; zssh => pssh_f 297 ELSE ; zssh => pssh(:,:,Kmm) 298 ENDIF 247 299 ! ! filtered "now" field 248 300 pssh(:,:,Kmm) = pssh(:,:,Kmm) + rn_atfp * ( pssh(:,:,Kbb) - 2 * pssh(:,:,Kmm) + pssh(:,:,Kaa) ) … … 266 318 END SUBROUTINE ssh_atf 267 319 320 268 321 SUBROUTINE wAimp( kt, Kmm ) 269 322 !!---------------------------------------------------------------------- … … 286 339 ! 287 340 INTEGER :: ji, jj, jk ! dummy loop indices 288 REAL(wp) :: zCu, zcff, z1_e3t 341 REAL(wp) :: zCu, zcff, z1_e3t, zdt ! local scalars 289 342 REAL(wp) , PARAMETER :: Cu_min = 0.15_wp ! local parameters 290 343 REAL(wp) , PARAMETER :: Cu_max = 0.30_wp ! local parameters … … 303 356 ! 304 357 ! Calculate Courant numbers 358 zdt = 2._wp * rn_Dt ! 2*rn_Dt and not rDt (for restartability) 305 359 IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN 306 DO_3D _00_00(1, jpkm1 )360 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 307 361 z1_e3t = 1._wp / e3t(ji,jj,jk,Kmm) 308 ! 2*rn_Dt and not rDt (for restartability) 309 Cu_adv(ji,jj,jk) = 2._wp * rn_Dt * ( ( MAX( ww(ji,jj,jk) , 0._wp ) - MIN( ww(ji,jj,jk+1) , 0._wp ) ) & 310 & + ( MAX( e2u(ji ,jj)*e3u(ji ,jj,jk,Kmm)*uu(ji ,jj,jk,Kmm) + un_td(ji ,jj,jk), 0._wp ) - & 311 & MIN( e2u(ji-1,jj)*e3u(ji-1,jj,jk,Kmm)*uu(ji-1,jj,jk,Kmm) + un_td(ji-1,jj,jk), 0._wp ) ) & 362 Cu_adv(ji,jj,jk) = zdt * & 363 & ( ( MAX( ww(ji,jj,jk) , 0._wp ) - MIN( ww(ji,jj,jk+1) , 0._wp ) ) & 364 & + ( MAX( e2u(ji ,jj) * e3u(ji ,jj,jk,Kmm) & 365 & * uu (ji ,jj,jk,Kmm) + un_td(ji ,jj,jk), 0._wp ) - & 366 & MIN( e2u(ji-1,jj) * e3u(ji-1,jj,jk,Kmm) & 367 & * uu (ji-1,jj,jk,Kmm) + un_td(ji-1,jj,jk), 0._wp ) ) & 312 368 & * r1_e1e2t(ji,jj) & 313 & + ( MAX( e1v(ji,jj )*e3v(ji,jj ,jk,Kmm)*vv(ji,jj ,jk,Kmm) + vn_td(ji,jj ,jk), 0._wp ) - & 314 & MIN( e1v(ji,jj-1)*e3v(ji,jj-1,jk,Kmm)*vv(ji,jj-1,jk,Kmm) + vn_td(ji,jj-1,jk), 0._wp ) ) & 369 & + ( MAX( e1v(ji,jj ) * e3v(ji,jj ,jk,Kmm) & 370 & * vv (ji,jj ,jk,Kmm) + vn_td(ji,jj ,jk), 0._wp ) - & 371 & MIN( e1v(ji,jj-1) * e3v(ji,jj-1,jk,Kmm) & 372 & * vv (ji,jj-1,jk,Kmm) + vn_td(ji,jj-1,jk), 0._wp ) ) & 315 373 & * r1_e1e2t(ji,jj) & 316 374 & ) * z1_e3t 317 375 END_3D 318 376 ELSE 319 DO_3D _00_00(1, jpkm1 )377 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 320 378 z1_e3t = 1._wp / e3t(ji,jj,jk,Kmm) 321 ! 2*rn_Dt and not rDt (for restartability)322 Cu_adv(ji,jj,jk) = 2._wp * rn_Dt * ( ( MAX( ww(ji,jj,jk) , 0._wp ) - MIN( ww(ji,jj,jk+1) , 0._wp ) ) &379 Cu_adv(ji,jj,jk) = zdt * & 380 & ( ( MAX( ww(ji,jj,jk) , 0._wp ) - MIN( ww(ji,jj,jk+1) , 0._wp ) ) & 323 381 & + ( MAX( e2u(ji ,jj)*e3u(ji ,jj,jk,Kmm)*uu(ji ,jj,jk,Kmm), 0._wp ) - & 324 382 & MIN( e2u(ji-1,jj)*e3u(ji-1,jj,jk,Kmm)*uu(ji-1,jj,jk,Kmm), 0._wp ) ) & … … 330 388 END_3D 331 389 ENDIF 332 CALL lbc_lnk( 'sshwzv', Cu_adv, 'T', 1. )390 CALL lbc_lnk( 'sshwzv', Cu_adv, 'T', 1.0_wp ) 333 391 ! 334 392 CALL iom_put("Courant",Cu_adv) 335 393 ! 336 394 IF( MAXVAL( Cu_adv(:,:,:) ) > Cu_min ) THEN ! Quick check if any breaches anywhere 337 DO_3DS _11_11( jpkm1, 2, -1 )395 DO_3DS( 1, 1, 1, 1, jpkm1, 2, -1 ) ! or scan Courant criterion and partition ! w where necessary 338 396 ! 339 397 zCu = MAX( Cu_adv(ji,jj,jk) , Cu_adv(ji,jj,jk-1) ) -
NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/OCE/DYN/wet_dry.F90
r12489 r13710 33 33 !! * Substitutions 34 34 # include "do_loop_substitute.h90" 35 # include "domzgr_substitute.h90" 35 36 !!---------------------------------------------------------------------- 36 37 !! critical depths,filters, limiters,and masks for Wetting and Drying … … 56 57 REAL(wp), PUBLIC :: ssh_ref !: height of z=0 with respect to the geoid; 57 58 58 LOGICAL, PUBLIC :: ll_wd !: Wetting/drying activation switch if either ln_wd_il or ln_wd_dl59 LOGICAL, PUBLIC :: ll_wd = .FALSE. !: Wetting/drying activation switch (ln_wd_il or ln_wd_dl) <- default def if wad_init not called 59 60 60 61 PUBLIC wad_init ! initialisation routine called by step.F90 … … 110 111 111 112 r_rn_wdmin1 = 1 / rn_wdmin1 112 ll_wd = .FALSE.113 113 IF( ln_wd_il .OR. ln_wd_dl ) THEN 114 114 ll_wd = .TRUE. … … 173 173 ! 174 174 wdmask(:,:) = 1._wp 175 DO_2D _01_01175 DO_2D( 0, 1, 0, 1 ) 176 176 ! 177 177 IF( tmask(ji,jj,1) < 0.5_wp ) CYCLE ! we don't care about land cells … … 197 197 wdramp(:,:) = min((ht_0(:,:) + psshb1(:,:) - rn_wdmin1)/(rn_wdmin0 - rn_wdmin1),1.0_wp) 198 198 !jth assume don't need a lbc_lnk here 199 DO_2D _10_10199 DO_2D( 1, 0, 1, 0 ) 200 200 wdrampu(ji,jj) = MIN( wdramp(ji,jj) , wdramp(ji+1,jj) ) 201 201 wdrampv(ji,jj) = MIN( wdramp(ji,jj) , wdramp(ji,jj+1) ) … … 210 210 jflag = 0 ! flag indicating if any further iterations are needed 211 211 ! 212 DO_2D _01_01212 DO_2D( 0, 1, 0, 1 ) 213 213 IF( tmask(ji, jj, 1) < 0.5_wp ) CYCLE 214 214 IF( ht_0(ji,jj) > zdepwd ) CYCLE … … 241 241 ENDIF 242 242 END_2D 243 CALL lbc_lnk_multi( 'wet_dry', zwdlmtu, 'U', 1. , zwdlmtv, 'V', 1.)243 CALL lbc_lnk_multi( 'wet_dry', zwdlmtu, 'U', 1.0_wp, zwdlmtv, 'V', 1.0_wp ) 244 244 ! 245 245 CALL mpp_max('wet_dry', jflag) !max over the global domain … … 257 257 ! 258 258 !!gm TO BE SUPPRESSED ? these lbc_lnk are useless since zwdlmtu and zwdlmtv are defined everywhere ! 259 CALL lbc_lnk_multi( 'wet_dry', puu(:,:,:,Kmm) , 'U', -1. , pvv(:,:,:,Kmm) , 'V', -1.)260 CALL lbc_lnk_multi( 'wet_dry', uu_b(:,:,Kmm), 'U', -1. , vv_b(:,:,Kmm), 'V', -1.)259 CALL lbc_lnk_multi( 'wet_dry', puu(:,:,:,Kmm) , 'U', -1.0_wp, pvv(:,:,:,Kmm) , 'V', -1.0_wp ) 260 CALL lbc_lnk_multi( 'wet_dry', uu_b(:,:,Kmm), 'U', -1.0_wp, vv_b(:,:,Kmm), 'V', -1.0_wp ) 261 261 !!gm 262 262 ! … … 306 306 zwdlmtv(:,:) = 1._wp 307 307 ! 308 DO_2D _01_01308 DO_2D( 0, 1, 0, 1 ) ! Horizontal Flux in u and v direction 309 309 ! 310 310 IF( tmask(ji, jj, 1 ) < 0.5_wp) CYCLE ! we don't care about land cells … … 332 332 jflag = 0 ! flag indicating if any further iterations are needed 333 333 ! 334 DO_2D _01_01334 DO_2D( 0, 1, 0, 1 ) 335 335 ! 336 336 IF( tmask(ji, jj, 1 ) < 0.5_wp ) CYCLE … … 366 366 END_2D 367 367 ! 368 CALL lbc_lnk_multi( 'wet_dry', zwdlmtu, 'U', 1. , zwdlmtv, 'V', 1.)368 CALL lbc_lnk_multi( 'wet_dry', zwdlmtu, 'U', 1.0_wp, zwdlmtv, 'V', 1.0_wp ) 369 369 ! 370 370 CALL mpp_max('wet_dry', jflag) !max over the global domain … … 378 378 ! 379 379 !!gm THIS lbc_lnk is useless since it is already done at the end of the jk1-loop 380 CALL lbc_lnk_multi( 'wet_dry', zflxu, 'U', -1. , zflxv, 'V', -1.)380 CALL lbc_lnk_multi( 'wet_dry', zflxu, 'U', -1.0_wp, zflxv, 'V', -1.0_wp ) 381 381 !!gm end 382 382 !
Note: See TracChangeset
for help on using the changeset viewer.
