Changeset 14834 for NEMO/trunk/src/OCE/DYN
- Timestamp:
- 2021-05-11T11:24:44+02:00 (3 years ago)
- Location:
- NEMO/trunk/src/OCE/DYN
- Files:
-
- 15 edited
- 2 copied
-
divhor.F90 (modified) (4 diffs)
-
dynadv_cen2.F90 (modified) (2 diffs)
-
dynadv_ubs.F90 (modified) (3 diffs)
-
dynatf.F90 (modified) (4 diffs)
-
dynatf_qco.F90 (modified) (4 diffs)
-
dynhpg.F90 (modified) (18 diffs)
-
dynkeg.F90 (modified) (4 diffs)
-
dynldf_iso.F90 (modified) (16 diffs)
-
dynldf_iso_lf.F90 (copied) (copied from NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/DYN/dynldf_iso_lf.F90)
-
dynldf_lap_blp.F90 (modified) (12 diffs)
-
dynldf_lap_blp_lf.F90 (copied) (copied from NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/DYN/dynldf_lap_blp_lf.F90)
-
dynspg_ts.F90 (modified) (1 diff)
-
dynvor.F90 (modified) (21 diffs)
-
dynzad.F90 (modified) (3 diffs)
-
dynzdf.F90 (modified) (3 diffs)
-
sshwzv.F90 (modified) (8 diffs)
-
wet_dry.F90 (modified) (1 diff)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/trunk/src/OCE/DYN/divhor.F90
r14820 r14834 64 64 ! 65 65 INTEGER :: ji, jj, jk ! dummy loop indices 66 REAL(wp) :: zraur, zdep ! local scalars67 REAL(wp), DIMENSION(jpi,jpj) :: ztmp68 66 !!---------------------------------------------------------------------- 69 67 ! … … 71 69 ! 72 70 IF( kt == nit000 ) THEN 73 IF(lwp) WRITE(numout,*) 74 IF(lwp) WRITE(numout,*) 'div_hor : horizontal velocity divergence ' 75 IF(lwp) WRITE(numout,*) '~~~~~~~ ' 76 hdiv(:,:,:) = 0._wp ! initialize hdiv for the halos at the first time step 71 IF( .NOT. l_istiled .OR. ntile == 1 ) THEN ! Do only on the first tile 72 IF(lwp) WRITE(numout,*) 73 IF(lwp) WRITE(numout,*) 'div_hor : horizontal velocity divergence ' 74 IF(lwp) WRITE(numout,*) '~~~~~~~ ' 75 ENDIF 76 DO_3D_OVR( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpk ) 77 hdiv(ji,jj,jk) = 0._wp ! initialize hdiv for the halos at the first time step 78 END_3D 77 79 ENDIF 78 80 ! 79 DO_3D ( 0, 0, 0, 0, 1, jpkm1 ) !== Horizontal divergence ==!81 DO_3D_OVR( nn_hls-1, nn_hls, nn_hls-1, nn_hls, 1, jpkm1 ) !== Horizontal divergence ==! 80 82 ! round brackets added to fix the order of floating point operations 81 83 ! needed to ensure halo 1 - halo 2 compatibility … … 84 86 & ) & ! bracket for halo 1 - halo 2 compatibility 85 87 & + ( e1v(ji,jj ) * e3v(ji,jj ,jk,Kmm) * vv(ji,jj ,jk,Kmm) & 86 & - e1v(ji,jj-1) * e3v(ji,jj-1,jk,Kmm) * vv(ji,jj-1,jk,Kmm) & 88 & - e1v(ji,jj-1) * e3v(ji,jj-1,jk,Kmm) * vv(ji,jj-1,jk,Kmm) & 87 89 & ) & ! bracket for halo 1 - halo 2 compatibility 88 90 & ) * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm) … … 95 97 ! 96 98 #endif 97 !98 99 IF( ln_isf ) CALL isf_hdiv( kt, Kmm, hdiv ) !== ice shelf ==! (update hdiv field) 99 100 ! 100 CALL lbc_lnk( 'divhor', hdiv, 'T', 1.0_wp ) ! (no sign change)101 IF (nn_hls==1) CALL lbc_lnk( 'divhor', hdiv, 'T', 1.0_wp ) ! (no sign change) 101 102 ! 102 103 IF( ln_timing ) CALL timing_stop('div_hor') -
NEMO/trunk/src/OCE/DYN/dynadv_cen2.F90
r13497 r14834 52 52 ! 53 53 INTEGER :: ji, jj, jk ! dummy loop indices 54 REAL(wp), DIMENSION( jpi,jpj,jpk) :: zfu_t, zfu_f, zfu_uw, zfu55 REAL(wp), DIMENSION( jpi,jpj,jpk) :: zfv_t, zfv_f, zfv_vw, zfv, zfw54 REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zfu_t, zfu_f, zfu_uw, zfu 55 REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zfv_t, zfv_f, zfv_vw, zfv, zfw 56 56 !!---------------------------------------------------------------------- 57 57 ! 58 IF( kt == nit000 .AND. lwp ) THEN 59 WRITE(numout,*) 60 WRITE(numout,*) 'dyn_adv_cen2 : 2nd order flux form momentum advection' 61 WRITE(numout,*) '~~~~~~~~~~~~' 58 IF( .NOT. l_istiled .OR. ntile == 1 ) THEN ! Do only on the first tile 59 IF( kt == nit000 .AND. lwp ) THEN 60 WRITE(numout,*) 61 WRITE(numout,*) 'dyn_adv_cen2 : 2nd order flux form momentum advection' 62 WRITE(numout,*) '~~~~~~~~~~~~' 63 ENDIF 62 64 ENDIF 63 65 ! … … 70 72 ! 71 73 DO jk = 1, jpkm1 ! horizontal transport 72 zfu(:,:,jk) = 0.25_wp * e2u(:,:) * e3u(:,:,jk,Kmm) * puu(:,:,jk,Kmm) 73 zfv(:,:,jk) = 0.25_wp * e1v(:,:) * e3v(:,:,jk,Kmm) * pvv(:,:,jk,Kmm) 74 DO_2D( 1, 1, 1, 1 ) 75 zfu(ji,jj,jk) = 0.25_wp * e2u(ji,jj) * e3u(ji,jj,jk,Kmm) * puu(ji,jj,jk,Kmm) 76 zfv(ji,jj,jk) = 0.25_wp * e1v(ji,jj) * e3v(ji,jj,jk,Kmm) * pvv(ji,jj,jk,Kmm) 77 END_2D 74 78 DO_2D( 1, 0, 1, 0 ) ! horizontal momentum fluxes (at T- and F-point) 75 79 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) ) -
NEMO/trunk/src/OCE/DYN/dynadv_ubs.F90
r14820 r14834 75 75 INTEGER :: ji, jj, jk ! dummy loop indices 76 76 REAL(wp) :: zui, zvj, zfuj, zfvi, zl_u, zl_v ! local scalars 77 REAL(wp), DIMENSION( jpi,jpj,jpk) :: zfu_t, zfu_f, zfu_uw, zfu78 REAL(wp), DIMENSION( jpi,jpj,jpk) :: zfv_t, zfv_f, zfv_vw, zfv, zfw79 REAL(wp), DIMENSION( jpi,jpj,jpk,2) :: zlu_uu, zlu_uv80 REAL(wp), DIMENSION( jpi,jpj,jpk,2) :: zlv_vv, zlv_vu77 REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zfu_t, zfu_f, zfu_uw, zfu 78 REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zfv_t, zfv_f, zfv_vw, zfv, zfw 79 REAL(wp), DIMENSION(A2D(nn_hls),jpk,2) :: zlu_uu, zlu_uv 80 REAL(wp), DIMENSION(A2D(nn_hls),jpk,2) :: zlv_vv, zlv_vu 81 81 !!---------------------------------------------------------------------- 82 82 ! 83 IF( kt == nit000 ) THEN 84 IF(lwp) WRITE(numout,*) 85 IF(lwp) WRITE(numout,*) 'dyn_adv_ubs : UBS flux form momentum advection' 86 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 83 IF( .NOT. l_istiled .OR. ntile == 1 ) THEN ! Do only on the first tile 84 IF( kt == nit000 ) THEN 85 IF(lwp) WRITE(numout,*) 86 IF(lwp) WRITE(numout,*) 'dyn_adv_ubs : UBS flux form momentum advection' 87 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 88 ENDIF 87 89 ENDIF 88 90 ! … … 105 107 ! ! =========================== ! 106 108 ! ! horizontal volume fluxes 107 zfu(:,:,jk) = e2u(:,:) * e3u(:,:,jk,Kmm) * puu(:,:,jk,Kmm) 108 zfv(:,:,jk) = e1v(:,:) * e3v(:,:,jk,Kmm) * pvv(:,:,jk,Kmm) 109 DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 110 zfu(ji,jj,jk) = e2u(ji,jj) * e3u(ji,jj,jk,Kmm) * puu(ji,jj,jk,Kmm) 111 zfv(ji,jj,jk) = e1v(ji,jj) * e3v(ji,jj,jk,Kmm) * pvv(ji,jj,jk,Kmm) 112 END_2D 109 113 ! 110 DO_2D( 0, 0, 0, 0) ! laplacian114 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) ! laplacian 111 115 ! round brackets added to fix the order of floating point operations 112 116 ! needed to ensure halo 1 - halo 2 compatibility 113 zlu_uu(ji,jj,jk,1) = ( (puu (ji+1,jj ,jk,Kbb) - puu (ji ,jj ,jk,Kbb)) + & 114 & (puu (ji-1,jj ,jk,Kbb) - puu (ji ,jj ,jk,Kbb)) ) * umask(ji ,jj ,jk) 115 zlv_vv(ji,jj,jk,1) = ( (pvv (ji ,jj+1,jk,Kbb) - pvv (ji ,jj ,jk,Kbb)) + & 116 & (pvv (ji ,jj-1,jk,Kbb) - pvv (ji ,jj ,jk,Kbb)) ) * vmask(ji ,jj ,jk) 117 zlu_uv(ji,jj,jk,1) = ( puu (ji ,jj+1,jk,Kbb) - puu (ji ,jj ,jk,Kbb) ) * fmask(ji ,jj ,jk) & 118 & - ( puu (ji ,jj ,jk,Kbb) - puu (ji ,jj-1,jk,Kbb) ) * fmask(ji ,jj-1,jk) 119 zlv_vu(ji,jj,jk,1) = ( pvv (ji+1,jj ,jk,Kbb) - pvv (ji ,jj ,jk,Kbb) ) * fmask(ji ,jj ,jk) & 120 & - ( pvv (ji ,jj ,jk,Kbb) - pvv (ji-1,jj ,jk,Kbb) ) * fmask(ji-1,jj ,jk) 121 ! 122 zlu_uu(ji,jj,jk,2) = ( (zfu(ji+1,jj ,jk) - zfu(ji ,jj ,jk)) + & 123 & (zfu(ji-1,jj ,jk) - zfu(ji ,jj ,jk)) ) * umask(ji ,jj ,jk) 124 zlv_vv(ji,jj,jk,2) = ( (zfv(ji ,jj+1,jk) - zfv(ji ,jj ,jk)) + & 125 & (zfv(ji ,jj-1,jk) - zfv(ji ,jj ,jk)) ) * vmask(ji ,jj ,jk) 126 zlu_uv(ji,jj,jk,2) = ( zfu(ji ,jj+1,jk) - zfu(ji ,jj ,jk) ) * fmask(ji ,jj ,jk) & 127 & - ( zfu(ji ,jj ,jk) - zfu(ji ,jj-1,jk) ) * fmask(ji ,jj-1,jk) 128 zlv_vu(ji,jj,jk,2) = ( zfv(ji+1,jj ,jk) - zfv(ji ,jj ,jk) ) * fmask(ji ,jj ,jk) & 129 & - ( zfv(ji ,jj ,jk) - zfv(ji-1,jj ,jk) ) * fmask(ji-1,jj ,jk) 117 zlu_uu(ji,jj,jk,1) = ( ( puu (ji+1,jj ,jk,Kbb) - puu (ji ,jj ,jk,Kbb) & 118 & ) & ! bracket for halo 1 - halo 2 compatibility 119 & + ( puu (ji-1,jj ,jk,Kbb) - puu (ji ,jj ,jk,Kbb) & 120 & ) & ! bracket for halo 1 - halo 2 compatibility 121 & ) * umask(ji ,jj ,jk) 122 zlv_vv(ji,jj,jk,1) = ( ( pvv (ji ,jj+1,jk,Kbb) - pvv (ji ,jj ,jk,Kbb) & 123 & ) & ! bracket for halo 1 - halo 2 compatibility 124 & + ( pvv (ji ,jj-1,jk,Kbb) - pvv (ji ,jj ,jk,Kbb) & 125 & ) & ! bracket for halo 1 - halo 2 compatibility 126 & ) * vmask(ji ,jj ,jk) 127 zlu_uv(ji,jj,jk,1) = ( puu (ji ,jj+1,jk,Kbb) - puu (ji ,jj ,jk,Kbb) ) * fmask(ji ,jj ,jk) & 128 & - ( puu (ji ,jj ,jk,Kbb) - puu (ji ,jj-1,jk,Kbb) ) * fmask(ji ,jj-1,jk) 129 zlv_vu(ji,jj,jk,1) = ( pvv (ji+1,jj ,jk,Kbb) - pvv (ji ,jj ,jk,Kbb) ) * fmask(ji ,jj ,jk) & 130 & - ( pvv (ji ,jj ,jk,Kbb) - pvv (ji-1,jj ,jk,Kbb) ) * fmask(ji-1,jj ,jk) 131 ! 132 ! round brackets added to fix the order of floating point operations 133 ! needed to ensure halo 1 - halo 2 compatibility 134 zlu_uu(ji,jj,jk,2) = ( ( zfu(ji+1,jj ,jk) - zfu(ji ,jj ,jk) & 135 & ) & ! bracket for halo 1 - halo 2 compatibility 136 & + ( zfu(ji-1,jj ,jk) - zfu(ji ,jj ,jk) & 137 & ) & ! bracket for halo 1 - halo 2 compatibility 138 & ) * umask(ji ,jj ,jk) 139 zlv_vv(ji,jj,jk,2) = ( ( zfv(ji ,jj+1,jk) - zfv(ji ,jj ,jk) & 140 & ) & ! bracket for halo 1 - halo 2 compatibility 141 & + ( zfv(ji ,jj-1,jk) - zfv(ji ,jj ,jk) & 142 & ) & ! bracket for halo 1 - halo 2 compatibility 143 & ) * vmask(ji ,jj ,jk) 144 zlu_uv(ji,jj,jk,2) = ( zfu(ji ,jj+1,jk) - zfu(ji ,jj ,jk) ) * fmask(ji ,jj ,jk) & 145 & - ( zfu(ji ,jj ,jk) - zfu(ji ,jj-1,jk) ) * fmask(ji ,jj-1,jk) 146 zlv_vu(ji,jj,jk,2) = ( zfv(ji+1,jj ,jk) - zfv(ji ,jj ,jk) ) * fmask(ji ,jj ,jk) & 147 & - ( zfv(ji ,jj ,jk) - zfv(ji-1,jj ,jk) ) * fmask(ji-1,jj ,jk) 130 148 END_2D 131 149 END DO 132 CALL lbc_lnk( 'dynadv_ubs', zlu_uu(:,:,:,1), 'U', -1.0_wp , zlu_uv(:,:,:,1), 'U', -1.0_wp, &133 &zlu_uu(:,:,:,2), 'U', -1.0_wp , zlu_uv(:,:,:,2), 'U', -1.0_wp, &134 &zlv_vv(:,:,:,1), 'V', -1.0_wp , zlv_vu(:,:,:,1), 'V', -1.0_wp, &135 &zlv_vv(:,:,:,2), 'V', -1.0_wp , zlv_vu(:,:,:,2), 'V', -1.0_wp )150 IF( nn_hls == 1 ) CALL lbc_lnk( 'dynadv_ubs', zlu_uu(:,:,:,1), 'U', -1.0_wp , zlu_uv(:,:,:,1), 'U', -1.0_wp, & 151 & zlu_uu(:,:,:,2), 'U', -1.0_wp , zlu_uv(:,:,:,2), 'U', -1.0_wp, & 152 & zlv_vv(:,:,:,1), 'V', -1.0_wp , zlv_vu(:,:,:,1), 'V', -1.0_wp, & 153 & zlv_vv(:,:,:,2), 'V', -1.0_wp , zlv_vu(:,:,:,2), 'V', -1.0_wp ) 136 154 ! 137 155 ! ! ====================== ! … … 139 157 DO jk = 1, jpkm1 ! ====================== ! 140 158 ! ! horizontal volume fluxes 141 zfu(:,:,jk) = 0.25_wp * e2u(:,:) * e3u(:,:,jk,Kmm) * puu(:,:,jk,Kmm) 142 zfv(:,:,jk) = 0.25_wp * e1v(:,:) * e3v(:,:,jk,Kmm) * pvv(:,:,jk,Kmm) 159 DO_2D( 1, 1, 1, 1 ) 160 zfu(ji,jj,jk) = 0.25_wp * e2u(ji,jj) * e3u(ji,jj,jk,Kmm) * puu(ji,jj,jk,Kmm) 161 zfv(ji,jj,jk) = 0.25_wp * e1v(ji,jj) * e3v(ji,jj,jk,Kmm) * pvv(ji,jj,jk,Kmm) 162 END_2D 143 163 ! 144 164 DO_2D( 1, 0, 1, 0 ) ! horizontal momentum fluxes at T- and F-point -
NEMO/trunk/src/OCE/DYN/dynatf.F90
r14433 r14834 201 201 IF( ln_linssh ) THEN ! Fixed volume ! 202 202 ! ! =============! 203 DO_3D( 1, 1, 1, 1, 1, jpkm1 )203 DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 ) 204 204 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) ) 205 205 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) ) … … 237 237 CALL dom_vvl_interpol( pe3t(:,:,:,Kmm), pe3u(:,:,:,Kmm), 'U' ) 238 238 CALL dom_vvl_interpol( pe3t(:,:,:,Kmm), pe3v(:,:,:,Kmm), 'V' ) 239 DO_3D( 1, 1, 1, 1, 1, jpkm1 )239 DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 ) 240 240 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) ) 241 241 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) ) … … 248 248 CALL dom_vvl_interpol( pe3t(:,:,:,Kmm), ze3u_f, 'U' ) 249 249 CALL dom_vvl_interpol( pe3t(:,:,:,Kmm), ze3v_f, 'V' ) 250 DO_3D( 1, 1, 1, 1, 1, jpkm1 )250 DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 ) 251 251 zue3a = pe3u(ji,jj,jk,Kaa) * puu(ji,jj,jk,Kaa) 252 252 zve3a = pe3v(ji,jj,jk,Kaa) * pvv(ji,jj,jk,Kaa) … … 285 285 ENDIF ! .NOT. l_1st_euler 286 286 ! 287 ! This is needed for dyn_ldf_blp to be restartable 288 IF( nn_hls == 2 ) CALL lbc_lnk( 'dynatf', puu(:,:,:,Kmm), 'U', -1.0_wp, pvv(:,:,:,Kmm), 'V', -1.0_wp ) 287 289 ! Set "now" and "before" barotropic velocities for next time step: 288 290 ! JC: Would be more clever to swap variables than to make a full vertical -
NEMO/trunk/src/OCE/DYN/dynatf_qco.F90
r14475 r14834 139 139 IF( ln_linssh ) THEN ! Fixed volume ! 140 140 ! ! =============! 141 DO_3D( 1, 1, 1, 1, 1, jpkm1 )141 DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 ) 142 142 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) ) 143 143 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) ) … … 149 149 IF( ln_dynadv_vec ) THEN ! Asselin filter applied on velocity 150 150 ! Before filtered scale factor at (u/v)-points 151 DO_3D( 1, 1, 1, 1, 1, jpkm1 )151 DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 ) 152 152 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) ) 153 153 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) ) … … 156 156 ELSE ! Asselin filter applied on thickness weighted velocity 157 157 ! 158 DO_3D( 1, 1, 1, 1, 1, jpkm1 )158 DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 ) 159 159 zue3a = ( 1._wp + r3u(ji,jj,Kaa) * umask(ji,jj,jk) ) * puu(ji,jj,jk,Kaa) 160 160 zve3a = ( 1._wp + r3v(ji,jj,Kaa) * vmask(ji,jj,jk) ) * pvv(ji,jj,jk,Kaa) … … 195 195 ENDIF ! .NOT. l_1st_euler 196 196 ! 197 ! This is needed for dyn_ldf_blp to be restartable 198 IF( nn_hls == 2 ) CALL lbc_lnk( 'dynatfqco', puu(:,:,:,Kmm), 'U', -1.0_wp, pvv(:,:,:,Kmm), 'V', -1.0_wp ) 199 197 200 ! Set "now" and "before" barotropic velocities for next time step: 198 201 ! JC: Would be more clever to swap variables than to make a full vertical -
NEMO/trunk/src/OCE/DYN/dynhpg.F90
r14433 r14834 266 266 INTEGER :: ji, jj, jk ! dummy loop indices 267 267 REAL(wp) :: zcoef0, zcoef1 ! temporary scalars 268 REAL(wp), DIMENSION(jpi,jpj) :: zhpi, zhpj 269 !!---------------------------------------------------------------------- 270 ! 271 IF( kt == nit000 ) THEN 272 IF(lwp) WRITE(numout,*) 273 IF(lwp) WRITE(numout,*) 'dyn:hpg_zco : hydrostatic pressure gradient trend' 274 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ z-coordinate case ' 268 REAL(wp), DIMENSION(A2D(nn_hls)) :: zhpi, zhpj 269 !!---------------------------------------------------------------------- 270 ! 271 IF( .NOT. l_istiled .OR. ntile == 1 ) THEN ! Do only on the first tile 272 IF( kt == nit000 ) THEN 273 IF(lwp) WRITE(numout,*) 274 IF(lwp) WRITE(numout,*) 'dyn:hpg_zco : hydrostatic pressure gradient trend' 275 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ z-coordinate case ' 276 ENDIF 275 277 ENDIF 276 278 ! … … 318 320 INTEGER :: iku, ikv ! temporary integers 319 321 REAL(wp) :: zcoef0, zcoef1, zcoef2, zcoef3 ! temporary scalars 320 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zhpi, zhpj 321 REAL(wp), DIMENSION(jpi,jpj,jpts) :: zgtsu, zgtsv 322 REAL(wp), DIMENSION(jpi,jpj) :: zgru, zgrv 323 !!---------------------------------------------------------------------- 324 ! 325 IF( kt == nit000 ) THEN 326 IF(lwp) WRITE(numout,*) 327 IF(lwp) WRITE(numout,*) 'dyn:hpg_zps : hydrostatic pressure gradient trend' 328 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ z-coordinate with partial steps - vector optimization' 322 REAL(wp), DIMENSION(A2D(nn_hls),jpk ) :: zhpi, zhpj 323 REAL(wp), DIMENSION(A2D(nn_hls),jpts) :: zgtsu, zgtsv 324 REAL(wp), DIMENSION(A2D(nn_hls) ) :: zgru, zgrv 325 !!---------------------------------------------------------------------- 326 ! 327 IF( .NOT. l_istiled .OR. ntile == 1 ) THEN ! Do only on the first tile 328 IF( kt == nit000 ) THEN 329 IF(lwp) WRITE(numout,*) 330 IF(lwp) WRITE(numout,*) 'dyn:hpg_zps : hydrostatic pressure gradient trend' 331 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ z-coordinate with partial steps - vector optimization' 332 ENDIF 329 333 ENDIF 330 334 … … 410 414 REAL(wp) :: zcoef0, zuap, zvap, ztmp ! local scalars 411 415 LOGICAL :: ll_tmp1, ll_tmp2 ! local logical variables 412 REAL(wp), DIMENSION( jpi,jpj,jpk):: zhpi, zhpj416 REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zhpi, zhpj 413 417 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: zcpx, zcpy !W/D pressure filter 414 418 !!---------------------------------------------------------------------- 415 419 ! 416 IF( ln_wd_il ) ALLOCATE(zcpx(jpi,jpj), zcpy(jpi,jpj)) 417 ! 418 IF( kt == nit000 ) THEN 419 IF(lwp) WRITE(numout,*) 420 IF(lwp) WRITE(numout,*) 'dyn:hpg_sco : hydrostatic pressure gradient trend' 421 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ s-coordinate case, OCE original scheme used' 420 IF( ln_wd_il ) ALLOCATE(zcpx(A2D(nn_hls)), zcpy(A2D(nn_hls))) 421 ! 422 IF( .NOT. l_istiled .OR. ntile == 1 ) THEN ! Do only on the first tile 423 IF( kt == nit000 ) THEN 424 IF(lwp) WRITE(numout,*) 425 IF(lwp) WRITE(numout,*) 'dyn:hpg_sco : hydrostatic pressure gradient trend' 426 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ s-coordinate case, OCE original scheme used' 427 ENDIF 422 428 ENDIF 423 429 ! … … 462 468 END IF 463 469 END_2D 464 CALL lbc_lnk( 'dynhpg', zcpx, 'U', 1.0_wp, zcpy, 'V', 1.0_wp )465 470 END IF 466 471 ! … … 548 553 REAL(wp) :: ze3w, ze3wi1, ze3wj1 ! local scalars 549 554 REAL(wp) :: zcoef0, zuap, zvap ! - - 550 REAL(wp), DIMENSION( jpi,jpj,jpk ) :: zhpi, zhpj551 REAL(wp), DIMENSION( jpi,jpj,jpts) :: zts_top552 REAL(wp), DIMENSION( jpi,jpj) :: zrhdtop_oce555 REAL(wp), DIMENSION(A2D(nn_hls),jpk ) :: zhpi, zhpj 556 REAL(wp), DIMENSION(A2D(nn_hls),jpts) :: zts_top 557 REAL(wp), DIMENSION(A2D(nn_hls)) :: zrhdtop_oce 553 558 !!---------------------------------------------------------------------- 554 559 ! … … 560 565 ! compute rhd at the ice/oce interface (ocean side) 561 566 ! usefull to reduce residual current in the test case ISOMIP with no melting 562 DO ji = 1, jpi 563 DO jj = 1, jpj 564 ikt = mikt(ji,jj) 565 zts_top(ji,jj,1) = ts(ji,jj,ikt,1,Kmm) 566 zts_top(ji,jj,2) = ts(ji,jj,ikt,2,Kmm) 567 END DO 568 END DO 567 DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 568 ikt = mikt(ji,jj) 569 zts_top(ji,jj,1) = ts(ji,jj,ikt,1,Kmm) 570 zts_top(ji,jj,2) = ts(ji,jj,ikt,2,Kmm) 571 END_2D 569 572 CALL eos( zts_top, risfdep, zrhdtop_oce ) 570 573 … … 636 639 INTEGER :: iktb, iktt ! jk indices at tracer points for top and bottom points 637 640 REAL(wp) :: zcoef0, zep, cffw ! temporary scalars 638 REAL(wp) :: z_grav_10, z1_12 641 REAL(wp) :: z_grav_10, z1_12, z1_cff 639 642 REAL(wp) :: cffu, cffx ! " " 640 643 REAL(wp) :: cffv, cffy ! " " 641 644 LOGICAL :: ll_tmp1, ll_tmp2 ! local logical variables 642 REAL(wp), DIMENSION( jpi,jpj,jpk) :: zhpi, zhpj643 644 REAL(wp), DIMENSION( jpi,jpj,jpk) :: zdzx, zdzy, zdzz ! Primitive grid differences ('delta_xyz')645 REAL(wp), DIMENSION( jpi,jpj,jpk) :: zdz_i, zdz_j, zdz_k ! Harmonic average of primitive grid differences ('d_xyz')646 REAL(wp), DIMENSION( jpi,jpj,jpk) :: zdrhox, zdrhoy, zdrhoz ! Primitive rho differences ('delta_rho')647 REAL(wp), DIMENSION( jpi,jpj,jpk) :: zdrho_i, zdrho_j, zdrho_k ! Harmonic average of primitive rho differences ('d_rho')648 REAL(wp), DIMENSION( jpi,jpj,jpk) :: z_rho_i, z_rho_j, z_rho_k ! Face intergrals649 REAL(wp), DIMENSION( jpi,jpj) :: zz_dz_i, zz_dz_j, zz_drho_i, zz_drho_j ! temporary arrays645 REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zhpi, zhpj 646 647 REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zdzx, zdzy, zdzz ! Primitive grid differences ('delta_xyz') 648 REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zdz_i, zdz_j, zdz_k ! Harmonic average of primitive grid differences ('d_xyz') 649 REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zdrhox, zdrhoy, zdrhoz ! Primitive rho differences ('delta_rho') 650 REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zdrho_i, zdrho_j, zdrho_k ! Harmonic average of primitive rho differences ('d_rho') 651 REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: z_rho_i, z_rho_j, z_rho_k ! Face intergrals 652 REAL(wp), DIMENSION(A2D(nn_hls)) :: zz_dz_i, zz_dz_j, zz_drho_i, zz_drho_j ! temporary arrays 650 653 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: zcpx, zcpy !W/D pressure filter 651 654 !!---------------------------------------------------------------------- 652 655 ! 653 656 IF( ln_wd_il ) THEN 654 ALLOCATE( zcpx( jpi,jpj) , zcpy(jpi,jpj) )657 ALLOCATE( zcpx(A2D(nn_hls)) , zcpy(A2D(nn_hls)) ) 655 658 DO_2D( 0, 0, 0, 0 ) 656 659 ll_tmp1 = MIN( ssh(ji,jj,Kmm) , ssh(ji+1,jj,Kmm) ) > & … … 689 692 END IF 690 693 END_2D 691 CALL lbc_lnk( 'dynhpg', zcpx, 'U', 1.0_wp, zcpy, 'V', 1.0_wp )692 694 END IF 693 695 694 IF( kt == nit000 ) THEN 695 IF(lwp) WRITE(numout,*) 696 IF(lwp) WRITE(numout,*) 'dyn:hpg_djc : hydrostatic pressure gradient trend' 697 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ s-coordinate case, density Jacobian with cubic polynomial scheme' 696 IF( .NOT. l_istiled .OR. ntile == 1 ) THEN ! Do only on the first tile 697 IF( kt == nit000 ) THEN 698 IF(lwp) WRITE(numout,*) 699 IF(lwp) WRITE(numout,*) 'dyn:hpg_djc : hydrostatic pressure gradient trend' 700 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ s-coordinate case, density Jacobian with cubic polynomial scheme' 701 ENDIF 698 702 ENDIF 699 703 … … 723 727 zdz_k (:,:,:) = 0._wp 724 728 725 DO_3D( 1, 1, 1, 1, 2, jpk-2 ) 726 cffw = 2._wp * zdrhoz(ji ,jj ,jk) * zdrhoz(ji,jj,jk+1) 727 IF( cffw > zep) THEN 728 zdrho_k(ji,jj,jk) = cffw / ( zdrhoz(ji,jj,jk) + zdrhoz(ji,jj,jk+1) ) 729 ENDIF 729 DO_3D( 1, 1, 1, 1, 2, jpk-2 ) 730 cffw = MAX( 2._wp * zdrhoz(ji,jj,jk) * zdrhoz(ji,jj,jk+1), 0._wp ) 731 z1_cff = zdrhoz(ji,jj,jk) + zdrhoz(ji,jj,jk+1) 732 zdrho_k(ji,jj,jk) = cffw / SIGN( MAX( ABS(z1_cff), zep ), z1_cff ) 730 733 zdz_k(ji,jj,jk) = 2._wp * zdzz(ji,jj,jk) * zdzz(ji,jj,jk+1) & 731 734 & / ( zdzz(ji,jj,jk) + zdzz(ji,jj,jk+1) ) … … 737 740 738 741 ! mb for sea-ice shelves we will need to re-write this upper boundary condition in the same form as the lower boundary condition 739 zdrho_k(:,:,1) = aco_bc_vrt * ( rhd (:,:,2) - rhd (:,:,1) ) - bco_bc_vrt * zdrho_k(:,:,2) 740 zdz_k (:,:,1) = aco_bc_vrt * (-gde3w(:,:,2) + gde3w(:,:,1) ) - bco_bc_vrt * zdz_k (:,:,2) 742 DO_2D( 1, 1, 1, 1 ) 743 zdrho_k(ji,jj,1) = aco_bc_vrt * ( rhd (ji,jj,2) - rhd (ji,jj,1) ) - bco_bc_vrt * zdrho_k(ji,jj,2) 744 zdz_k (ji,jj,1) = aco_bc_vrt * (-gde3w(ji,jj,2) + gde3w(ji,jj,1) ) - bco_bc_vrt * zdz_k (ji,jj,2) 745 END_2D 741 746 742 747 DO_2D( 1, 1, 1, 1 ) … … 785 790 ! 5. compute and store elementary horizontal differences in provisional arrays 786 791 !---------------------------------------------------------------------------------------- 787 788 DO_3D( 1, 0, 1, 0, 1, jpkm1 ) 789 zdrhox(ji,jj,jk) = rhd (ji+1,jj ,jk) - rhd (ji,jj,jk ) 790 zdzx (ji,jj,jk) = - gde3w(ji+1,jj ,jk) + gde3w(ji,jj,jk ) 791 zdrhoy(ji,jj,jk) = rhd (ji ,jj+1,jk) - rhd (ji,jj,jk ) 792 zdzy (ji,jj,jk) = - gde3w(ji ,jj+1,jk) + gde3w(ji,jj,jk ) 793 END_3D 794 795 CALL lbc_lnk( 'dynhpg', zdrhox, 'U', 1., zdzx, 'U', 1., zdrhoy, 'V', 1., zdzy, 'V', 1. ) 792 zdrhox(:,:,:) = 0._wp 793 zdzx (:,:,:) = 0._wp 794 zdrhoy(:,:,:) = 0._wp 795 zdzy (:,:,:) = 0._wp 796 797 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 ) 798 zdrhox(ji,jj,jk) = rhd (ji+1,jj ,jk) - rhd (ji ,jj ,jk) 799 zdzx (ji,jj,jk) = gde3w(ji ,jj ,jk) - gde3w(ji+1,jj ,jk) 800 zdrhoy(ji,jj,jk) = rhd (ji ,jj+1,jk) - rhd (ji ,jj ,jk) 801 zdzy (ji,jj,jk) = gde3w(ji ,jj ,jk) - gde3w(ji ,jj+1,jk) 802 END_3D 803 804 IF( nn_hls == 1 ) CALL lbc_lnk( 'dynhpg', zdrhox, 'U', -1._wp, zdzx, 'U', -1._wp, zdrhoy, 'V', -1._wp, zdzy, 'V', -1._wp ) 796 805 797 806 !------------------------------------------------------------------------- … … 800 809 801 810 DO_3D( 0, 1, 0, 1, 1, jpkm1 ) 802 cffu = 2._wp * zdrhox(ji-1,jj ,jk) * zdrhox(ji,jj,jk ) 803 IF( cffu > zep ) THEN 804 zdrho_i(ji,jj,jk) = cffu / ( zdrhox(ji-1,jj,jk) + zdrhox(ji,jj,jk) ) 805 ELSE 806 zdrho_i(ji,jj,jk ) = 0._wp 807 ENDIF 808 809 cffx = 2._wp * zdzx (ji-1,jj ,jk) * zdzx (ji,jj,jk ) 810 IF( cffx > zep ) THEN 811 zdz_i(ji,jj,jk) = cffx / ( zdzx(ji-1,jj,jk) + zdzx(ji,jj,jk) ) 812 ELSE 813 zdz_i(ji,jj,jk) = 0._wp 814 ENDIF 815 816 cffv = 2._wp * zdrhoy(ji ,jj-1,jk) * zdrhoy(ji,jj,jk ) 817 IF( cffv > zep ) THEN 818 zdrho_j(ji,jj,jk) = cffv / ( zdrhoy(ji,jj-1,jk) + zdrhoy(ji,jj,jk) ) 819 ELSE 820 zdrho_j(ji,jj,jk) = 0._wp 821 ENDIF 822 823 cffy = 2._wp * zdzy (ji ,jj-1,jk) * zdzy (ji,jj,jk ) 824 IF( cffy > zep ) THEN 825 zdz_j(ji,jj,jk) = cffy / ( zdzy(ji,jj-1,jk) + zdzy(ji,jj,jk) ) 826 ELSE 827 zdz_j(ji,jj,jk) = 0._wp 828 ENDIF 811 cffu = MAX( 2._wp * zdrhox(ji-1,jj,jk) * zdrhox(ji,jj,jk), 0._wp ) 812 z1_cff = zdrhox(ji-1,jj,jk) + zdrhox(ji,jj,jk) 813 zdrho_i(ji,jj,jk) = cffu / SIGN( MAX( ABS(z1_cff), zep ), z1_cff ) 814 815 cffx = MAX( 2._wp * zdzx(ji-1,jj,jk) * zdzx(ji,jj,jk), 0._wp ) 816 z1_cff = zdzx(ji-1,jj,jk) + zdzx(ji,jj,jk) 817 zdz_i(ji,jj,jk) = cffx / SIGN( MAX( ABS(z1_cff), zep ), z1_cff ) 818 819 cffv = MAX( 2._wp * zdrhoy(ji,jj-1,jk) * zdrhoy(ji,jj,jk), 0._wp ) 820 z1_cff = zdrhoy(ji,jj-1,jk) + zdrhoy(ji,jj,jk) 821 zdrho_j(ji,jj,jk) = cffv / SIGN( MAX( ABS(z1_cff), zep ), z1_cff ) 822 823 cffy = MAX( 2._wp * zdzy(ji,jj-1,jk) * zdzy(ji,jj,jk), 0._wp ) 824 z1_cff = zdzy(ji,jj-1,jk) + zdzy(ji,jj,jk) 825 zdz_j(ji,jj,jk) = cffy / SIGN( MAX( ABS(z1_cff), zep ), z1_cff ) 829 826 END_3D 830 827 … … 840 837 zz_drho_j(:,:) = zdrho_j(:,:,jk) 841 838 zz_dz_j (:,:) = zdz_j (:,:,jk) 842 DO_2D( 0, 1, 0, 1) 843 ! Walls coming from left: should check from 2 to jpi-1 (and jpj=2-jpj) 844 IF (ji < jpi) THEN 845 IF ( umask(ji,jj,jk) > 0.5_wp .AND. umask(ji-1,jj,jk) < 0.5_wp .AND. umask(ji+1,jj,jk) > 0.5_wp) THEN 846 zz_drho_i(ji,jj) = aco_bc_hor * ( rhd (ji+1,jj,jk) - rhd (ji,jj,jk) ) - bco_bc_hor * zdrho_i(ji+1,jj,jk) 847 zz_dz_i (ji,jj) = aco_bc_hor * (-gde3w(ji+1,jj,jk) + gde3w(ji,jj,jk) ) - bco_bc_hor * zdz_i (ji+1,jj,jk) 848 END IF 839 ! Walls coming from left: should check from 2 to jpi-1 (and jpj=2-jpj) 840 DO_2D( 0, 0, 0, 1 ) 841 IF ( umask(ji,jj,jk) > 0.5_wp .AND. umask(ji-1,jj,jk) < 0.5_wp .AND. umask(ji+1,jj,jk) > 0.5_wp) THEN 842 zz_drho_i(ji,jj) = aco_bc_hor * ( rhd (ji+1,jj,jk) - rhd (ji,jj,jk) ) - bco_bc_hor * zdrho_i(ji+1,jj,jk) 843 zz_dz_i (ji,jj) = aco_bc_hor * (-gde3w(ji+1,jj,jk) + gde3w(ji,jj,jk) ) - bco_bc_hor * zdz_i (ji+1,jj,jk) 849 844 END IF 850 ! Walls coming from right: should check from 3 to jpi (and jpj=2-jpj)851 IF (ji > 2) THEN852 IF ( umask(ji,jj,jk) < 0.5_wp .AND. umask(ji-1,jj,jk) > 0.5_wp .AND. umask(ji-2,jj,jk) > 0.5_wp) THEN853 zz_drho_i(ji,jj) = aco_bc_hor * ( rhd (ji,jj,jk) - rhd (ji-1,jj,jk) ) - bco_bc_hor * zdrho_i(ji-1,jj,jk)854 zz_dz_i (ji,jj) = aco_bc_hor * (-gde3w(ji,jj,jk) + gde3w(ji-1,jj,jk) ) - bco_bc_hor * zdz_i(ji-1,jj,jk)855 END IF845 END_2D 846 ! Walls coming from right: should check from 3 to jpi (and jpj=2-jpj) 847 DO_2D( -1, 1, 0, 1 ) 848 IF ( umask(ji,jj,jk) < 0.5_wp .AND. umask(ji-1,jj,jk) > 0.5_wp .AND. umask(ji-2,jj,jk) > 0.5_wp) THEN 849 zz_drho_i(ji,jj) = aco_bc_hor * ( rhd (ji,jj,jk) - rhd (ji-1,jj,jk) ) - bco_bc_hor * zdrho_i(ji-1,jj,jk) 850 zz_dz_i (ji,jj) = aco_bc_hor * (-gde3w(ji,jj,jk) + gde3w(ji-1,jj,jk) ) - bco_bc_hor * zdz_i (ji-1,jj,jk) 856 851 END IF 857 ! Walls coming from left: should check from 2 to jpj-1 (and jpi=2-jpi)858 IF (jj < jpj) THEN859 IF ( vmask(ji,jj,jk) > 0.5_wp .AND. vmask(ji,jj-1,jk) < 0.5_wp .AND. vmask(ji,jj+1,jk) > 0.5_wp) THEN860 zz_drho_j(ji,jj) = aco_bc_hor * ( rhd (ji,jj+1,jk) - rhd (ji,jj,jk) ) - bco_bc_hor * zdrho_j(ji,jj+1,jk)861 zz_dz_j (ji,jj) = aco_bc_hor * (-gde3w(ji,jj+1,jk) + gde3w(ji,jj,jk) ) - bco_bc_hor * zdz_j(ji,jj+1,jk)862 END IF863 END IF 864 ! Walls coming from right: should check from 3 to jpj (and jpi=2-jpi)865 IF (jj > 2) THEN866 IF ( vmask(ji,jj,jk) < 0.5_wp .AND. vmask(ji,jj-1,jk) > 0.5_wp .AND. vmask(ji,jj-2,jk) > 0.5_wp) THEN867 zz_drho_j(ji,jj) = aco_bc_hor * ( rhd (ji,jj,jk) - rhd (ji,jj-1,jk) ) - bco_bc_hor * zdrho_j(ji,jj-1,jk)868 zz_dz_j (ji,jj) = aco_bc_hor * (-gde3w(ji,jj,jk) + gde3w(ji,jj-1,jk) ) - bco_bc_hor * zdz_j(ji,jj-1,jk)869 END IF852 END_2D 853 ! Walls coming from left: should check from 2 to jpj-1 (and jpi=2-jpi) 854 DO_2D( 0, 1, 0, 0 ) 855 IF ( vmask(ji,jj,jk) > 0.5_wp .AND. vmask(ji,jj-1,jk) < 0.5_wp .AND. vmask(ji,jj+1,jk) > 0.5_wp) THEN 856 zz_drho_j(ji,jj) = aco_bc_hor * ( rhd (ji,jj+1,jk) - rhd (ji,jj,jk) ) - bco_bc_hor * zdrho_j(ji,jj+1,jk) 857 zz_dz_j (ji,jj) = aco_bc_hor * (-gde3w(ji,jj+1,jk) + gde3w(ji,jj,jk) ) - bco_bc_hor * zdz_j (ji,jj+1,jk) 858 END IF 859 END_2D 860 ! Walls coming from right: should check from 3 to jpj (and jpi=2-jpi) 861 DO_2D( 0, 1, -1, 1 ) 862 IF ( vmask(ji,jj,jk) < 0.5_wp .AND. vmask(ji,jj-1,jk) > 0.5_wp .AND. vmask(ji,jj-2,jk) > 0.5_wp) THEN 863 zz_drho_j(ji,jj) = aco_bc_hor * ( rhd (ji,jj,jk) - rhd (ji,jj-1,jk) ) - bco_bc_hor * zdrho_j(ji,jj-1,jk) 864 zz_dz_j (ji,jj) = aco_bc_hor * (-gde3w(ji,jj,jk) + gde3w(ji,jj-1,jk) ) - bco_bc_hor * zdz_j (ji,jj-1,jk) 870 865 END IF 871 866 END_2D … … 974 969 REAL(wp) :: zrhdt1 975 970 REAL(wp) :: zdpdx1, zdpdx2, zdpdy1, zdpdy2 976 REAL(wp), DIMENSION( jpi,jpj) :: zpgu, zpgv ! 2D workspace977 REAL(wp), DIMENSION( jpi,jpj) :: zsshu_n, zsshv_n978 REAL(wp), DIMENSION( jpi,jpj,jpk) :: zdept, zrhh979 REAL(wp), DIMENSION( jpi,jpj,jpk) :: zhpi, zu, zv, fsp, xsp, asp, bsp, csp, dsp971 REAL(wp), DIMENSION(A2D(nn_hls)) :: zpgu, zpgv ! 2D workspace 972 REAL(wp), DIMENSION(A2D(nn_hls)) :: zsshu_n, zsshv_n 973 REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zdept, zrhh 974 REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zhpi, zu, zv, fsp, xsp, asp, bsp, csp, dsp 980 975 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: zcpx, zcpy !W/D pressure filter 981 976 !!---------------------------------------------------------------------- 982 977 ! 983 IF( kt == nit000 ) THEN 984 IF(lwp) WRITE(numout,*) 985 IF(lwp) WRITE(numout,*) 'dyn:hpg_prj : hydrostatic pressure gradient trend' 986 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ s-coordinate case, cubic spline pressure Jacobian' 978 IF( .NOT. l_istiled .OR. ntile == 1 ) THEN ! Do only on the first tile 979 IF( kt == nit000 ) THEN 980 IF(lwp) WRITE(numout,*) 981 IF(lwp) WRITE(numout,*) 'dyn:hpg_prj : hydrostatic pressure gradient trend' 982 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ s-coordinate case, cubic spline pressure Jacobian' 983 ENDIF 987 984 ENDIF 988 985 … … 1001 998 ! 1002 999 IF( ln_wd_il ) THEN 1003 ALLOCATE( zcpx( jpi,jpj) , zcpy(jpi,jpj) )1000 ALLOCATE( zcpx(A2D(nn_hls)) , zcpy(A2D(nn_hls)) ) 1004 1001 DO_2D( 0, 0, 0, 0 ) 1005 ll_tmp1 = MIN( ssh(ji,jj,Kmm) , ssh(ji+1,jj,Kmm) ) > & 1006 & MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) .AND. & 1007 & MAX( ssh(ji,jj,Kmm) + ht_0(ji,jj), ssh(ji+1,jj,Kmm) + ht_0(ji+1,jj) ) & 1008 & > rn_wdmin1 + rn_wdmin2 1009 ll_tmp2 = ( ABS( ssh(ji,jj,Kmm) - ssh(ji+1,jj,Kmm) ) > 1.E-12 ) .AND. ( & 1010 & MAX( ssh(ji,jj,Kmm) , ssh(ji+1,jj,Kmm) ) > & 1011 & MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 ) 1012 1013 IF(ll_tmp1) THEN 1014 zcpx(ji,jj) = 1.0_wp 1015 ELSE IF(ll_tmp2) THEN 1016 ! no worries about ssh(ji+1,jj,Kmm) - ssh(ji ,jj,Kmm) = 0, it won't happen ! here 1017 zcpx(ji,jj) = ABS( (ssh(ji+1,jj,Kmm) + ht_0(ji+1,jj) - ssh(ji,jj,Kmm) - ht_0(ji,jj)) & 1018 & / (ssh(ji+1,jj,Kmm) - ssh(ji ,jj,Kmm)) ) 1019 1020 zcpx(ji,jj) = max(min( zcpx(ji,jj) , 1.0_wp),0.0_wp) 1021 ELSE 1022 zcpx(ji,jj) = 0._wp 1023 END IF 1024 1025 ll_tmp1 = MIN( ssh(ji,jj,Kmm) , ssh(ji,jj+1,Kmm) ) > & 1026 & MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) .AND. & 1027 & MAX( ssh(ji,jj,Kmm) + ht_0(ji,jj), ssh(ji,jj+1,Kmm) + ht_0(ji,jj+1) ) & 1028 & > rn_wdmin1 + rn_wdmin2 1029 ll_tmp2 = ( ABS( ssh(ji,jj,Kmm) - ssh(ji,jj+1,Kmm) ) > 1.E-12 ) .AND. ( & 1030 & MAX( ssh(ji,jj,Kmm) , ssh(ji,jj+1,Kmm) ) > & 1031 & MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 ) 1032 1033 IF(ll_tmp1) THEN 1034 zcpy(ji,jj) = 1.0_wp 1035 ELSE IF(ll_tmp2) THEN 1036 ! no worries about ssh(ji,jj+1,Kmm) - ssh(ji,jj ,Kmm) = 0, it won't happen ! here 1037 zcpy(ji,jj) = ABS( (ssh(ji,jj+1,Kmm) + ht_0(ji,jj+1) - ssh(ji,jj,Kmm) - ht_0(ji,jj)) & 1038 & / (ssh(ji,jj+1,Kmm) - ssh(ji,jj ,Kmm)) ) 1039 zcpy(ji,jj) = max(min( zcpy(ji,jj) , 1.0_wp),0.0_wp) 1040 1002 ll_tmp1 = MIN( ssh(ji,jj,Kmm) , ssh(ji+1,jj,Kmm) ) > & 1003 & MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) .AND. & 1004 & MAX( ssh(ji,jj,Kmm) + ht_0(ji,jj), ssh(ji+1,jj,Kmm) + ht_0(ji+1,jj) ) > & 1005 & rn_wdmin1 + rn_wdmin2 1006 ll_tmp2 = ( ABS( ssh(ji,jj,Kmm) - ssh(ji+1,jj,Kmm) ) > 1.E-12 ) .AND. & 1007 & ( MAX( ssh(ji,jj,Kmm) , ssh(ji+1,jj,Kmm) ) > & 1008 & MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 ) 1009 1010 IF(ll_tmp1) THEN 1011 zcpx(ji,jj) = 1.0_wp 1012 ELSE IF(ll_tmp2) THEN 1013 ! no worries about ssh(ji+1,jj,Kmm) - ssh(ji ,jj,Kmm) = 0, it won't happen ! here 1014 zcpx(ji,jj) = ABS( (ssh(ji+1,jj,Kmm) + ht_0(ji+1,jj) - ssh(ji,jj,Kmm) - ht_0(ji,jj)) & 1015 & / (ssh(ji+1,jj,Kmm) - ssh(ji ,jj,Kmm)) ) 1016 zcpx(ji,jj) = MAX(MIN( zcpx(ji,jj) , 1.0_wp),0.0_wp) 1017 ELSE 1018 zcpx(ji,jj) = 0._wp 1019 END IF 1020 1021 ll_tmp1 = MIN( ssh(ji,jj,Kmm) , ssh(ji,jj+1,Kmm) ) > & 1022 & MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) .AND. & 1023 & MAX( ssh(ji,jj,Kmm) + ht_0(ji,jj), ssh(ji,jj+1,Kmm) + ht_0(ji,jj+1) ) > & 1024 & rn_wdmin1 + rn_wdmin2 1025 ll_tmp2 = ( ABS( ssh(ji,jj,Kmm) - ssh(ji,jj+1,Kmm) ) > 1.E-12 ) .AND. & 1026 & ( MAX( ssh(ji,jj,Kmm) , ssh(ji,jj+1,Kmm) ) > & 1027 & MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 ) 1028 1029 IF(ll_tmp1) THEN 1030 zcpy(ji,jj) = 1.0_wp 1031 ELSE IF(ll_tmp2) THEN 1032 ! no worries about ssh(ji,jj+1,Kmm) - ssh(ji,jj ,Kmm) = 0, it won't happen ! here 1033 zcpy(ji,jj) = ABS( (ssh(ji,jj+1,Kmm) + ht_0(ji,jj+1) - ssh(ji,jj,Kmm) - ht_0(ji,jj)) & 1034 & / (ssh(ji,jj+1,Kmm) - ssh(ji,jj ,Kmm)) ) 1035 zcpy(ji,jj) = MAX(MIN( zcpy(ji,jj) , 1.0_wp),0.0_wp) 1041 1036 ELSE 1042 1037 zcpy(ji,jj) = 0._wp 1043 1038 ENDIF 1044 1039 END_2D 1045 CALL lbc_lnk( 'dynhpg', zcpx, 'U', 1.0_wp, zcpy, 'V', 1.0_wp )1046 1040 ENDIF 1047 1041 1048 1042 ! Clean 3-D work arrays 1049 1043 zhpi(:,:,:) = 0._wp 1050 zrhh(:,:,:) = rhd( :,:,:)1044 zrhh(:,:,:) = rhd(A2D(nn_hls),:) 1051 1045 1052 1046 ! Preparing vertical density profile "zrhh(:,:,:)" for hybrid-sco coordinate 1053 1047 DO_2D( 1, 1, 1, 1 ) 1054 jk = mbkt(ji,jj)1055 IF( jk <= 1 ) THEN ; zrhh(ji,jj, : ) = 0._wp1056 ELSEIF( jk == 2 ) THEN ; zrhh(ji,jj,jk+1:jpk) = rhd(ji,jj,jk)1057 ELSEIF( jk < jpkm1 ) THEN1058 DO jkk = jk+1, jpk1059 zrhh(ji,jj,jkk) = interp1(gde3w(ji,jj,jkk ), gde3w(ji,jj,jkk-1), &1060 & gde3w(ji,jj,jkk-2), zrhh (ji,jj,jkk-1), zrhh(ji,jj,jkk-2))1061 END DO1062 ENDIF1048 jk = mbkt(ji,jj) 1049 IF( jk <= 1 ) THEN ; zrhh(ji,jj, : ) = 0._wp 1050 ELSEIF( jk == 2 ) THEN ; zrhh(ji,jj,jk+1:jpk) = rhd(ji,jj,jk) 1051 ELSEIF( jk < jpkm1 ) THEN 1052 DO jkk = jk+1, jpk 1053 zrhh(ji,jj,jkk) = interp1(gde3w(ji,jj,jkk ), gde3w(ji,jj,jkk-1), & 1054 & gde3w(ji,jj,jkk-2), zrhh (ji,jj,jkk-1), zrhh(ji,jj,jkk-2)) 1055 END DO 1056 ENDIF 1063 1057 END_2D 1064 1058 … … 1082 1076 ! Integrate the hydrostatic pressure "zhpi(:,:,:)" at "T(ji,jj,1)" 1083 1077 DO_2D( 0, 1, 0, 1 ) 1084 zrhdt1 = zrhh(ji,jj,1) - interp3( zdept(ji,jj,1), asp(ji,jj,1), bsp(ji,jj,1), &1085 & csp(ji,jj,1), dsp(ji,jj,1) ) * 0.25_wp * e3w(ji,jj,1,Kmm)1086 1087 ! assuming linear profile across the top half surface layer1088 zhpi(ji,jj,1) = 0.5_wp * e3w(ji,jj,1,Kmm) * zrhdt11078 zrhdt1 = zrhh(ji,jj,1) - interp3( zdept(ji,jj,1), asp(ji,jj,1), bsp(ji,jj,1), & 1079 & csp(ji,jj,1), dsp(ji,jj,1) ) * 0.25_wp * e3w(ji,jj,1,Kmm) 1080 1081 ! assuming linear profile across the top half surface layer 1082 zhpi(ji,jj,1) = 0.5_wp * e3w(ji,jj,1,Kmm) * zrhdt1 1089 1083 END_2D 1090 1084 1091 1085 ! Calculate the pressure "zhpi(:,:,:)" at "T(ji,jj,2:jpkm1)" 1092 1086 DO_3D( 0, 1, 0, 1, 2, jpkm1 ) 1093 zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + &1094 & integ_spline( zdept(ji,jj,jk-1), zdept(ji,jj,jk), &1095 & asp (ji,jj,jk-1), bsp (ji,jj,jk-1), &1096 & csp (ji,jj,jk-1), dsp (ji,jj,jk-1) )1087 zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + & 1088 & integ_spline( zdept(ji,jj,jk-1), zdept(ji,jj,jk), & 1089 & asp (ji,jj,jk-1), bsp (ji,jj,jk-1), & 1090 & csp (ji,jj,jk-1), dsp (ji,jj,jk-1) ) 1097 1091 END_3D 1098 1092 … … 1107 1101 ! & r1_e1e2v(ji,jj) * vmask(ji,jj,1) * 0.5_wp 1108 1102 !!gm not this: 1109 zsshu_n(ji,jj) = (e1e2u(ji,jj) * ssh(ji,jj,Kmm) + e1e2u(ji+1, jj) * ssh(ji+1,jj,Kmm)) * & 1110 & r1_e1e2u(ji,jj) * umask(ji,jj,1) * 0.5_wp 1111 zsshv_n(ji,jj) = (e1e2v(ji,jj) * ssh(ji,jj,Kmm) + e1e2v(ji+1, jj) * ssh(ji,jj+1,Kmm)) * & 1112 & r1_e1e2v(ji,jj) * vmask(ji,jj,1) * 0.5_wp 1113 END_2D 1114 1115 CALL lbc_lnk ('dynhpg', zsshu_n, 'U', 1.0_wp, zsshv_n, 'V', 1.0_wp ) 1103 zsshu_n(ji,jj) = (e1e2u(ji,jj) * ssh(ji,jj,Kmm) + e1e2u(ji+1, jj) * ssh(ji+1,jj,Kmm)) * & 1104 & r1_e1e2u(ji,jj) * umask(ji,jj,1) * 0.5_wp 1105 zsshv_n(ji,jj) = (e1e2v(ji,jj) * ssh(ji,jj,Kmm) + e1e2v(ji+1, jj) * ssh(ji,jj+1,Kmm)) * & 1106 & r1_e1e2v(ji,jj) * vmask(ji,jj,1) * 0.5_wp 1107 END_2D 1116 1108 1117 1109 DO_2D( 0, 0, 0, 0 ) 1118 zu(ji,jj,1) = - ( e3u(ji,jj,1,Kmm) - zsshu_n(ji,jj) )1119 zv(ji,jj,1) = - ( e3v(ji,jj,1,Kmm) - zsshv_n(ji,jj) )1110 zu(ji,jj,1) = - ( e3u(ji,jj,1,Kmm) - zsshu_n(ji,jj) ) 1111 zv(ji,jj,1) = - ( e3v(ji,jj,1,Kmm) - zsshv_n(ji,jj) ) 1120 1112 END_2D 1121 1113 1122 1114 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 1123 zu(ji,jj,jk) = zu(ji,jj,jk-1) - e3u(ji,jj,jk,Kmm)1124 zv(ji,jj,jk) = zv(ji,jj,jk-1) - e3v(ji,jj,jk,Kmm)1115 zu(ji,jj,jk) = zu(ji,jj,jk-1) - e3u(ji,jj,jk,Kmm) 1116 zv(ji,jj,jk) = zv(ji,jj,jk-1) - e3v(ji,jj,jk,Kmm) 1125 1117 END_3D 1126 1118 1127 1119 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 1128 zu(ji,jj,jk) = zu(ji,jj,jk) + 0.5_wp * e3u(ji,jj,jk,Kmm)1129 zv(ji,jj,jk) = zv(ji,jj,jk) + 0.5_wp * e3v(ji,jj,jk,Kmm)1120 zu(ji,jj,jk) = zu(ji,jj,jk) + 0.5_wp * e3u(ji,jj,jk,Kmm) 1121 zv(ji,jj,jk) = zv(ji,jj,jk) + 0.5_wp * e3v(ji,jj,jk,Kmm) 1130 1122 END_3D 1131 1123 1132 1124 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 1133 zu(ji,jj,jk) = MIN( zu(ji,jj,jk) , MAX( -zdept(ji,jj,jk) , -zdept(ji+1,jj,jk) ) )1134 zu(ji,jj,jk) = MAX( zu(ji,jj,jk) , MIN( -zdept(ji,jj,jk) , -zdept(ji+1,jj,jk) ) )1135 zv(ji,jj,jk) = MIN( zv(ji,jj,jk) , MAX( -zdept(ji,jj,jk) , -zdept(ji,jj+1,jk) ) )1136 zv(ji,jj,jk) = MAX( zv(ji,jj,jk) , MIN( -zdept(ji,jj,jk) , -zdept(ji,jj+1,jk) ) )1125 zu(ji,jj,jk) = MIN( zu(ji,jj,jk) , MAX( -zdept(ji,jj,jk) , -zdept(ji+1,jj,jk) ) ) 1126 zu(ji,jj,jk) = MAX( zu(ji,jj,jk) , MIN( -zdept(ji,jj,jk) , -zdept(ji+1,jj,jk) ) ) 1127 zv(ji,jj,jk) = MIN( zv(ji,jj,jk) , MAX( -zdept(ji,jj,jk) , -zdept(ji,jj+1,jk) ) ) 1128 zv(ji,jj,jk) = MAX( zv(ji,jj,jk) , MIN( -zdept(ji,jj,jk) , -zdept(ji,jj+1,jk) ) ) 1137 1129 END_3D 1138 1130 1139 1131 1140 1132 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 1141 zpwes = 0._wp; zpwed = 0._wp 1142 zpnss = 0._wp; zpnsd = 0._wp 1143 zuijk = zu(ji,jj,jk) 1144 zvijk = zv(ji,jj,jk) 1145 1146 !!!!! for u equation 1147 IF( jk <= mbku(ji,jj) ) THEN 1148 IF( -zdept(ji+1,jj,jk) >= -zdept(ji,jj,jk) ) THEN 1149 jis = ji + 1; jid = ji 1150 ELSE 1151 jis = ji; jid = ji +1 1133 zpwes = 0._wp; zpwed = 0._wp 1134 zpnss = 0._wp; zpnsd = 0._wp 1135 zuijk = zu(ji,jj,jk) 1136 zvijk = zv(ji,jj,jk) 1137 1138 !!!!! for u equation 1139 IF( jk <= mbku(ji,jj) ) THEN 1140 IF( -zdept(ji+1,jj,jk) >= -zdept(ji,jj,jk) ) THEN 1141 jis = ji + 1; jid = ji 1142 ELSE 1143 jis = ji; jid = ji +1 1144 ENDIF 1145 1146 ! integrate the pressure on the shallow side 1147 jk1 = jk 1148 DO WHILE ( -zdept(jis,jj,jk1) > zuijk ) 1149 IF( jk1 == mbku(ji,jj) ) THEN 1150 zuijk = -zdept(jis,jj,jk1) 1151 EXIT 1152 ENDIF 1153 zdeps = MIN(zdept(jis,jj,jk1+1), -zuijk) 1154 zpwes = zpwes + & 1155 integ_spline(zdept(jis,jj,jk1), zdeps, & 1156 asp(jis,jj,jk1), bsp(jis,jj,jk1), & 1157 csp(jis,jj,jk1), dsp(jis,jj,jk1)) 1158 jk1 = jk1 + 1 1159 END DO 1160 1161 ! integrate the pressure on the deep side 1162 jk1 = jk 1163 DO WHILE ( -zdept(jid,jj,jk1) < zuijk ) 1164 IF( jk1 == 1 ) THEN 1165 zdeps = zdept(jid,jj,1) + MIN(zuijk, ssh(jid,jj,Kmm)*znad) 1166 zrhdt1 = zrhh(jid,jj,1) - interp3(zdept(jid,jj,1), asp(jid,jj,1), & 1167 bsp(jid,jj,1) , csp(jid,jj,1), & 1168 dsp(jid,jj,1)) * zdeps 1169 zpwed = zpwed + 0.5_wp * (zrhh(jid,jj,1) + zrhdt1) * zdeps 1170 EXIT 1171 ENDIF 1172 zdeps = MAX(zdept(jid,jj,jk1-1), -zuijk) 1173 zpwed = zpwed + & 1174 integ_spline(zdeps, zdept(jid,jj,jk1), & 1175 asp(jid,jj,jk1-1), bsp(jid,jj,jk1-1), & 1176 csp(jid,jj,jk1-1), dsp(jid,jj,jk1-1) ) 1177 jk1 = jk1 - 1 1178 END DO 1179 1180 ! update the momentum trends in u direction 1181 zdpdx1 = zcoef0 * r1_e1u(ji,jj) * ( zhpi(ji+1,jj,jk) - zhpi(ji,jj,jk) ) 1182 IF( .NOT.ln_linssh ) THEN 1183 zdpdx2 = zcoef0 * r1_e1u(ji,jj) * & 1184 & ( REAL(jis-jid, wp) * (zpwes + zpwed) + (ssh(ji+1,jj,Kmm)-ssh(ji,jj,Kmm)) ) 1185 ELSE 1186 zdpdx2 = zcoef0 * r1_e1u(ji,jj) * REAL(jis-jid, wp) * (zpwes + zpwed) 1187 ENDIF 1188 IF( ln_wd_il ) THEN 1189 zdpdx1 = zdpdx1 * zcpx(ji,jj) * wdrampu(ji,jj) 1190 zdpdx2 = zdpdx2 * zcpx(ji,jj) * wdrampu(ji,jj) 1191 ENDIF 1192 puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + (zdpdx1 + zdpdx2 - zpgu(ji,jj)) * umask(ji,jj,jk) 1152 1193 ENDIF 1153 1194 1154 ! integrate the pressure on the shallow side 1155 jk1 = jk 1156 DO WHILE ( -zdept(jis,jj,jk1) > zuijk ) 1157 IF( jk1 == mbku(ji,jj) ) THEN 1158 zuijk = -zdept(jis,jj,jk1) 1159 EXIT 1160 ENDIF 1161 zdeps = MIN(zdept(jis,jj,jk1+1), -zuijk) 1162 zpwes = zpwes + & 1163 integ_spline(zdept(jis,jj,jk1), zdeps, & 1164 asp(jis,jj,jk1), bsp(jis,jj,jk1), & 1165 csp(jis,jj,jk1), dsp(jis,jj,jk1)) 1166 jk1 = jk1 + 1 1167 END DO 1168 1169 ! integrate the pressure on the deep side 1170 jk1 = jk 1171 DO WHILE ( -zdept(jid,jj,jk1) < zuijk ) 1172 IF( jk1 == 1 ) THEN 1173 zdeps = zdept(jid,jj,1) + MIN(zuijk, ssh(jid,jj,Kmm)*znad) 1174 zrhdt1 = zrhh(jid,jj,1) - interp3(zdept(jid,jj,1), asp(jid,jj,1), & 1175 bsp(jid,jj,1), csp(jid,jj,1), & 1176 dsp(jid,jj,1)) * zdeps 1177 zpwed = zpwed + 0.5_wp * (zrhh(jid,jj,1) + zrhdt1) * zdeps 1178 EXIT 1179 ENDIF 1180 zdeps = MAX(zdept(jid,jj,jk1-1), -zuijk) 1181 zpwed = zpwed + & 1182 integ_spline(zdeps, zdept(jid,jj,jk1), & 1183 asp(jid,jj,jk1-1), bsp(jid,jj,jk1-1), & 1184 csp(jid,jj,jk1-1), dsp(jid,jj,jk1-1) ) 1185 jk1 = jk1 - 1 1186 END DO 1187 1188 ! update the momentum trends in u direction 1189 1190 zdpdx1 = zcoef0 * r1_e1u(ji,jj) * ( zhpi(ji+1,jj,jk) - zhpi(ji,jj,jk) ) 1191 IF( .NOT.ln_linssh ) THEN 1192 zdpdx2 = zcoef0 * r1_e1u(ji,jj) * & 1193 & ( REAL(jis-jid, wp) * (zpwes + zpwed) + (ssh(ji+1,jj,Kmm)-ssh(ji,jj,Kmm)) ) 1194 ELSE 1195 zdpdx2 = zcoef0 * r1_e1u(ji,jj) * REAL(jis-jid, wp) * (zpwes + zpwed) 1195 !!!!! for v equation 1196 IF( jk <= mbkv(ji,jj) ) THEN 1197 IF( -zdept(ji,jj+1,jk) >= -zdept(ji,jj,jk) ) THEN 1198 jjs = jj + 1; jjd = jj 1199 ELSE 1200 jjs = jj ; jjd = jj + 1 1201 ENDIF 1202 1203 ! integrate the pressure on the shallow side 1204 jk1 = jk 1205 DO WHILE ( -zdept(ji,jjs,jk1) > zvijk ) 1206 IF( jk1 == mbkv(ji,jj) ) THEN 1207 zvijk = -zdept(ji,jjs,jk1) 1208 EXIT 1209 ENDIF 1210 zdeps = MIN(zdept(ji,jjs,jk1+1), -zvijk) 1211 zpnss = zpnss + & 1212 integ_spline(zdept(ji,jjs,jk1), zdeps, & 1213 asp(ji,jjs,jk1), bsp(ji,jjs,jk1), & 1214 csp(ji,jjs,jk1), dsp(ji,jjs,jk1) ) 1215 jk1 = jk1 + 1 1216 END DO 1217 1218 ! integrate the pressure on the deep side 1219 jk1 = jk 1220 DO WHILE ( -zdept(ji,jjd,jk1) < zvijk ) 1221 IF( jk1 == 1 ) THEN 1222 zdeps = zdept(ji,jjd,1) + MIN(zvijk, ssh(ji,jjd,Kmm)*znad) 1223 zrhdt1 = zrhh(ji,jjd,1) - interp3(zdept(ji,jjd,1), asp(ji,jjd,1), & 1224 bsp(ji,jjd,1) , csp(ji,jjd,1), & 1225 dsp(ji,jjd,1) ) * zdeps 1226 zpnsd = zpnsd + 0.5_wp * (zrhh(ji,jjd,1) + zrhdt1) * zdeps 1227 EXIT 1228 ENDIF 1229 zdeps = MAX(zdept(ji,jjd,jk1-1), -zvijk) 1230 zpnsd = zpnsd + & 1231 integ_spline(zdeps, zdept(ji,jjd,jk1), & 1232 asp(ji,jjd,jk1-1), bsp(ji,jjd,jk1-1), & 1233 csp(ji,jjd,jk1-1), dsp(ji,jjd,jk1-1) ) 1234 jk1 = jk1 - 1 1235 END DO 1236 1237 ! update the momentum trends in v direction 1238 zdpdy1 = zcoef0 * r1_e2v(ji,jj) * ( zhpi(ji,jj+1,jk) - zhpi(ji,jj,jk) ) 1239 IF( .NOT.ln_linssh ) THEN 1240 zdpdy2 = zcoef0 * r1_e2v(ji,jj) * & 1241 ( REAL(jjs-jjd, wp) * (zpnss + zpnsd) + (ssh(ji,jj+1,Kmm)-ssh(ji,jj,Kmm)) ) 1242 ELSE 1243 zdpdy2 = zcoef0 * r1_e2v(ji,jj) * REAL(jjs-jjd, wp) * (zpnss + zpnsd ) 1244 ENDIF 1245 IF( ln_wd_il ) THEN 1246 zdpdy1 = zdpdy1 * zcpy(ji,jj) * wdrampv(ji,jj) 1247 zdpdy2 = zdpdy2 * zcpy(ji,jj) * wdrampv(ji,jj) 1248 ENDIF 1249 1250 pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + (zdpdy1 + zdpdy2 - zpgv(ji,jj)) * vmask(ji,jj,jk) 1196 1251 ENDIF 1197 IF( ln_wd_il ) THEN1198 zdpdx1 = zdpdx1 * zcpx(ji,jj) * wdrampu(ji,jj)1199 zdpdx2 = zdpdx2 * zcpx(ji,jj) * wdrampu(ji,jj)1200 ENDIF1201 puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + (zdpdx1 + zdpdx2 - zpgu(ji,jj)) * umask(ji,jj,jk)1202 ENDIF1203 1204 !!!!! for v equation1205 IF( jk <= mbkv(ji,jj) ) THEN1206 IF( -zdept(ji,jj+1,jk) >= -zdept(ji,jj,jk) ) THEN1207 jjs = jj + 1; jjd = jj1208 ELSE1209 jjs = jj ; jjd = jj + 11210 ENDIF1211 1212 ! integrate the pressure on the shallow side1213 jk1 = jk1214 DO WHILE ( -zdept(ji,jjs,jk1) > zvijk )1215 IF( jk1 == mbkv(ji,jj) ) THEN1216 zvijk = -zdept(ji,jjs,jk1)1217 EXIT1218 ENDIF1219 zdeps = MIN(zdept(ji,jjs,jk1+1), -zvijk)1220 zpnss = zpnss + &1221 integ_spline(zdept(ji,jjs,jk1), zdeps, &1222 asp(ji,jjs,jk1), bsp(ji,jjs,jk1), &1223 csp(ji,jjs,jk1), dsp(ji,jjs,jk1) )1224 jk1 = jk1 + 11225 END DO1226 1227 ! integrate the pressure on the deep side1228 jk1 = jk1229 DO WHILE ( -zdept(ji,jjd,jk1) < zvijk )1230 IF( jk1 == 1 ) THEN1231 zdeps = zdept(ji,jjd,1) + MIN(zvijk, ssh(ji,jjd,Kmm)*znad)1232 zrhdt1 = zrhh(ji,jjd,1) - interp3(zdept(ji,jjd,1), asp(ji,jjd,1), &1233 bsp(ji,jjd,1), csp(ji,jjd,1), &1234 dsp(ji,jjd,1) ) * zdeps1235 zpnsd = zpnsd + 0.5_wp * (zrhh(ji,jjd,1) + zrhdt1) * zdeps1236 EXIT1237 ENDIF1238 zdeps = MAX(zdept(ji,jjd,jk1-1), -zvijk)1239 zpnsd = zpnsd + &1240 integ_spline(zdeps, zdept(ji,jjd,jk1), &1241 asp(ji,jjd,jk1-1), bsp(ji,jjd,jk1-1), &1242 csp(ji,jjd,jk1-1), dsp(ji,jjd,jk1-1) )1243 jk1 = jk1 - 11244 END DO1245 1246 1247 ! update the momentum trends in v direction1248 1249 zdpdy1 = zcoef0 * r1_e2v(ji,jj) * ( zhpi(ji,jj+1,jk) - zhpi(ji,jj,jk) )1250 IF( .NOT.ln_linssh ) THEN1251 zdpdy2 = zcoef0 * r1_e2v(ji,jj) * &1252 ( REAL(jjs-jjd, wp) * (zpnss + zpnsd) + (ssh(ji,jj+1,Kmm)-ssh(ji,jj,Kmm)) )1253 ELSE1254 zdpdy2 = zcoef0 * r1_e2v(ji,jj) * REAL(jjs-jjd, wp) * (zpnss + zpnsd )1255 ENDIF1256 IF( ln_wd_il ) THEN1257 zdpdy1 = zdpdy1 * zcpy(ji,jj) * wdrampv(ji,jj)1258 zdpdy2 = zdpdy2 * zcpy(ji,jj) * wdrampv(ji,jj)1259 ENDIF1260 1261 pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + (zdpdy1 + zdpdy2 - zpgv(ji,jj)) * vmask(ji,jj,jk)1262 ENDIF1263 1252 ! 1264 1253 END_3D … … 1279 1268 !! Reference: CJC Kruger, Constrained Cubic Spline Interpoltation 1280 1269 !!---------------------------------------------------------------------- 1281 REAL(wp), DIMENSION( :,:,:), INTENT(in ) :: fsp, xsp ! value and coordinate1282 REAL(wp), DIMENSION( :,:,:), INTENT( out) :: asp, bsp, csp, dsp ! coefficients of the interpoated function1283 INTEGER , INTENT(in ) :: polynomial_type ! 1: cubic spline ; 2: Linear1270 REAL(wp), DIMENSION(A2D(nn_hls),jpk), INTENT(in ) :: fsp, xsp ! value and coordinate 1271 REAL(wp), DIMENSION(A2D(nn_hls),jpk), INTENT( out) :: asp, bsp, csp, dsp ! coefficients of the interpoated function 1272 INTEGER , INTENT(in ) :: polynomial_type ! 1: cubic spline ; 2: Linear 1284 1273 ! 1285 1274 INTEGER :: ji, jj, jk ! dummy loop indices 1286 INTEGER :: jpi, jpj, jpkm11287 1275 REAL(wp) :: zdf1, zdf2, zddf1, zddf2, ztmp1, ztmp2, zdxtmp 1288 1276 REAL(wp) :: zdxtmp1, zdxtmp2, zalpha 1289 REAL(wp) :: zdf(size(fsp,3)) 1290 !!---------------------------------------------------------------------- 1291 ! 1292 !!gm WHAT !!!!! THIS IS VERY DANGEROUS !!!!! 1293 jpi = size(fsp,1) 1294 jpj = size(fsp,2) 1295 jpkm1 = MAX( 1, size(fsp,3) - 1 ) 1277 REAL(wp) :: zdf(jpk) 1278 !!---------------------------------------------------------------------- 1296 1279 ! 1297 1280 IF (polynomial_type == 1) THEN ! Constrained Cubic Spline 1298 DO ji = 1, jpi 1299 DO jj = 1, jpj 1300 !!Fritsch&Butland's method, 1984 (preferred, but more computation) 1301 ! DO jk = 2, jpkm1-1 1302 ! zdxtmp1 = xsp(ji,jj,jk) - xsp(ji,jj,jk-1) 1303 ! zdxtmp2 = xsp(ji,jj,jk+1) - xsp(ji,jj,jk) 1304 ! zdf1 = ( fsp(ji,jj,jk) - fsp(ji,jj,jk-1) ) / zdxtmp1 1305 ! zdf2 = ( fsp(ji,jj,jk+1) - fsp(ji,jj,jk) ) / zdxtmp2 1306 ! 1307 ! zalpha = ( zdxtmp1 + 2._wp * zdxtmp2 ) / ( zdxtmp1 + zdxtmp2 ) / 3._wp 1308 ! 1309 ! IF(zdf1 * zdf2 <= 0._wp) THEN 1310 ! zdf(jk) = 0._wp 1311 ! ELSE 1312 ! zdf(jk) = zdf1 * zdf2 / ( ( 1._wp - zalpha ) * zdf1 + zalpha * zdf2 ) 1313 ! ENDIF 1314 ! END DO 1315 1316 !!Simply geometric average 1317 DO jk = 2, jpkm1-1 1318 zdf1 = (fsp(ji,jj,jk ) - fsp(ji,jj,jk-1)) / (xsp(ji,jj,jk ) - xsp(ji,jj,jk-1)) 1319 zdf2 = (fsp(ji,jj,jk+1) - fsp(ji,jj,jk )) / (xsp(ji,jj,jk+1) - xsp(ji,jj,jk )) 1320 1321 IF(zdf1 * zdf2 <= 0._wp) THEN 1322 zdf(jk) = 0._wp 1323 ELSE 1324 zdf(jk) = 2._wp * zdf1 * zdf2 / (zdf1 + zdf2) 1325 ENDIF 1326 END DO 1327 1328 zdf(1) = 1.5_wp * ( fsp(ji,jj,2) - fsp(ji,jj,1) ) / & 1329 & ( xsp(ji,jj,2) - xsp(ji,jj,1) ) - 0.5_wp * zdf(2) 1330 zdf(jpkm1) = 1.5_wp * ( fsp(ji,jj,jpkm1) - fsp(ji,jj,jpkm1-1) ) / & 1331 & ( xsp(ji,jj,jpkm1) - xsp(ji,jj,jpkm1-1) ) - 0.5_wp * zdf(jpkm1 - 1) 1332 1333 DO jk = 1, jpkm1 - 1 1334 zdxtmp = xsp(ji,jj,jk+1) - xsp(ji,jj,jk) 1335 ztmp1 = (zdf(jk+1) + 2._wp * zdf(jk)) / zdxtmp 1336 ztmp2 = 6._wp * (fsp(ji,jj,jk+1) - fsp(ji,jj,jk)) / zdxtmp / zdxtmp 1337 zddf1 = -2._wp * ztmp1 + ztmp2 1338 ztmp1 = (2._wp * zdf(jk+1) + zdf(jk)) / zdxtmp 1339 zddf2 = 2._wp * ztmp1 - ztmp2 1340 1341 dsp(ji,jj,jk) = (zddf2 - zddf1) / 6._wp / zdxtmp 1342 csp(ji,jj,jk) = ( xsp(ji,jj,jk+1) * zddf1 - xsp(ji,jj,jk)*zddf2 ) / 2._wp / zdxtmp 1343 bsp(ji,jj,jk) = ( fsp(ji,jj,jk+1) - fsp(ji,jj,jk) ) / zdxtmp - & 1344 & csp(ji,jj,jk) * ( xsp(ji,jj,jk+1) + xsp(ji,jj,jk) ) - & 1345 & dsp(ji,jj,jk) * ((xsp(ji,jj,jk+1) + xsp(ji,jj,jk))**2 - & 1346 & xsp(ji,jj,jk+1) * xsp(ji,jj,jk)) 1347 asp(ji,jj,jk) = fsp(ji,jj,jk) - xsp(ji,jj,jk) * (bsp(ji,jj,jk) + & 1348 & (xsp(ji,jj,jk) * (csp(ji,jj,jk) + & 1349 & dsp(ji,jj,jk) * xsp(ji,jj,jk)))) 1350 END DO 1281 DO_2D( 1, 1, 1, 1 ) 1282 !!Fritsch&Butland's method, 1984 (preferred, but more computation) 1283 ! DO jk = 2, jpkm1-1 1284 ! zdxtmp1 = xsp(ji,jj,jk) - xsp(ji,jj,jk-1) 1285 ! zdxtmp2 = xsp(ji,jj,jk+1) - xsp(ji,jj,jk) 1286 ! zdf1 = ( fsp(ji,jj,jk) - fsp(ji,jj,jk-1) ) / zdxtmp1 1287 ! zdf2 = ( fsp(ji,jj,jk+1) - fsp(ji,jj,jk) ) / zdxtmp2 1288 ! 1289 ! zalpha = ( zdxtmp1 + 2._wp * zdxtmp2 ) / ( zdxtmp1 + zdxtmp2 ) / 3._wp 1290 ! 1291 ! IF(zdf1 * zdf2 <= 0._wp) THEN 1292 ! zdf(jk) = 0._wp 1293 ! ELSE 1294 ! zdf(jk) = zdf1 * zdf2 / ( ( 1._wp - zalpha ) * zdf1 + zalpha * zdf2 ) 1295 ! ENDIF 1296 ! END DO 1297 1298 !!Simply geometric average 1299 DO jk = 2, jpk-2 1300 zdf1 = (fsp(ji,jj,jk ) - fsp(ji,jj,jk-1)) / (xsp(ji,jj,jk ) - xsp(ji,jj,jk-1)) 1301 zdf2 = (fsp(ji,jj,jk+1) - fsp(ji,jj,jk )) / (xsp(ji,jj,jk+1) - xsp(ji,jj,jk )) 1302 1303 IF(zdf1 * zdf2 <= 0._wp) THEN 1304 zdf(jk) = 0._wp 1305 ELSE 1306 zdf(jk) = 2._wp * zdf1 * zdf2 / (zdf1 + zdf2) 1307 ENDIF 1351 1308 END DO 1352 END DO 1309 1310 zdf(1) = 1.5_wp * ( fsp(ji,jj,2) - fsp(ji,jj,1) ) / & 1311 & ( xsp(ji,jj,2) - xsp(ji,jj,1) ) - 0.5_wp * zdf(2) 1312 zdf(jpkm1) = 1.5_wp * ( fsp(ji,jj,jpkm1) - fsp(ji,jj,jpkm1-1) ) / & 1313 & ( xsp(ji,jj,jpkm1) - xsp(ji,jj,jpkm1-1) ) - 0.5_wp * zdf(jpk - 2) 1314 1315 DO jk = 1, jpk-2 1316 zdxtmp = xsp(ji,jj,jk+1) - xsp(ji,jj,jk) 1317 ztmp1 = (zdf(jk+1) + 2._wp * zdf(jk)) / zdxtmp 1318 ztmp2 = 6._wp * (fsp(ji,jj,jk+1) - fsp(ji,jj,jk)) / zdxtmp / zdxtmp 1319 zddf1 = -2._wp * ztmp1 + ztmp2 1320 ztmp1 = (2._wp * zdf(jk+1) + zdf(jk)) / zdxtmp 1321 zddf2 = 2._wp * ztmp1 - ztmp2 1322 1323 dsp(ji,jj,jk) = (zddf2 - zddf1) / 6._wp / zdxtmp 1324 csp(ji,jj,jk) = ( xsp(ji,jj,jk+1) * zddf1 - xsp(ji,jj,jk)*zddf2 ) / 2._wp / zdxtmp 1325 bsp(ji,jj,jk) = ( fsp(ji,jj,jk+1) - fsp(ji,jj,jk) ) / zdxtmp - & 1326 & csp(ji,jj,jk) * ( xsp(ji,jj,jk+1) + xsp(ji,jj,jk) ) - & 1327 & dsp(ji,jj,jk) * ((xsp(ji,jj,jk+1) + xsp(ji,jj,jk))**2 - & 1328 & xsp(ji,jj,jk+1) * xsp(ji,jj,jk)) 1329 asp(ji,jj,jk) = fsp(ji,jj,jk) - xsp(ji,jj,jk) * (bsp(ji,jj,jk) + & 1330 & (xsp(ji,jj,jk) * (csp(ji,jj,jk) + & 1331 & dsp(ji,jj,jk) * xsp(ji,jj,jk)))) 1332 END DO 1333 END_2D 1353 1334 1354 1335 ELSEIF ( polynomial_type == 2 ) THEN ! Linear 1355 DO ji = 1, jpi 1356 DO jj = 1, jpj 1357 DO jk = 1, jpkm1-1 1358 zdxtmp =xsp(ji,jj,jk+1) - xsp(ji,jj,jk) 1359 ztmp1 = fsp(ji,jj,jk+1) - fsp(ji,jj,jk) 1360 1361 dsp(ji,jj,jk) = 0._wp 1362 csp(ji,jj,jk) = 0._wp 1363 bsp(ji,jj,jk) = ztmp1 / zdxtmp 1364 asp(ji,jj,jk) = fsp(ji,jj,jk) - bsp(ji,jj,jk) * xsp(ji,jj,jk) 1365 END DO 1366 END DO 1367 END DO 1336 DO_3D( 1, 1, 1, 1, 1, jpk-2 ) 1337 zdxtmp =xsp(ji,jj,jk+1) - xsp(ji,jj,jk) 1338 ztmp1 = fsp(ji,jj,jk+1) - fsp(ji,jj,jk) 1339 1340 dsp(ji,jj,jk) = 0._wp 1341 csp(ji,jj,jk) = 0._wp 1342 bsp(ji,jj,jk) = ztmp1 / zdxtmp 1343 asp(ji,jj,jk) = fsp(ji,jj,jk) - bsp(ji,jj,jk) * xsp(ji,jj,jk) 1344 END_3D 1368 1345 ! 1369 1346 ELSE -
NEMO/trunk/src/OCE/DYN/dynkeg.F90
r14820 r14834 78 78 INTEGER :: ji, jj, jk ! dummy loop indices 79 79 REAL(wp) :: zu, zv ! local scalars 80 REAL(wp), DIMENSION( jpi,jpj,jpk):: zhke80 REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zhke 81 81 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: ztrdu, ztrdv 82 82 !!---------------------------------------------------------------------- … … 84 84 IF( ln_timing ) CALL timing_start('dyn_keg') 85 85 ! 86 IF( kt == nit000 ) THEN 87 IF(lwp) WRITE(numout,*) 88 IF(lwp) WRITE(numout,*) 'dyn_keg : kinetic energy gradient trend, scheme number=', kscheme 89 IF(lwp) WRITE(numout,*) '~~~~~~~' 86 IF( .NOT. l_istiled .OR. ntile == 1 ) THEN ! Do only on the first tile 87 IF( kt == nit000 ) THEN 88 IF(lwp) WRITE(numout,*) 89 IF(lwp) WRITE(numout,*) 'dyn_keg : kinetic energy gradient trend, scheme number=', kscheme 90 IF(lwp) WRITE(numout,*) '~~~~~~~' 91 ENDIF 90 92 ENDIF 91 93 … … 109 111 END_3D 110 112 CASE ( nkeg_HW ) !-- Hollingsworth scheme --! 111 DO_3D( 0, 0, 0, 0, 1, jpkm1 )113 DO_3D( 0, nn_hls-1, 0, nn_hls-1, 1, jpkm1 ) 112 114 ! round brackets added to fix the order of floating point operations 113 115 ! needed to ensure halo 1 - halo 2 compatibility … … 121 123 & + pvv(ji ,jj ,jk,Kmm) * pvv(ji ,jj ,jk,Kmm) ) & 122 124 & + ( ( pvv(ji-1,jj-1,jk,Kmm) + pvv(ji+1,jj-1,jk,Kmm) ) * ( pvv(ji-1,jj-1,jk,Kmm) + pvv(ji+1,jj-1,jk,Kmm) ) & 123 & + ( pvv(ji-1,jj ,jk,Kmm) + pvv(ji+1,jj ,jk,Kmm) ) * ( pvv(ji-1,jj ,jk,Kmm) + pvv(ji+1,jj ,jk,Kmm) ) & 125 & + ( pvv(ji-1,jj ,jk,Kmm) + pvv(ji+1,jj ,jk,Kmm) ) * ( pvv(ji-1,jj ,jk,Kmm) + pvv(ji+1,jj ,jk,Kmm) ) & 124 126 & ) ! bracket for halo 1 - halo 2 compatibility 125 127 zhke(ji,jj,jk) = r1_48 * ( zv + zu ) 126 128 END_3D 127 CALL lbc_lnk( 'dynkeg', zhke, 'T', 1.0_wp )129 IF (nn_hls==1) CALL lbc_lnk( 'dynkeg', zhke, 'T', 1.0_wp ) 128 130 ! 129 131 END SELECT -
NEMO/trunk/src/OCE/DYN/dynldf_iso.F90
r14433 r14834 28 28 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 29 29 USE prtctl ! Print control 30 #if defined key_loop_fusion 31 USE dynldf_iso_lf, ONLY: dyn_ldf_iso_lf ! lateral mixing - loop fusion version (dyn_ldf_iso routine ) 32 #endif 30 33 31 34 IMPLICIT NONE … … 36 39 37 40 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: akzu, akzv !: vertical component of rotated lateral viscosity 38 39 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: zfuw, zdiu, zdju, zdj1u ! 2D workspace (dyn_ldf_iso)40 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: zfvw, zdiv, zdjv, zdj1v ! - -41 41 42 42 !! * Substitutions … … 54 54 !! *** ROUTINE dyn_ldf_iso_alloc *** 55 55 !!---------------------------------------------------------------------- 56 ALLOCATE( akzu(jpi,jpj,jpk) , zfuw(jpi,jpk) , zdiu(jpi,jpk) , zdju(jpi,jpk) , zdj1u(jpi,jpk) , & 57 & akzv(jpi,jpj,jpk) , zfvw(jpi,jpk) , zdiv(jpi,jpk) , zdjv(jpi,jpk) , zdj1v(jpi,jpk) , STAT=dyn_ldf_iso_alloc ) 58 ! 59 IF( dyn_ldf_iso_alloc /= 0 ) CALL ctl_warn('dyn_ldf_iso_alloc: array allocate failed.') 56 dyn_ldf_iso_alloc = 0 57 IF( .NOT. ALLOCATED( akzu ) ) THEN 58 ALLOCATE( akzu(jpi,jpj,jpk), akzv(jpi,jpj,jpk), STAT=dyn_ldf_iso_alloc ) 59 ! 60 IF( dyn_ldf_iso_alloc /= 0 ) CALL ctl_warn('dyn_ldf_iso_alloc: array allocate failed.') 61 ENDIF 60 62 END FUNCTION dyn_ldf_iso_alloc 61 63 … … 112 114 REAL(wp) :: zabe2, zmskf, zmkf, zvav, zvwslpi, zvwslpj ! - - 113 115 REAL(wp) :: zcof0, zcof1, zcof2, zcof3, zcof4, zaht_0 ! - - 114 REAL(wp), DIMENSION(jpi,jpj) :: ziut, zivf, zdku, zdk1u ! 2D workspace 115 REAL(wp), DIMENSION(jpi,jpj) :: zjuf, zjvt, zdkv, zdk1v ! - - 116 REAL(wp), DIMENSION(A2D(nn_hls)) :: ziut, zivf, zdku, zdk1u ! 2D workspace 117 REAL(wp), DIMENSION(A2D(nn_hls)) :: zjuf, zjvt, zdkv, zdk1v ! - - 118 REAL(wp), DIMENSION(A1Di(nn_hls),jpk) :: zfuw, zdiu, zdju, zdj1u ! - - 119 REAL(wp), DIMENSION(A1Di(nn_hls),jpk) :: zfvw, zdiv, zdjv, zdj1v ! - - 116 120 !!---------------------------------------------------------------------- 117 121 ! 118 IF( kt == nit000 ) THEN 119 IF(lwp) WRITE(numout,*) 120 IF(lwp) WRITE(numout,*) 'dyn_ldf_iso : iso-neutral laplacian diffusive operator or ' 121 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ s-coordinate horizontal diffusive operator' 122 ! ! allocate dyn_ldf_bilap arrays 123 IF( dyn_ldf_iso_alloc() /= 0 ) CALL ctl_stop('STOP', 'dyn_ldf_iso: failed to allocate arrays') 122 #if defined key_loop_fusion 123 CALL dyn_ldf_iso_lf( kt, Kbb, Kmm, puu, pvv, Krhs ) 124 #else 125 126 IF( .NOT. l_istiled .OR. ntile == 1 ) THEN ! Do only on the first tile 127 IF( kt == nit000 ) THEN 128 IF(lwp) WRITE(numout,*) 129 IF(lwp) WRITE(numout,*) 'dyn_ldf_iso : iso-neutral laplacian diffusive operator or ' 130 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ s-coordinate horizontal diffusive operator' 131 ! ! allocate dyn_ldf_iso arrays 132 IF( dyn_ldf_iso_alloc() /= 0 ) CALL ctl_stop('STOP', 'dyn_ldf_iso: failed to allocate arrays') 133 ENDIF 124 134 ENDIF 125 135 … … 128 138 IF( ln_dynldf_hor .AND. ln_traldf_iso ) THEN 129 139 ! 130 DO_3D ( 0, 0, 0, 0, 1, jpk ) ! set the slopes of iso-level140 DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpk ) ! set the slopes of iso-level 131 141 uslp (ji,jj,jk) = - ( gdept(ji+1,jj,jk,Kbb) - gdept(ji ,jj ,jk,Kbb) ) * r1_e1u(ji,jj) * umask(ji,jj,jk) 132 142 vslp (ji,jj,jk) = - ( gdept(ji,jj+1,jk,Kbb) - gdept(ji ,jj ,jk,Kbb) ) * r1_e2v(ji,jj) * vmask(ji,jj,jk) … … 135 145 END_3D 136 146 ! Lateral boundary conditions on the slopes 137 CALL lbc_lnk( 'dynldf_iso', uslp , 'U', -1.0_wp, vslp , 'V', -1.0_wp, wslpi, 'W', -1.0_wp, wslpj, 'W', -1.0_wp )147 IF (nn_hls == 1) CALL lbc_lnk( 'dynldf_iso', uslp , 'U', -1.0_wp, vslp , 'V', -1.0_wp, wslpi, 'W', -1.0_wp, wslpj, 'W', -1.0_wp ) 138 148 ! 139 149 ENDIF 140 150 141 151 zaht_0 = 0.5_wp * rn_Ud * rn_Ld ! aht_0 from namtra_ldf = zaht_max … … 150 160 ! zdkv(jk=1)=zdkv(jk=2) 151 161 152 zdk1u(:,:) = ( puu(:,:,jk,Kbb) -puu(:,:,jk+1,Kbb) ) * umask(:,:,jk+1) 153 zdk1v(:,:) = ( pvv(:,:,jk,Kbb) -pvv(:,:,jk+1,Kbb) ) * vmask(:,:,jk+1) 162 DO_2D( 1, 1, 1, 1 ) 163 zdk1u(ji,jj) = ( puu(ji,jj,jk,Kbb) -puu(ji,jj,jk+1,Kbb) ) * umask(ji,jj,jk+1) 164 zdk1v(ji,jj) = ( pvv(ji,jj,jk,Kbb) -pvv(ji,jj,jk+1,Kbb) ) * vmask(ji,jj,jk+1) 165 END_2D 154 166 155 167 IF( jk == 1 ) THEN … … 157 169 zdkv(:,:) = zdk1v(:,:) 158 170 ELSE 159 zdku(:,:) = ( puu(:,:,jk-1,Kbb) - puu(:,:,jk,Kbb) ) * umask(:,:,jk) 160 zdkv(:,:) = ( pvv(:,:,jk-1,Kbb) - pvv(:,:,jk,Kbb) ) * vmask(:,:,jk) 171 DO_2D( 1, 1, 1, 1 ) 172 zdku(ji,jj) = ( puu(ji,jj,jk-1,Kbb) - puu(ji,jj,jk,Kbb) ) * umask(ji,jj,jk) 173 zdkv(ji,jj) = ( pvv(ji,jj,jk-1,Kbb) - pvv(ji,jj,jk,Kbb) ) * vmask(ji,jj,jk) 174 END_2D 161 175 ENDIF 162 176 … … 286 300 287 301 ! ! =============== 288 DO jj = 2, jpjm1! Vertical slab302 DO jj = ntsj, ntej ! Vertical slab 289 303 ! ! =============== 290 304 … … 299 313 300 314 DO jk = 1, jpk 301 DO ji = 2, jpi315 DO ji = ntsi, ntei + nn_hls 302 316 ! i-gradient of u at jj 303 317 zdiu (ji,jk) = tmask(ji,jj ,jk) * ( puu(ji,jj ,jk,Kbb) - puu(ji-1,jj ,jk,Kbb) ) … … 311 325 END DO 312 326 DO jk = 1, jpk 313 DO ji = 1, jpim1327 DO ji = ntsi - nn_hls, ntei 314 328 ! i-gradient of v at jj 315 329 zdiv (ji,jk) = fmask(ji,jj ,jk) * ( pvv(ji+1,jj,jk,Kbb) - pvv(ji ,jj ,jk,Kbb) ) … … 322 336 323 337 ! Surface and bottom vertical fluxes set to zero 324 DO ji = 1, jpi338 DO ji = ntsi - nn_hls, ntei + nn_hls 325 339 zfuw(ji, 1 ) = 0.e0 326 340 zfvw(ji, 1 ) = 0.e0 … … 331 345 ! interior (2=<jk=<jpk-1) on U field 332 346 DO jk = 2, jpkm1 333 DO ji = 2, jpim1347 DO ji = ntsi, ntei 334 348 zcof0 = 0.5_wp * zaht_0 * umask(ji,jj,jk) 335 349 ! … … 357 371 ! interior (2=<jk=<jpk-1) on V field 358 372 DO jk = 2, jpkm1 359 DO ji = 2, jpim1373 DO ji = ntsi, ntei 360 374 zcof0 = 0.5_wp * zaht_0 * vmask(ji,jj,jk) 361 375 ! … … 385 399 ! ------------------------------------------------------------------- 386 400 DO jk = 1, jpkm1 387 DO ji = 2, jpim1401 DO ji = ntsi, ntei 388 402 puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + ( zfuw(ji,jk) - zfuw(ji,jk+1) ) * r1_e1e2u(ji,jj) & 389 403 & / e3u(ji,jj,jk,Kmm) … … 395 409 END DO ! End of slab 396 410 ! ! =============== 411 #endif 397 412 END SUBROUTINE dyn_ldf_iso 398 413 -
NEMO/trunk/src/OCE/DYN/dynldf_lap_blp.F90
r14433 r14834 14 14 USE oce ! ocean dynamics and tracers 15 15 USE dom_oce ! ocean space and time domain 16 USE domutl, ONLY : is_tile 16 17 USE ldfdyn ! lateral diffusion: eddy viscosity coef. 17 18 USE ldfslp ! iso-neutral slopes … … 21 22 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 22 23 USE lib_mpp 23 24 #if defined key_loop_fusion 25 USE dynldf_lap_blp_lf 26 #endif 27 24 28 IMPLICIT NONE 25 29 PRIVATE … … 39 43 40 44 SUBROUTINE dyn_ldf_lap( kt, Kbb, Kmm, pu, pv, pu_rhs, pv_rhs, kpass ) 45 !! 46 INTEGER , INTENT(in ) :: kt ! ocean time-step index 47 INTEGER , INTENT(in ) :: Kbb, Kmm ! ocean time level indices 48 INTEGER , INTENT(in ) :: kpass ! =1/2 first or second passage 49 REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: pu, pv ! before velocity [m/s] 50 REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: pu_rhs, pv_rhs ! velocity trend [m/s2] 51 !! 52 #if defined key_loop_fusion 53 CALL dyn_ldf_lap_lf( kt, Kbb, Kmm, pu, pv, pu_rhs, pv_rhs, kpass ) 54 #else 55 CALL dyn_ldf_lap_t( kt, Kbb, Kmm, pu, pv, is_tile(pu), pu_rhs, pv_rhs, is_tile(pu_rhs), kpass ) 56 #endif 57 58 END SUBROUTINE dyn_ldf_lap 59 60 61 SUBROUTINE dyn_ldf_lap_t( kt, Kbb, Kmm, pu, pv, ktuv, pu_rhs, pv_rhs, ktuv_rhs, kpass ) 41 62 !!---------------------------------------------------------------------- 42 63 !! *** ROUTINE dyn_ldf_lap *** … … 52 73 !! Reference : S.Griffies, R.Hallberg 2000 Mon.Wea.Rev., DOI:/ 53 74 !!---------------------------------------------------------------------- 54 INTEGER , INTENT(in ) :: kt ! ocean time-step index 55 INTEGER , INTENT(in ) :: Kbb, Kmm ! ocean time level indices 56 INTEGER , INTENT(in ) :: kpass ! =1/2 first or second passage 57 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in ) :: pu, pv ! before velocity [m/s] 58 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pu_rhs, pv_rhs ! velocity trend [m/s2] 75 INTEGER , INTENT(in ) :: kt ! ocean time-step index 76 INTEGER , INTENT(in ) :: Kbb, Kmm ! ocean time level indices 77 INTEGER , INTENT(in ) :: kpass ! =1/2 first or second passage 78 INTEGER , INTENT(in ) :: ktuv, ktuv_rhs 79 REAL(wp), DIMENSION(A2D_T(ktuv) ,JPK), INTENT(in ) :: pu, pv ! before velocity [m/s] 80 REAL(wp), DIMENSION(A2D_T(ktuv_rhs),JPK), INTENT(inout) :: pu_rhs, pv_rhs ! velocity trend [m/s2] 59 81 ! 60 82 INTEGER :: ji, jj, jk ! dummy loop indices 83 INTEGER :: iij 61 84 REAL(wp) :: zsign ! local scalars 62 85 REAL(wp) :: zua, zva ! local scalars … … 65 88 !!---------------------------------------------------------------------- 66 89 ! 67 IF( kt == nit000 .AND. lwp ) THEN 68 WRITE(numout,*) 69 WRITE(numout,*) 'dyn_ldf : iso-level harmonic (laplacian) operator, pass=', kpass 70 WRITE(numout,*) '~~~~~~~ ' 90 IF( .NOT. l_istiled .OR. ntile == 1 ) THEN ! Do only on the first tile 91 IF( kt == nit000 .AND. lwp ) THEN 92 WRITE(numout,*) 93 WRITE(numout,*) 'dyn_ldf : iso-level harmonic (laplacian) operator, pass=', kpass 94 WRITE(numout,*) '~~~~~~~ ' 95 ENDIF 96 ENDIF 97 ! 98 ! Define pu_rhs/pv_rhs halo points for multi-point haloes in bilaplacian case 99 IF( nldf_dyn == np_blp .AND. kpass == 1 ) THEN ; iij = nn_hls 100 ELSE ; iij = 1 71 101 ENDIF 72 102 ! … … 79 109 CASE ( np_typ_rot ) !== Vorticity-Divergence operator ==! 80 110 ! 81 ALLOCATE( zcur( jpi,jpj) , zdiv(jpi,jpj) )111 ALLOCATE( zcur(A2D(nn_hls)) , zdiv(A2D(nn_hls)) ) 82 112 ! 83 113 DO jk = 1, jpkm1 ! Horizontal slab 84 114 ! 85 DO_2D( 0, 1, 0, 1)115 DO_2D( iij-1, iij, iij-1, iij ) 86 116 ! ! ahm * e3 * curl (computed from 1 to jpim1/jpjm1) 87 117 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 … … 94 124 END_2D 95 125 ! 96 DO_2D( 0, 0, 0, 0 )! - curl( curl) + grad( div )126 DO_2D( iij-1, iij-1, iij-1, iij-1 ) ! - curl( curl) + grad( div ) 97 127 pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + zsign * umask(ji,jj,jk) * ( & ! * by umask is mandatory for dyn_ldf_blp use 98 128 & - ( zcur(ji ,jj) - zcur(ji,jj-1) ) * r1_e2u(ji,jj) / e3u(ji,jj,jk,Kmm) & … … 110 140 CASE ( np_typ_sym ) !== Symmetric operator ==! 111 141 ! 112 ALLOCATE( zten( jpi,jpj) , zshe(jpi,jpj) )142 ALLOCATE( zten(A2D(nn_hls)) , zshe(A2D(nn_hls)) ) 113 143 ! 114 144 DO jk = 1, jpkm1 ! Horizontal slab 115 145 ! 116 DO_2D( 0, 1, 0, 1)146 DO_2D( iij-1, iij, iij-1, iij ) 117 147 ! ! shearing stress component (F-point) NB : ahmf has already been multiplied by fmask 118 148 zshe(ji-1,jj-1) = ahmf(ji-1,jj-1,jk) & … … 129 159 END_2D 130 160 ! 131 DO_2D( 0, 0, 0, 0)161 DO_2D( iij-1, iij-1, iij-1, iij-1 ) 132 162 pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + zsign * r1_e1e2u(ji,jj) / e3u(ji,jj,jk,Kmm) & 133 163 & * ( ( zten(ji+1,jj ) * e2t(ji+1,jj )*e2t(ji+1,jj ) * e3t(ji+1,jj ,jk,Kmm) & … … 150 180 END SELECT 151 181 ! 152 END SUBROUTINE dyn_ldf_lap 182 END SUBROUTINE dyn_ldf_lap_t 153 183 154 184 … … 171 201 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pu_rhs, pv_rhs ! momentum trend 172 202 ! 173 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zulap, zvlap ! laplacian at u- and v-point 174 !!---------------------------------------------------------------------- 175 ! 176 IF( kt == nit000 ) THEN 177 IF(lwp) WRITE(numout,*) 178 IF(lwp) WRITE(numout,*) 'dyn_ldf_blp : bilaplacian operator momentum ' 179 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~' 203 REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zulap, zvlap ! laplacian at u- and v-point 204 !!---------------------------------------------------------------------- 205 ! 206 #if defined key_loop_fusion 207 CALL dyn_ldf_blp_lf( kt, Kbb, Kmm, pu, pv, pu_rhs, pv_rhs ) 208 #else 209 IF( .NOT. l_istiled .OR. ntile == 1 ) THEN ! Do only on the first tile 210 IF( kt == nit000 ) THEN 211 IF(lwp) WRITE(numout,*) 212 IF(lwp) WRITE(numout,*) 'dyn_ldf_blp : bilaplacian operator momentum ' 213 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~' 214 ENDIF 180 215 ENDIF 181 216 ! … … 185 220 CALL dyn_ldf_lap( kt, Kbb, Kmm, pu, pv, zulap, zvlap, 1 ) ! rotated laplacian applied to pt (output in zlap,Kbb) 186 221 ! 187 CALL lbc_lnk( 'dynldf_lap_blp', zulap, 'U', -1.0_wp, zvlap, 'V', -1.0_wp ) ! Lateral boundary conditions222 IF (nn_hls==1) CALL lbc_lnk( 'dynldf_lap_blp', zulap, 'U', -1.0_wp, zvlap, 'V', -1.0_wp ) ! Lateral boundary conditions 188 223 ! 189 224 CALL dyn_ldf_lap( kt, Kbb, Kmm, zulap, zvlap, pu_rhs, pv_rhs, 2 ) ! rotated laplacian applied to zlap (output in pt(:,:,:,:,Krhs)) 190 225 ! 226 #endif 191 227 END SUBROUTINE dyn_ldf_blp 192 228 -
NEMO/trunk/src/OCE/DYN/dynspg_ts.F90
r14433 r14834 730 730 IF (ln_bt_fw) THEN 731 731 IF( .NOT.( kt == nit000 .AND. l_1st_euler ) ) THEN 732 DO_2D( 1, 1, 1, 1)732 DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 733 733 zun_save = un_adv(ji,jj) 734 734 zvn_save = vn_adv(ji,jj) -
NEMO/trunk/src/OCE/DYN/dynvor.F90
r14820 r14834 240 240 INTEGER :: ji, jj, jk ! dummy loop indices 241 241 REAL(wp) :: zx1, zy1, zx2, zy2 ! local scalars 242 REAL(wp), DIMENSION(jpi,jpj) :: zwx, zwy, zwt ! 2D workspace 243 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zwz ! 3D workspace, jpkm1 -> avoid lbc_lnk on jpk that is not defined 244 !!---------------------------------------------------------------------- 245 ! 246 IF( kt == nit000 ) THEN 247 IF(lwp) WRITE(numout,*) 248 IF(lwp) WRITE(numout,*) 'dyn:vor_enT : vorticity term: t-point energy conserving scheme' 249 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 242 REAL(wp), DIMENSION(A2D(nn_hls)) :: zwx, zwy, zwt ! 2D workspace 243 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zwz ! 3D workspace, jpkm1 -> avoid lbc_lnk on jpk that is not defined 244 !!---------------------------------------------------------------------- 245 ! 246 IF( .NOT. l_istiled .OR. ntile == 1 ) THEN ! Do only on the first tile 247 IF( kt == nit000 ) THEN 248 IF(lwp) WRITE(numout,*) 249 IF(lwp) WRITE(numout,*) 'dyn:vor_enT : vorticity term: t-point energy conserving scheme' 250 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 251 ENDIF 250 252 ENDIF 251 253 ! … … 254 256 ! 255 257 CASE ( np_RVO , np_CRV ) !* relative vorticity at f-point is used 256 ALLOCATE( zwz( jpi,jpj,jpk) )258 ALLOCATE( zwz(A2D(nn_hls),jpk) ) 257 259 DO jk = 1, jpkm1 ! Horizontal slab 258 DO_2D( 1, 0, 1, 0)260 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 259 261 zwz(ji,jj,jk) = ( e2v(ji+1,jj) * pv(ji+1,jj,jk) - e2v(ji,jj) * pv(ji,jj,jk) & 260 262 & - e1u(ji,jj+1) * pu(ji,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk) ) * r1_e1e2f(ji,jj) 261 263 END_2D 262 264 IF( ln_dynvor_msk ) THEN ! mask relative vorticity 263 DO_2D( 1, 0, 1, 0)265 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 264 266 zwz(ji,jj,jk) = zwz(ji,jj,jk) * fmask(ji,jj,jk) 265 267 END_2D 266 268 ENDIF 267 269 END DO 268 CALL lbc_lnk( 'dynvor', zwz, 'F', 1.0_wp )270 IF (nn_hls==1) CALL lbc_lnk( 'dynvor', zwz, 'F', 1.0_wp ) 269 271 ! 270 272 END SELECT … … 277 279 ! 278 280 CASE ( np_COR ) !* Coriolis (planetary vorticity) 279 zwt(:,:) = ff_t(:,:) * e1e2t(:,:)*e3t(:,:,jk,Kmm) 281 DO_2D( 0, 1, 0, 1 ) 282 zwt(ji,jj) = ff_t(ji,jj) * e1e2t(ji,jj)*e3t(ji,jj,jk,Kmm) 283 END_2D 280 284 CASE ( np_RVO ) !* relative vorticity 281 285 DO_2D( 0, 1, 0, 1 ) … … 356 360 INTEGER :: ji, jj, jk ! dummy loop indices 357 361 REAL(wp) :: zx1, zy1, zx2, zy2, ze3f, zmsk ! local scalars 358 REAL(wp), DIMENSION(jpi,jpj) :: zwx, zwy, zwz ! 2D workspace 359 !!---------------------------------------------------------------------- 360 ! 361 IF( kt == nit000 ) THEN 362 IF(lwp) WRITE(numout,*) 363 IF(lwp) WRITE(numout,*) 'dyn:vor_ene : vorticity term: energy conserving scheme' 364 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 362 REAL(wp), DIMENSION(A2D(nn_hls)) :: zwx, zwy, zwz ! 2D workspace 363 !!---------------------------------------------------------------------- 364 ! 365 IF( .NOT. l_istiled .OR. ntile == 1 ) THEN ! Do only on the first tile 366 IF( kt == nit000 ) THEN 367 IF(lwp) WRITE(numout,*) 368 IF(lwp) WRITE(numout,*) 'dyn:vor_ene : vorticity term: energy conserving scheme' 369 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 370 ENDIF 365 371 ENDIF 366 372 ! … … 371 377 SELECT CASE( kvor ) !== vorticity considered ==! 372 378 CASE ( np_COR ) !* Coriolis (planetary vorticity) 373 zwz(:,:) = ff_f(:,:) 379 DO_2D( 1, 0, 1, 0 ) 380 zwz(ji,jj) = ff_f(ji,jj) 381 END_2D 374 382 CASE ( np_RVO ) !* relative vorticity 375 383 DO_2D( 1, 0, 1, 0 ) … … 437 445 #endif 438 446 ! !== horizontal fluxes ==! 439 zwx(:,:) = e2u(:,:) * e3u(:,:,jk,Kmm) * pu(:,:,jk) 440 zwy(:,:) = e1v(:,:) * e3v(:,:,jk,Kmm) * pv(:,:,jk) 447 DO_2D( 1, 1, 1, 1 ) 448 zwx(ji,jj) = e2u(ji,jj) * e3u(ji,jj,jk,Kmm) * pu(ji,jj,jk) 449 zwy(ji,jj) = e1v(ji,jj) * e3v(ji,jj,jk,Kmm) * pv(ji,jj,jk) 450 END_2D 441 451 ! 442 452 ! !== compute and add the vorticity term trend =! … … 483 493 INTEGER :: ji, jj, jk ! dummy loop indices 484 494 REAL(wp) :: zuav, zvau, ze3f, zmsk ! local scalars 485 REAL(wp), DIMENSION(jpi,jpj) :: zwx, zwy, zwz, zww ! 2D workspace 486 !!---------------------------------------------------------------------- 487 ! 488 IF( kt == nit000 ) THEN 489 IF(lwp) WRITE(numout,*) 490 IF(lwp) WRITE(numout,*) 'dyn:vor_ens : vorticity term: enstrophy conserving scheme' 491 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 495 REAL(wp), DIMENSION(A2D(nn_hls)) :: zwx, zwy, zwz ! 2D workspace 496 !!---------------------------------------------------------------------- 497 ! 498 IF( .NOT. l_istiled .OR. ntile == 1 ) THEN ! Do only on the first tile 499 IF( kt == nit000 ) THEN 500 IF(lwp) WRITE(numout,*) 501 IF(lwp) WRITE(numout,*) 'dyn:vor_ens : vorticity term: enstrophy conserving scheme' 502 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 503 ENDIF 492 504 ENDIF 493 505 ! ! =============== … … 497 509 SELECT CASE( kvor ) !== vorticity considered ==! 498 510 CASE ( np_COR ) !* Coriolis (planetary vorticity) 499 zwz(:,:) = ff_f(:,:) 511 DO_2D( 1, 0, 1, 0 ) 512 zwz(ji,jj) = ff_f(ji,jj) 513 END_2D 500 514 CASE ( np_RVO ) !* relative vorticity 501 515 DO_2D( 1, 0, 1, 0 ) … … 564 578 #endif 565 579 ! !== horizontal fluxes ==! 566 zwx(:,:) = e2u(:,:) * e3u(:,:,jk,Kmm) * pu(:,:,jk) 567 zwy(:,:) = e1v(:,:) * e3v(:,:,jk,Kmm) * pv(:,:,jk) 580 DO_2D( 1, 1, 1, 1 ) 581 zwx(ji,jj) = e2u(ji,jj) * e3u(ji,jj,jk,Kmm) * pu(ji,jj,jk) 582 zwy(ji,jj) = e1v(ji,jj) * e3v(ji,jj,jk,Kmm) * pv(ji,jj,jk) 583 END_2D 568 584 ! 569 585 ! !== compute and add the vorticity term trend =! … … 609 625 REAL(wp) :: zua, zva ! local scalars 610 626 REAL(wp) :: zmsk, ze3f ! local scalars 611 REAL(wp), DIMENSION(jpi,jpj) :: zwx , zwy , z1_e3f 612 REAL(wp), DIMENSION(jpi,jpj) :: ztnw, ztne, ztsw, ztse 613 REAL(wp), DIMENSION(jpi,jpj,jpkm1) :: zwz ! 3D workspace, jpkm1 -> jpkm1 -> avoid lbc_lnk on jpk that is not defined 614 !!---------------------------------------------------------------------- 615 ! 616 IF( kt == nit000 ) THEN 617 IF(lwp) WRITE(numout,*) 618 IF(lwp) WRITE(numout,*) 'dyn:vor_een : vorticity term: energy and enstrophy conserving scheme' 619 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 627 REAL(wp), DIMENSION(A2D(nn_hls)) :: z1_e3f 628 #if defined key_loop_fusion 629 REAL(wp) :: ztne, ztnw, ztnw_ip1, ztse, ztse_jp1, ztsw_jp1, ztsw_ip1 630 REAL(wp) :: zwx, zwx_im1, zwx_jp1, zwx_im1_jp1 631 REAL(wp) :: zwy, zwy_ip1, zwy_jm1, zwy_ip1_jm1 632 #else 633 REAL(wp), DIMENSION(A2D(nn_hls)) :: zwx , zwy 634 REAL(wp), DIMENSION(A2D(nn_hls)) :: ztnw, ztne, ztsw, ztse 635 #endif 636 REAL(wp), DIMENSION(A2D(nn_hls),jpkm1) :: zwz ! 3D workspace, jpkm1 -> jpkm1 -> avoid lbc_lnk on jpk that is not defined 637 !!---------------------------------------------------------------------- 638 ! 639 IF( .NOT. l_istiled .OR. ntile == 1 ) THEN ! Do only on the first tile 640 IF( kt == nit000 ) THEN 641 IF(lwp) WRITE(numout,*) 642 IF(lwp) WRITE(numout,*) 'dyn:vor_een : vorticity term: energy and enstrophy conserving scheme' 643 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 644 ENDIF 620 645 ENDIF 621 646 ! … … 625 650 ! 626 651 #if defined key_qco || defined key_linssh 627 DO_2D( 1, 0, 1, 0) ! == reciprocal of e3 at F-point (key_qco)652 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) ! == reciprocal of e3 at F-point (key_qco) 628 653 z1_e3f(ji,jj) = 1._wp / e3f_vor(ji,jj,jk) 629 654 END_2D … … 631 656 SELECT CASE( nn_e3f_typ ) ! == reciprocal of e3 at F-point 632 657 CASE ( 0 ) ! original formulation (masked averaging of e3t divided by 4) 633 DO_2D( 1, 0, 1, 0)658 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 634 659 ! round brackets added to fix the order of floating point operations 635 660 ! needed to ensure halo 1 - halo 2 compatibility … … 643 668 END_2D 644 669 CASE ( 1 ) ! new formulation (masked averaging of e3t divided by the sum of mask) 645 DO_2D( 1, 0, 1, 0)670 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 646 671 ! round brackets added to fix the order of floating point operations 647 672 ! needed to ensure halo 1 - halo 2 compatibility … … 662 687 ! 663 688 CASE ( np_COR ) !* Coriolis (planetary vorticity) 664 DO_2D( 1, 0, 1, 0)689 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 665 690 zwz(ji,jj,jk) = ff_f(ji,jj) * z1_e3f(ji,jj) 666 691 END_2D 667 692 CASE ( np_RVO ) !* relative vorticity 668 DO_2D( 1, 0, 1, 0)693 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 669 694 zwz(ji,jj,jk) = ( e2v(ji+1,jj ) * pv(ji+1,jj,jk) - e2v(ji,jj) * pv(ji,jj,jk) & 670 695 & - e1u(ji ,jj+1) * pu(ji,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk) ) * r1_e1e2f(ji,jj)*z1_e3f(ji,jj) 671 696 END_2D 672 697 IF( ln_dynvor_msk ) THEN ! mask the relative vorticity 673 DO_2D( 1, 0, 1, 0)698 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 674 699 zwz(ji,jj,jk) = zwz(ji,jj,jk) * fmask(ji,jj,jk) 675 700 END_2D 676 701 ENDIF 677 702 CASE ( np_MET ) !* metric term 678 DO_2D( 1, 0, 1, 0)703 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 679 704 zwz(ji,jj,jk) = ( ( pv(ji+1,jj,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj) & 680 705 & - ( pu(ji,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj) ) * z1_e3f(ji,jj) 681 706 END_2D 682 707 CASE ( np_CRV ) !* Coriolis + relative vorticity 683 DO_2D( 1, 0, 1, 0)684 685 708 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 709 ! round brackets added to fix the order of floating point operations 710 ! needed to ensure halo 1 - halo 2 compatibility 686 711 zwz(ji,jj,jk) = ( ff_f(ji,jj) + ( ( e2v(ji+1,jj ) * pv(ji+1,jj,jk) - e2v(ji,jj) * pv(ji,jj,jk) & 687 712 & ) & ! bracket for halo 1 - halo 2 compatibility 688 & - ( e1u(ji ,jj+1) * pu(ji,jj+1,jk) - e1u(ji,jj) * pu(ji,jj,jk) & 713 & - ( e1u(ji ,jj+1) * pu(ji,jj+1,jk) - e1u(ji,jj) * pu(ji,jj,jk) & 689 714 & ) & ! bracket for halo 1 - halo 2 compatibility 690 715 & ) * r1_e1e2f(ji,jj) ) * z1_e3f(ji,jj) 691 716 END_2D 692 717 IF( ln_dynvor_msk ) THEN ! mask the relative vorticity 693 DO_2D( 1, 0, 1, 0)718 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 694 719 zwz(ji,jj,jk) = ( zwz(ji,jj,jk) - ff_f(ji,jj) ) * fmask(ji,jj,jk) + ff_f(ji,jj) 695 720 END_2D 696 721 ENDIF 697 722 CASE ( np_CME ) !* Coriolis + metric 698 DO_2D( 1, 0, 1, 0)723 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 699 724 zwz(ji,jj,jk) = ( ff_f(ji,jj) + ( pv(ji+1,jj ,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj) & 700 725 & - ( pu(ji ,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj) ) * z1_e3f(ji,jj) … … 707 732 ! ! =============== 708 733 ! 709 CALL lbc_lnk( 'dynvor', zwz, 'F', 1.0_wp ) 710 ! 711 ! ! =============== 712 DO jk = 1, jpkm1 ! Horizontal slab 713 ! ! =============== 714 ! 734 IF (nn_hls==1) CALL lbc_lnk( 'dynvor', zwz, 'F', 1.0_wp ) 735 ! 736 ! ! =============== 737 ! ! Horizontal slab 738 ! ! =============== 739 #if defined key_loop_fusion 740 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 715 741 ! !== horizontal fluxes ==! 716 zwx(:,:) = e2u(:,:) * e3u(:,:,jk,Kmm) * pu(:,:,jk) 717 zwy(:,:) = e1v(:,:) * e3v(:,:,jk,Kmm) * pv(:,:,jk) 742 zwx = e2u(ji ,jj ) * e3u(ji ,jj ,jk,Kmm) * pu(ji ,jj ,jk) 743 zwx_im1 = e2u(ji-1,jj ) * e3u(ji-1,jj ,jk,Kmm) * pu(ji-1,jj ,jk) 744 zwx_jp1 = e2u(ji ,jj+1) * e3u(ji ,jj+1,jk,Kmm) * pu(ji ,jj+1,jk) 745 zwx_im1_jp1 = e2u(ji-1,jj+1) * e3u(ji-1,jj+1,jk,Kmm) * pu(ji-1,jj+1,jk) 746 zwy = e1v(ji ,jj ) * e3v(ji ,jj ,jk,Kmm) * pv(ji ,jj ,jk) 747 zwy_ip1 = e1v(ji+1,jj ) * e3v(ji+1,jj ,jk,Kmm) * pv(ji+1,jj ,jk) 748 zwy_jm1 = e1v(ji ,jj-1) * e3v(ji ,jj-1,jk,Kmm) * pv(ji ,jj-1,jk) 749 zwy_ip1_jm1 = e1v(ji+1,jj-1) * e3v(ji+1,jj-1,jk,Kmm) * pv(ji+1,jj-1,jk) 750 ! !== compute and add the vorticity term trend =! 751 ztne = zwz(ji-1,jj ,jk) + zwz(ji ,jj ,jk) + zwz(ji ,jj-1,jk) 752 ztnw = zwz(ji-1,jj-1,jk) + zwz(ji-1,jj ,jk) + zwz(ji ,jj ,jk) 753 ztnw_ip1 = zwz(ji ,jj-1,jk) + zwz(ji ,jj ,jk) + zwz(ji+1,jj ,jk) 754 ztse = zwz(ji ,jj ,jk) + zwz(ji ,jj-1,jk) + zwz(ji-1,jj-1,jk) 755 ztse_jp1 = zwz(ji ,jj+1,jk) + zwz(ji ,jj ,jk) + zwz(ji-1,jj ,jk) 756 ztsw_jp1 = zwz(ji ,jj ,jk) + zwz(ji-1,jj ,jk) + zwz(ji-1,jj+1,jk) 757 ztsw_ip1 = zwz(ji+1,jj-1,jk) + zwz(ji ,jj-1,jk) + zwz(ji ,jj ,jk) 758 ! 759 zua = + r1_12 * r1_e1u(ji,jj) * ( ztne * zwy + ztnw_ip1 * zwy_ip1 & 760 & + ztse * zwy_jm1 + ztsw_ip1 * zwy_ip1_jm1 ) 761 zva = - r1_12 * r1_e2v(ji,jj) * ( ztsw_jp1 * zwx_im1_jp1 + ztse_jp1 * zwx_jp1 & 762 & + ztnw * zwx_im1 + ztne * zwx ) 763 pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + zua 764 pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) + zva 765 END_3D 766 #else 767 DO jk = 1, jpkm1 768 ! 769 ! !== horizontal fluxes ==! 770 DO_2D( 1, 1, 1, 1 ) 771 zwx(ji,jj) = e2u(ji,jj) * e3u(ji,jj,jk,Kmm) * pu(ji,jj,jk) 772 zwy(ji,jj) = e1v(ji,jj) * e3v(ji,jj,jk,Kmm) * pv(ji,jj,jk) 773 END_2D 718 774 ! 719 775 ! !== compute and add the vorticity term trend =! … … 733 789 pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) + zva 734 790 END_2D 735 ! ! =============== 736 END DO ! End of slab 791 END DO 792 #endif 793 ! ! =============== 794 ! ! End of slab 737 795 ! ! =============== 738 796 END SUBROUTINE vor_een … … 766 824 REAL(wp) :: zua, zva ! local scalars 767 825 REAL(wp) :: zmsk, z1_e3t ! local scalars 768 REAL(wp), DIMENSION(jpi,jpj) :: zwx , zwy 769 REAL(wp), DIMENSION(jpi,jpj) :: ztnw, ztne, ztsw, ztse 770 REAL(wp), DIMENSION(jpi,jpj,jpkm1) :: zwz ! 3D workspace, avoid lbc_lnk on jpk that is not defined 771 !!---------------------------------------------------------------------- 772 ! 773 IF( kt == nit000 ) THEN 774 IF(lwp) WRITE(numout,*) 775 IF(lwp) WRITE(numout,*) 'dyn:vor_eeT : vorticity term: energy and enstrophy conserving scheme' 776 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 826 REAL(wp), DIMENSION(A2D(nn_hls)) :: zwx , zwy 827 REAL(wp), DIMENSION(A2D(nn_hls)) :: ztnw, ztne, ztsw, ztse 828 REAL(wp), DIMENSION(A2D(nn_hls),jpkm1) :: zwz ! 3D workspace, avoid lbc_lnk on jpk that is not defined 829 !!---------------------------------------------------------------------- 830 ! 831 IF( .NOT. l_istiled .OR. ntile == 1 ) THEN ! Do only on the first tile 832 IF( kt == nit000 ) THEN 833 IF(lwp) WRITE(numout,*) 834 IF(lwp) WRITE(numout,*) 'dyn:vor_eeT : vorticity term: energy and enstrophy conserving scheme' 835 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 836 ENDIF 777 837 ENDIF 778 838 ! … … 784 844 SELECT CASE( kvor ) !== vorticity considered ==! 785 845 CASE ( np_COR ) !* Coriolis (planetary vorticity) 786 DO_2D( 1, 0, 1, 0)846 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 787 847 zwz(ji,jj,jk) = ff_f(ji,jj) 788 848 END_2D 789 849 CASE ( np_RVO ) !* relative vorticity 790 DO_2D( 1, 0, 1, 0)850 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 791 851 ! round brackets added to fix the order of floating point operations 792 852 ! needed to ensure halo 1 - halo 2 compatibility … … 796 856 END_2D 797 857 IF( ln_dynvor_msk ) THEN ! mask the relative vorticity 798 DO_2D( 1, 0, 1, 0)858 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 799 859 zwz(ji,jj,jk) = zwz(ji,jj,jk) * fmask(ji,jj,jk) 800 860 END_2D 801 861 ENDIF 802 862 CASE ( np_MET ) !* metric term 803 DO_2D( 1, 0, 1, 0)863 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 804 864 zwz(ji,jj,jk) = ( pv(ji+1,jj ,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj) & 805 865 & - ( pu(ji ,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj) 806 866 END_2D 807 867 CASE ( np_CRV ) !* Coriolis + relative vorticity 808 DO_2D( 1, 0, 1, 0)868 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 809 869 ! round brackets added to fix the order of floating point operations 810 870 ! needed to ensure halo 1 - halo 2 compatibility 811 871 zwz(ji,jj,jk) = ( ff_f(ji,jj) + ( (e2v(ji+1,jj ) * pv(ji+1,jj ,jk) - e2v(ji,jj) * pv(ji,jj,jk)) & 812 872 & - (e1u(ji ,jj+1) * pu(ji ,jj+1,jk) - e1u(ji,jj) * pu(ji,jj,jk)) ) & 813 & * r1_e1e2f(ji,jj) )873 & * r1_e1e2f(ji,jj) ) 814 874 END_2D 815 875 IF( ln_dynvor_msk ) THEN ! mask the relative vorticity 816 DO_2D( 1, 0, 1, 0)876 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 817 877 zwz(ji,jj,jk) = ( zwz(ji,jj,jk) - ff_f(ji,jj) ) * fmask(ji,jj,jk) + ff_f(ji,jj) 818 878 END_2D 819 879 ENDIF 820 880 CASE ( np_CME ) !* Coriolis + metric 821 DO_2D( 1, 0, 1, 0)881 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 822 882 zwz(ji,jj,jk) = ff_f(ji,jj) + ( pv(ji+1,jj ,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj) & 823 883 & - ( pu(ji ,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj) … … 831 891 ! ! =============== 832 892 ! 833 CALL lbc_lnk( 'dynvor', zwz, 'F', 1.0_wp )893 IF (nn_hls==1) CALL lbc_lnk( 'dynvor', zwz, 'F', 1.0_wp ) 834 894 ! 835 895 ! ! =============== … … 838 898 ! 839 899 ! !== horizontal fluxes ==! 840 zwx(:,:) = e2u(:,:) * e3u(:,:,jk,Kmm) * pu(:,:,jk) 841 zwy(:,:) = e1v(:,:) * e3v(:,:,jk,Kmm) * pv(:,:,jk) 900 DO_2D( 1, 1, 1, 1 ) 901 zwx(ji,jj) = e2u(ji,jj) * e3u(ji,jj,jk,Kmm) * pu(ji,jj,jk) 902 zwy(ji,jj) = e1v(ji,jj) * e3v(ji,jj,jk,Kmm) * pv(ji,jj,jk) 903 END_2D 842 904 ! 843 905 ! !== compute and add the vorticity term trend =! -
NEMO/trunk/src/OCE/DYN/dynzad.F90
r14072 r14834 60 60 INTEGER :: ji, jj, jk ! dummy loop indices 61 61 REAL(wp) :: zua, zva ! local scalars 62 REAL(wp), DIMENSION( jpi,jpj) :: zww63 REAL(wp), DIMENSION( jpi,jpj,jpk) :: zwuw, zwvw62 REAL(wp), DIMENSION(A2D(nn_hls)) :: zww 63 REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zwuw, zwvw 64 64 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ztrdu, ztrdv 65 65 !!---------------------------------------------------------------------- … … 67 67 IF( ln_timing ) CALL timing_start('dyn_zad') 68 68 ! 69 IF( kt == nit000 ) THEN 70 IF(lwp) WRITE(numout,*) 71 IF(lwp) WRITE(numout,*) 'dyn_zad : 2nd order vertical advection scheme' 69 IF( .NOT. l_istiled .OR. ntile == 1 ) THEN ! Do only on the first tile 70 IF( kt == nit000 ) THEN 71 IF(lwp) WRITE(numout,*) 72 IF(lwp) WRITE(numout,*) 'dyn_zad : 2nd order vertical advection scheme' 73 ENDIF 72 74 ENDIF 73 75 … … 79 81 80 82 DO jk = 2, jpkm1 ! Vertical momentum advection at level w and u- and v- vertical 81 DO_2D( 0, 1, 0, 1 ) ! vertical fluxes 82 IF( ln_vortex_force ) THEN 83 zww(ji,jj) = 0.25_wp * e1e2t(ji,jj) * ( ww(ji,jj,jk) + wsd(ji,jj,jk) ) 84 ELSE 85 zww(ji,jj) = 0.25_wp * e1e2t(ji,jj) * ww(ji,jj,jk) 86 ENDIF 87 END_2D 83 IF( ln_vortex_force ) THEN ! vertical fluxes 84 DO_2D( 0, 1, 0, 1 ) 85 zww(ji,jj) = 0.25_wp * e1e2t(ji,jj) * ( ww(ji,jj,jk) + wsd(ji,jj,jk) ) 86 END_2D 87 ELSE 88 DO_2D( 0, 1, 0, 1 ) 89 zww(ji,jj) = 0.25_wp * e1e2t(ji,jj) * ww(ji,jj,jk) 90 END_2D 91 ENDIF 88 92 DO_2D( 0, 0, 0, 0 ) ! vertical momentum advection at w-point 89 93 zwuw(ji,jj,jk) = ( zww(ji+1,jj ) + zww(ji,jj) ) * ( puu(ji,jj,jk-1,Kmm) - puu(ji,jj,jk,Kmm) ) -
NEMO/trunk/src/OCE/DYN/dynzdf.F90
r13497 r14834 19 19 USE zdfdrg ! vertical physics: top/bottom drag coef. 20 20 USE dynadv ,ONLY: ln_dynadv_vec ! dynamics: advection form 21 #if defined key_loop_fusion 22 USE dynldf_iso_lf,ONLY: akzu, akzv ! dynamics: vertical component of rotated lateral mixing 23 #else 21 24 USE dynldf_iso,ONLY: akzu, akzv ! dynamics: vertical component of rotated lateral mixing 25 #endif 22 26 USE ldfdyn ! lateral diffusion: eddy viscosity coef. and type of operator 23 27 USE trd_oce ! trends: ocean variables … … 78 82 REAL(wp) :: zWui, zWvi ! - - 79 83 REAL(wp) :: zWus, zWvs ! - - 80 REAL(wp), DIMENSION( jpi,jpj,jpk) :: zwi, zwd, zws ! 3D workspace84 REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zwi, zwd, zws ! 3D workspace 81 85 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ztrdu, ztrdv ! - - 82 86 !!--------------------------------------------------------------------- … … 84 88 IF( ln_timing ) CALL timing_start('dyn_zdf') 85 89 ! 86 IF( kt == nit000 ) THEN !* initialization 87 IF(lwp) WRITE(numout,*) 88 IF(lwp) WRITE(numout,*) 'dyn_zdf_imp : vertical momentum diffusion implicit operator' 89 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ ' 90 ! 91 If( ln_linssh ) THEN ; r_vvl = 0._wp ! non-linear free surface indicator 92 ELSE ; r_vvl = 1._wp 90 IF( .NOT. l_istiled .OR. ntile == 1 ) THEN ! Do only on the first tile 91 IF( kt == nit000 ) THEN !* initialization 92 IF(lwp) WRITE(numout,*) 93 IF(lwp) WRITE(numout,*) 'dyn_zdf_imp : vertical momentum diffusion implicit operator' 94 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ ' 95 ! 96 If( ln_linssh ) THEN ; r_vvl = 0._wp ! non-linear free surface indicator 97 ELSE ; r_vvl = 1._wp 98 ENDIF 93 99 ENDIF 94 100 ENDIF -
NEMO/trunk/src/OCE/DYN/sshwzv.F90
r14205 r14834 78 78 REAL(wp), DIMENSION(jpi,jpj,jpt), INTENT(inout) :: pssh ! sea-surface height 79 79 ! 80 INTEGER :: j k ! dummy loop index80 INTEGER :: ji, jj, jk ! dummy loop index 81 81 REAL(wp) :: zcoef ! local scalar 82 82 REAL(wp), DIMENSION(jpi,jpj) :: zhdiv ! 2D workspace … … 103 103 ! 104 104 zhdiv(:,:) = 0._wp 105 DO jk = 1, jpkm1! Horizontal divergence of barotropic transports106 zhdiv( :,:) = zhdiv(:,:) + e3t(:,:,jk,Kmm) * hdiv(:,:,jk)107 END DO105 DO_3D( 1, nn_hls, 1, nn_hls, 1, jpkm1 ) ! Horizontal divergence of barotropic transports 106 zhdiv(ji,jj) = zhdiv(ji,jj) + e3t(ji,jj,jk,Kmm) * hdiv(ji,jj,jk) 107 END_3D 108 108 ! ! Sea surface elevation time stepping 109 109 ! In time-split case we need a first guess of the ssh after (using the baroclinic timestep) in order to 110 110 ! compute the vertical velocity which can be used to compute the non-linear terms of the momentum equations. 111 111 ! 112 pssh(:,:,Kaa) = ( pssh(:,:,Kbb) - rDt * ( zcoef * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) ) ) * ssmask(:,:) 112 DO_2D_OVR( 1, nn_hls, 1, nn_hls ) ! Loop bounds limited by hdiv definition in div_hor 113 pssh(ji,jj,Kaa) = ( pssh(ji,jj,Kbb) - rDt * ( zcoef * ( emp_b(ji,jj) + emp(ji,jj) ) + zhdiv(ji,jj) ) ) * ssmask(ji,jj) 114 END_2D 115 ! pssh must be defined everywhere (true for dyn_spg_ts, not for dyn_spg_exp) 116 IF ( .NOT. ln_dynspg_ts .AND. nn_hls == 2 ) CALL lbc_lnk( 'sshwzv', pssh(:,:,Kaa), 'T', 1.0_wp ) 113 117 ! 114 118 #if defined key_agrif … … 119 123 IF ( .NOT.ln_dynspg_ts ) THEN 120 124 IF( ln_bdy ) THEN 121 CALL lbc_lnk( 'sshwzv', pssh(:,:,Kaa), 'T', 1.0_wp ) ! Not sure that's necessary125 IF (nn_hls==1) CALL lbc_lnk( 'sshwzv', pssh(:,:,Kaa), 'T', 1.0_wp ) ! Not sure that's necessary 122 126 CALL bdy_ssh( pssh(:,:,Kaa) ) ! Duplicate sea level across open boundaries 123 127 ENDIF … … 178 182 ! horizontal divergence of thickness diffusion transport ( velocity multiplied by e3t) 179 183 ! - ML - note: computation already done in dom_vvl_sf_nxt. Could be optimized (not critical and clearer this way) 180 DO_2D( 0, 0, 0, 0)184 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 181 185 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) ) 182 186 END_2D 183 187 END DO 184 CALL lbc_lnk('sshwzv', zhdiv, 'T', 1.0_wp) ! - ML - Perhaps not necessary: not used for horizontal "connexions"188 IF (nn_hls==1) CALL lbc_lnk('sshwzv', zhdiv, 'T', 1.0_wp) ! - ML - Perhaps not necessary: not used for horizontal "connexions" 185 189 ! ! Is it problematic to have a wrong vertical velocity in boundary cells? 186 190 ! ! Same question holds for hdiv. Perhaps just for security 187 DO jk = jpkm1, 1, -1! integrate from the bottom the hor. divergence191 DO_3DS( 1, 1, 1, 1, jpkm1, 1, -1 ) ! integrate from the bottom the hor. divergence 188 192 ! computation of w 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)193 END DO193 pww(ji,jj,jk) = pww(ji,jj,jk+1) - ( e3t(ji,jj,jk,Kmm) * hdiv(ji,jj,jk) & 194 & + zhdiv(ji,jj,jk) & 195 & + r1_Dt * ( e3t(ji,jj,jk,Kaa) & 196 & - e3t(ji,jj,jk,Kbb) ) ) * tmask(ji,jj,jk) 197 END_3D 194 198 ! IF( ln_vvl_layer ) pww(:,:,:) = 0.e0 195 199 DEALLOCATE( zhdiv ) … … 197 201 ELSEIF( ln_linssh ) THEN !== linear free surface cases ==! 198 202 ! !=================================! 199 DO jk = jpkm1, 1, -1! integrate from the bottom the hor. divergence200 pww( :,:,jk) = pww(:,:,jk+1) - ( e3t(:,:,jk,Kmm) * hdiv(:,:,jk) ) * tmask(:,:,jk)201 END DO203 DO_3DS( 1, 1, 1, 1, jpkm1, 1, -1 ) ! integrate from the bottom the hor. divergence 204 pww(ji,jj,jk) = pww(ji,jj,jk+1) - ( e3t(ji,jj,jk,Kmm) * hdiv(ji,jj,jk) ) * tmask(ji,jj,jk) 205 END_3D 202 206 ! !==========================================! 203 207 ELSE !== Quasi-Eulerian vertical coordinate ==! ('key_qco') 204 208 ! !==========================================! 205 DO jk = jpkm1, 1, -1! integrate from the bottom the hor. divergence206 pww( :,:,jk) = pww(:,:,jk+1) - ( e3t(:,:,jk,Kmm) * hdiv(:,:,jk)&207 & + r1_Dt * ( e3t(:,:,jk,Kaa) &208 & - e3t(:,:,jk,Kbb) ) ) * tmask(:,:,jk)209 END DO209 DO_3DS( 1, 1, 1, 1, jpkm1, 1, -1 ) ! integrate from the bottom the hor. divergence 210 pww(ji,jj,jk) = pww(ji,jj,jk+1) - ( e3t(ji,jj,jk,Kmm) * hdiv(ji,jj,jk) & 211 & + r1_Dt * ( e3t(ji,jj,jk,Kaa) & 212 & - e3t(ji,jj,jk,Kbb) ) ) * tmask(ji,jj,jk) 213 END_3D 210 214 ENDIF 211 215 … … 357 361 zdt = 2._wp * rn_Dt ! 2*rn_Dt and not rDt (for restartability) 358 362 IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN 359 DO_3D( 0, 0, 0, 0, 1, jpkm1 )363 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 ) 360 364 z1_e3t = 1._wp / e3t(ji,jj,jk,Kmm) 361 365 Cu_adv(ji,jj,jk) = zdt * & … … 374 378 END_3D 375 379 ELSE 376 DO_3D( 0, 0, 0, 0, 1, jpkm1 )380 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 ) 377 381 z1_e3t = 1._wp / e3t(ji,jj,jk,Kmm) 378 382 Cu_adv(ji,jj,jk) = zdt * & … … 387 391 END_3D 388 392 ENDIF 389 CALL lbc_lnk( 'sshwzv', Cu_adv, 'T', 1.0_wp )393 IF (nn_hls==1) CALL lbc_lnk( 'sshwzv', Cu_adv, 'T', 1.0_wp ) 390 394 ! 391 395 CALL iom_put("Courant",Cu_adv) 392 396 ! 393 397 IF( MAXVAL( Cu_adv(:,:,:) ) > Cu_min ) THEN ! Quick check if any breaches anywhere 394 DO_3DS( 1, 1, 1, 1, jpkm1, 2, -1 ) ! or scan Courant criterion and partition ! w where necessary398 DO_3DS( nn_hls, nn_hls, nn_hls, nn_hls, jpkm1, 2, -1 ) ! or scan Courant criterion and partition ! w where necessary 395 399 ! 396 400 zCu = MAX( Cu_adv(ji,jj,jk) , Cu_adv(ji,jj,jk-1) ) -
NEMO/trunk/src/OCE/DYN/wet_dry.F90
r14433 r14834 117 117 IF( ierr /= 0 ) CALL ctl_stop('STOP', 'wad_init : Array allocation error') 118 118 ENDIF 119 120 IF( ln_tile .AND. ln_wd_il ) CALL ctl_warn('Tiling has not been tested with ln_wd_il = T') 119 121 ! 120 122 END SUBROUTINE wad_init
Note: See TracChangeset
for help on using the changeset viewer.
