Changeset 5965 for branches/2014/dev_r4650_UKMO14.5_SST_BIAS_CORRECTION/NEMOGCM/NEMO/OPA_SRC/TRA/zpshde.F90
- Timestamp:
- 2015-12-01T16:35:30+01:00 (8 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2014/dev_r4650_UKMO14.5_SST_BIAS_CORRECTION/NEMOGCM/NEMO/OPA_SRC/TRA/zpshde.F90
r3294 r5965 8 8 !! - ! 2004-03 (C. Ethe) adapted for passive tracers 9 9 !! 3.3 ! 2010-05 (C. Ethe, G. Madec) merge TRC-TRA 10 !! 3.6 ! 2014-11 (P. Mathiot) Add zps_hde_isf (needed to open a cavity) 10 11 !!====================================================================== 11 12 … … 27 28 PRIVATE 28 29 29 PUBLIC zps_hde ! routine called by step.F90 30 PUBLIC zps_hde ! routine called by step.F90 31 PUBLIC zps_hde_isf ! routine called by step.F90 30 32 31 33 !! * Substitutions … … 40 42 41 43 SUBROUTINE zps_hde( kt, kjpt, pta, pgtu, pgtv, & 42 44 & prd, pgru, pgrv ) 43 45 !!---------------------------------------------------------------------- 44 46 !! *** ROUTINE zps_hde *** … … 74 76 !! Idem for di(s) and dj(s) 75 77 !! 76 !! For rho, we call eos _insitu_2d which will compute rd~(t~,s~) at77 !! the good depth zh from interpolated T and S for the different78 !! formulationof the equation of state (eos).78 !! For rho, we call eos which will compute rd~(t~,s~) at the right 79 !! depth zh from interpolated T and S for the different formulations 80 !! of the equation of state (eos). 79 81 !! Gradient formulation for rho : 80 !! di(rho) = rd~ - rd(i,j,k) orrd(i+1,j,k) - rd~82 !! di(rho) = rd~ - rd(i,j,k) or rd(i+1,j,k) - rd~ 81 83 !! 82 !! ** Action : - pgtu, pgtv: horizontal gradient of tracer at u- & v-points 83 !! - pgru, pgrv: horizontal gradient of rho (if present) at u- & v-points 84 !! ** Action : compute for top interfaces 85 !! - pgtu, pgtv: horizontal gradient of tracer at u- & v-points 86 !! - pgru, pgrv: horizontal gradient of rho (if present) at u- & v-points 84 87 !!---------------------------------------------------------------------- 85 !86 88 INTEGER , INTENT(in ) :: kt ! ocean time-step index 87 89 INTEGER , INTENT(in ) :: kjpt ! number of tracers … … 89 91 REAL(wp), DIMENSION(jpi,jpj, kjpt), INTENT( out) :: pgtu, pgtv ! hor. grad. of ptra at u- & v-pts 90 92 REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ), OPTIONAL :: prd ! 3D density anomaly fields 91 REAL(wp), DIMENSION(jpi,jpj ), INTENT( out), OPTIONAL :: pgru, pgrv ! hor. grad . of prd at u- & v-pts93 REAL(wp), DIMENSION(jpi,jpj ), INTENT( out), OPTIONAL :: pgru, pgrv ! hor. grad of prd at u- & v-pts (bottom) 92 94 ! 93 95 INTEGER :: ji, jj, jn ! Dummy loop indices 94 96 INTEGER :: iku, ikv, ikum1, ikvm1 ! partial step level (ocean bottom level) at u- and v-points 95 97 REAL(wp) :: ze3wu, ze3wv, zmaxu, zmaxv ! temporary scalars 96 REAL(wp), POINTER, DIMENSION(:,: ) :: zri, zrj, zhi, zhj97 REAL(wp), POINTER, DIMENSION(:,:,:) :: zti, ztj ! interpolated value of tracer98 REAL(wp), DIMENSION(jpi,jpj) :: zri, zrj, zhi, zhj ! NB: 3rd dim=1 to use eos 99 REAL(wp), DIMENSION(jpi,jpj,kjpt) :: zti, ztj ! 98 100 !!---------------------------------------------------------------------- 99 101 ! 100 102 IF( nn_timing == 1 ) CALL timing_start( 'zps_hde') 101 103 ! 102 CALL wrk_alloc( jpi, jpj, zri, zrj, zhi, zhj ) 103 CALL wrk_alloc( jpi, jpj, kjpt, zti, ztj ) 104 pgtu(:,:,:)=0.0_wp ; pgtv(:,:,:)=0.0_wp ; 105 zti (:,:,:)=0.0_wp ; ztj (:,:,:)=0.0_wp ; 106 zhi (:,: )=0.0_wp ; zhj (:,: )=0.0_wp ; 104 107 ! 105 108 DO jn = 1, kjpt !== Interpolation of tracers at the last ocean level ==! 106 109 ! 107 # if defined key_vectopt_loop 108 jj = 1 109 DO ji = 1, jpij-jpi ! vector opt. (forced unrolled) 110 # else 111 DO jj = 1, jpjm1 112 DO ji = 1, jpim1 113 # endif 110 DO jj = 1, jpjm1 111 DO ji = 1, jpim1 114 112 iku = mbku(ji,jj) ; ikum1 = MAX( iku - 1 , 1 ) ! last and before last ocean level at u- & v-points 115 113 ikv = mbkv(ji,jj) ; ikvm1 = MAX( ikv - 1 , 1 ) ! if level first is a p-step, ik.m1=1 … … 121 119 zmaxu = ze3wu / fse3w(ji+1,jj,iku) 122 120 ! interpolated values of tracers 123 zti (ji,jj,jn) = pta(ji+1,jj,iku,jn) + zmaxu * ( pta(ji+1,jj,ikum1,jn) - pta(ji+1,jj,iku,jn) )121 zti (ji,jj,jn) = pta(ji+1,jj,iku,jn) + zmaxu * ( pta(ji+1,jj,ikum1,jn) - pta(ji+1,jj,iku,jn) ) 124 122 ! gradient of tracers 125 123 pgtu(ji,jj,jn) = umask(ji,jj,1) * ( zti(ji,jj,jn) - pta(ji,jj,iku,jn) ) … … 127 125 zmaxu = -ze3wu / fse3w(ji,jj,iku) 128 126 ! interpolated values of tracers 129 zti (ji,jj,jn) = pta(ji,jj,iku,jn) + zmaxu * ( pta(ji,jj,ikum1,jn) - pta(ji,jj,iku,jn) )127 zti (ji,jj,jn) = pta(ji,jj,iku,jn) + zmaxu * ( pta(ji,jj,ikum1,jn) - pta(ji,jj,iku,jn) ) 130 128 ! gradient of tracers 131 129 pgtu(ji,jj,jn) = umask(ji,jj,1) * ( pta(ji+1,jj,iku,jn) - zti(ji,jj,jn) ) … … 136 134 zmaxv = ze3wv / fse3w(ji,jj+1,ikv) 137 135 ! interpolated values of tracers 138 ztj (ji,jj,jn) = pta(ji,jj+1,ikv,jn) + zmaxv * ( pta(ji,jj+1,ikvm1,jn) - pta(ji,jj+1,ikv,jn) )136 ztj (ji,jj,jn) = pta(ji,jj+1,ikv,jn) + zmaxv * ( pta(ji,jj+1,ikvm1,jn) - pta(ji,jj+1,ikv,jn) ) 139 137 ! gradient of tracers 140 138 pgtv(ji,jj,jn) = vmask(ji,jj,1) * ( ztj(ji,jj,jn) - pta(ji,jj,ikv,jn) ) … … 142 140 zmaxv = -ze3wv / fse3w(ji,jj,ikv) 143 141 ! interpolated values of tracers 144 ztj (ji,jj,jn) = pta(ji,jj,ikv,jn) + zmaxv * ( pta(ji,jj,ikvm1,jn) - pta(ji,jj,ikv,jn) )142 ztj (ji,jj,jn) = pta(ji,jj,ikv,jn) + zmaxv * ( pta(ji,jj,ikvm1,jn) - pta(ji,jj,ikv,jn) ) 145 143 ! gradient of tracers 146 144 pgtv(ji,jj,jn) = vmask(ji,jj,1) * ( pta(ji,jj+1,ikv,jn) - ztj(ji,jj,jn) ) 147 145 ENDIF 148 # if ! defined key_vectopt_loop 149 END DO 150 # endif 146 END DO 151 147 END DO 152 148 CALL lbc_lnk( pgtu(:,:,jn), 'U', -1. ) ; CALL lbc_lnk( pgtv(:,:,jn), 'V', -1. ) ! Lateral boundary cond. … … 156 152 ! horizontal derivative of density anomalies (rd) 157 153 IF( PRESENT( prd ) ) THEN ! depth of the partial step level 158 # if defined key_vectopt_loop 159 jj = 1 160 DO ji = 1, jpij-jpi ! vector opt. (forced unrolled) 161 # else 162 DO jj = 1, jpjm1 163 DO ji = 1, jpim1 164 # endif 154 pgru(:,:)=0.0_wp ; pgrv(:,:)=0.0_wp ; 155 DO jj = 1, jpjm1 156 DO ji = 1, jpim1 165 157 iku = mbku(ji,jj) 166 158 ikv = mbkv(ji,jj) … … 173 165 ELSE ; zhj(ji,jj) = fsdept(ji,jj+1,ikv) ! - - case 2 174 166 ENDIF 175 # if ! defined key_vectopt_loop 176 END DO 177 # endif 167 END DO 178 168 END DO 179 169 … … 184 174 185 175 ! Gradient of density at the last level 186 # if defined key_vectopt_loop 187 jj = 1 188 DO ji = 1, jpij-jpi ! vector opt. (forced unrolled) 189 # else 190 DO jj = 1, jpjm1 191 DO ji = 1, jpim1 192 # endif 176 DO jj = 1, jpjm1 177 DO ji = 1, jpim1 193 178 iku = mbku(ji,jj) 194 179 ikv = mbkv(ji,jj) 195 180 ze3wu = fse3w(ji+1,jj ,iku) - fse3w(ji,jj,iku) 196 181 ze3wv = fse3w(ji ,jj+1,ikv) - fse3w(ji,jj,ikv) 197 IF( ze3wu >= 0._wp ) THEN ; pgru(ji,jj) = umask(ji,jj,1) * ( zri(ji ,jj) - prd(ji,jj,iku) ) ! i: 1 198 ELSE ; pgru(ji,jj) = umask(ji,jj,1) * ( prd(ji+1,jj,iku) - zri(ji,jj) ) ! i: 2 199 ENDIF 200 IF( ze3wv >= 0._wp ) THEN ; pgrv(ji,jj) = vmask(ji,jj,1) * ( zrj(ji,jj ) - prd(ji,jj,ikv) ) ! j: 1 201 ELSE ; pgrv(ji,jj) = vmask(ji,jj,1) * ( prd(ji,jj+1,ikv) - zrj(ji,jj) ) ! j: 2 202 ENDIF 203 # if ! defined key_vectopt_loop 204 END DO 205 # endif 182 IF( ze3wu >= 0._wp ) THEN ; pgru(ji,jj) = umask(ji,jj,1) * ( zri(ji ,jj ) - prd(ji,jj,iku) ) ! i: 1 183 ELSE ; pgru(ji,jj) = umask(ji,jj,1) * ( prd(ji+1,jj,iku) - zri(ji,jj ) ) ! i: 2 184 ENDIF 185 IF( ze3wv >= 0._wp ) THEN ; pgrv(ji,jj) = vmask(ji,jj,1) * ( zrj(ji,jj ) - prd(ji,jj,ikv) ) ! j: 1 186 ELSE ; pgrv(ji,jj) = vmask(ji,jj,1) * ( prd(ji,jj+1,ikv) - zrj(ji,jj ) ) ! j: 2 187 ENDIF 188 END DO 206 189 END DO 207 190 CALL lbc_lnk( pgru , 'U', -1. ) ; CALL lbc_lnk( pgrv , 'V', -1. ) ! Lateral boundary conditions … … 209 192 END IF 210 193 ! 211 CALL wrk_dealloc( jpi, jpj, zri, zrj, zhi, zhj )212 CALL wrk_dealloc( jpi, jpj, kjpt, zti, ztj )213 !214 194 IF( nn_timing == 1 ) CALL timing_stop( 'zps_hde') 215 195 ! 216 196 END SUBROUTINE zps_hde 217 197 ! 198 SUBROUTINE zps_hde_isf( kt, kjpt, pta, pgtu, pgtv, & 199 & prd, pgru, pgrv, pmru, pmrv, pgzu, pgzv, pge3ru, pge3rv, & 200 & pgtui, pgtvi, pgrui, pgrvi, pmrui, pmrvi, pgzui, pgzvi, pge3rui, pge3rvi ) 201 !!---------------------------------------------------------------------- 202 !! *** ROUTINE zps_hde *** 203 !! 204 !! ** Purpose : Compute the horizontal derivative of T, S and rho 205 !! at u- and v-points with a linear interpolation for z-coordinate 206 !! with partial steps. 207 !! 208 !! ** Method : In z-coord with partial steps, scale factors on last 209 !! levels are different for each grid point, so that T, S and rd 210 !! points are not at the same depth as in z-coord. To have horizontal 211 !! gradients again, we interpolate T and S at the good depth : 212 !! Linear interpolation of T, S 213 !! Computation of di(tb) and dj(tb) by vertical interpolation: 214 !! di(t) = t~ - t(i,j,k) or t(i+1,j,k) - t~ 215 !! dj(t) = t~ - t(i,j,k) or t(i,j+1,k) - t~ 216 !! This formulation computes the two cases: 217 !! CASE 1 CASE 2 218 !! k-1 ___ ___________ k-1 ___ ___________ 219 !! Ti T~ T~ Ti+1 220 !! _____ _____ 221 !! k | |Ti+1 k Ti | | 222 !! | |____ ____| | 223 !! ___ | | | ___ | | | 224 !! 225 !! case 1-> e3w(i+1) >= e3w(i) ( and e3w(j+1) >= e3w(j) ) then 226 !! t~ = t(i+1,j ,k) + (e3w(i+1) - e3w(i)) * dk(Ti+1)/e3w(i+1) 227 !! ( t~ = t(i ,j+1,k) + (e3w(j+1) - e3w(j)) * dk(Tj+1)/e3w(j+1) ) 228 !! or 229 !! case 2-> e3w(i+1) <= e3w(i) ( and e3w(j+1) <= e3w(j) ) then 230 !! t~ = t(i,j,k) + (e3w(i) - e3w(i+1)) * dk(Ti)/e3w(i ) 231 !! ( t~ = t(i,j,k) + (e3w(j) - e3w(j+1)) * dk(Tj)/e3w(j ) ) 232 !! Idem for di(s) and dj(s) 233 !! 234 !! For rho, we call eos which will compute rd~(t~,s~) at the right 235 !! depth zh from interpolated T and S for the different formulations 236 !! of the equation of state (eos). 237 !! Gradient formulation for rho : 238 !! di(rho) = rd~ - rd(i,j,k) or rd(i+1,j,k) - rd~ 239 !! 240 !! ** Action : compute for top and bottom interfaces 241 !! - pgtu, pgtv, pgtui, pgtvi: horizontal gradient of tracer at u- & v-points 242 !! - pgru, pgrv, pgrui, pgtvi: horizontal gradient of rho (if present) at u- & v-points 243 !! - pmru, pmrv, pmrui, pmrvi: horizontal sum of rho at u- & v- point (used in dynhpg with vvl) 244 !! - pgzu, pgzv, pgzui, pgzvi: horizontal gradient of z at u- and v- point (used in dynhpg with vvl) 245 !! - pge3ru, pge3rv, pge3rui, pge3rvi: horizontal gradient of rho weighted by local e3w at u- & v-points 246 !!---------------------------------------------------------------------- 247 INTEGER , INTENT(in ) :: kt ! ocean time-step index 248 INTEGER , INTENT(in ) :: kjpt ! number of tracers 249 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: pta ! 4D tracers fields 250 REAL(wp), DIMENSION(jpi,jpj, kjpt), INTENT( out) :: pgtu, pgtv ! hor. grad. of ptra at u- & v-pts 251 REAL(wp), DIMENSION(jpi,jpj, kjpt), INTENT( out) :: pgtui, pgtvi ! hor. grad. of stra at u- & v-pts (ISF) 252 REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ), OPTIONAL :: prd ! 3D density anomaly fields 253 REAL(wp), DIMENSION(jpi,jpj ), INTENT( out), OPTIONAL :: pgru, pgrv ! hor. grad of prd at u- & v-pts (bottom) 254 REAL(wp), DIMENSION(jpi,jpj ), INTENT( out), OPTIONAL :: pmru, pmrv ! hor. sum of prd at u- & v-pts (bottom) 255 REAL(wp), DIMENSION(jpi,jpj ), INTENT( out), OPTIONAL :: pgzu, pgzv ! hor. grad of z at u- & v-pts (bottom) 256 REAL(wp), DIMENSION(jpi,jpj ), INTENT( out), OPTIONAL :: pge3ru, pge3rv ! hor. grad of prd weighted by local e3w at u- & v-pts (bottom) 257 REAL(wp), DIMENSION(jpi,jpj ), INTENT( out), OPTIONAL :: pgrui, pgrvi ! hor. grad of prd at u- & v-pts (top) 258 REAL(wp), DIMENSION(jpi,jpj ), INTENT( out), OPTIONAL :: pmrui, pmrvi ! hor. sum of prd at u- & v-pts (top) 259 REAL(wp), DIMENSION(jpi,jpj ), INTENT( out), OPTIONAL :: pgzui, pgzvi ! hor. grad of z at u- & v-pts (top) 260 REAL(wp), DIMENSION(jpi,jpj ), INTENT( out), OPTIONAL :: pge3rui, pge3rvi ! hor. grad of prd weighted by local e3w at u- & v-pts (top) 261 ! 262 INTEGER :: ji, jj, jn ! Dummy loop indices 263 INTEGER :: iku, ikv, ikum1, ikvm1,ikup1, ikvp1 ! partial step level (ocean bottom level) at u- and v-points 264 REAL(wp) :: ze3wu, ze3wv, zmaxu, zmaxv, zdzwu, zdzwv, zdzwuip1, zdzwvjp1 ! temporary scalars 265 REAL(wp), DIMENSION(jpi,jpj) :: zri, zrj, zhi, zhj ! NB: 3rd dim=1 to use eos 266 REAL(wp), DIMENSION(jpi,jpj,kjpt) :: zti, ztj ! 267 !!---------------------------------------------------------------------- 268 ! 269 IF( nn_timing == 1 ) CALL timing_start( 'zps_hde_isf') 270 ! 271 pgtu(:,:,:)=0.0_wp ; pgtv(:,:,:)=0.0_wp ; 272 pgtui(:,:,:)=0.0_wp ; pgtvi(:,:,:)=0.0_wp ; 273 zti (:,:,:)=0.0_wp ; ztj (:,:,:)=0.0_wp ; 274 zhi (:,: )=0.0_wp ; zhj (:,: )=0.0_wp ; 275 ! 276 DO jn = 1, kjpt !== Interpolation of tracers at the last ocean level ==! 277 ! 278 DO jj = 1, jpjm1 279 DO ji = 1, jpim1 280 iku = mbku(ji,jj) ; ikum1 = MAX( iku - 1 , 1 ) ! last and before last ocean level at u- & v-points 281 ikv = mbkv(ji,jj) ; ikvm1 = MAX( ikv - 1 , 1 ) ! if level first is a p-step, ik.m1=1 282 ! (ISF) case partial step top and bottom in adjacent cell in vertical 283 ! cannot used e3w because if 2 cell water column, we have ps at top and bottom 284 ! in this case e3w(i,j) - e3w(i,j+1) is not the distance between Tj~ and Tj 285 ! the only common depth between cells (i,j) and (i,j+1) is gdepw_0 286 ze3wu = (gdept_0(ji+1,jj,iku) - gdepw_0(ji+1,jj,iku)) - (gdept_0(ji,jj,iku) - gdepw_0(ji,jj,iku)) 287 ze3wv = (gdept_0(ji,jj+1,ikv) - gdepw_0(ji,jj+1,ikv)) - (gdept_0(ji,jj,ikv) - gdepw_0(ji,jj,ikv)) 288 ! 289 ! i- direction 290 IF( ze3wu >= 0._wp ) THEN ! case 1 291 zmaxu = ze3wu / fse3w(ji+1,jj,iku) 292 ! interpolated values of tracers 293 zti (ji,jj,jn) = pta(ji+1,jj,iku,jn) + zmaxu * ( pta(ji+1,jj,ikum1,jn) - pta(ji+1,jj,iku,jn) ) 294 ! gradient of tracers 295 pgtu(ji,jj,jn) = umask(ji,jj,iku) * ( zti(ji,jj,jn) - pta(ji,jj,iku,jn) ) 296 ELSE ! case 2 297 zmaxu = -ze3wu / fse3w(ji,jj,iku) 298 ! interpolated values of tracers 299 zti (ji,jj,jn) = pta(ji,jj,iku,jn) + zmaxu * ( pta(ji,jj,ikum1,jn) - pta(ji,jj,iku,jn) ) 300 ! gradient of tracers 301 pgtu(ji,jj,jn) = umask(ji,jj,iku) * ( pta(ji+1,jj,iku,jn) - zti(ji,jj,jn) ) 302 ENDIF 303 ! 304 ! j- direction 305 IF( ze3wv >= 0._wp ) THEN ! case 1 306 zmaxv = ze3wv / fse3w(ji,jj+1,ikv) 307 ! interpolated values of tracers 308 ztj (ji,jj,jn) = pta(ji,jj+1,ikv,jn) + zmaxv * ( pta(ji,jj+1,ikvm1,jn) - pta(ji,jj+1,ikv,jn) ) 309 ! gradient of tracers 310 pgtv(ji,jj,jn) = vmask(ji,jj,ikv) * ( ztj(ji,jj,jn) - pta(ji,jj,ikv,jn) ) 311 ELSE ! case 2 312 zmaxv = -ze3wv / fse3w(ji,jj,ikv) 313 ! interpolated values of tracers 314 ztj (ji,jj,jn) = pta(ji,jj,ikv,jn) + zmaxv * ( pta(ji,jj,ikvm1,jn) - pta(ji,jj,ikv,jn) ) 315 ! gradient of tracers 316 pgtv(ji,jj,jn) = vmask(ji,jj,ikv) * ( pta(ji,jj+1,ikv,jn) - ztj(ji,jj,jn) ) 317 ENDIF 318 END DO 319 END DO 320 CALL lbc_lnk( pgtu(:,:,jn), 'U', -1. ) ; CALL lbc_lnk( pgtv(:,:,jn), 'V', -1. ) ! Lateral boundary cond. 321 ! 322 END DO 323 324 ! horizontal derivative of density anomalies (rd) 325 IF( PRESENT( prd ) ) THEN ! depth of the partial step level 326 pgru(:,:)=0.0_wp ; pgrv(:,:)=0.0_wp ; 327 pgzu(:,:)=0.0_wp ; pgzv(:,:)=0.0_wp ; 328 pmru(:,:)=0.0_wp ; pmru(:,:)=0.0_wp ; 329 pge3ru(:,:)=0.0_wp ; pge3rv(:,:)=0.0_wp ; 330 DO jj = 1, jpjm1 331 DO ji = 1, jpim1 332 iku = mbku(ji,jj) 333 ikv = mbkv(ji,jj) 334 ze3wu = (gdept_0(ji+1,jj,iku) - gdepw_0(ji+1,jj,iku)) - (gdept_0(ji,jj,iku) - gdepw_0(ji,jj,iku)) 335 ze3wv = (gdept_0(ji,jj+1,ikv) - gdepw_0(ji,jj+1,ikv)) - (gdept_0(ji,jj,ikv) - gdepw_0(ji,jj,ikv)) 336 337 IF( ze3wu >= 0._wp ) THEN ; zhi(ji,jj) = fsdept(ji+1,jj,iku) - ze3wu ! i-direction: case 1 338 ELSE ; zhi(ji,jj) = fsdept(ji ,jj,iku) + ze3wu ! - - case 2 339 ENDIF 340 IF( ze3wv >= 0._wp ) THEN ; zhj(ji,jj) = fsdept(ji,jj+1,ikv) - ze3wv ! j-direction: case 1 341 ELSE ; zhj(ji,jj) = fsdept(ji,jj ,ikv) + ze3wv ! - - case 2 342 ENDIF 343 END DO 344 END DO 345 346 ! Compute interpolated rd from zti, ztj for the 2 cases at the depth of the partial 347 ! step and store it in zri, zrj for each case 348 CALL eos( zti, zhi, zri ) 349 CALL eos( ztj, zhj, zrj ) 350 351 ! Gradient of density at the last level 352 DO jj = 1, jpjm1 353 DO ji = 1, jpim1 354 iku = mbku(ji,jj) ; ikum1 = MAX( iku - 1 , 1 ) ! last and before last ocean level at u- & v-points 355 ikv = mbkv(ji,jj) ; ikvm1 = MAX( ikv - 1 , 1 ) ! last and before last ocean level at u- & v-points 356 ze3wu = (gdept_0(ji+1,jj,iku) - gdepw_0(ji+1,jj,iku)) - (gdept_0(ji,jj,iku) - gdepw_0(ji,jj,iku)) 357 ze3wv = (gdept_0(ji,jj+1,ikv) - gdepw_0(ji,jj+1,ikv)) - (gdept_0(ji,jj,ikv) - gdepw_0(ji,jj,ikv)) 358 IF( ze3wu >= 0._wp ) THEN 359 pgzu(ji,jj) = (fsde3w(ji+1,jj,iku) - ze3wu) - fsde3w(ji,jj,iku) 360 pgru(ji,jj) = umask(ji,jj,iku) * ( zri(ji ,jj) - prd(ji,jj,iku) ) ! i: 1 361 pmru(ji,jj) = umask(ji,jj,iku) * ( zri(ji ,jj) + prd(ji,jj,iku) ) ! i: 1 362 pge3ru(ji,jj) = umask(ji,jj,iku) & 363 * ( (fse3w(ji+1,jj,iku) - ze3wu )* ( zri(ji ,jj ) + prd(ji+1,jj,ikum1) + 2._wp) & 364 - fse3w(ji ,jj,iku) * ( prd(ji ,jj,iku) + prd(ji ,jj,ikum1) + 2._wp) ) ! j: 2 365 ELSE 366 pgzu(ji,jj) = fsde3w(ji+1,jj,iku) - (fsde3w(ji,jj,iku) + ze3wu) 367 pgru(ji,jj) = umask(ji,jj,iku) * ( prd(ji+1,jj,iku) - zri(ji,jj) ) ! i: 2 368 pmru(ji,jj) = umask(ji,jj,iku) * ( prd(ji+1,jj,iku) + zri(ji,jj) ) ! i: 2 369 pge3ru(ji,jj) = umask(ji,jj,iku) & 370 * ( fse3w(ji+1,jj,iku) * ( prd(ji+1,jj,iku) + prd(ji+1,jj,ikum1) + 2._wp) & 371 -(fse3w(ji ,jj,iku) + ze3wu) * ( zri(ji ,jj ) + prd(ji ,jj,ikum1) + 2._wp) ) ! j: 2 372 ENDIF 373 IF( ze3wv >= 0._wp ) THEN 374 pgzv(ji,jj) = (fsde3w(ji,jj+1,ikv) - ze3wv) - fsde3w(ji,jj,ikv) 375 pgrv(ji,jj) = vmask(ji,jj,ikv) * ( zrj(ji,jj ) - prd(ji,jj,ikv) ) ! j: 1 376 pmrv(ji,jj) = vmask(ji,jj,ikv) * ( zrj(ji,jj ) + prd(ji,jj,ikv) ) ! j: 1 377 pge3rv(ji,jj) = vmask(ji,jj,ikv) & 378 * ( (fse3w(ji,jj+1,ikv) - ze3wv )* ( zrj(ji,jj ) + prd(ji,jj+1,ikvm1) + 2._wp) & 379 - fse3w(ji,jj ,ikv) * ( prd(ji,jj ,ikv) + prd(ji,jj ,ikvm1) + 2._wp) ) ! j: 2 380 ELSE 381 pgzv(ji,jj) = fsde3w(ji,jj+1,ikv) - (fsde3w(ji,jj,ikv) + ze3wv) 382 pgrv(ji,jj) = vmask(ji,jj,ikv) * ( prd(ji,jj+1,ikv) - zrj(ji,jj) ) ! j: 2 383 pmrv(ji,jj) = vmask(ji,jj,ikv) * ( prd(ji,jj+1,ikv) + zrj(ji,jj) ) ! j: 2 384 pge3rv(ji,jj) = vmask(ji,jj,ikv) & 385 * ( fse3w(ji,jj+1,ikv) * ( prd(ji,jj+1,ikv) + prd(ji,jj+1,ikvm1) + 2._wp) & 386 -(fse3w(ji,jj ,ikv) + ze3wv) * ( zrj(ji,jj ) + prd(ji,jj ,ikvm1) + 2._wp) ) ! j: 2 387 ENDIF 388 END DO 389 END DO 390 CALL lbc_lnk( pgru , 'U', -1. ) ; CALL lbc_lnk( pgrv , 'V', -1. ) ! Lateral boundary conditions 391 CALL lbc_lnk( pmru , 'U', 1. ) ; CALL lbc_lnk( pmrv , 'V', 1. ) ! Lateral boundary conditions 392 CALL lbc_lnk( pgzu , 'U', -1. ) ; CALL lbc_lnk( pgzv , 'V', -1. ) ! Lateral boundary conditions 393 CALL lbc_lnk( pge3ru , 'U', -1. ) ; CALL lbc_lnk( pge3rv , 'V', -1. ) ! Lateral boundary conditions 394 ! 395 END IF 396 ! (ISH) compute grui and gruvi 397 DO jn = 1, kjpt !== Interpolation of tracers at the last ocean level ==! ! 398 DO jj = 1, jpjm1 399 DO ji = 1, jpim1 400 iku = miku(ji,jj) ; ikup1 = miku(ji,jj) + 1 401 ikv = mikv(ji,jj) ; ikvp1 = mikv(ji,jj) + 1 402 ! 403 ! (ISF) case partial step top and bottom in adjacent cell in vertical 404 ! cannot used e3w because if 2 cell water column, we have ps at top and bottom 405 ! in this case e3w(i,j) - e3w(i,j+1) is not the distance between Tj~ and Tj 406 ! the only common depth between cells (i,j) and (i,j+1) is gdepw_0 407 ze3wu = (gdepw_0(ji+1,jj,iku+1) - gdept_0(ji+1,jj,iku)) - (gdepw_0(ji,jj,iku+1) - gdept_0(ji,jj,iku)) 408 ze3wv = (gdepw_0(ji,jj+1,ikv+1) - gdept_0(ji,jj+1,ikv)) - (gdepw_0(ji,jj,ikv+1) - gdept_0(ji,jj,ikv)) 409 ! i- direction 410 IF( ze3wu >= 0._wp ) THEN ! case 1 411 zmaxu = ze3wu / fse3w(ji+1,jj,iku+1) 412 ! interpolated values of tracers 413 zti(ji,jj,jn) = pta(ji+1,jj,iku,jn) + zmaxu * ( pta(ji+1,jj,iku+1,jn) - pta(ji+1,jj,iku,jn) ) 414 ! gradient of tracers 415 pgtui(ji,jj,jn) = umask(ji,jj,iku) * ( zti(ji,jj,jn) - pta(ji,jj,iku,jn) ) 416 ELSE ! case 2 417 zmaxu = - ze3wu / fse3w(ji,jj,iku+1) 418 ! interpolated values of tracers 419 zti(ji,jj,jn) = pta(ji,jj,iku,jn) + zmaxu * ( pta(ji,jj,iku+1,jn) - pta(ji,jj,iku,jn) ) 420 ! gradient of tracers 421 pgtui(ji,jj,jn) = umask(ji,jj,iku) * ( pta(ji+1,jj,iku,jn) - zti(ji,jj,jn) ) 422 ENDIF 423 ! 424 ! j- direction 425 IF( ze3wv >= 0._wp ) THEN ! case 1 426 zmaxv = ze3wv / fse3w(ji,jj+1,ikv+1) 427 ! interpolated values of tracers 428 ztj(ji,jj,jn) = pta(ji,jj+1,ikv,jn) + zmaxv * ( pta(ji,jj+1,ikv+1,jn) - pta(ji,jj+1,ikv,jn) ) 429 ! gradient of tracers 430 pgtvi(ji,jj,jn) = vmask(ji,jj,ikv) * ( ztj(ji,jj,jn) - pta(ji,jj,ikv,jn) ) 431 ELSE ! case 2 432 zmaxv = - ze3wv / fse3w(ji,jj,ikv+1) 433 ! interpolated values of tracers 434 ztj(ji,jj,jn) = pta(ji,jj,ikv,jn) + zmaxv * ( pta(ji,jj,ikv+1,jn) - pta(ji,jj,ikv,jn) ) 435 ! gradient of tracers 436 pgtvi(ji,jj,jn) = vmask(ji,jj,ikv) * ( pta(ji,jj+1,ikv,jn) - ztj(ji,jj,jn) ) 437 ENDIF 438 END DO!! 439 END DO!! 440 CALL lbc_lnk( pgtui(:,:,jn), 'U', -1. ) ; CALL lbc_lnk( pgtvi(:,:,jn), 'V', -1. ) ! Lateral boundary cond. 441 ! 442 END DO 443 444 ! horizontal derivative of density anomalies (rd) 445 IF( PRESENT( prd ) ) THEN ! depth of the partial step level 446 pgrui(:,:) =0.0_wp ; pgrvi(:,:) =0.0_wp ; 447 pgzui(:,:) =0.0_wp ; pgzvi(:,:) =0.0_wp ; 448 pmrui(:,:) =0.0_wp ; pmrui(:,:) =0.0_wp ; 449 pge3rui(:,:)=0.0_wp ; pge3rvi(:,:)=0.0_wp ; 450 451 DO jj = 1, jpjm1 452 DO ji = 1, jpim1 453 iku = miku(ji,jj) 454 ikv = mikv(ji,jj) 455 ze3wu = (gdepw_0(ji+1,jj,iku+1) - gdept_0(ji+1,jj,iku)) - (gdepw_0(ji,jj,iku+1) - gdept_0(ji,jj,iku)) 456 ze3wv = (gdepw_0(ji,jj+1,ikv+1) - gdept_0(ji,jj+1,ikv)) - (gdepw_0(ji,jj,ikv+1) - gdept_0(ji,jj,ikv)) 457 458 IF( ze3wu >= 0._wp ) THEN ; zhi(ji,jj) = fsdept(ji+1,jj,iku) + ze3wu ! i-direction: case 1 459 ELSE ; zhi(ji,jj) = fsdept(ji ,jj,iku) - ze3wu ! - - case 2 460 ENDIF 461 IF( ze3wv >= 0._wp ) THEN ; zhj(ji,jj) = fsdept(ji,jj+1,ikv) + ze3wv ! j-direction: case 1 462 ELSE ; zhj(ji,jj) = fsdept(ji,jj ,ikv) - ze3wv ! - - case 2 463 ENDIF 464 END DO 465 END DO 466 467 ! Compute interpolated rd from zti, ztj for the 2 cases at the depth of the partial 468 ! step and store it in zri, zrj for each case 469 CALL eos( zti, zhi, zri ) 470 CALL eos( ztj, zhj, zrj ) 471 472 ! Gradient of density at the last level 473 DO jj = 1, jpjm1 474 DO ji = 1, jpim1 475 iku = miku(ji,jj) ; ikup1 = miku(ji,jj) + 1 476 ikv = mikv(ji,jj) ; ikvp1 = mikv(ji,jj) + 1 477 ze3wu = (gdepw_0(ji+1,jj,iku+1) - gdept_0(ji+1,jj,iku)) - (gdepw_0(ji,jj,iku+1) - gdept_0(ji,jj,iku)) 478 ze3wv = (gdepw_0(ji,jj+1,ikv+1) - gdept_0(ji,jj+1,ikv)) - (gdepw_0(ji,jj,ikv+1) - gdept_0(ji,jj,ikv)) 479 IF( ze3wu >= 0._wp ) THEN 480 pgzui (ji,jj) = (fsde3w(ji+1,jj,iku) + ze3wu) - fsde3w(ji,jj,iku) 481 pgrui (ji,jj) = umask(ji,jj,iku) * ( zri(ji,jj) - prd(ji,jj,iku) ) ! i: 1 482 pmrui (ji,jj) = umask(ji,jj,iku) * ( zri(ji,jj) + prd(ji,jj,iku) ) ! i: 1 483 pge3rui(ji,jj) = umask(ji,jj,iku+1) & 484 * ( (fse3w(ji+1,jj,iku+1) - ze3wu) * (zri(ji,jj ) + prd(ji+1,jj,iku+1) + 2._wp) & 485 - fse3w(ji ,jj,iku+1) * (prd(ji,jj,iku) + prd(ji ,jj,iku+1) + 2._wp) ) ! i: 1 486 ELSE 487 pgzui (ji,jj) = fsde3w(ji+1,jj,iku) - (fsde3w(ji,jj,iku) - ze3wu) 488 pgrui (ji,jj) = umask(ji,jj,iku) * ( prd(ji+1,jj,iku) - zri(ji,jj) ) ! i: 2 489 pmrui (ji,jj) = umask(ji,jj,iku) * ( prd(ji+1,jj,iku) + zri(ji,jj) ) ! i: 2 490 pge3rui(ji,jj) = umask(ji,jj,iku+1) & 491 * ( fse3w(ji+1,jj,iku+1) * (prd(ji+1,jj,iku) + prd(ji+1,jj,iku+1) + 2._wp) & 492 -(fse3w(ji ,jj,iku+1) + ze3wu) * (zri(ji,jj ) + prd(ji ,jj,iku+1) + 2._wp) ) ! i: 2 493 ENDIF 494 IF( ze3wv >= 0._wp ) THEN 495 pgzvi (ji,jj) = (fsde3w(ji,jj+1,ikv) + ze3wv) - fsde3w(ji,jj,ikv) 496 pgrvi (ji,jj) = vmask(ji,jj,ikv) * ( zrj(ji,jj ) - prd(ji,jj,ikv) ) ! j: 1 497 pmrvi (ji,jj) = vmask(ji,jj,ikv) * ( zrj(ji,jj ) + prd(ji,jj,ikv) ) ! j: 1 498 pge3rvi(ji,jj) = vmask(ji,jj,ikv+1) & 499 * ( (fse3w(ji,jj+1,ikv+1) - ze3wv) * ( zrj(ji,jj ) + prd(ji,jj+1,ikv+1) + 2._wp) & 500 - fse3w(ji,jj ,ikv+1) * ( prd(ji,jj,ikv) + prd(ji,jj ,ikv+1) + 2._wp) ) ! j: 1 501 ! + 2 due to the formulation in density and not in anomalie in hpg sco 502 ELSE 503 pgzvi (ji,jj) = fsde3w(ji,jj+1,ikv) - (fsde3w(ji,jj,ikv) - ze3wv) 504 pgrvi (ji,jj) = vmask(ji,jj,ikv) * ( prd(ji,jj+1,ikv) - zrj(ji,jj) ) ! j: 2 505 pmrvi (ji,jj) = vmask(ji,jj,ikv) * ( prd(ji,jj+1,ikv) + zrj(ji,jj) ) ! j: 2 506 pge3rvi(ji,jj) = vmask(ji,jj,ikv+1) & 507 * ( fse3w(ji,jj+1,ikv+1) * ( prd(ji,jj+1,ikv) + prd(ji,jj+1,ikv+1) + 2._wp) & 508 -(fse3w(ji,jj ,ikv+1) + ze3wv) * ( zrj(ji,jj ) + prd(ji,jj ,ikv+1) + 2._wp) ) ! j: 2 509 ENDIF 510 END DO 511 END DO 512 CALL lbc_lnk( pgrui , 'U', -1. ) ; CALL lbc_lnk( pgrvi , 'V', -1. ) ! Lateral boundary conditions 513 CALL lbc_lnk( pmrui , 'U', 1. ) ; CALL lbc_lnk( pmrvi , 'V', 1. ) ! Lateral boundary conditions 514 CALL lbc_lnk( pgzui , 'U', -1. ) ; CALL lbc_lnk( pgzvi , 'V', -1. ) ! Lateral boundary conditions 515 CALL lbc_lnk( pge3rui , 'U', -1. ) ; CALL lbc_lnk( pge3rvi , 'V', -1. ) ! Lateral boundary conditions 516 ! 517 END IF 518 ! 519 IF( nn_timing == 1 ) CALL timing_stop( 'zps_hde_isf') 520 ! 521 END SUBROUTINE zps_hde_isf 218 522 !!====================================================================== 219 523 END MODULE zpshde
Note: See TracChangeset
for help on using the changeset viewer.
