Changeset 12377 for NEMO/trunk/src/NST/agrif_oce_update.F90
- Timestamp:
- 2020-02-12T15:39:06+01:00 (4 years ago)
- Location:
- NEMO/trunk
- Files:
-
- 2 edited
-
. (modified) (1 prop)
-
src/NST/agrif_oce_update.F90 (modified) (71 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/trunk
- Property svn:externals
-
old new 3 3 ^/utils/build/mk@HEAD mk 4 4 ^/utils/tools@HEAD tools 5 ^/vendors/AGRIF/dev @HEAD ext/AGRIF5 ^/vendors/AGRIF/dev_r11615_ENHANCE-04_namelists_as_internalfiles_agrif@HEAD ext/AGRIF 6 6 ^/vendors/FCM@HEAD ext/FCM 7 7 ^/vendors/IOIPSL@HEAD ext/IOIPSL
-
- Property svn:externals
-
NEMO/trunk/src/NST/agrif_oce_update.F90
r10068 r12377 1 # define TWO_WAY /* TWO WAY NESTING*/2 #undef DECAL_FEEDBACK /* SEPARATION of INTERFACES*/3 #undef VOL_REFLUX /* VOLUME REFLUXING*/1 #undef DECAL_FEEDBACK /* SEPARATION of INTERFACES */ 2 #undef DECAL_FEEDBACK_2D /* SEPARATION of INTERFACES (Barotropic mode) */ 3 #undef VOL_REFLUX /* VOLUME REFLUXING*/ 4 4 5 5 MODULE agrif_oce_update … … 25 25 USE lib_mpp ! MPP library 26 26 USE domvvl ! Need interpolation routines 27 USE vremap ! Vertical remapping 27 28 28 29 IMPLICIT NONE … … 46 47 IF (Agrif_Root()) RETURN 47 48 ! 48 #if defined TWO_WAY49 49 IF (lwp.AND.lk_agrif_debug) Write(*,*) 'Update tracers from grid Number',Agrif_Fixed() 50 50 51 #if defined key_vertical 52 ! Effect of this has to be carrefully checked 53 ! depending on what the nesting tools ensure for 54 ! volume conservation: 55 Agrif_UseSpecialValueInUpdate = .FALSE. 56 #else 51 57 Agrif_UseSpecialValueInUpdate = .TRUE. 58 #endif 52 59 Agrif_SpecialValueFineGrid = 0._wp 53 60 ! … … 64 71 Agrif_UseSpecialValueInUpdate = .FALSE. 65 72 ! 66 #endif67 73 ! 68 74 END SUBROUTINE Agrif_Update_Tra … … 75 81 IF (Agrif_Root()) RETURN 76 82 ! 77 #if defined TWO_WAY78 83 IF (lwp.AND.lk_agrif_debug) Write(*,*) 'Update momentum from grid Number',Agrif_Fixed() 79 84 … … 95 100 # endif 96 101 97 # if ! defined DECAL_FEEDBACK 102 # if ! defined DECAL_FEEDBACK_2D 98 103 CALL Agrif_Update_Variable(e1u_id,procname = updateU2d) 99 104 CALL Agrif_Update_Variable(e2v_id,procname = updateV2d) … … 103 108 # endif 104 109 ! 105 # if ! defined DECAL_FEEDBACK 110 # if ! defined DECAL_FEEDBACK_2D 106 111 ! Account for updated thicknesses at boundary edges 107 112 IF (.NOT.ln_linssh) THEN … … 113 118 IF ( ln_dynspg_ts .AND. ln_bt_fw ) THEN 114 119 ! Update time integrated transports 115 # if ! defined DECAL_FEEDBACK 120 # if ! defined DECAL_FEEDBACK_2D 116 121 CALL Agrif_Update_Variable(ub2b_update_id,procname = updateub2b) 117 122 CALL Agrif_Update_Variable(vb2b_update_id,procname = updatevb2b) … … 121 126 # endif 122 127 END IF 123 #endif124 128 ! 125 129 END SUBROUTINE Agrif_Update_Dyn … … 131 135 ! 132 136 IF (Agrif_Root()) RETURN 133 !134 #if defined TWO_WAY135 137 ! 136 138 Agrif_UseSpecialValueInUpdate = .TRUE. 137 139 Agrif_SpecialValueFineGrid = 0. 138 # if ! defined DECAL_FEEDBACK 140 # if ! defined DECAL_FEEDBACK_2D 139 141 CALL Agrif_Update_Variable(sshn_id,procname = updateSSH) 140 142 # else … … 147 149 IF ( ln_dynspg_ts.AND.ln_bt_fw ) THEN 148 150 ! Refluxing on ssh: 149 # if defined DECAL_FEEDBACK 151 # if defined DECAL_FEEDBACK_2D 150 152 CALL Agrif_Update_Variable(ub2b_update_id,locupdate1=(/0, 0/),locupdate2=(/1, 1/),procname = reflux_sshu) 151 153 CALL Agrif_Update_Variable(vb2b_update_id,locupdate1=(/1, 1/),locupdate2=(/0, 0/),procname = reflux_sshv) … … 157 159 # endif 158 160 ! 159 #endif160 !161 161 END SUBROUTINE Agrif_Update_ssh 162 162 … … 170 170 IF (Agrif_Root()) RETURN 171 171 ! 172 # if defined TWO_WAY173 174 172 Agrif_UseSpecialValueInUpdate = .TRUE. 175 173 Agrif_SpecialValueFineGrid = 0. … … 180 178 181 179 Agrif_UseSpecialValueInUpdate = .FALSE. 182 183 # endif184 180 185 181 END SUBROUTINE Agrif_Update_Tke … … 192 188 ! 193 189 IF (Agrif_Root()) RETURN 194 !195 #if defined TWO_WAY196 190 ! 197 191 IF (lwp.AND.lk_agrif_debug) Write(*,*) 'Update e3 from grid Number',Agrif_Fixed(), 'Step', Agrif_Nb_Step() … … 210 204 CALL Agrif_ParentGrid_To_ChildGrid() 211 205 ! 212 #endif213 !214 206 END SUBROUTINE Agrif_Update_vvl 215 207 … … 230 222 ! ----------------------- 231 223 ! 232 e3u _a(:,:,:) = e3u_n(:,:,:)233 e3v _a(:,:,:) = e3v_n(:,:,:)234 ! u a(:,:,:) = e3u_b(:,:,:)235 ! v a(:,:,:) = e3v_b(:,:,:)236 hu _a(:,:) = hu_n(:,:)237 hv _a(:,:) = hv_n(:,:)224 e3u(:,:,:,Krhs_a) = e3u(:,:,:,Kmm_a) 225 e3v(:,:,:,Krhs_a) = e3v(:,:,:,Kmm_a) 226 ! uu(:,:,:,Krhs_a) = e3u(:,:,:,Kbb_a) 227 ! vv(:,:,:,Krhs_a) = e3v(:,:,:,Kbb_a) 228 hu(:,:,Krhs_a) = hu(:,:,Kmm_a) 229 hv(:,:,Krhs_a) = hv(:,:,Kmm_a) 238 230 239 231 ! 1) NOW fields … … 242 234 ! Vertical scale factor interpolations 243 235 ! ------------------------------------ 244 CALL dom_vvl_interpol( e3t _n(:,:,:), e3u_n(:,:,:) , 'U' )245 CALL dom_vvl_interpol( e3t _n(:,:,:), e3v_n(:,:,:) , 'V' )246 CALL dom_vvl_interpol( e3u _n(:,:,:), e3f_n(:,:,:) , 'F' )247 248 CALL dom_vvl_interpol( e3u _n(:,:,:), e3uw_n(:,:,:), 'UW' )249 CALL dom_vvl_interpol( e3v _n(:,:,:), e3vw_n(:,:,:), 'VW' )236 CALL dom_vvl_interpol( e3t(:,:,:,Kmm_a), e3u(:,:,:,Kmm_a) , 'U' ) 237 CALL dom_vvl_interpol( e3t(:,:,:,Kmm_a), e3v(:,:,:,Kmm_a) , 'V' ) 238 CALL dom_vvl_interpol( e3u(:,:,:,Kmm_a), e3f(:,:,:) , 'F' ) 239 240 CALL dom_vvl_interpol( e3u(:,:,:,Kmm_a), e3uw(:,:,:,Kmm_a), 'UW' ) 241 CALL dom_vvl_interpol( e3v(:,:,:,Kmm_a), e3vw(:,:,:,Kmm_a), 'VW' ) 250 242 251 243 ! Update total depths: 252 244 ! -------------------- 253 hu _n(:,:) = 0._wp ! Ocean depth at U-points254 hv _n(:,:) = 0._wp ! Ocean depth at V-points245 hu(:,:,Kmm_a) = 0._wp ! Ocean depth at U-points 246 hv(:,:,Kmm_a) = 0._wp ! Ocean depth at V-points 255 247 DO jk = 1, jpkm1 256 hu _n(:,:) = hu_n(:,:) + e3u_n(:,:,jk) * umask(:,:,jk)257 hv _n(:,:) = hv_n(:,:) + e3v_n(:,:,jk) * vmask(:,:,jk)248 hu(:,:,Kmm_a) = hu(:,:,Kmm_a) + e3u(:,:,jk,Kmm_a) * umask(:,:,jk) 249 hv(:,:,Kmm_a) = hv(:,:,Kmm_a) + e3v(:,:,jk,Kmm_a) * vmask(:,:,jk) 258 250 END DO 259 251 ! ! Inverse of the local depth 260 r1_hu _n(:,:) = ssumask(:,:) / ( hu_n(:,:) + 1._wp - ssumask(:,:) )261 r1_hv _n(:,:) = ssvmask(:,:) / ( hv_n(:,:) + 1._wp - ssvmask(:,:) )252 r1_hu(:,:,Kmm_a) = ssumask(:,:) / ( hu(:,:,Kmm_a) + 1._wp - ssumask(:,:) ) 253 r1_hv(:,:,Kmm_a) = ssvmask(:,:) / ( hv(:,:,Kmm_a) + 1._wp - ssvmask(:,:) ) 262 254 263 255 … … 268 260 ! Vertical scale factor interpolations 269 261 ! ------------------------------------ 270 CALL dom_vvl_interpol( e3t _b(:,:,:), e3u_b(:,:,:), 'U' )271 CALL dom_vvl_interpol( e3t _b(:,:,:), e3v_b(:,:,:), 'V' )272 273 CALL dom_vvl_interpol( e3u _b(:,:,:), e3uw_b(:,:,:), 'UW' )274 CALL dom_vvl_interpol( e3v _b(:,:,:), e3vw_b(:,:,:), 'VW' )262 CALL dom_vvl_interpol( e3t(:,:,:,Kbb_a), e3u(:,:,:,Kbb_a), 'U' ) 263 CALL dom_vvl_interpol( e3t(:,:,:,Kbb_a), e3v(:,:,:,Kbb_a), 'V' ) 264 265 CALL dom_vvl_interpol( e3u(:,:,:,Kbb_a), e3uw(:,:,:,Kbb_a), 'UW' ) 266 CALL dom_vvl_interpol( e3v(:,:,:,Kbb_a), e3vw(:,:,:,Kbb_a), 'VW' ) 275 267 276 268 ! Update total depths: 277 269 ! -------------------- 278 hu _b(:,:) = 0._wp ! Ocean depth at U-points279 hv _b(:,:) = 0._wp ! Ocean depth at V-points270 hu(:,:,Kbb_a) = 0._wp ! Ocean depth at U-points 271 hv(:,:,Kbb_a) = 0._wp ! Ocean depth at V-points 280 272 DO jk = 1, jpkm1 281 hu _b(:,:) = hu_b(:,:) + e3u_b(:,:,jk) * umask(:,:,jk)282 hv _b(:,:) = hv_b(:,:) + e3v_b(:,:,jk) * vmask(:,:,jk)273 hu(:,:,Kbb_a) = hu(:,:,Kbb_a) + e3u(:,:,jk,Kbb_a) * umask(:,:,jk) 274 hv(:,:,Kbb_a) = hv(:,:,Kbb_a) + e3v(:,:,jk,Kbb_a) * vmask(:,:,jk) 283 275 END DO 284 276 ! ! Inverse of the local depth 285 r1_hu _b(:,:) = ssumask(:,:) / ( hu_b(:,:) + 1._wp - ssumask(:,:) )286 r1_hv _b(:,:) = ssvmask(:,:) / ( hv_b(:,:) + 1._wp - ssvmask(:,:) )277 r1_hu(:,:,Kbb_a) = ssumask(:,:) / ( hu(:,:,Kbb_a) + 1._wp - ssumask(:,:) ) 278 r1_hv(:,:,Kbb_a) = ssvmask(:,:) / ( hv(:,:,Kbb_a) + 1._wp - ssvmask(:,:) ) 287 279 ENDIF 288 280 ! … … 300 292 !! 301 293 INTEGER :: ji,jj,jk,jn 302 REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk,n1:n2) :: tabres_child 294 INTEGER :: N_in, N_out 295 REAL(wp) :: ztb, ztnu, ztno 303 296 REAL(wp) :: h_in(k1:k2) 304 297 REAL(wp) :: h_out(1:jpk) 305 INTEGER :: N_in, N_out 306 REAL(wp) :: zrho_xy, h_diff 307 REAL(wp) :: tabin(k1:k2,n1:n2) 298 REAL(wp) :: tabin(k1:k2,1:jpts) 299 REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk,1:jpts) :: tabres_child 308 300 !!--------------------------------------------- 309 301 ! 310 302 IF (before) THEN 311 AGRIF_SpecialValue = -999._wp 312 zrho_xy = Agrif_rhox() * Agrif_rhoy() 303 !jc_alt 304 ! AGRIF_SpecialValue = -999._wp 313 305 DO jn = n1,n2-1 314 306 DO jk=k1,k2 315 307 DO jj=j1,j2 316 308 DO ji=i1,i2 317 tabres(ji,jj,jk,jn) = (tsn(ji,jj,jk,jn) * e3t_n(ji,jj,jk) ) & 318 * tmask(ji,jj,jk) + (tmask(ji,jj,jk)-1)*999._wp 309 !jc_alt 310 ! tabres(ji,jj,jk,jn) = (ts(ji,jj,jk,jn,Kmm_a) * e3t(ji,jj,jk,Kmm_a) ) & 311 ! & * tmask(ji,jj,jk) + (tmask(ji,jj,jk)-1._wp) * 999._wp 312 tabres(ji,jj,jk,jn) = ts(ji,jj,jk,jn,Kmm_a) * e3t(ji,jj,jk,Kmm_a) 319 313 END DO 320 314 END DO … … 324 318 DO jj=j1,j2 325 319 DO ji=i1,i2 326 tabres(ji,jj,jk,n2) = tmask(ji,jj,jk) * e3t_n(ji,jj,jk) & 327 + (tmask(ji,jj,jk)-1)*999._wp 320 !jc_alt 321 ! tabres(ji,jj,jk,n2) = tmask(ji,jj,jk) * e3t(ji,jj,jk,Kmm_a) & 322 ! & + (tmask(ji,jj,jk) - 1._wp) * 999._wp 323 tabres(ji,jj,jk,n2) = tmask(ji,jj,jk) * e3t(ji,jj,jk,Kmm_a) 328 324 END DO 329 325 END DO 330 326 END DO 331 327 ELSE 332 tabres_child(:,:,:,:) = 0. 328 tabres_child(:,:,:,:) = 0._wp 333 329 AGRIF_SpecialValue = 0._wp 334 330 DO jj=j1,j2 … … 336 332 N_in = 0 337 333 DO jk=k1,k2 !k2 = jpk of child grid 338 IF (tabres(ji,jj,jk,n2) == 0 ) EXIT 334 ! jc_alt 335 ! IF (tabres(ji,jj,jk,n2) < -900._wp ) EXIT 336 IF (tabres(ji,jj,jk,n2) == 0._wp ) EXIT 339 337 N_in = N_in + 1 340 338 tabin(jk,:) = tabres(ji,jj,jk,n1:n2-1)/tabres(ji,jj,jk,n2) … … 343 341 N_out = 0 344 342 DO jk=1,jpk ! jpk of parent grid 345 IF (tmask(ji,jj,jk) < -900) EXIT ! TODO: Will not work with ISF343 IF (tmask(ji,jj,jk) == 0 ) EXIT ! TODO: Will not work with ISF 346 344 N_out = N_out + 1 347 h_out(N_out) = e3t _n(ji,jj,jk)345 h_out(N_out) = e3t(ji,jj,jk,Kmm_a) 348 346 ENDDO 349 IF (N_in > 0) THEN !Remove this? 350 h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in)) 351 IF (h_diff < -1.e-4) THEN 352 print *,'CHECK YOUR bathy T points ...',ji,jj,h_diff,sum(h_in(1:N_in)),sum(h_out(1:N_out)) 353 print *,h_in(1:N_in) 354 print *,h_out(1:N_out) 355 STOP 356 ENDIF 357 DO jn=n1,n2-1 358 CALL reconstructandremap(tabin(1:N_in,jn),h_in(1:N_in),tabres_child(ji,jj,1:N_out,jn),h_out(1:N_out),N_in,N_out) 359 ENDDO 347 IF (N_in*N_out > 0) THEN !Remove this? 348 CALL reconstructandremap(tabin(1:N_in,1:jpts),h_in(1:N_in),tabres_child(ji,jj,1:N_out,1:jpts),h_out(1:N_out),N_in,N_out,jpts) 360 349 ENDIF 361 350 ENDDO … … 364 353 IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN 365 354 ! Add asselin part 366 DO jn = n1,n2-1 367 DO jk=1,jpk 368 DO jj=j1,j2 369 DO ji=i1,i2 370 IF( tabres_child(ji,jj,jk,jn) .NE. 0. ) THEN 371 tsb(ji,jj,jk,jn) = tsb(ji,jj,jk,jn) & 372 & + atfp * ( tabres_child(ji,jj,jk,jn) & 373 & - tsn(ji,jj,jk,jn) ) * tmask(ji,jj,jk) 355 DO jn = 1,jpts 356 DO jk = 1, jpkm1 357 DO jj = j1, j2 358 DO ji = i1, i2 359 IF( tabres_child(ji,jj,jk,jn) /= 0._wp ) THEN 360 ztb = ts(ji,jj,jk,jn,Kbb_a) * e3t(ji,jj,jk,Kbb_a) ! fse3t_b prior update should be used 361 ztnu = tabres_child(ji,jj,jk,jn) * e3t(ji,jj,jk,Kmm_a) 362 ztno = ts(ji,jj,jk,jn,Kmm_a) * e3t(ji,jj,jk,Krhs_a) 363 ts(ji,jj,jk,jn,Kbb_a) = ( ztb + atfp * ( ztnu - ztno) ) & 364 & * tmask(ji,jj,jk) / e3t(ji,jj,jk,Kbb_a) 374 365 ENDIF 375 END DO376 END DO377 END DO378 END DO379 ENDIF 380 DO jn = n1,n2-1381 DO jk =1,jpk382 DO jj =j1,j2383 DO ji =i1,i2384 IF( tabres_child(ji,jj,jk,jn) .NE. 0.) THEN385 ts n(ji,jj,jk,jn) = tabres_child(ji,jj,jk,jn) * tmask(ji,jj,jk)366 END DO 367 END DO 368 END DO 369 END DO 370 ENDIF 371 DO jn = 1,jpts 372 DO jk = 1, jpkm1 373 DO jj = j1, j2 374 DO ji = i1, i2 375 IF( tabres_child(ji,jj,jk,jn) /= 0._wp ) THEN 376 ts(ji,jj,jk,jn,Kmm_a) = tabres_child(ji,jj,jk,jn) 386 377 END IF 387 378 END DO … … 389 380 END DO 390 381 END DO 382 ! 383 IF ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN 384 ts(i1:i2,j1:j2,1:jpkm1,1:jpts,Kbb_a) = ts(i1:i2,j1:j2,1:jpkm1,1:jpts,Kmm_a) 385 ENDIF 391 386 ENDIF 392 387 ! … … 413 408 DO ji=i1,i2 414 409 !> jc tmp 415 tabres(ji,jj,jk,jn) = ts n(ji,jj,jk,jn) * e3t_n(ji,jj,jk) / e3t_0(ji,jj,jk)416 ! tabres(ji,jj,jk,jn) = ts n(ji,jj,jk,jn) * e3t_n(ji,jj,jk)410 tabres(ji,jj,jk,jn) = ts(ji,jj,jk,jn,Kmm_a) * e3t(ji,jj,jk,Kmm_a) / e3t_0(ji,jj,jk) 411 ! tabres(ji,jj,jk,jn) = ts(ji,jj,jk,jn,Kmm_a) * e3t(ji,jj,jk,Kmm_a) 417 412 !< jc tmp 418 413 END DO … … 434 429 DO ji = i1, i2 435 430 IF( tabres(ji,jj,jk,jn) /= 0._wp ) THEN 436 ztb = ts b(ji,jj,jk,jn) * e3t_b(ji,jj,jk) ! fse3t_b prior update should be used431 ztb = ts(ji,jj,jk,jn,Kbb_a) * e3t(ji,jj,jk,Kbb_a) ! fse3t_b prior update should be used 437 432 ztnu = tabres(ji,jj,jk,jn) 438 ztno = ts n(ji,jj,jk,jn) * e3t_a(ji,jj,jk)439 ts b(ji,jj,jk,jn) = ( ztb + atfp * ( ztnu - ztno) ) &440 & * tmask(ji,jj,jk) / e3t _b(ji,jj,jk)433 ztno = ts(ji,jj,jk,jn,Kmm_a) * e3t(ji,jj,jk,Krhs_a) 434 ts(ji,jj,jk,jn,Kbb_a) = ( ztb + atfp * ( ztnu - ztno) ) & 435 & * tmask(ji,jj,jk) / e3t(ji,jj,jk,Kbb_a) 441 436 ENDIF 442 437 END DO … … 450 445 DO ji=i1,i2 451 446 IF( tabres(ji,jj,jk,jn) /= 0._wp ) THEN 452 ts n(ji,jj,jk,jn) = tabres(ji,jj,jk,jn) / e3t_n(ji,jj,jk)447 ts(ji,jj,jk,jn,Kmm_a) = tabres(ji,jj,jk,jn) / e3t(ji,jj,jk,Kmm_a) 453 448 END IF 454 449 END DO … … 458 453 ! 459 454 IF ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN 460 ts b(i1:i2,j1:j2,k1:k2,1:jpts) = tsn(i1:i2,j1:j2,k1:k2,1:jpts)455 ts(i1:i2,j1:j2,k1:k2,1:jpts,Kbb_a) = ts(i1:i2,j1:j2,k1:k2,1:jpts,Kmm_a) 461 456 ENDIF 462 457 ! … … 478 473 ! 479 474 INTEGER :: ji, jj, jk 480 REAL(wp):: zrhoy 475 REAL(wp):: zrhoy, zub, zunu, zuno 481 476 ! VERTICAL REFINEMENT BEGIN 482 477 REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: tabres_child … … 491 486 IF( before ) THEN 492 487 zrhoy = Agrif_Rhoy() 493 AGRIF_SpecialValue = -999._wp 488 !jc_alt 489 ! AGRIF_SpecialValue = -999._wp 494 490 DO jk=k1,k2 495 491 DO jj=j1,j2 496 492 DO ji=i1,i2 497 tabres(ji,jj,jk,1) = zrhoy * e2u(ji,jj) * e3u_n(ji,jj,jk) * umask(ji,jj,jk) * un(ji,jj,jk) & 498 + (umask(ji,jj,jk)-1)*999._wp 499 tabres(ji,jj,jk,2) = zrhoy * umask(ji,jj,jk) * e2u(ji,jj) * e3u_n(ji,jj,jk) & 500 + (umask(ji,jj,jk)-1)*999._wp 493 !jc_alt 494 ! tabres(ji,jj,jk,1) = zrhoy * e2u(ji,jj) * e3u(ji,jj,jk,Kmm_a) * umask(ji,jj,jk) * uu(ji,jj,jk,Kmm_a) & 495 ! & + (umask(ji,jj,jk)-1._wp)*999._wp 496 tabres(ji,jj,jk,1) = zrhoy * e2u(ji,jj) * e3u(ji,jj,jk,Kmm_a) * umask(ji,jj,jk) * uu(ji,jj,jk,Kmm_a) 497 !jc_alt 498 ! tabres(ji,jj,jk,2) = zrhoy * umask(ji,jj,jk) * e2u(ji,jj) * e3u(ji,jj,jk,Kmm_a) & 499 ! & + (umask(ji,jj,jk)-1._wp)*999._wp 500 tabres(ji,jj,jk,2) = zrhoy * umask(ji,jj,jk) * e2u(ji,jj) * e3u(ji,jj,jk,Kmm_a) 501 501 END DO 502 502 END DO … … 511 511 tabin(:) = 0._wp 512 512 DO jk=k1,k2 !k2=jpk of child grid 513 IF( tabres(ji,jj,jk,2) < -900) EXIT 513 !jc_alt 514 ! IF( tabres(ji,jj,jk,2) < -900._wp) EXIT 515 IF( tabres(ji,jj,jk,2) == 0.) EXIT 514 516 N_in = N_in + 1 515 517 tabin(jk) = tabres(ji,jj,jk,1)/tabres(ji,jj,jk,2) … … 520 522 IF (umask(ji,jj,jk) == 0) EXIT 521 523 N_out = N_out + 1 522 h_out(N_out) = e3u _n(ji,jj,jk)524 h_out(N_out) = e3u(ji,jj,jk,Kmm_a) 523 525 ENDDO 524 526 IF (N_in * N_out > 0) THEN 525 527 h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in)) 528 excess = 0._wp 526 529 IF (h_diff < -1.e-4) THEN 527 530 !Even if bathy at T points match it's possible for the U points to be deeper in the child grid. 528 531 !In this case we need to move transport from the child grid cells below bed of parent grid into the bottom cell. 529 excess = 0._wp530 532 DO jk=N_in,1,-1 531 533 thick = MIN(-1*h_diff, h_in(jk)) … … 540 542 ENDDO 541 543 ENDIF 542 CALL reconstructandremap(tabin(1:N_in),h_in(1:N_in),tabres_child(ji,jj,1:N_out),h_out(1:N_out),N_in,N_out )544 CALL reconstructandremap(tabin(1:N_in),h_in(1:N_in),tabres_child(ji,jj,1:N_out),h_out(1:N_out),N_in,N_out,1) 543 545 tabres_child(ji,jj,N_out) = tabres_child(ji,jj,N_out) + excess/(e2u(ji,jj)*h_out(N_out)) 544 546 ENDIF 545 547 ENDDO 546 548 ENDDO 547 549 ! 548 550 DO jk=1,jpk 549 551 DO jj=j1,j2 550 552 DO ji=i1,i2 551 553 IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part 552 ub(ji,jj,jk) = ub(ji,jj,jk) & 553 & + atfp * ( tabres_child(ji,jj,jk) - un(ji,jj,jk) ) * umask(ji,jj,jk) 554 zub = uu(ji,jj,jk,Kbb_a) * e3u(ji,jj,jk,Kbb_a) ! fse3t_b prior update should be used 555 zuno = uu(ji,jj,jk,Kmm_a) * e3u(ji,jj,jk,Krhs_a) 556 zunu = tabres_child(ji,jj,jk) * e3u(ji,jj,jk,Kmm_a) 557 uu(ji,jj,jk,Kbb_a) = ( zub + atfp * ( zunu - zuno) ) & 558 & * umask(ji,jj,jk) / e3u(ji,jj,jk,Kbb_a) 554 559 ENDIF 555 560 ! 556 un(ji,jj,jk) = tabres_child(ji,jj,jk) * umask(ji,jj,jk) 557 END DO 558 END DO 559 END DO 561 uu(ji,jj,jk,Kmm_a) = tabres_child(ji,jj,jk) * umask(ji,jj,jk) 562 END DO 563 END DO 564 END DO 565 ! 566 IF ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN 567 uu(i1:i2,j1:j2,1:jpkm1,Kbb_a) = uu(i1:i2,j1:j2,1:jpkm1,Kmm_a) 568 ENDIF 569 ! 560 570 ENDIF 561 571 ! … … 579 589 zrhoy = Agrif_Rhoy() 580 590 DO jk = k1, k2 581 tabres(i1:i2,j1:j2,jk,1) = zrhoy * e2u(i1:i2,j1:j2) * e3u _n(i1:i2,j1:j2,jk) * un(i1:i2,j1:j2,jk)591 tabres(i1:i2,j1:j2,jk,1) = zrhoy * e2u(i1:i2,j1:j2) * e3u(i1:i2,j1:j2,jk,Kmm_a) * uu(i1:i2,j1:j2,jk,Kmm_a) 582 592 END DO 583 593 ELSE … … 588 598 ! 589 599 IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part 590 zub = u b(ji,jj,jk) * e3u_b(ji,jj,jk) ! fse3t_b prior update should be used591 zuno = u n(ji,jj,jk) * e3u_a(ji,jj,jk)600 zub = uu(ji,jj,jk,Kbb_a) * e3u(ji,jj,jk,Kbb_a) ! fse3t_b prior update should be used 601 zuno = uu(ji,jj,jk,Kmm_a) * e3u(ji,jj,jk,Krhs_a) 592 602 zunu = tabres(ji,jj,jk,1) 593 u b(ji,jj,jk) = ( zub + atfp * ( zunu - zuno) ) &594 & * umask(ji,jj,jk) / e3u _b(ji,jj,jk)603 uu(ji,jj,jk,Kbb_a) = ( zub + atfp * ( zunu - zuno) ) & 604 & * umask(ji,jj,jk) / e3u(ji,jj,jk,Kbb_a) 595 605 ENDIF 596 606 ! 597 u n(ji,jj,jk) = tabres(ji,jj,jk,1) * umask(ji,jj,jk) / e3u_n(ji,jj,jk)607 uu(ji,jj,jk,Kmm_a) = tabres(ji,jj,jk,1) * umask(ji,jj,jk) / e3u(ji,jj,jk,Kmm_a) 598 608 END DO 599 609 END DO … … 601 611 ! 602 612 IF ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN 603 u b(i1:i2,j1:j2,k1:k2) = un(i1:i2,j1:j2,k1:k2)613 uu(i1:i2,j1:j2,k1:k2,Kbb_a) = uu(i1:i2,j1:j2,k1:k2,Kmm_a) 604 614 ENDIF 605 615 ! … … 632 642 IF (western_side) THEN 633 643 DO jj=j1,j2 634 zcor = u n_b(i1-1,jj) * hu_a(i1-1,jj) * r1_hu_n(i1-1,jj) - un_b(i1-1,jj)635 u n_b(i1-1,jj) = un_b(i1-1,jj) + zcor644 zcor = uu_b(i1-1,jj,Kmm_a) * hu(i1-1,jj,Krhs_a) * r1_hu(i1-1,jj,Kmm_a) - uu_b(i1-1,jj,Kmm_a) 645 uu_b(i1-1,jj,Kmm_a) = uu_b(i1-1,jj,Kmm_a) + zcor 636 646 DO jk=1,jpkm1 637 u n(i1-1,jj,jk) = un(i1-1,jj,jk) + zcor * umask(i1-1,jj,jk)647 uu(i1-1,jj,jk,Kmm_a) = uu(i1-1,jj,jk,Kmm_a) + zcor * umask(i1-1,jj,jk) 638 648 END DO 639 649 END DO … … 642 652 IF (eastern_side) THEN 643 653 DO jj=j1,j2 644 zcor = u n_b(i2+1,jj) * hu_a(i2+1,jj) * r1_hu_n(i2+1,jj) - un_b(i2+1,jj)645 u n_b(i2+1,jj) = un_b(i2+1,jj) + zcor654 zcor = uu_b(i2+1,jj,Kmm_a) * hu(i2+1,jj,Krhs_a) * r1_hu(i2+1,jj,Kmm_a) - uu_b(i2+1,jj,Kmm_a) 655 uu_b(i2+1,jj,Kmm_a) = uu_b(i2+1,jj,Kmm_a) + zcor 646 656 DO jk=1,jpkm1 647 u n(i2+1,jj,jk) = un(i2+1,jj,jk) + zcor * umask(i2+1,jj,jk)657 uu(i2+1,jj,jk,Kmm_a) = uu(i2+1,jj,jk,Kmm_a) + zcor * umask(i2+1,jj,jk) 648 658 END DO 649 659 END DO … … 665 675 ! 666 676 INTEGER :: ji, jj, jk 667 REAL(wp) :: zrhox 677 REAL(wp) :: zrhox, zvb, zvnu, zvno 668 678 ! VERTICAL REFINEMENT BEGIN 669 679 REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: tabres_child … … 678 688 IF( before ) THEN 679 689 zrhox = Agrif_Rhox() 680 AGRIF_SpecialValue = -999._wp 690 !jc_alt 691 ! AGRIF_SpecialValue = -999._wp 681 692 DO jk=k1,k2 682 693 DO jj=j1,j2 683 694 DO ji=i1,i2 684 tabres(ji,jj,jk,1) = zrhox * e1v(ji,jj) * e3v_n(ji,jj,jk) * vmask(ji,jj,jk) * vn(ji,jj,jk) & 685 + (vmask(ji,jj,jk)-1)*999._wp 686 tabres(ji,jj,jk,2) = vmask(ji,jj,jk) * zrhox * e1v(ji,jj) * e3v_n(ji,jj,jk) & 687 + (vmask(ji,jj,jk)-1)*999._wp 695 !jc_alt 696 ! tabres(ji,jj,jk,1) = zrhox * e1v(ji,jj) * e3v(ji,jj,jk,Kmm_a) * vmask(ji,jj,jk) * vv(ji,jj,jk,Kmm_a) & 697 ! & + (vmask(ji,jj,jk)-1._wp) * 999._wp 698 tabres(ji,jj,jk,1) = zrhox * e1v(ji,jj) * e3v(ji,jj,jk,Kmm_a) * vmask(ji,jj,jk) * vv(ji,jj,jk,Kmm_a) 699 !jc_alt 700 ! tabres(ji,jj,jk,2) = zrhox * e1v(ji,jj) * e3v(ji,jj,jk,Kmm_a) * vmask(ji,jj,jk) & 701 ! & + (vmask(ji,jj,jk)-1._wp) * 999._wp 702 tabres(ji,jj,jk,2) = zrhox * e1v(ji,jj) * e3v(ji,jj,jk,Kmm_a) * vmask(ji,jj,jk) 688 703 END DO 689 704 END DO … … 696 711 N_in = 0 697 712 DO jk=k1,k2 698 IF (tabres(ji,jj,jk,2) < -900) EXIT 713 !jc_alt 714 ! IF (tabres(ji,jj,jk,2) < -900._wp) EXIT 715 IF (tabres(ji,jj,jk,2) == 0) EXIT 699 716 N_in = N_in + 1 700 717 tabin(jk) = tabres(ji,jj,jk,1)/tabres(ji,jj,jk,2) … … 705 722 IF (vmask(ji,jj,jk) == 0) EXIT 706 723 N_out = N_out + 1 707 h_out(N_out) = e3v _n(ji,jj,jk)724 h_out(N_out) = e3v(ji,jj,jk,Kmm_a) 708 725 ENDDO 709 726 IF (N_in * N_out > 0) THEN 710 727 h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in)) 728 excess = 0._wp 711 729 IF (h_diff < -1.e-4) then 712 !Even if bathy at T points match it's possible for the Upoints to be deeper in the child grid.730 !Even if bathy at T points match it's possible for the V points to be deeper in the child grid. 713 731 !In this case we need to move transport from the child grid cells below bed of parent grid into the bottom cell. 714 excess = 0._wp715 732 DO jk=N_in,1,-1 716 733 thick = MIN(-1*h_diff, h_in(jk)) … … 725 742 ENDDO 726 743 ENDIF 727 CALL reconstructandremap(tabin(1:N_in),h_in(1:N_in),tabres_child(ji,jj,1:N_out),h_out(1:N_out),N_in,N_out )744 CALL reconstructandremap(tabin(1:N_in),h_in(1:N_in),tabres_child(ji,jj,1:N_out),h_out(1:N_out),N_in,N_out,1) 728 745 tabres_child(ji,jj,N_out) = tabres_child(ji,jj,N_out) + excess/(e1v(ji,jj)*h_out(N_out)) 729 746 ENDIF 730 747 ENDDO 731 748 ENDDO 732 733 DO jk=1,jpk 749 ! 750 DO jk=1,jpkm1 734 751 DO jj=j1,j2 735 752 DO ji=i1,i2 736 ! 737 IF( .NOT.(lk_agrif_fstep.AND.(neuler==0)) ) THEN ! Add asselin part 738 vb(ji,jj,jk) = vb(ji,jj,jk) & 739 & + atfp * ( tabres_child(ji,jj,jk) - vn(ji,jj,jk) ) * vmask(ji,jj,jk) 753 IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part 754 zvb = vv(ji,jj,jk,Kbb_a) * e3v(ji,jj,jk,Kbb_a) ! fse3t_b prior update should be used 755 zvno = vv(ji,jj,jk,Kmm_a) * e3v(ji,jj,jk,Krhs_a) 756 zvnu = tabres_child(ji,jj,jk) * e3v(ji,jj,jk,Kmm_a) 757 vv(ji,jj,jk,Kbb_a) = ( zvb + atfp * ( zvnu - zvno) ) & 758 & * vmask(ji,jj,jk) / e3v(ji,jj,jk,Kbb_a) 740 759 ENDIF 741 760 ! 742 vn(ji,jj,jk) = tabres_child(ji,jj,jk) * vmask(ji,jj,jk) 743 END DO 744 END DO 745 END DO 761 vv(ji,jj,jk,Kmm_a) = tabres_child(ji,jj,jk) * vmask(ji,jj,jk) 762 END DO 763 END DO 764 END DO 765 ! 766 IF ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN 767 vv(i1:i2,j1:j2,1:jpkm1,Kbb_a) = vv(i1:i2,j1:j2,1:jpkm1,Kmm_a) 768 ENDIF 769 ! 746 770 ENDIF 747 771 ! … … 767 791 DO jj=j1,j2 768 792 DO ji=i1,i2 769 tabres(ji,jj,jk,1) = zrhox * e1v(ji,jj) * e3v _n(ji,jj,jk) * vn(ji,jj,jk)793 tabres(ji,jj,jk,1) = zrhox * e1v(ji,jj) * e3v(ji,jj,jk,Kmm_a) * vv(ji,jj,jk,Kmm_a) 770 794 END DO 771 795 END DO … … 778 802 ! 779 803 IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part 780 zvb = v b(ji,jj,jk) * e3v_b(ji,jj,jk) ! fse3t_b prior update should be used781 zvno = v n(ji,jj,jk) * e3v_a(ji,jj,jk)804 zvb = vv(ji,jj,jk,Kbb_a) * e3v(ji,jj,jk,Kbb_a) ! fse3t_b prior update should be used 805 zvno = vv(ji,jj,jk,Kmm_a) * e3v(ji,jj,jk,Krhs_a) 782 806 zvnu = tabres(ji,jj,jk,1) 783 v b(ji,jj,jk) = ( zvb + atfp * ( zvnu - zvno) ) &784 & * vmask(ji,jj,jk) / e3v _b(ji,jj,jk)807 vv(ji,jj,jk,Kbb_a) = ( zvb + atfp * ( zvnu - zvno) ) & 808 & * vmask(ji,jj,jk) / e3v(ji,jj,jk,Kbb_a) 785 809 ENDIF 786 810 ! 787 v n(ji,jj,jk) = tabres(ji,jj,jk,1) * vmask(ji,jj,jk) / e3v_n(ji,jj,jk)811 vv(ji,jj,jk,Kmm_a) = tabres(ji,jj,jk,1) * vmask(ji,jj,jk) / e3v(ji,jj,jk,Kmm_a) 788 812 END DO 789 813 END DO … … 791 815 ! 792 816 IF ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN 793 v b(i1:i2,j1:j2,k1:k2) = vn(i1:i2,j1:j2,k1:k2)817 vv(i1:i2,j1:j2,k1:k2,Kbb_a) = vv(i1:i2,j1:j2,k1:k2,Kmm_a) 794 818 ENDIF 795 819 ! … … 822 846 IF (southern_side) THEN 823 847 DO ji=i1,i2 824 zcor = v n_b(ji,j1-1) * hv_a(ji,j1-1) * r1_hv_n(ji,j1-1) - vn_b(ji,j1-1)825 v n_b(ji,j1-1) = vn_b(ji,j1-1) + zcor848 zcor = vv_b(ji,j1-1,Kmm_a) * hv(ji,j1-1,Krhs_a) * r1_hv(ji,j1-1,Kmm_a) - vv_b(ji,j1-1,Kmm_a) 849 vv_b(ji,j1-1,Kmm_a) = vv_b(ji,j1-1,Kmm_a) + zcor 826 850 DO jk=1,jpkm1 827 v n(ji,j1-1,jk) = vn(ji,j1-1,jk) + zcor * vmask(ji,j1-1,jk)851 vv(ji,j1-1,jk,Kmm_a) = vv(ji,j1-1,jk,Kmm_a) + zcor * vmask(ji,j1-1,jk) 828 852 END DO 829 853 END DO … … 832 856 IF (northern_side) THEN 833 857 DO ji=i1,i2 834 zcor = v n_b(ji,j2+1) * hv_a(ji,j2+1) * r1_hv_n(ji,j2+1) - vn_b(ji,j2+1)835 v n_b(ji,j2+1) = vn_b(ji,j2+1) + zcor858 zcor = vv_b(ji,j2+1,Kmm_a) * hv(ji,j2+1,Krhs_a) * r1_hv(ji,j2+1,Kmm_a) - vv_b(ji,j2+1,Kmm_a) 859 vv_b(ji,j2+1,Kmm_a) = vv_b(ji,j2+1,Kmm_a) + zcor 836 860 DO jk=1,jpkm1 837 v n(ji,j2+1,jk) = vn(ji,j2+1,jk) + zcor * vmask(ji,j2+1,jk)861 vv(ji,j2+1,jk,Kmm_a) = vv(ji,j2+1,jk,Kmm_a) + zcor * vmask(ji,j2+1,jk) 838 862 END DO 839 863 END DO … … 862 886 DO jj=j1,j2 863 887 DO ji=i1,i2 864 tabres(ji,jj) = zrhoy * u n_b(ji,jj) * hu_n(ji,jj) * e2u(ji,jj)888 tabres(ji,jj) = zrhoy * uu_b(ji,jj,Kmm_a) * hu(ji,jj,Kmm_a) * e2u(ji,jj) 865 889 END DO 866 890 END DO … … 873 897 spgu(ji,jj) = 0._wp 874 898 DO jk=1,jpkm1 875 spgu(ji,jj) = spgu(ji,jj) + e3u _n(ji,jj,jk) * un(ji,jj,jk)899 spgu(ji,jj) = spgu(ji,jj) + e3u(ji,jj,jk,Kmm_a) * uu(ji,jj,jk,Kmm_a) 876 900 END DO 877 901 ! 878 zcorr = (tabres(ji,jj) - spgu(ji,jj)) * r1_hu _n(ji,jj)902 zcorr = (tabres(ji,jj) - spgu(ji,jj)) * r1_hu(ji,jj,Kmm_a) 879 903 DO jk=1,jpkm1 880 u n(ji,jj,jk) = un(ji,jj,jk) + zcorr * umask(ji,jj,jk)904 uu(ji,jj,jk,Kmm_a) = uu(ji,jj,jk,Kmm_a) + zcorr * umask(ji,jj,jk) 881 905 END DO 882 906 ! … … 884 908 IF ( .NOT.ln_dynspg_ts .OR. (ln_dynspg_ts.AND.(.NOT.ln_bt_fw)) ) THEN 885 909 IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part 886 zcorr = (tabres(ji,jj) - u n_b(ji,jj) * hu_a(ji,jj)) * r1_hu_b(ji,jj)887 u b_b(ji,jj) = ub_b(ji,jj) + atfp * zcorr * umask(ji,jj,1)910 zcorr = (tabres(ji,jj) - uu_b(ji,jj,Kmm_a) * hu(ji,jj,Krhs_a)) * r1_hu(ji,jj,Kbb_a) 911 uu_b(ji,jj,Kbb_a) = uu_b(ji,jj,Kbb_a) + atfp * zcorr * umask(ji,jj,1) 888 912 END IF 889 913 ENDIF 890 u n_b(ji,jj) = tabres(ji,jj) * r1_hu_n(ji,jj) * umask(ji,jj,1)914 uu_b(ji,jj,Kmm_a) = tabres(ji,jj) * r1_hu(ji,jj,Kmm_a) * umask(ji,jj,1) 891 915 ! 892 916 ! Correct "before" velocities to hold correct bt component: 893 917 spgu(ji,jj) = 0.e0 894 918 DO jk=1,jpkm1 895 spgu(ji,jj) = spgu(ji,jj) + e3u _b(ji,jj,jk) * ub(ji,jj,jk)919 spgu(ji,jj) = spgu(ji,jj) + e3u(ji,jj,jk,Kbb_a) * uu(ji,jj,jk,Kbb_a) 896 920 END DO 897 921 ! 898 zcorr = u b_b(ji,jj) - spgu(ji,jj) * r1_hu_b(ji,jj)922 zcorr = uu_b(ji,jj,Kbb_a) - spgu(ji,jj) * r1_hu(ji,jj,Kbb_a) 899 923 DO jk=1,jpkm1 900 u b(ji,jj,jk) = ub(ji,jj,jk) + zcorr * umask(ji,jj,jk)924 uu(ji,jj,jk,Kbb_a) = uu(ji,jj,jk,Kbb_a) + zcorr * umask(ji,jj,jk) 901 925 END DO 902 926 ! … … 905 929 ! 906 930 IF ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN 907 u b_b(i1:i2,j1:j2) = un_b(i1:i2,j1:j2)931 uu_b(i1:i2,j1:j2,Kbb_a) = uu_b(i1:i2,j1:j2,Kmm_a) 908 932 ENDIF 909 933 ENDIF … … 928 952 DO jj=j1,j2 929 953 DO ji=i1,i2 930 tabres(ji,jj) = zrhox * v n_b(ji,jj) * hv_n(ji,jj) * e1v(ji,jj)954 tabres(ji,jj) = zrhox * vv_b(ji,jj,Kmm_a) * hv(ji,jj,Kmm_a) * e1v(ji,jj) 931 955 END DO 932 956 END DO … … 939 963 spgv(ji,jj) = 0.e0 940 964 DO jk=1,jpkm1 941 spgv(ji,jj) = spgv(ji,jj) + e3v _n(ji,jj,jk) * vn(ji,jj,jk)965 spgv(ji,jj) = spgv(ji,jj) + e3v(ji,jj,jk,Kmm_a) * vv(ji,jj,jk,Kmm_a) 942 966 END DO 943 967 ! 944 zcorr = (tabres(ji,jj) - spgv(ji,jj)) * r1_hv _n(ji,jj)968 zcorr = (tabres(ji,jj) - spgv(ji,jj)) * r1_hv(ji,jj,Kmm_a) 945 969 DO jk=1,jpkm1 946 v n(ji,jj,jk) = vn(ji,jj,jk) + zcorr * vmask(ji,jj,jk)970 vv(ji,jj,jk,Kmm_a) = vv(ji,jj,jk,Kmm_a) + zcorr * vmask(ji,jj,jk) 947 971 END DO 948 972 ! … … 950 974 IF ( .NOT.ln_dynspg_ts .OR. (ln_dynspg_ts.AND.(.NOT.ln_bt_fw)) ) THEN 951 975 IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part 952 zcorr = (tabres(ji,jj) - v n_b(ji,jj) * hv_a(ji,jj)) * r1_hv_b(ji,jj)953 v b_b(ji,jj) = vb_b(ji,jj) + atfp * zcorr * vmask(ji,jj,1)976 zcorr = (tabres(ji,jj) - vv_b(ji,jj,Kmm_a) * hv(ji,jj,Krhs_a)) * r1_hv(ji,jj,Kbb_a) 977 vv_b(ji,jj,Kbb_a) = vv_b(ji,jj,Kbb_a) + atfp * zcorr * vmask(ji,jj,1) 954 978 END IF 955 979 ENDIF 956 v n_b(ji,jj) = tabres(ji,jj) * r1_hv_n(ji,jj) * vmask(ji,jj,1)980 vv_b(ji,jj,Kmm_a) = tabres(ji,jj) * r1_hv(ji,jj,Kmm_a) * vmask(ji,jj,1) 957 981 ! 958 982 ! Correct "before" velocities to hold correct bt component: 959 983 spgv(ji,jj) = 0.e0 960 984 DO jk=1,jpkm1 961 spgv(ji,jj) = spgv(ji,jj) + e3v _b(ji,jj,jk) * vb(ji,jj,jk)985 spgv(ji,jj) = spgv(ji,jj) + e3v(ji,jj,jk,Kbb_a) * vv(ji,jj,jk,Kbb_a) 962 986 END DO 963 987 ! 964 zcorr = v b_b(ji,jj) - spgv(ji,jj) * r1_hv_b(ji,jj)988 zcorr = vv_b(ji,jj,Kbb_a) - spgv(ji,jj) * r1_hv(ji,jj,Kbb_a) 965 989 DO jk=1,jpkm1 966 v b(ji,jj,jk) = vb(ji,jj,jk) + zcorr * vmask(ji,jj,jk)990 vv(ji,jj,jk,Kbb_a) = vv(ji,jj,jk,Kbb_a) + zcorr * vmask(ji,jj,jk) 967 991 END DO 968 992 ! … … 971 995 ! 972 996 IF ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN 973 v b_b(i1:i2,j1:j2) = vn_b(i1:i2,j1:j2)997 vv_b(i1:i2,j1:j2,Kbb_a) = vv_b(i1:i2,j1:j2,Kmm_a) 974 998 ENDIF 975 999 ! … … 993 1017 DO jj=j1,j2 994 1018 DO ji=i1,i2 995 tabres(ji,jj) = ssh n(ji,jj)1019 tabres(ji,jj) = ssh(ji,jj,Kmm_a) 996 1020 END DO 997 1021 END DO … … 1000 1024 DO jj=j1,j2 1001 1025 DO ji=i1,i2 1002 ssh b(ji,jj) = sshb(ji,jj) &1003 & + atfp * ( tabres(ji,jj) - ssh n(ji,jj) ) * tmask(ji,jj,1)1026 ssh(ji,jj,Kbb_a) = ssh(ji,jj,Kbb_a) & 1027 & + atfp * ( tabres(ji,jj) - ssh(ji,jj,Kmm_a) ) * tmask(ji,jj,1) 1004 1028 END DO 1005 1029 END DO … … 1008 1032 DO jj=j1,j2 1009 1033 DO ji=i1,i2 1010 ssh n(ji,jj) = tabres(ji,jj) * tmask(ji,jj,1)1034 ssh(ji,jj,Kmm_a) = tabres(ji,jj) * tmask(ji,jj,1) 1011 1035 END DO 1012 1036 END DO 1013 1037 ! 1014 1038 IF ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN 1015 ssh b(i1:i2,j1:j2) = sshn(i1:i2,j1:j2)1039 ssh(i1:i2,j1:j2,Kbb_a) = ssh(i1:i2,j1:j2,Kmm_a) 1016 1040 ENDIF 1017 1041 ! … … 1094 1118 DO jj=j1,j2 1095 1119 zcor = rdt * r1_e1e2t(i1 ,jj) * e2u(i1,jj) * (ub2_b(i1,jj)-tabres(i1,jj)) 1096 ssh n(i1 ,jj) = sshn(i1 ,jj) + zcor1097 IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) ssh b(i1 ,jj) = sshb(i1 ,jj) + atfp * zcor1120 ssh(i1 ,jj,Kmm_a) = ssh(i1 ,jj,Kmm_a) + zcor 1121 IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) ssh(i1 ,jj,Kbb_a) = ssh(i1 ,jj,Kbb_a) + atfp * zcor 1098 1122 END DO 1099 1123 ENDIF … … 1101 1125 DO jj=j1,j2 1102 1126 zcor = - rdt * r1_e1e2t(i2+1,jj) * e2u(i2,jj) * (ub2_b(i2,jj)-tabres(i2,jj)) 1103 ssh n(i2+1,jj) = sshn(i2+1,jj) + zcor1104 IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) ssh b(i2+1,jj) = sshb(i2+1,jj) + atfp * zcor1127 ssh(i2+1,jj,Kmm_a) = ssh(i2+1,jj,Kmm_a) + zcor 1128 IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) ssh(i2+1,jj,Kbb_a) = ssh(i2+1,jj,Kbb_a) + atfp * zcor 1105 1129 END DO 1106 1130 ENDIF … … 1182 1206 DO ji=i1,i2 1183 1207 zcor = rdt * r1_e1e2t(ji,j1 ) * e1v(ji,j1 ) * (vb2_b(ji,j1)-tabres(ji,j1)) 1184 ssh n(ji,j1 ) = sshn(ji,j1) + zcor1185 IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) ssh b(ji,j1 ) = sshb(ji,j1) + atfp * zcor1208 ssh(ji,j1 ,Kmm_a) = ssh(ji,j1 ,Kmm_a) + zcor 1209 IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) ssh(ji,j1 ,Kbb_a) = ssh(ji,j1,Kbb_a) + atfp * zcor 1186 1210 END DO 1187 1211 ENDIF … … 1189 1213 DO ji=i1,i2 1190 1214 zcor = - rdt * r1_e1e2t(ji,j2+1) * e1v(ji,j2 ) * (vb2_b(ji,j2)-tabres(ji,j2)) 1191 ssh n(ji,j2+1) = sshn(ji,j2+1) + zcor1192 IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) ssh b(ji,j2+1) = sshb(ji,j2+1) + atfp * zcor1215 ssh(ji,j2+1,Kmm_a) = ssh(ji,j2+1,Kmm_a) + zcor 1216 IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) ssh(ji,j2+1,Kbb_a) = ssh(ji,j2+1,Kbb_a) + atfp * zcor 1193 1217 END DO 1194 1218 ENDIF … … 1319 1343 DO jj=j1,j2 1320 1344 DO ji=i1,i2 1321 ptab(ji,jj,jk) = e3t_0(ji,jj,jk) * (1._wp + ssh n(ji,jj) &1345 ptab(ji,jj,jk) = e3t_0(ji,jj,jk) * (1._wp + ssh(ji,jj,Kmm_a) & 1322 1346 & *ssmask(ji,jj)/(ht_0(ji,jj)-1._wp + ssmask(ji,jj))) 1323 1347 END DO … … 1330 1354 ! Save "old" scale factor (prior update) for subsequent asselin correction 1331 1355 ! of prognostic variables 1332 e3t _a(i1:i2,j1:j2,1:jpkm1) = e3t_n(i1:i2,j1:j2,1:jpkm1)1333 1334 ! One should also save e3t _b, but lacking of workspace...1335 ! hdiv n(i1:i2,j1:j2,1:jpkm1) = e3t_b(i1:i2,j1:j2,1:jpkm1)1356 e3t(i1:i2,j1:j2,1:jpkm1,Krhs_a) = e3t(i1:i2,j1:j2,1:jpkm1,Kmm_a) 1357 1358 ! One should also save e3t(:,:,:,Kbb_a), but lacking of workspace... 1359 ! hdiv(i1:i2,j1:j2,1:jpkm1) = e3t(i1:i2,j1:j2,1:jpkm1,Kbb_a) 1336 1360 1337 1361 IF (.NOT.(lk_agrif_fstep.AND.(neuler==0) )) THEN … … 1339 1363 DO jj=j1,j2 1340 1364 DO ji=i1,i2 1341 e3t _b(ji,jj,jk) = e3t_b(ji,jj,jk) &1342 & + atfp * ( ptab(ji,jj,jk) - e3t _n(ji,jj,jk) )1365 e3t(ji,jj,jk,Kbb_a) = e3t(ji,jj,jk,Kbb_a) & 1366 & + atfp * ( ptab(ji,jj,jk) - e3t(ji,jj,jk,Kmm_a) ) 1343 1367 END DO 1344 1368 END DO 1345 1369 END DO 1346 1370 ! 1347 e3w _b (i1:i2,j1:j2,1) = e3w_0(i1:i2,j1:j2,1) + e3t_b(i1:i2,j1:j2,1) - e3t_0(i1:i2,j1:j2,1)1348 gdepw _b(i1:i2,j1:j2,1) = 0.0_wp1349 gdept _b(i1:i2,j1:j2,1) = 0.5_wp * e3w_b(i1:i2,j1:j2,1)1371 e3w (i1:i2,j1:j2,1,Kbb_a) = e3w_0(i1:i2,j1:j2,1) + e3t(i1:i2,j1:j2,1,Kbb_a) - e3t_0(i1:i2,j1:j2,1) 1372 gdepw(i1:i2,j1:j2,1,Kbb_a) = 0.0_wp 1373 gdept(i1:i2,j1:j2,1,Kbb_a) = 0.5_wp * e3w(i1:i2,j1:j2,1,Kbb_a) 1350 1374 ! 1351 1375 DO jk = 2, jpk … … 1353 1377 DO ji = i1,i2 1354 1378 zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk)) 1355 e3w _b(ji,jj,jk) = e3w_0(ji,jj,jk) + ( 1.0_wp - 0.5_wp * tmask(ji,jj,jk) ) * &1356 & ( e3t _b(ji,jj,jk-1) - e3t_0(ji,jj,jk-1) ) &1379 e3w(ji,jj,jk,Kbb_a) = e3w_0(ji,jj,jk) + ( 1.0_wp - 0.5_wp * tmask(ji,jj,jk) ) * & 1380 & ( e3t(ji,jj,jk-1,Kbb_a) - e3t_0(ji,jj,jk-1) ) & 1357 1381 & + 0.5_wp * tmask(ji,jj,jk) * & 1358 & ( e3t _b(ji,jj,jk) - e3t_0(ji,jj,jk ) )1359 gdepw _b(ji,jj,jk) = gdepw_b(ji,jj,jk-1) + e3t_b(ji,jj,jk-1)1360 gdept _b(ji,jj,jk) = zcoef * ( gdepw_b(ji,jj,jk ) + 0.5 * e3w_b(ji,jj,jk)) &1361 & + (1-zcoef) * ( gdept _b(ji,jj,jk-1) + e3w_b(ji,jj,jk))1382 & ( e3t(ji,jj,jk ,Kbb_a) - e3t_0(ji,jj,jk ) ) 1383 gdepw(ji,jj,jk,Kbb_a) = gdepw(ji,jj,jk-1,Kbb_a) + e3t(ji,jj,jk-1,Kbb_a) 1384 gdept(ji,jj,jk,Kbb_a) = zcoef * ( gdepw(ji,jj,jk ,Kbb_a) + 0.5 * e3w(ji,jj,jk,Kbb_a)) & 1385 & + (1-zcoef) * ( gdept(ji,jj,jk-1,Kbb_a) + e3w(ji,jj,jk,Kbb_a)) 1362 1386 END DO 1363 1387 END DO … … 1370 1394 ! 1371 1395 ! Update vertical scale factor at T-points: 1372 e3t _n(i1:i2,j1:j2,1:jpkm1) = ptab(i1:i2,j1:j2,1:jpkm1)1396 e3t(i1:i2,j1:j2,1:jpkm1,Kmm_a) = ptab(i1:i2,j1:j2,1:jpkm1) 1373 1397 ! 1374 1398 ! Update total depth: 1375 ht _n(i1:i2,j1:j2) = 0._wp1399 ht(i1:i2,j1:j2) = 0._wp 1376 1400 DO jk = 1, jpkm1 1377 ht _n(i1:i2,j1:j2) = ht_n(i1:i2,j1:j2) + e3t_n(i1:i2,j1:j2,jk) * tmask(i1:i2,j1:j2,jk)1401 ht(i1:i2,j1:j2) = ht(i1:i2,j1:j2) + e3t(i1:i2,j1:j2,jk,Kmm_a) * tmask(i1:i2,j1:j2,jk) 1378 1402 END DO 1379 1403 ! 1380 1404 ! Update vertical scale factor at W-points and depths: 1381 e3w _n (i1:i2,j1:j2,1) = e3w_0(i1:i2,j1:j2,1) + e3t_n(i1:i2,j1:j2,1) - e3t_0(i1:i2,j1:j2,1)1382 gdept _n(i1:i2,j1:j2,1) = 0.5_wp * e3w_n(i1:i2,j1:j2,1)1383 gdepw _n(i1:i2,j1:j2,1) = 0.0_wp1384 gde3w _n(i1:i2,j1:j2,1) = gdept_n(i1:i2,j1:j2,1) - (ht_n(i1:i2,j1:j2)-ht_0(i1:i2,j1:j2)) ! Last term in the rhs is ssh1405 e3w (i1:i2,j1:j2,1,Kmm_a) = e3w_0(i1:i2,j1:j2,1) + e3t(i1:i2,j1:j2,1,Kmm_a) - e3t_0(i1:i2,j1:j2,1) 1406 gdept(i1:i2,j1:j2,1,Kmm_a) = 0.5_wp * e3w(i1:i2,j1:j2,1,Kmm_a) 1407 gdepw(i1:i2,j1:j2,1,Kmm_a) = 0.0_wp 1408 gde3w(i1:i2,j1:j2,1) = gdept(i1:i2,j1:j2,1,Kmm_a) - (ht(i1:i2,j1:j2)-ht_0(i1:i2,j1:j2)) ! Last term in the rhs is ssh 1385 1409 ! 1386 1410 DO jk = 2, jpk … … 1388 1412 DO ji = i1,i2 1389 1413 zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk)) 1390 e3w _n(ji,jj,jk) = e3w_0(ji,jj,jk) + ( 1.0_wp - 0.5_wp * tmask(ji,jj,jk) ) * ( e3t_n(ji,jj,jk-1) - e3t_0(ji,jj,jk-1) ) &1391 & + 0.5_wp * tmask(ji,jj,jk) * ( e3t _n(ji,jj,jk) - e3t_0(ji,jj,jk ) )1392 gdepw _n(ji,jj,jk) = gdepw_n(ji,jj,jk-1) + e3t_n(ji,jj,jk-1)1393 gdept _n(ji,jj,jk) = zcoef * ( gdepw_n(ji,jj,jk ) + 0.5 * e3w_n(ji,jj,jk)) &1394 & + (1-zcoef) * ( gdept _n(ji,jj,jk-1) + e3w_n(ji,jj,jk))1395 gde3w _n(ji,jj,jk) = gdept_n(ji,jj,jk) - (ht_n(ji,jj)-ht_0(ji,jj)) ! Last term in the rhs is ssh1414 e3w(ji,jj,jk,Kmm_a) = e3w_0(ji,jj,jk) + ( 1.0_wp - 0.5_wp * tmask(ji,jj,jk) ) * ( e3t(ji,jj,jk-1,Kmm_a) - e3t_0(ji,jj,jk-1) ) & 1415 & + 0.5_wp * tmask(ji,jj,jk) * ( e3t(ji,jj,jk ,Kmm_a) - e3t_0(ji,jj,jk ) ) 1416 gdepw(ji,jj,jk,Kmm_a) = gdepw(ji,jj,jk-1,Kmm_a) + e3t(ji,jj,jk-1,Kmm_a) 1417 gdept(ji,jj,jk,Kmm_a) = zcoef * ( gdepw(ji,jj,jk ,Kmm_a) + 0.5 * e3w(ji,jj,jk,Kmm_a)) & 1418 & + (1-zcoef) * ( gdept(ji,jj,jk-1,Kmm_a) + e3w(ji,jj,jk,Kmm_a)) 1419 gde3w(ji,jj,jk) = gdept(ji,jj,jk,Kmm_a) - (ht(ji,jj)-ht_0(ji,jj)) ! Last term in the rhs is ssh 1396 1420 END DO 1397 1421 END DO … … 1399 1423 ! 1400 1424 IF ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN 1401 e3t _b (i1:i2,j1:j2,1:jpk) = e3t_n (i1:i2,j1:j2,1:jpk)1402 e3w _b (i1:i2,j1:j2,1:jpk) = e3w_n (i1:i2,j1:j2,1:jpk)1403 gdepw _b(i1:i2,j1:j2,1:jpk) = gdepw_n(i1:i2,j1:j2,1:jpk)1404 gdept _b(i1:i2,j1:j2,1:jpk) = gdept_n(i1:i2,j1:j2,1:jpk)1425 e3t (i1:i2,j1:j2,1:jpk,Kbb_a) = e3t (i1:i2,j1:j2,1:jpk,Kmm_a) 1426 e3w (i1:i2,j1:j2,1:jpk,Kbb_a) = e3w (i1:i2,j1:j2,1:jpk,Kmm_a) 1427 gdepw(i1:i2,j1:j2,1:jpk,Kbb_a) = gdepw(i1:i2,j1:j2,1:jpk,Kmm_a) 1428 gdept(i1:i2,j1:j2,1:jpk,Kbb_a) = gdept(i1:i2,j1:j2,1:jpk,Kmm_a) 1405 1429 ENDIF 1406 1430 !
Note: See TracChangeset
for help on using the changeset viewer.
