- Timestamp:
- 2020-09-14T17:40:34+02:00 (4 years ago)
- Location:
- NEMO/branches/2019/dev_r11351_fldread_with_XIOS
- Files:
-
- 2 edited
-
. (modified) (1 prop)
-
src/NST/agrif_oce_update.F90 (modified) (62 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_r11351_fldread_with_XIOS
- Property svn:externals
-
old new 3 3 ^/utils/build/mk@HEAD mk 4 4 ^/utils/tools@HEAD tools 5 ^/vendors/AGRIF/dev @HEADext/AGRIF5 ^/vendors/AGRIF/dev_r12970_AGRIF_CMEMS ext/AGRIF 6 6 ^/vendors/FCM@HEAD ext/FCM 7 7 ^/vendors/IOIPSL@HEAD ext/IOIPSL 8 9 # SETTE 10 ^/utils/CI/sette@13382 sette
-
- Property svn:externals
-
NEMO/branches/2019/dev_r11351_fldread_with_XIOS/src/NST/agrif_oce_update.F90
r10068 r13463 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 28 USE lbclnk 27 29 28 30 IMPLICIT NONE … … 46 48 IF (Agrif_Root()) RETURN 47 49 ! 48 #if defined TWO_WAY49 50 IF (lwp.AND.lk_agrif_debug) Write(*,*) 'Update tracers from grid Number',Agrif_Fixed() 50 51 52 #if defined key_vertical 53 ! Effect of this has to be carrefully checked 54 ! depending on what the nesting tools ensure for 55 ! volume conservation: 56 Agrif_UseSpecialValueInUpdate = .FALSE. 57 #else 51 58 Agrif_UseSpecialValueInUpdate = .TRUE. 59 #endif 52 60 Agrif_SpecialValueFineGrid = 0._wp 53 61 ! … … 64 72 Agrif_UseSpecialValueInUpdate = .FALSE. 65 73 ! 66 #endif67 74 ! 68 75 END SUBROUTINE Agrif_Update_Tra … … 75 82 IF (Agrif_Root()) RETURN 76 83 ! 77 #if defined TWO_WAY78 84 IF (lwp.AND.lk_agrif_debug) Write(*,*) 'Update momentum from grid Number',Agrif_Fixed() 79 85 80 86 Agrif_UseSpecialValueInUpdate = .FALSE. 81 Agrif_SpecialValueFineGrid = 0. 87 Agrif_SpecialValueFineGrid = 0._wp 88 89 use_sign_north = .TRUE. 90 sign_north = -1._wp 91 82 92 ! 83 93 # if ! defined DECAL_FEEDBACK … … 95 105 # endif 96 106 97 # if ! defined DECAL_FEEDBACK 107 # if ! defined DECAL_FEEDBACK_2D 98 108 CALL Agrif_Update_Variable(e1u_id,procname = updateU2d) 99 109 CALL Agrif_Update_Variable(e2v_id,procname = updateV2d) … … 103 113 # endif 104 114 ! 105 # if ! defined DECAL_FEEDBACK 115 # if ! defined DECAL_FEEDBACK_2D 106 116 ! Account for updated thicknesses at boundary edges 107 117 IF (.NOT.ln_linssh) THEN … … 113 123 IF ( ln_dynspg_ts .AND. ln_bt_fw ) THEN 114 124 ! Update time integrated transports 115 # if ! defined DECAL_FEEDBACK 125 # if ! defined DECAL_FEEDBACK_2D 116 126 CALL Agrif_Update_Variable(ub2b_update_id,procname = updateub2b) 117 127 CALL Agrif_Update_Variable(vb2b_update_id,procname = updatevb2b) … … 121 131 # endif 122 132 END IF 123 #endif 133 ! 134 use_sign_north = .FALSE. 124 135 ! 125 136 END SUBROUTINE Agrif_Update_Dyn … … 132 143 IF (Agrif_Root()) RETURN 133 144 ! 134 #if defined TWO_WAY135 !136 145 Agrif_UseSpecialValueInUpdate = .TRUE. 137 Agrif_SpecialValueFineGrid = 0. 138 # if ! defined DECAL_FEEDBACK 146 Agrif_SpecialValueFineGrid = 0._wp 147 # if ! defined DECAL_FEEDBACK_2D 139 148 CALL Agrif_Update_Variable(sshn_id,procname = updateSSH) 140 149 # else … … 146 155 # if defined VOL_REFLUX 147 156 IF ( ln_dynspg_ts.AND.ln_bt_fw ) THEN 157 use_sign_north = .TRUE. 158 sign_north = -1._wp 148 159 ! Refluxing on ssh: 149 # if defined DECAL_FEEDBACK 160 # if defined DECAL_FEEDBACK_2D 150 161 CALL Agrif_Update_Variable(ub2b_update_id,locupdate1=(/0, 0/),locupdate2=(/1, 1/),procname = reflux_sshu) 151 162 CALL Agrif_Update_Variable(vb2b_update_id,locupdate1=(/1, 1/),locupdate2=(/0, 0/),procname = reflux_sshv) … … 154 165 CALL Agrif_Update_Variable(vb2b_update_id,locupdate1=(/ 0, 0/),locupdate2=(/-1,-1/),procname = reflux_sshv) 155 166 # endif 167 use_sign_north = .FALSE. 156 168 END IF 157 169 # endif 158 170 ! 159 #endif160 !161 171 END SUBROUTINE Agrif_Update_ssh 162 172 … … 170 180 IF (Agrif_Root()) RETURN 171 181 ! 172 # if defined TWO_WAY173 174 182 Agrif_UseSpecialValueInUpdate = .TRUE. 175 183 Agrif_SpecialValueFineGrid = 0. … … 180 188 181 189 Agrif_UseSpecialValueInUpdate = .FALSE. 182 183 # endif184 190 185 191 END SUBROUTINE Agrif_Update_Tke … … 192 198 ! 193 199 IF (Agrif_Root()) RETURN 194 !195 #if defined TWO_WAY196 200 ! 197 201 IF (lwp.AND.lk_agrif_debug) Write(*,*) 'Update e3 from grid Number',Agrif_Fixed(), 'Step', Agrif_Nb_Step() … … 210 214 CALL Agrif_ParentGrid_To_ChildGrid() 211 215 ! 212 #endif213 !214 216 END SUBROUTINE Agrif_Update_vvl 215 217 … … 230 232 ! ----------------------- 231 233 ! 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(:,:)234 e3u(:,:,:,Krhs_a) = e3u(:,:,:,Kmm_a) 235 e3v(:,:,:,Krhs_a) = e3v(:,:,:,Kmm_a) 236 ! uu(:,:,:,Krhs_a) = e3u(:,:,:,Kbb_a) 237 ! vv(:,:,:,Krhs_a) = e3v(:,:,:,Kbb_a) 238 hu(:,:,Krhs_a) = hu(:,:,Kmm_a) 239 hv(:,:,Krhs_a) = hv(:,:,Kmm_a) 238 240 239 241 ! 1) NOW fields … … 242 244 ! Vertical scale factor interpolations 243 245 ! ------------------------------------ 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' )246 CALL dom_vvl_interpol( e3t(:,:,:,Kmm_a), e3u(:,:,:,Kmm_a) , 'U' ) 247 CALL dom_vvl_interpol( e3t(:,:,:,Kmm_a), e3v(:,:,:,Kmm_a) , 'V' ) 248 CALL dom_vvl_interpol( e3u(:,:,:,Kmm_a), e3f(:,:,:) , 'F' ) 249 250 CALL dom_vvl_interpol( e3u(:,:,:,Kmm_a), e3uw(:,:,:,Kmm_a), 'UW' ) 251 CALL dom_vvl_interpol( e3v(:,:,:,Kmm_a), e3vw(:,:,:,Kmm_a), 'VW' ) 250 252 251 253 ! Update total depths: 252 254 ! -------------------- 253 hu _n(:,:) = 0._wp ! Ocean depth at U-points254 hv _n(:,:) = 0._wp ! Ocean depth at V-points255 hu(:,:,Kmm_a) = 0._wp ! Ocean depth at U-points 256 hv(:,:,Kmm_a) = 0._wp ! Ocean depth at V-points 255 257 DO jk = 1, jpkm1 256 hu _n(:,:) = hu_n(:,:) + e3u_n(:,:,jk) * umask(:,:,jk)257 hv _n(:,:) = hv_n(:,:) + e3v_n(:,:,jk) * vmask(:,:,jk)258 hu(:,:,Kmm_a) = hu(:,:,Kmm_a) + e3u(:,:,jk,Kmm_a) * umask(:,:,jk) 259 hv(:,:,Kmm_a) = hv(:,:,Kmm_a) + e3v(:,:,jk,Kmm_a) * vmask(:,:,jk) 258 260 END DO 259 261 ! ! 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(:,:) )262 r1_hu(:,:,Kmm_a) = ssumask(:,:) / ( hu(:,:,Kmm_a) + 1._wp - ssumask(:,:) ) 263 r1_hv(:,:,Kmm_a) = ssvmask(:,:) / ( hv(:,:,Kmm_a) + 1._wp - ssvmask(:,:) ) 262 264 263 265 264 266 ! 2) BEFORE fields: 265 267 !------------------ 266 IF (.NOT.(lk_agrif_fstep.AND.( neuler==0) )) THEN268 IF (.NOT.(lk_agrif_fstep.AND.(l_1st_euler) )) THEN 267 269 ! 268 270 ! Vertical scale factor interpolations 269 271 ! ------------------------------------ 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' )272 CALL dom_vvl_interpol( e3t(:,:,:,Kbb_a), e3u(:,:,:,Kbb_a), 'U' ) 273 CALL dom_vvl_interpol( e3t(:,:,:,Kbb_a), e3v(:,:,:,Kbb_a), 'V' ) 274 275 CALL dom_vvl_interpol( e3u(:,:,:,Kbb_a), e3uw(:,:,:,Kbb_a), 'UW' ) 276 CALL dom_vvl_interpol( e3v(:,:,:,Kbb_a), e3vw(:,:,:,Kbb_a), 'VW' ) 275 277 276 278 ! Update total depths: 277 279 ! -------------------- 278 hu _b(:,:) = 0._wp ! Ocean depth at U-points279 hv _b(:,:) = 0._wp ! Ocean depth at V-points280 hu(:,:,Kbb_a) = 0._wp ! Ocean depth at U-points 281 hv(:,:,Kbb_a) = 0._wp ! Ocean depth at V-points 280 282 DO jk = 1, jpkm1 281 hu _b(:,:) = hu_b(:,:) + e3u_b(:,:,jk) * umask(:,:,jk)282 hv _b(:,:) = hv_b(:,:) + e3v_b(:,:,jk) * vmask(:,:,jk)283 hu(:,:,Kbb_a) = hu(:,:,Kbb_a) + e3u(:,:,jk,Kbb_a) * umask(:,:,jk) 284 hv(:,:,Kbb_a) = hv(:,:,Kbb_a) + e3v(:,:,jk,Kbb_a) * vmask(:,:,jk) 283 285 END DO 284 286 ! ! 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(:,:) )287 r1_hu(:,:,Kbb_a) = ssumask(:,:) / ( hu(:,:,Kbb_a) + 1._wp - ssumask(:,:) ) 288 r1_hv(:,:,Kbb_a) = ssvmask(:,:) / ( hv(:,:,Kbb_a) + 1._wp - ssvmask(:,:) ) 287 289 ENDIF 288 290 ! … … 300 302 !! 301 303 INTEGER :: ji,jj,jk,jn 302 REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk,n1:n2) :: tabres_child 304 INTEGER :: N_in, N_out 305 REAL(wp) :: ztb, ztnu, ztno 303 306 REAL(wp) :: h_in(k1:k2) 304 307 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) 308 REAL(wp) :: tabin(k1:k2,1:jpts) 309 REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk,1:jpts) :: tabres_child 308 310 !!--------------------------------------------- 309 311 ! 310 312 IF (before) THEN 311 AGRIF_SpecialValue = -999._wp 312 zrho_xy = Agrif_rhox() * Agrif_rhoy() 313 !jc_alt 314 ! AGRIF_SpecialValue = -999._wp 313 315 DO jn = n1,n2-1 314 316 DO jk=k1,k2 315 317 DO jj=j1,j2 316 318 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 319 !jc_alt 320 ! tabres(ji,jj,jk,jn) = (ts(ji,jj,jk,jn,Kmm_a) * e3t(ji,jj,jk,Kmm_a) ) & 321 ! & * tmask(ji,jj,jk) + (tmask(ji,jj,jk)-1._wp) * 999._wp 322 tabres(ji,jj,jk,jn) = ts(ji,jj,jk,jn,Kmm_a) * e3t(ji,jj,jk,Kmm_a) 319 323 END DO 320 324 END DO … … 324 328 DO jj=j1,j2 325 329 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 330 !jc_alt 331 ! tabres(ji,jj,jk,n2) = tmask(ji,jj,jk) * e3t(ji,jj,jk,Kmm_a) & 332 ! & + (tmask(ji,jj,jk) - 1._wp) * 999._wp 333 tabres(ji,jj,jk,n2) = tmask(ji,jj,jk) * e3t(ji,jj,jk,Kmm_a) 328 334 END DO 329 335 END DO 330 336 END DO 331 337 ELSE 332 tabres_child(:,:,:,:) = 0. 338 tabres_child(:,:,:,:) = 0._wp 333 339 AGRIF_SpecialValue = 0._wp 334 340 DO jj=j1,j2 … … 336 342 N_in = 0 337 343 DO jk=k1,k2 !k2 = jpk of child grid 338 IF (tabres(ji,jj,jk,n2) == 0 ) EXIT 344 ! jc_alt 345 ! IF (tabres(ji,jj,jk,n2) < -900._wp ) EXIT 346 IF (tabres(ji,jj,jk,n2) == 0._wp ) EXIT 339 347 N_in = N_in + 1 340 348 tabin(jk,:) = tabres(ji,jj,jk,n1:n2-1)/tabres(ji,jj,jk,n2) … … 343 351 N_out = 0 344 352 DO jk=1,jpk ! jpk of parent grid 345 IF (tmask(ji,jj,jk) < -900) EXIT ! TODO: Will not work with ISF353 IF (tmask(ji,jj,jk) == 0 ) EXIT ! TODO: Will not work with ISF 346 354 N_out = N_out + 1 347 h_out(N_out) = e3t _n(ji,jj,jk)355 h_out(N_out) = e3t(ji,jj,jk,Kmm_a) 348 356 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 357 IF (N_in*N_out > 0) THEN !Remove this? 358 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 359 ENDIF 361 360 ENDDO 362 361 ENDDO 363 362 364 IF (.NOT.(lk_agrif_fstep.AND.( neuler==0))) THEN363 IF (.NOT.(lk_agrif_fstep.AND.(l_1st_euler))) THEN 365 364 ! 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) 365 DO jn = 1,jpts 366 DO jk = 1, jpkm1 367 DO jj = j1, j2 368 DO ji = i1, i2 369 IF( tabres_child(ji,jj,jk,jn) /= 0._wp ) THEN 370 ztb = ts(ji,jj,jk,jn,Kbb_a) * e3t(ji,jj,jk,Kbb_a) ! fse3t_b prior update should be used 371 ztnu = tabres_child(ji,jj,jk,jn) * e3t(ji,jj,jk,Kmm_a) 372 ztno = ts(ji,jj,jk,jn,Kmm_a) * e3t(ji,jj,jk,Krhs_a) 373 ts(ji,jj,jk,jn,Kbb_a) = ( ztb + rn_atfp * ( ztnu - ztno) ) & 374 & * tmask(ji,jj,jk) / e3t(ji,jj,jk,Kbb_a) 374 375 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)376 END DO 377 END DO 378 END DO 379 END DO 380 ENDIF 381 DO jn = 1,jpts 382 DO jk = 1, jpkm1 383 DO jj = j1, j2 384 DO ji = i1, i2 385 IF( tabres_child(ji,jj,jk,jn) /= 0._wp ) THEN 386 ts(ji,jj,jk,jn,Kmm_a) = tabres_child(ji,jj,jk,jn) 386 387 END IF 387 388 END DO … … 389 390 END DO 390 391 END DO 392 ! 393 IF ((l_1st_euler).AND.(Agrif_Nb_Step()==0) ) THEN 394 ts(i1:i2,j1:j2,1:jpkm1,1:jpts,Kbb_a) = ts(i1:i2,j1:j2,1:jpkm1,1:jpts,Kmm_a) 395 ENDIF 391 396 ENDIF 392 397 ! … … 413 418 DO ji=i1,i2 414 419 !> 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)420 tabres(ji,jj,jk,jn) = ts(ji,jj,jk,jn,Kmm_a) * e3t(ji,jj,jk,Kmm_a) / e3t_0(ji,jj,jk) 421 ! tabres(ji,jj,jk,jn) = ts(ji,jj,jk,jn,Kmm_a) * e3t(ji,jj,jk,Kmm_a) 417 422 !< jc tmp 418 423 END DO … … 427 432 ENDDO 428 433 !< jc tmp 429 IF (.NOT.(lk_agrif_fstep.AND.( neuler==0))) THEN434 IF (.NOT.(lk_agrif_fstep.AND.(l_1st_euler))) THEN 430 435 ! Add asselin part 431 436 DO jn = 1,jpts … … 434 439 DO ji = i1, i2 435 440 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 used441 ztb = ts(ji,jj,jk,jn,Kbb_a) * e3t(ji,jj,jk,Kbb_a) ! fse3t_b prior update should be used 437 442 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)443 ztno = ts(ji,jj,jk,jn,Kmm_a) * e3t(ji,jj,jk,Krhs_a) 444 ts(ji,jj,jk,jn,Kbb_a) = ( ztb + rn_atfp * ( ztnu - ztno) ) & 445 & * tmask(ji,jj,jk) / e3t(ji,jj,jk,Kbb_a) 441 446 ENDIF 442 447 END DO … … 450 455 DO ji=i1,i2 451 456 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)457 ts(ji,jj,jk,jn,Kmm_a) = tabres(ji,jj,jk,jn) / e3t(ji,jj,jk,Kmm_a) 453 458 END IF 454 459 END DO … … 457 462 END DO 458 463 ! 459 IF (( neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN460 ts b(i1:i2,j1:j2,k1:k2,1:jpts) = tsn(i1:i2,j1:j2,k1:k2,1:jpts)464 IF ((l_1st_euler).AND.(Agrif_Nb_Step()==0) ) THEN 465 ts(i1:i2,j1:j2,k1:k2,1:jpts,Kbb_a) = ts(i1:i2,j1:j2,k1:k2,1:jpts,Kmm_a) 461 466 ENDIF 462 467 ! … … 478 483 ! 479 484 INTEGER :: ji, jj, jk 480 REAL(wp):: zrhoy 485 REAL(wp):: zrhoy, zub, zunu, zuno 481 486 ! VERTICAL REFINEMENT BEGIN 482 487 REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: tabres_child … … 491 496 IF( before ) THEN 492 497 zrhoy = Agrif_Rhoy() 493 AGRIF_SpecialValue = -999._wp 498 !jc_alt 499 ! AGRIF_SpecialValue = -999._wp 494 500 DO jk=k1,k2 495 501 DO jj=j1,j2 496 502 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 503 !jc_alt 504 ! 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) & 505 ! & + (umask(ji,jj,jk)-1._wp)*999._wp 506 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) 507 !jc_alt 508 ! tabres(ji,jj,jk,2) = zrhoy * umask(ji,jj,jk) * e2u(ji,jj) * e3u(ji,jj,jk,Kmm_a) & 509 ! & + (umask(ji,jj,jk)-1._wp)*999._wp 510 tabres(ji,jj,jk,2) = zrhoy * umask(ji,jj,jk) * e2u(ji,jj) * e3u(ji,jj,jk,Kmm_a) 501 511 END DO 502 512 END DO … … 511 521 tabin(:) = 0._wp 512 522 DO jk=k1,k2 !k2=jpk of child grid 513 IF( tabres(ji,jj,jk,2) < -900) EXIT 523 !jc_alt 524 ! IF( tabres(ji,jj,jk,2) < -900._wp) EXIT 525 IF( tabres(ji,jj,jk,2) == 0.) EXIT 514 526 N_in = N_in + 1 515 527 tabin(jk) = tabres(ji,jj,jk,1)/tabres(ji,jj,jk,2) … … 520 532 IF (umask(ji,jj,jk) == 0) EXIT 521 533 N_out = N_out + 1 522 h_out(N_out) = e3u _n(ji,jj,jk)534 h_out(N_out) = e3u(ji,jj,jk,Kmm_a) 523 535 ENDDO 524 536 IF (N_in * N_out > 0) THEN 525 537 h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in)) 538 excess = 0._wp 526 539 IF (h_diff < -1.e-4) THEN 527 540 !Even if bathy at T points match it's possible for the U points to be deeper in the child grid. 528 541 !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 542 DO jk=N_in,1,-1 531 543 thick = MIN(-1*h_diff, h_in(jk)) … … 540 552 ENDDO 541 553 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 )554 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 555 tabres_child(ji,jj,N_out) = tabres_child(ji,jj,N_out) + excess/(e2u(ji,jj)*h_out(N_out)) 544 556 ENDIF 545 557 ENDDO 546 558 ENDDO 547 559 ! 548 560 DO jk=1,jpk 549 561 DO jj=j1,j2 550 562 DO ji=i1,i2 551 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) 563 IF (.NOT.(lk_agrif_fstep.AND.(l_1st_euler))) THEN ! Add asselin part 564 zub = uu(ji,jj,jk,Kbb_a) * e3u(ji,jj,jk,Kbb_a) ! fse3t_b prior update should be used 565 zuno = uu(ji,jj,jk,Kmm_a) * e3u(ji,jj,jk,Krhs_a) 566 zunu = tabres_child(ji,jj,jk) * e3u(ji,jj,jk,Kmm_a) 567 uu(ji,jj,jk,Kbb_a) = ( zub + rn_atfp * ( zunu - zuno) ) & 568 & * umask(ji,jj,jk) / e3u(ji,jj,jk,Kbb_a) 554 569 ENDIF 555 570 ! 556 un(ji,jj,jk) = tabres_child(ji,jj,jk) * umask(ji,jj,jk) 557 END DO 558 END DO 559 END DO 571 uu(ji,jj,jk,Kmm_a) = tabres_child(ji,jj,jk) * umask(ji,jj,jk) 572 END DO 573 END DO 574 END DO 575 ! 576 IF ((l_1st_euler).AND.(Agrif_Nb_Step()==0) ) THEN 577 uu(i1:i2,j1:j2,1:jpkm1,Kbb_a) = uu(i1:i2,j1:j2,1:jpkm1,Kmm_a) 578 ENDIF 579 ! 560 580 ENDIF 561 581 ! … … 579 599 zrhoy = Agrif_Rhoy() 580 600 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)601 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 602 END DO 583 603 ELSE … … 587 607 tabres(ji,jj,jk,1) = tabres(ji,jj,jk,1) * r1_e2u(ji,jj) 588 608 ! 589 IF (.NOT.(lk_agrif_fstep.AND.( neuler==0))) THEN ! Add asselin part590 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)609 IF (.NOT.(lk_agrif_fstep.AND.(l_1st_euler))) THEN ! Add asselin part 610 zub = uu(ji,jj,jk,Kbb_a) * e3u(ji,jj,jk,Kbb_a) ! fse3t_b prior update should be used 611 zuno = uu(ji,jj,jk,Kmm_a) * e3u(ji,jj,jk,Krhs_a) 592 612 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)613 uu(ji,jj,jk,Kbb_a) = ( zub + rn_atfp * ( zunu - zuno) ) & 614 & * umask(ji,jj,jk) / e3u(ji,jj,jk,Kbb_a) 595 615 ENDIF 596 616 ! 597 u n(ji,jj,jk) = tabres(ji,jj,jk,1) * umask(ji,jj,jk) / e3u_n(ji,jj,jk)598 END DO 599 END DO 600 END DO 601 ! 602 IF (( neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN603 u b(i1:i2,j1:j2,k1:k2) = un(i1:i2,j1:j2,k1:k2)617 uu(ji,jj,jk,Kmm_a) = tabres(ji,jj,jk,1) * umask(ji,jj,jk) / e3u(ji,jj,jk,Kmm_a) 618 END DO 619 END DO 620 END DO 621 ! 622 IF ((l_1st_euler).AND.(Agrif_Nb_Step()==0) ) THEN 623 uu(i1:i2,j1:j2,k1:k2,Kbb_a) = uu(i1:i2,j1:j2,k1:k2,Kmm_a) 604 624 ENDIF 605 625 ! … … 632 652 IF (western_side) THEN 633 653 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) + zcor654 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) 655 uu_b(i1-1,jj,Kmm_a) = uu_b(i1-1,jj,Kmm_a) + zcor 636 656 DO jk=1,jpkm1 637 u n(i1-1,jj,jk) = un(i1-1,jj,jk) + zcor * umask(i1-1,jj,jk)657 uu(i1-1,jj,jk,Kmm_a) = uu(i1-1,jj,jk,Kmm_a) + zcor * umask(i1-1,jj,jk) 638 658 END DO 639 659 END DO … … 642 662 IF (eastern_side) THEN 643 663 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) + zcor664 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) 665 uu_b(i2+1,jj,Kmm_a) = uu_b(i2+1,jj,Kmm_a) + zcor 646 666 DO jk=1,jpkm1 647 u n(i2+1,jj,jk) = un(i2+1,jj,jk) + zcor * umask(i2+1,jj,jk)667 uu(i2+1,jj,jk,Kmm_a) = uu(i2+1,jj,jk,Kmm_a) + zcor * umask(i2+1,jj,jk) 648 668 END DO 649 669 END DO … … 665 685 ! 666 686 INTEGER :: ji, jj, jk 667 REAL(wp) :: zrhox 687 REAL(wp) :: zrhox, zvb, zvnu, zvno 668 688 ! VERTICAL REFINEMENT BEGIN 669 689 REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: tabres_child … … 678 698 IF( before ) THEN 679 699 zrhox = Agrif_Rhox() 680 AGRIF_SpecialValue = -999._wp 700 !jc_alt 701 ! AGRIF_SpecialValue = -999._wp 681 702 DO jk=k1,k2 682 703 DO jj=j1,j2 683 704 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 705 !jc_alt 706 ! 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) & 707 ! & + (vmask(ji,jj,jk)-1._wp) * 999._wp 708 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) 709 !jc_alt 710 ! tabres(ji,jj,jk,2) = zrhox * e1v(ji,jj) * e3v(ji,jj,jk,Kmm_a) * vmask(ji,jj,jk) & 711 ! & + (vmask(ji,jj,jk)-1._wp) * 999._wp 712 tabres(ji,jj,jk,2) = zrhox * e1v(ji,jj) * e3v(ji,jj,jk,Kmm_a) * vmask(ji,jj,jk) 688 713 END DO 689 714 END DO … … 696 721 N_in = 0 697 722 DO jk=k1,k2 698 IF (tabres(ji,jj,jk,2) < -900) EXIT 723 !jc_alt 724 ! IF (tabres(ji,jj,jk,2) < -900._wp) EXIT 725 IF (tabres(ji,jj,jk,2) == 0) EXIT 699 726 N_in = N_in + 1 700 727 tabin(jk) = tabres(ji,jj,jk,1)/tabres(ji,jj,jk,2) … … 705 732 IF (vmask(ji,jj,jk) == 0) EXIT 706 733 N_out = N_out + 1 707 h_out(N_out) = e3v _n(ji,jj,jk)734 h_out(N_out) = e3v(ji,jj,jk,Kmm_a) 708 735 ENDDO 709 736 IF (N_in * N_out > 0) THEN 710 737 h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in)) 738 excess = 0._wp 711 739 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.740 !Even if bathy at T points match it's possible for the V points to be deeper in the child grid. 713 741 !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 742 DO jk=N_in,1,-1 716 743 thick = MIN(-1*h_diff, h_in(jk)) … … 725 752 ENDDO 726 753 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 )754 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 755 tabres_child(ji,jj,N_out) = tabres_child(ji,jj,N_out) + excess/(e1v(ji,jj)*h_out(N_out)) 729 756 ENDIF 730 757 ENDDO 731 758 ENDDO 732 733 DO jk=1,jpk 759 ! 760 DO jk=1,jpkm1 734 761 DO jj=j1,j2 735 762 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) 763 IF (.NOT.(lk_agrif_fstep.AND.(l_1st_euler))) THEN ! Add asselin part 764 zvb = vv(ji,jj,jk,Kbb_a) * e3v(ji,jj,jk,Kbb_a) ! fse3t_b prior update should be used 765 zvno = vv(ji,jj,jk,Kmm_a) * e3v(ji,jj,jk,Krhs_a) 766 zvnu = tabres_child(ji,jj,jk) * e3v(ji,jj,jk,Kmm_a) 767 vv(ji,jj,jk,Kbb_a) = ( zvb + rn_atfp * ( zvnu - zvno) ) & 768 & * vmask(ji,jj,jk) / e3v(ji,jj,jk,Kbb_a) 740 769 ENDIF 741 770 ! 742 vn(ji,jj,jk) = tabres_child(ji,jj,jk) * vmask(ji,jj,jk) 743 END DO 744 END DO 745 END DO 771 vv(ji,jj,jk,Kmm_a) = tabres_child(ji,jj,jk) * vmask(ji,jj,jk) 772 END DO 773 END DO 774 END DO 775 ! 776 IF ((l_1st_euler).AND.(Agrif_Nb_Step()==0) ) THEN 777 vv(i1:i2,j1:j2,1:jpkm1,Kbb_a) = vv(i1:i2,j1:j2,1:jpkm1,Kmm_a) 778 ENDIF 779 ! 746 780 ENDIF 747 781 ! … … 767 801 DO jj=j1,j2 768 802 DO ji=i1,i2 769 tabres(ji,jj,jk,1) = zrhox * e1v(ji,jj) * e3v _n(ji,jj,jk) * vn(ji,jj,jk)803 tabres(ji,jj,jk,1) = zrhox * e1v(ji,jj) * e3v(ji,jj,jk,Kmm_a) * vv(ji,jj,jk,Kmm_a) 770 804 END DO 771 805 END DO … … 777 811 tabres(ji,jj,jk,1) = tabres(ji,jj,jk,1) * r1_e1v(ji,jj) 778 812 ! 779 IF (.NOT.(lk_agrif_fstep.AND.( neuler==0))) THEN ! Add asselin part780 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)813 IF (.NOT.(lk_agrif_fstep.AND.(l_1st_euler))) THEN ! Add asselin part 814 zvb = vv(ji,jj,jk,Kbb_a) * e3v(ji,jj,jk,Kbb_a) ! fse3t_b prior update should be used 815 zvno = vv(ji,jj,jk,Kmm_a) * e3v(ji,jj,jk,Krhs_a) 782 816 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)817 vv(ji,jj,jk,Kbb_a) = ( zvb + rn_atfp * ( zvnu - zvno) ) & 818 & * vmask(ji,jj,jk) / e3v(ji,jj,jk,Kbb_a) 785 819 ENDIF 786 820 ! 787 v n(ji,jj,jk) = tabres(ji,jj,jk,1) * vmask(ji,jj,jk) / e3v_n(ji,jj,jk)788 END DO 789 END DO 790 END DO 791 ! 792 IF (( neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN793 v b(i1:i2,j1:j2,k1:k2) = vn(i1:i2,j1:j2,k1:k2)821 vv(ji,jj,jk,Kmm_a) = tabres(ji,jj,jk,1) * vmask(ji,jj,jk) / e3v(ji,jj,jk,Kmm_a) 822 END DO 823 END DO 824 END DO 825 ! 826 IF ((l_1st_euler).AND.(Agrif_Nb_Step()==0) ) THEN 827 vv(i1:i2,j1:j2,k1:k2,Kbb_a) = vv(i1:i2,j1:j2,k1:k2,Kmm_a) 794 828 ENDIF 795 829 ! … … 802 836 SUBROUTINE correct_v_bdy( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before, nb, ndir ) 803 837 !!--------------------------------------------- 804 !! *** ROUTINE correct_ u_bdy ***838 !! *** ROUTINE correct_v_bdy *** 805 839 !!--------------------------------------------- 806 840 INTEGER , INTENT(in ) :: i1, i2, j1, j2, k1, k2, n1, n2 … … 822 856 IF (southern_side) THEN 823 857 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) + zcor858 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) 859 vv_b(ji,j1-1,Kmm_a) = vv_b(ji,j1-1,Kmm_a) + zcor 826 860 DO jk=1,jpkm1 827 v n(ji,j1-1,jk) = vn(ji,j1-1,jk) + zcor * vmask(ji,j1-1,jk)861 vv(ji,j1-1,jk,Kmm_a) = vv(ji,j1-1,jk,Kmm_a) + zcor * vmask(ji,j1-1,jk) 828 862 END DO 829 863 END DO … … 832 866 IF (northern_side) THEN 833 867 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) + zcor868 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) 869 vv_b(ji,j2+1,Kmm_a) = vv_b(ji,j2+1,Kmm_a) + zcor 836 870 DO jk=1,jpkm1 837 v n(ji,j2+1,jk) = vn(ji,j2+1,jk) + zcor * vmask(ji,j2+1,jk)871 vv(ji,j2+1,jk,Kmm_a) = vv(ji,j2+1,jk,Kmm_a) + zcor * vmask(ji,j2+1,jk) 838 872 END DO 839 873 END DO … … 862 896 DO jj=j1,j2 863 897 DO ji=i1,i2 864 tabres(ji,jj) = zrhoy * u n_b(ji,jj) * hu_n(ji,jj) * e2u(ji,jj)898 tabres(ji,jj) = zrhoy * uu_b(ji,jj,Kmm_a) * hu(ji,jj,Kmm_a) * e2u(ji,jj) 865 899 END DO 866 900 END DO … … 873 907 spgu(ji,jj) = 0._wp 874 908 DO jk=1,jpkm1 875 spgu(ji,jj) = spgu(ji,jj) + e3u _n(ji,jj,jk) * un(ji,jj,jk)909 spgu(ji,jj) = spgu(ji,jj) + e3u(ji,jj,jk,Kmm_a) * uu(ji,jj,jk,Kmm_a) 876 910 END DO 877 911 ! 878 zcorr = (tabres(ji,jj) - spgu(ji,jj)) * r1_hu _n(ji,jj)912 zcorr = (tabres(ji,jj) - spgu(ji,jj)) * r1_hu(ji,jj,Kmm_a) 879 913 DO jk=1,jpkm1 880 u n(ji,jj,jk) = un(ji,jj,jk) + zcorr * umask(ji,jj,jk)914 uu(ji,jj,jk,Kmm_a) = uu(ji,jj,jk,Kmm_a) + zcorr * umask(ji,jj,jk) 881 915 END DO 882 916 ! 883 917 ! Update barotropic velocities: 884 918 IF ( .NOT.ln_dynspg_ts .OR. (ln_dynspg_ts.AND.(.NOT.ln_bt_fw)) ) THEN 885 IF (.NOT.(lk_agrif_fstep.AND.( neuler==0))) THEN ! Add asselin part886 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)919 IF (.NOT.(lk_agrif_fstep.AND.(l_1st_euler))) THEN ! Add asselin part 920 zcorr = (tabres(ji,jj) - uu_b(ji,jj,Kmm_a) * hu(ji,jj,Krhs_a)) * r1_hu(ji,jj,Kbb_a) 921 uu_b(ji,jj,Kbb_a) = uu_b(ji,jj,Kbb_a) + rn_atfp * zcorr * umask(ji,jj,1) 888 922 END IF 889 923 ENDIF 890 u n_b(ji,jj) = tabres(ji,jj) * r1_hu_n(ji,jj) * umask(ji,jj,1)924 uu_b(ji,jj,Kmm_a) = tabres(ji,jj) * r1_hu(ji,jj,Kmm_a) * umask(ji,jj,1) 891 925 ! 892 926 ! Correct "before" velocities to hold correct bt component: 893 927 spgu(ji,jj) = 0.e0 894 928 DO jk=1,jpkm1 895 spgu(ji,jj) = spgu(ji,jj) + e3u _b(ji,jj,jk) * ub(ji,jj,jk)929 spgu(ji,jj) = spgu(ji,jj) + e3u(ji,jj,jk,Kbb_a) * uu(ji,jj,jk,Kbb_a) 896 930 END DO 897 931 ! 898 zcorr = u b_b(ji,jj) - spgu(ji,jj) * r1_hu_b(ji,jj)932 zcorr = uu_b(ji,jj,Kbb_a) - spgu(ji,jj) * r1_hu(ji,jj,Kbb_a) 899 933 DO jk=1,jpkm1 900 u b(ji,jj,jk) = ub(ji,jj,jk) + zcorr * umask(ji,jj,jk)934 uu(ji,jj,jk,Kbb_a) = uu(ji,jj,jk,Kbb_a) + zcorr * umask(ji,jj,jk) 901 935 END DO 902 936 ! … … 904 938 END DO 905 939 ! 906 IF (( neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN907 u b_b(i1:i2,j1:j2) = un_b(i1:i2,j1:j2)940 IF ((l_1st_euler).AND.(Agrif_Nb_Step()==0) ) THEN 941 uu_b(i1:i2,j1:j2,Kbb_a) = uu_b(i1:i2,j1:j2,Kmm_a) 908 942 ENDIF 909 943 ENDIF … … 928 962 DO jj=j1,j2 929 963 DO ji=i1,i2 930 tabres(ji,jj) = zrhox * v n_b(ji,jj) * hv_n(ji,jj) * e1v(ji,jj)964 tabres(ji,jj) = zrhox * vv_b(ji,jj,Kmm_a) * hv(ji,jj,Kmm_a) * e1v(ji,jj) 931 965 END DO 932 966 END DO … … 939 973 spgv(ji,jj) = 0.e0 940 974 DO jk=1,jpkm1 941 spgv(ji,jj) = spgv(ji,jj) + e3v _n(ji,jj,jk) * vn(ji,jj,jk)975 spgv(ji,jj) = spgv(ji,jj) + e3v(ji,jj,jk,Kmm_a) * vv(ji,jj,jk,Kmm_a) 942 976 END DO 943 977 ! 944 zcorr = (tabres(ji,jj) - spgv(ji,jj)) * r1_hv _n(ji,jj)978 zcorr = (tabres(ji,jj) - spgv(ji,jj)) * r1_hv(ji,jj,Kmm_a) 945 979 DO jk=1,jpkm1 946 v n(ji,jj,jk) = vn(ji,jj,jk) + zcorr * vmask(ji,jj,jk)980 vv(ji,jj,jk,Kmm_a) = vv(ji,jj,jk,Kmm_a) + zcorr * vmask(ji,jj,jk) 947 981 END DO 948 982 ! 949 983 ! Update barotropic velocities: 950 984 IF ( .NOT.ln_dynspg_ts .OR. (ln_dynspg_ts.AND.(.NOT.ln_bt_fw)) ) THEN 951 IF (.NOT.(lk_agrif_fstep.AND.( neuler==0))) THEN ! Add asselin part952 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)985 IF (.NOT.(lk_agrif_fstep.AND.(l_1st_euler))) THEN ! Add asselin part 986 zcorr = (tabres(ji,jj) - vv_b(ji,jj,Kmm_a) * hv(ji,jj,Krhs_a)) * r1_hv(ji,jj,Kbb_a) 987 vv_b(ji,jj,Kbb_a) = vv_b(ji,jj,Kbb_a) + rn_atfp * zcorr * vmask(ji,jj,1) 954 988 END IF 955 989 ENDIF 956 v n_b(ji,jj) = tabres(ji,jj) * r1_hv_n(ji,jj) * vmask(ji,jj,1)990 vv_b(ji,jj,Kmm_a) = tabres(ji,jj) * r1_hv(ji,jj,Kmm_a) * vmask(ji,jj,1) 957 991 ! 958 992 ! Correct "before" velocities to hold correct bt component: 959 993 spgv(ji,jj) = 0.e0 960 994 DO jk=1,jpkm1 961 spgv(ji,jj) = spgv(ji,jj) + e3v _b(ji,jj,jk) * vb(ji,jj,jk)995 spgv(ji,jj) = spgv(ji,jj) + e3v(ji,jj,jk,Kbb_a) * vv(ji,jj,jk,Kbb_a) 962 996 END DO 963 997 ! 964 zcorr = v b_b(ji,jj) - spgv(ji,jj) * r1_hv_b(ji,jj)998 zcorr = vv_b(ji,jj,Kbb_a) - spgv(ji,jj) * r1_hv(ji,jj,Kbb_a) 965 999 DO jk=1,jpkm1 966 v b(ji,jj,jk) = vb(ji,jj,jk) + zcorr * vmask(ji,jj,jk)1000 vv(ji,jj,jk,Kbb_a) = vv(ji,jj,jk,Kbb_a) + zcorr * vmask(ji,jj,jk) 967 1001 END DO 968 1002 ! … … 970 1004 END DO 971 1005 ! 972 IF (( neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN973 v b_b(i1:i2,j1:j2) = vn_b(i1:i2,j1:j2)1006 IF ((l_1st_euler).AND.(Agrif_Nb_Step()==0) ) THEN 1007 vv_b(i1:i2,j1:j2,Kbb_a) = vv_b(i1:i2,j1:j2,Kmm_a) 974 1008 ENDIF 975 1009 ! … … 993 1027 DO jj=j1,j2 994 1028 DO ji=i1,i2 995 tabres(ji,jj) = ssh n(ji,jj)1029 tabres(ji,jj) = ssh(ji,jj,Kmm_a) 996 1030 END DO 997 1031 END DO 998 1032 ELSE 999 IF (.NOT.(lk_agrif_fstep.AND.( neuler==0))) THEN1033 IF (.NOT.(lk_agrif_fstep.AND.(l_1st_euler))) THEN 1000 1034 DO jj=j1,j2 1001 1035 DO ji=i1,i2 1002 ssh b(ji,jj) = sshb(ji,jj) &1003 & + atfp * ( tabres(ji,jj) - sshn(ji,jj) ) * tmask(ji,jj,1)1036 ssh(ji,jj,Kbb_a) = ssh(ji,jj,Kbb_a) & 1037 & + rn_atfp * ( tabres(ji,jj) - ssh(ji,jj,Kmm_a) ) * tmask(ji,jj,1) 1004 1038 END DO 1005 1039 END DO … … 1008 1042 DO jj=j1,j2 1009 1043 DO ji=i1,i2 1010 ssh n(ji,jj) = tabres(ji,jj) * tmask(ji,jj,1)1011 END DO 1012 END DO 1013 ! 1014 IF (( neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN1015 ssh b(i1:i2,j1:j2) = sshn(i1:i2,j1:j2)1044 ssh(ji,jj,Kmm_a) = tabres(ji,jj) * tmask(ji,jj,1) 1045 END DO 1046 END DO 1047 ! 1048 IF ((l_1st_euler).AND.(Agrif_Nb_Step()==0) ) THEN 1049 ssh(i1:i2,j1:j2,Kbb_a) = ssh(i1:i2,j1:j2,Kmm_a) 1016 1050 ENDIF 1017 1051 ! … … 1093 1127 IF (western_side) THEN 1094 1128 DO jj=j1,j2 1095 zcor = r dt * 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))) sshb(i1 ,jj) = sshb(i1 ,jj) +atfp * zcor1129 zcor = rn_Dt * r1_e1e2t(i1 ,jj) * e2u(i1,jj) * (ub2_b(i1,jj)-tabres(i1,jj)) 1130 ssh(i1 ,jj,Kmm_a) = ssh(i1 ,jj,Kmm_a) + zcor 1131 IF (.NOT.(lk_agrif_fstep.AND.(l_1st_euler))) ssh(i1 ,jj,Kbb_a) = ssh(i1 ,jj,Kbb_a) + rn_atfp * zcor 1098 1132 END DO 1099 1133 ENDIF 1100 1134 IF (eastern_side) THEN 1101 1135 DO jj=j1,j2 1102 zcor = - r dt * 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))) sshb(i2+1,jj) = sshb(i2+1,jj) +atfp * zcor1136 zcor = - rn_Dt * r1_e1e2t(i2+1,jj) * e2u(i2,jj) * (ub2_b(i2,jj)-tabres(i2,jj)) 1137 ssh(i2+1,jj,Kmm_a) = ssh(i2+1,jj,Kmm_a) + zcor 1138 IF (.NOT.(lk_agrif_fstep.AND.(l_1st_euler))) ssh(i2+1,jj,Kbb_a) = ssh(i2+1,jj,Kbb_a) + rn_atfp * zcor 1105 1139 END DO 1106 1140 ENDIF … … 1181 1215 IF (southern_side) THEN 1182 1216 DO ji=i1,i2 1183 zcor = r dt * 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))) sshb(ji,j1 ) = sshb(ji,j1) +atfp * zcor1217 zcor = rn_Dt * r1_e1e2t(ji,j1 ) * e1v(ji,j1 ) * (vb2_b(ji,j1)-tabres(ji,j1)) 1218 ssh(ji,j1 ,Kmm_a) = ssh(ji,j1 ,Kmm_a) + zcor 1219 IF (.NOT.(lk_agrif_fstep.AND.(l_1st_euler))) ssh(ji,j1 ,Kbb_a) = ssh(ji,j1,Kbb_a) + rn_atfp * zcor 1186 1220 END DO 1187 1221 ENDIF 1188 1222 IF (northern_side) THEN 1189 1223 DO ji=i1,i2 1190 zcor = - r dt * 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))) sshb(ji,j2+1) = sshb(ji,j2+1) +atfp * zcor1224 zcor = - rn_Dt * r1_e1e2t(ji,j2+1) * e1v(ji,j2 ) * (vb2_b(ji,j2)-tabres(ji,j2)) 1225 ssh(ji,j2+1,Kmm_a) = ssh(ji,j2+1,Kmm_a) + zcor 1226 IF (.NOT.(lk_agrif_fstep.AND.(l_1st_euler))) ssh(ji,j2+1,Kbb_a) = ssh(ji,j2+1,Kbb_a) + rn_atfp * zcor 1193 1227 END DO 1194 1228 ENDIF … … 1319 1353 DO jj=j1,j2 1320 1354 DO ji=i1,i2 1321 ptab(ji,jj,jk) = e3t_0(ji,jj,jk) * (1._wp + ssh n(ji,jj) &1355 ptab(ji,jj,jk) = e3t_0(ji,jj,jk) * (1._wp + ssh(ji,jj,Kmm_a) & 1322 1356 & *ssmask(ji,jj)/(ht_0(ji,jj)-1._wp + ssmask(ji,jj))) 1323 1357 END DO … … 1330 1364 ! Save "old" scale factor (prior update) for subsequent asselin correction 1331 1365 ! 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)1336 1337 IF (.NOT.(lk_agrif_fstep.AND.( neuler==0) )) THEN1366 e3t(i1:i2,j1:j2,1:jpkm1,Krhs_a) = e3t(i1:i2,j1:j2,1:jpkm1,Kmm_a) 1367 1368 ! One should also save e3t(:,:,:,Kbb_a), but lacking of workspace... 1369 ! hdiv(i1:i2,j1:j2,1:jpkm1) = e3t(i1:i2,j1:j2,1:jpkm1,Kbb_a) 1370 1371 IF (.NOT.(lk_agrif_fstep.AND.(l_1st_euler) )) THEN 1338 1372 DO jk = 1, jpkm1 1339 1373 DO jj=j1,j2 1340 1374 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) )1375 e3t(ji,jj,jk,Kbb_a) = e3t(ji,jj,jk,Kbb_a) & 1376 & + rn_atfp * ( ptab(ji,jj,jk) - e3t(ji,jj,jk,Kmm_a) ) 1343 1377 END DO 1344 1378 END DO 1345 1379 END DO 1346 1380 ! 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)1381 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) 1382 gdepw(i1:i2,j1:j2,1,Kbb_a) = 0.0_wp 1383 gdept(i1:i2,j1:j2,1,Kbb_a) = 0.5_wp * e3w(i1:i2,j1:j2,1,Kbb_a) 1350 1384 ! 1351 1385 DO jk = 2, jpk … … 1353 1387 DO ji = i1,i2 1354 1388 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) ) &1389 e3w(ji,jj,jk,Kbb_a) = e3w_0(ji,jj,jk) + ( 1.0_wp - 0.5_wp * tmask(ji,jj,jk) ) * & 1390 & ( e3t(ji,jj,jk-1,Kbb_a) - e3t_0(ji,jj,jk-1) ) & 1357 1391 & + 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))1392 & ( e3t(ji,jj,jk ,Kbb_a) - e3t_0(ji,jj,jk ) ) 1393 gdepw(ji,jj,jk,Kbb_a) = gdepw(ji,jj,jk-1,Kbb_a) + e3t(ji,jj,jk-1,Kbb_a) 1394 gdept(ji,jj,jk,Kbb_a) = zcoef * ( gdepw(ji,jj,jk ,Kbb_a) + 0.5 * e3w(ji,jj,jk,Kbb_a)) & 1395 & + (1-zcoef) * ( gdept(ji,jj,jk-1,Kbb_a) + e3w(ji,jj,jk,Kbb_a)) 1362 1396 END DO 1363 1397 END DO … … 1370 1404 ! 1371 1405 ! Update vertical scale factor at T-points: 1372 e3t _n(i1:i2,j1:j2,1:jpkm1) = ptab(i1:i2,j1:j2,1:jpkm1)1406 e3t(i1:i2,j1:j2,1:jpkm1,Kmm_a) = ptab(i1:i2,j1:j2,1:jpkm1) 1373 1407 ! 1374 1408 ! Update total depth: 1375 ht _n(i1:i2,j1:j2) = 0._wp1409 ht(i1:i2,j1:j2) = 0._wp 1376 1410 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)1411 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 1412 END DO 1379 1413 ! 1380 1414 ! 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 ssh1415 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) 1416 gdept(i1:i2,j1:j2,1,Kmm_a) = 0.5_wp * e3w(i1:i2,j1:j2,1,Kmm_a) 1417 gdepw(i1:i2,j1:j2,1,Kmm_a) = 0.0_wp 1418 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 1419 ! 1386 1420 DO jk = 2, jpk … … 1388 1422 DO ji = i1,i2 1389 1423 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 ssh1396 END DO 1397 END DO 1398 END DO 1399 ! 1400 IF (( neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN1401 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)1424 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) ) & 1425 & + 0.5_wp * tmask(ji,jj,jk) * ( e3t(ji,jj,jk ,Kmm_a) - e3t_0(ji,jj,jk ) ) 1426 gdepw(ji,jj,jk,Kmm_a) = gdepw(ji,jj,jk-1,Kmm_a) + e3t(ji,jj,jk-1,Kmm_a) 1427 gdept(ji,jj,jk,Kmm_a) = zcoef * ( gdepw(ji,jj,jk ,Kmm_a) + 0.5 * e3w(ji,jj,jk,Kmm_a)) & 1428 & + (1-zcoef) * ( gdept(ji,jj,jk-1,Kmm_a) + e3w(ji,jj,jk,Kmm_a)) 1429 gde3w(ji,jj,jk) = gdept(ji,jj,jk,Kmm_a) - (ht(ji,jj)-ht_0(ji,jj)) ! Last term in the rhs is ssh 1430 END DO 1431 END DO 1432 END DO 1433 ! 1434 IF ((l_1st_euler).AND.(Agrif_Nb_Step()==0) ) THEN 1435 e3t (i1:i2,j1:j2,1:jpk,Kbb_a) = e3t (i1:i2,j1:j2,1:jpk,Kmm_a) 1436 e3w (i1:i2,j1:j2,1:jpk,Kbb_a) = e3w (i1:i2,j1:j2,1:jpk,Kmm_a) 1437 gdepw(i1:i2,j1:j2,1:jpk,Kbb_a) = gdepw(i1:i2,j1:j2,1:jpk,Kmm_a) 1438 gdept(i1:i2,j1:j2,1:jpk,Kbb_a) = gdept(i1:i2,j1:j2,1:jpk,Kmm_a) 1405 1439 ENDIF 1406 1440 !
Note: See TracChangeset
for help on using the changeset viewer.
