- Timestamp:
- 2017-11-17T17:19:55+01:00 (6 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2017/dev_r8624_AGRIF3_VVL/NEMOGCM/NEMO/NST_SRC/agrif_opa_update.F90
r7646 r8741 12 12 USE wrk_nemo 13 13 USE zdf_oce ! vertical physics: ocean variables 14 USE domvvl ! Need interpolation routines 14 15 15 16 IMPLICIT NONE 16 17 PRIVATE 17 18 18 PUBLIC Agrif_Update_Tra, Agrif_Update_Dyn,Update_Scales 19 PUBLIC Agrif_Update_Tra, Agrif_Update_Dyn, Update_Scales, Agrif_Update_vvl 20 19 21 # if defined key_zdftke 20 22 PUBLIC Agrif_Update_Tke … … 27 29 CONTAINS 28 30 29 RECURSIVESUBROUTINE Agrif_Update_Tra( )31 SUBROUTINE Agrif_Update_Tra( ) 30 32 !!--------------------------------------------- 31 33 !! *** ROUTINE Agrif_Update_Tra *** … … 56 58 Agrif_UseSpecialValueInUpdate = .FALSE. 57 59 ! 58 IF ( lk_agrif_doupd ) THEN ! Initialisation: do recursive update:59 CALL Agrif_ChildGrid_To_ParentGrid()60 CALL Agrif_Update_Tra()61 CALL Agrif_ParentGrid_To_ChildGrid()62 ENDIF63 !64 60 #endif 65 61 ! 66 62 END SUBROUTINE Agrif_Update_Tra 67 63 68 69 RECURSIVE SUBROUTINE Agrif_Update_Dyn( ) 64 SUBROUTINE Agrif_Update_Dyn( ) 70 65 !!--------------------------------------------- 71 66 !! *** ROUTINE Agrif_Update_Dyn *** … … 140 135 #endif 141 136 ! 142 ! Do recursive update:143 IF ( lk_agrif_doupd ) THEN ! Initialisation: do recursive update:144 CALL Agrif_ChildGrid_To_ParentGrid()145 CALL Agrif_Update_Dyn()146 CALL Agrif_ParentGrid_To_ChildGrid()147 ENDIF148 !149 137 END SUBROUTINE Agrif_Update_Dyn 150 138 151 139 # if defined key_zdftke 152 140 153 SUBROUTINE Agrif_Update_Tke( kt)141 SUBROUTINE Agrif_Update_Tke( ) 154 142 !!--------------------------------------------- 155 143 !! *** ROUTINE Agrif_Update_Tke *** 156 144 !!--------------------------------------------- 157 145 !! 158 INTEGER, INTENT(in) :: kt 146 ! 147 IF (Agrif_Root()) RETURN 159 148 ! 160 149 IF( ( Agrif_NbStepint() .NE. 0 ) .AND. (kt /= 0) ) RETURN … … 176 165 # endif /* key_zdftke */ 177 166 178 SUBROUTINE updateTS( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before ) 167 SUBROUTINE Agrif_Update_vvl( ) 168 !!--------------------------------------------- 169 !! *** ROUTINE Agrif_Update_vvl *** 170 !!--------------------------------------------- 171 ! 172 IF (Agrif_Root()) RETURN 173 ! 174 #if defined TWO_WAY 175 ! 176 IF (lwp.AND.lk_agrif_debug) Write(*,*) 'Update e3 from grid Number',Agrif_Fixed(), 'Step', Agrif_Nb_Step() 177 ! 178 Agrif_UseSpecialValueInUpdate = .TRUE. 179 Agrif_SpecialValueFineGrid = 0. 180 ! 181 # if ! defined DECAL_FEEDBACK 182 CALL Agrif_Update_Variable(e3t_id, procname=updatee3t) 183 # else 184 CALL Agrif_Update_Variable(e3t_id, locupdate=(/1,0/), procname=updatee3t) 185 # endif 186 ! 187 Agrif_UseSpecialValueInUpdate = .FALSE. 188 ! 189 CALL Agrif_ChildGrid_To_ParentGrid() 190 CALL dom_vvl_update_UVF 191 CALL Agrif_ParentGrid_To_ChildGrid() 192 ! 193 #endif 194 ! 195 END SUBROUTINE Agrif_Update_vvl 196 197 SUBROUTINE dom_vvl_update_UVF 198 !!--------------------------------------------- 199 !! *** ROUTINE dom_vvl_update_UVF *** 200 !!--------------------------------------------- 201 !! 202 INTEGER :: jk 203 REAL(wp):: zcoef 204 !!--------------------------------------------- 205 206 IF (lwp.AND.lk_agrif_debug) Write(*,*) 'Finalize e3 on grid Number', & 207 & Agrif_Fixed(), 'Step', Agrif_Nb_Step() 208 209 ! Save "old" scale factor (prior update) for subsequent asselin correction 210 ! of prognostic variables (needed to update initial state only) 211 ! ------------------------------------------------------------- 212 ! 213 e3u_a(:,:,:) = e3u_n(:,:,:) 214 e3v_a(:,:,:) = e3v_n(:,:,:) 215 ! ua(:,:,:) = e3u_b(:,:,:) 216 ! va(:,:,:) = e3v_b(:,:,:) 217 hu_a(:,:) = hu_n(:,:) 218 hv_a(:,:) = hv_n(:,:) 219 220 ! 1) NOW fields 221 !-------------- 222 223 ! Vertical scale factor interpolations 224 ! ------------------------------------ 225 CALL dom_vvl_interpol( e3t_n(:,:,:), e3u_n(:,:,:) , 'U' ) 226 CALL dom_vvl_interpol( e3t_n(:,:,:), e3v_n(:,:,:) , 'V' ) 227 CALL dom_vvl_interpol( e3u_n(:,:,:), e3f_n(:,:,:) , 'F' ) 228 229 CALL dom_vvl_interpol( e3u_n(:,:,:), e3uw_n(:,:,:), 'UW' ) 230 CALL dom_vvl_interpol( e3v_n(:,:,:), e3vw_n(:,:,:), 'VW' ) 231 232 ! Update total depths: 233 ! -------------------- 234 hu_n(:,:) = 0._wp ! Ocean depth at U-points 235 hv_n(:,:) = 0._wp ! Ocean depth at V-points 236 DO jk = 1, jpkm1 237 hu_n(:,:) = hu_n(:,:) + e3u_n(:,:,jk) * umask(:,:,jk) 238 hv_n(:,:) = hv_n(:,:) + e3v_n(:,:,jk) * vmask(:,:,jk) 239 END DO 240 ! ! Inverse of the local depth 241 r1_hu_n(:,:) = ssumask(:,:) / ( hu_n(:,:) + 1._wp - ssumask(:,:) ) 242 r1_hv_n(:,:) = ssvmask(:,:) / ( hv_n(:,:) + 1._wp - ssvmask(:,:) ) 243 244 245 ! 2) BEFORE fields: 246 !------------------ 247 IF ( (.NOT.(lk_agrif_fstep.AND.(neuler==0)).AND.(ln_dynspg_exp)) & 248 & .OR.(.NOT.(lk_agrif_fstep.AND.(neuler==0)).AND.(ln_dynspg_ts & 249 & .AND.(.NOT.ln_bt_fw)))) THEN 250 ! 251 ! Vertical scale factor interpolations 252 ! ------------------------------------ 253 CALL dom_vvl_interpol( e3t_b(:,:,:), e3u_b(:,:,:), 'U' ) 254 CALL dom_vvl_interpol( e3t_b(:,:,:), e3v_b(:,:,:), 'V' ) 255 256 CALL dom_vvl_interpol( e3u_b(:,:,:), e3uw_b(:,:,:), 'UW' ) 257 CALL dom_vvl_interpol( e3v_b(:,:,:), e3vw_b(:,:,:), 'VW' ) 258 259 ! Update total depths: 260 ! -------------------- 261 hu_b(:,:) = 0._wp ! Ocean depth at U-points 262 hv_b(:,:) = 0._wp ! Ocean depth at V-points 263 DO jk = 1, jpkm1 264 hu_b(:,:) = hu_b(:,:) + e3u_b(:,:,jk) * umask(:,:,jk) 265 hv_b(:,:) = hv_b(:,:) + e3v_b(:,:,jk) * vmask(:,:,jk) 266 END DO 267 ! ! Inverse of the local depth 268 r1_hu_b(:,:) = ssumask(:,:) / ( hu_b(:,:) + 1._wp - ssumask(:,:) ) 269 r1_hv_b(:,:) = ssvmask(:,:) / ( hv_b(:,:) + 1._wp - ssvmask(:,:) ) 270 ENDIF 271 ! 272 END SUBROUTINE dom_vvl_update_UVF 273 274 SUBROUTINE updateTS( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before, nb, ndir ) 179 275 !!--------------------------------------------- 180 276 !! *** ROUTINE updateT *** … … 183 279 REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres 184 280 LOGICAL, INTENT(in) :: before 281 INTEGER, INTENT(in) :: nb, ndir 185 282 !! 283 LOGICAL :: western_side, eastern_side, southern_side, northern_side 186 284 INTEGER :: ji,jj,jk,jn 285 REAL(wp) :: ztb, ztnu, ztno 187 286 !!--------------------------------------------- 188 287 ! … … 192 291 DO jj=j1,j2 193 292 DO ji=i1,i2 194 tabres(ji,jj,jk,jn) = tsn(ji,jj,jk,jn) 293 !> jc tmp 294 tabres(ji,jj,jk,jn) = tsn(ji,jj,jk,jn) * e3t_n(ji,jj,jk) / e3t_0(ji,jj,jk) 295 ! tabres(ji,jj,jk,jn) = tsn(ji,jj,jk,jn) * e3t_n(ji,jj,jk) 296 !< jc tmp 195 297 END DO 196 298 END DO … … 198 300 END DO 199 301 ELSE 302 !> jc tmp 303 DO jn = n1,n2 304 tabres(i1:i2,j1:j2,k1:k2,jn) = tabres(i1:i2,j1:j2,k1:k2,jn) * e3t_0(i1:i2,j1:j2,k1:k2) & 305 & * tmask(i1:i2,j1:j2,k1:k2) 306 ENDDO 307 !< jc tmp 200 308 IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN 201 309 ! Add asselin part … … 205 313 DO ji=i1,i2 206 314 IF( tabres(ji,jj,jk,jn) .NE. 0. ) THEN 207 tsb(ji,jj,jk,jn) = tsb(ji,jj,jk,jn) & 208 & + atfp * ( tabres(ji,jj,jk,jn) & 209 & - tsn(ji,jj,jk,jn) ) * tmask(ji,jj,jk) 315 ztb = tsb(ji,jj,jk,jn) * e3t_b(ji,jj,jk) ! fse3t_b prior update should be used 316 ztnu = tabres(ji,jj,jk,jn) 317 ztno = tsn(ji,jj,jk,jn) * e3t_a(ji,jj,jk) 318 tsb(ji,jj,jk,jn) = ( ztb + atfp * ( ztnu - ztno) ) & 319 & * tmask(ji,jj,jk) / e3t_b(ji,jj,jk) 210 320 ENDIF 211 321 ENDDO … … 219 329 DO ji=i1,i2 220 330 IF( tabres(ji,jj,jk,jn) .NE. 0. ) THEN 221 tsn(ji,jj,jk,jn) = tabres(ji,jj,jk,jn) * tmask(ji,jj,jk)331 tsn(ji,jj,jk,jn) = tabres(ji,jj,jk,jn) / e3t_n(ji,jj,jk) 222 332 END IF 223 333 END DO … … 225 335 END DO 226 336 END DO 337 ! 338 IF ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN 339 tsb(i1:i2,j1:j2,k1:k2,n1:n2) = tsn(i1:i2,j1:j2,k1:k2,n1:n2) 340 ENDIF 341 ! 342 ! 343 # if defined DECAL_FEEDBACK 344 IF (.NOT.ln_linssh) THEN 345 western_side = (nb == 1).AND.(ndir == 1) 346 eastern_side = (nb == 1).AND.(ndir == 2) 347 southern_side = (nb == 2).AND.(ndir == 1) 348 northern_side = (nb == 2).AND.(ndir == 2) 349 ! 350 ! Asselin correction 351 IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN 352 IF (southern_side) THEN 353 DO jn = n1,n2 354 DO jk=k1,k2 355 DO ji=i1,i2 356 ztb = tsb(ji,j1-1,jk,jn) * e3t_b(ji,j1-1,jk) ! fse3t_b prior update should be used 357 ztnu = tsn(ji,j1-1,jk,jn) * e3t_n(ji,j1-1,jk) 358 ztno = tsn(ji,j1-1,jk,jn) * e3t_a(ji,j1-1,jk) 359 tsb(ji,j1-1,jk,jn) = ( ztb + atfp * ( ztnu - ztno) ) & 360 & * tmask(ji,j1-1,jk) / e3t_b(ji,j1-1,jk) 361 END DO 362 ENDDO 363 ENDDO 364 ENDIF 365 IF (northern_side) THEN 366 DO jn = n1,n2 367 DO jk=k1,k2 368 DO ji=i1,i2 369 ztb = tsb(ji,j2+1,jk,jn) * e3t_b(ji,j2+1,jk) ! fse3t_b prior update should be used 370 ztnu = tsn(ji,j2+1,jk,jn) * e3t_n(ji,j2+1,jk) 371 ztno = tsn(ji,j2+1,jk,jn) * e3t_a(ji,j2+1,jk) 372 tsb(ji,j2+1,jk,jn) = ( ztb + atfp * ( ztnu - ztno) ) & 373 & * tmask(ji,j2+1,jk) / e3t_b(ji,j2+1,jk) 374 END DO 375 ENDDO 376 ENDDO 377 ENDIF 378 IF (western_side) THEN 379 DO jn = n1,n2 380 DO jk=k1,k2 381 DO jj=j1,j2 382 ztb = tsb(i1-1,jj,jk,jn) * e3t_b(i1-1,jj,jk) ! fse3t_b prior update should be used 383 ztnu = tsn(i1-1,jj,jk,jn) * e3t_n(i1-1,jj,jk) 384 ztno = tsn(i1-1,jj,jk,jn) * e3t_a(i1-1,jj,jk) 385 tsb(i1-1,jj,jk,jn) = ( ztb + atfp * ( ztnu - ztno) ) & 386 & * tmask(i1-1,jj,jk) / e3t_b(i1-1,jj,jk) 387 END DO 388 ENDDO 389 ENDDO 390 ENDIF 391 IF (eastern_side) THEN 392 DO jn = n1,n2 393 DO jk=k1,k2 394 DO jj=j1,j2 395 ztb = tsb(i2+1,jj,jk,jn) * e3t_b(i2+1,jj,jk) ! fse3t_b prior update should be used 396 ztnu = tsn(i2+1,jj,jk,jn) * e3t_n(i2+1,jj,jk) 397 ztno = tsn(i2+1,jj,jk,jn) * e3t_a(i2+1,jj,jk) 398 tsb(i2+1,jj,jk,jn) = ( ztb + atfp * ( ztnu - ztno) ) & 399 & * tmask(i2+1,jj,jk) / e3t_b(i2+1,jj,jk) 400 END DO 401 ENDDO 402 ENDDO 403 ENDIF 404 ENDIF ! Asselin correction 405 406 IF (southern_side) THEN 407 DO jn = n1,n2 408 DO jk=k1,k2 409 DO ji=i1,i2 410 tsn(ji,j1-1,jk,jn) = tsn(ji,j1-1,jk,jn) * e3t_a(ji,j1-1,jk) / e3t_n(ji,j1-1,jk) 411 END DO 412 ENDDO 413 ENDDO 414 ENDIF 415 IF (northern_side) THEN 416 DO jn = n1,n2 417 DO jk=k1,k2 418 DO ji=i1,i2 419 tsn(ji,j2+1,jk,jn) = tsn(ji,j2+1,jk,jn) * e3t_a(ji,j2+1,jk) / e3t_n(ji,j2+1,jk) 420 END DO 421 ENDDO 422 ENDDO 423 ENDIF 424 IF (western_side) THEN 425 DO jn = n1,n2 426 DO jk=k1,k2 427 DO jj=j1,j2 428 tsn(i1-1,jj,jk,jn) = tsn(i1-1,jj,jk,jn) * e3t_a(i1-1,jj,jk) / e3t_n(i1-1,jj,jk) 429 END DO 430 ENDDO 431 ENDDO 432 ENDIF 433 IF (eastern_side) THEN 434 DO jn = n1,n2 435 DO jk=k1,k2 436 DO jj=j1,j2 437 tsn(i2+1,jj,jk,jn) = tsn(i2+1,jj,jk,jn) * e3t_a(i2+1,jj,jk) / e3t_n(i2+1,jj,jk) 438 END DO 439 ENDDO 440 ENDDO 441 ENDIF 442 ENDIF 443 #endif 227 444 ENDIF 228 445 ! … … 238 455 LOGICAL , INTENT(in ) :: before 239 456 ! 240 INTEGER :: 241 REAL(wp) :: zrhoy457 INTEGER :: ji, jj, jk 458 REAL(wp) :: zrhoy, zub, zunu, zuno 242 459 !!--------------------------------------------- 243 460 ! … … 251 468 DO jj=j1,j2 252 469 DO ji=i1,i2 253 tabres(ji,jj,jk) = tabres(ji,jj,jk) * r1_e2u(ji,jj) / e3u_n(ji,jj,jk)470 tabres(ji,jj,jk) = tabres(ji,jj,jk) * r1_e2u(ji,jj) 254 471 ! 255 472 IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part 256 ub(ji,jj,jk) = ub(ji,jj,jk) & 257 & + atfp * ( tabres(ji,jj,jk) - un(ji,jj,jk) ) * umask(ji,jj,jk) 473 zub = ub(ji,jj,jk) * e3u_b(ji,jj,jk) 474 zuno = un(ji,jj,jk) * e3u_a(ji,jj,jk) 475 zunu = tabres(ji,jj,jk) 476 ub(ji,jj,jk) = ( zub + atfp * ( zunu - zuno) ) & 477 & * umask(ji,jj,jk) / e3u_b(ji,jj,jk) 258 478 ENDIF 259 479 ! 260 un(ji,jj,jk) = tabres(ji,jj,jk) * umask(ji,jj,jk) 261 END DO 262 END DO 263 END DO 480 un(ji,jj,jk) = tabres(ji,jj,jk) * umask(ji,jj,jk) / e3u_n(ji,jj,jk) 481 END DO 482 END DO 483 END DO 484 ! 485 IF ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN 486 ub(i1:i2,j1:j2,k1:k2) = un(i1:i2,j1:j2,k1:k2) 487 ENDIF 488 ! 264 489 ENDIF 265 490 ! … … 267 492 268 493 269 SUBROUTINE updatev( tabres, i1, i2, j1, j2, k1, k2, before 494 SUBROUTINE updatev( tabres, i1, i2, j1, j2, k1, k2, before) 270 495 !!--------------------------------------------- 271 496 !! *** ROUTINE updatev *** 272 497 !!--------------------------------------------- 273 INTEGER :: i1,i2,j1,j2,k1,k2274 INTEGER :: ji,jj,jk275 REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2) :: tabres276 LOGICAL :: before277 !!278 REAL(wp) :: zrhox 498 INTEGER , INTENT(in ) :: i1, i2, j1, j2, k1, k2 499 REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: tabres 500 LOGICAL , INTENT(in ) :: before 501 ! 502 INTEGER :: ji, jj, jk 503 REAL(wp) :: zrhox, zvb, zvnu, zvno 279 504 !!--------------------------------------------- 280 505 ! … … 292 517 DO jj=j1,j2 293 518 DO ji=i1,i2 294 tabres(ji,jj,jk) = tabres(ji,jj,jk) * r1_e1v(ji,jj) / e3v_n(ji,jj,jk)519 tabres(ji,jj,jk) = tabres(ji,jj,jk) * r1_e1v(ji,jj) 295 520 ! 296 521 IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part 297 vb(ji,jj,jk) = vb(ji,jj,jk) & 298 & + atfp * ( tabres(ji,jj,jk) - vn(ji,jj,jk) ) * vmask(ji,jj,jk) 522 zvb = vb(ji,jj,jk) * e3v_b(ji,jj,jk) 523 zvno = vn(ji,jj,jk) * e3v_a(ji,jj,jk) 524 zvnu = tabres(ji,jj,jk) 525 vb(ji,jj,jk) = ( zvb + atfp * ( zvnu - zvno) ) & 526 & * vmask(ji,jj,jk) / e3v_b(ji,jj,jk) 299 527 ENDIF 300 528 ! 301 vn(ji,jj,jk) = tabres(ji,jj,jk) * vmask(ji,jj,jk) 302 END DO 303 END DO 304 END DO 529 vn(ji,jj,jk) = tabres(ji,jj,jk) * vmask(ji,jj,jk) / e3v_n(ji,jj,jk) 530 END DO 531 END DO 532 END DO 533 ! 534 IF ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN 535 vb(i1:i2,j1:j2,k1:k2) = vn(i1:i2,j1:j2,k1:k2) 536 ENDIF 537 ! 305 538 ENDIF 306 539 ! … … 316 549 LOGICAL, INTENT(in) :: before 317 550 !! 318 INTEGER :: ji, jj, jk551 INTEGER :: ji, jj, jk 319 552 REAL(wp) :: zrhoy 320 553 REAL(wp) :: zcorr … … 331 564 DO jj=j1,j2 332 565 DO ji=i1,i2 333 tabres(ji,jj) = tabres(ji,jj) * r1_ hu_n(ji,jj) * r1_e2u(ji,jj)566 tabres(ji,jj) = tabres(ji,jj) * r1_e2u(ji,jj) 334 567 ! 335 568 ! Update "now" 3d velocities: … … 338 571 spgu(ji,jj) = spgu(ji,jj) + e3u_n(ji,jj,jk) * un(ji,jj,jk) 339 572 END DO 340 spgu(ji,jj) = spgu(ji,jj) * r1_hu_n(ji,jj)341 573 ! 342 zcorr = tabres(ji,jj) - spgu(ji,jj)574 zcorr = (tabres(ji,jj) - spgu(ji,jj)) * r1_hu_n(ji,jj) 343 575 DO jk=1,jpkm1 344 576 un(ji,jj,jk) = un(ji,jj,jk) + zcorr * umask(ji,jj,jk) … … 348 580 IF ( .NOT.ln_dynspg_ts .OR. (ln_dynspg_ts.AND.(.NOT.ln_bt_fw)) ) THEN 349 581 IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part 350 zcorr = tabres(ji,jj) - un_b(ji,jj)582 zcorr = (tabres(ji,jj) - un_b(ji,jj) * hu_a(ji,jj)) * r1_hu_b(ji,jj) 351 583 ub_b(ji,jj) = ub_b(ji,jj) + atfp * zcorr * umask(ji,jj,1) 352 584 END IF 353 ENDIF 354 un_b(ji,jj) = tabres(ji,jj) * umask(ji,jj,1)585 ENDIF 586 un_b(ji,jj) = tabres(ji,jj) * r1_hu_n(ji,jj) * umask(ji,jj,1) 355 587 ! 356 588 ! Correct "before" velocities to hold correct bt component: … … 359 591 spgu(ji,jj) = spgu(ji,jj) + e3u_b(ji,jj,jk) * ub(ji,jj,jk) 360 592 END DO 361 spgu(ji,jj) = spgu(ji,jj) * r1_hu_b(ji,jj)362 593 ! 363 zcorr = ub_b(ji,jj) - spgu(ji,jj) 594 zcorr = ub_b(ji,jj) - spgu(ji,jj) * r1_hu_b(ji,jj) 364 595 DO jk=1,jpkm1 365 596 ub(ji,jj,jk) = ub(ji,jj,jk) + zcorr * umask(ji,jj,jk) … … 368 599 END DO 369 600 END DO 601 ! 602 IF ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN 603 ub_b(i1:i2,j1:j2) = un_b(i1:i2,j1:j2) 604 ENDIF 370 605 ENDIF 371 606 ! … … 381 616 LOGICAL, INTENT(in) :: before 382 617 !! 383 INTEGER :: ji, jj, jk618 INTEGER :: ji, jj, jk 384 619 REAL(wp) :: zrhox 385 620 REAL(wp) :: zcorr … … 396 631 DO jj=j1,j2 397 632 DO ji=i1,i2 398 tabres(ji,jj) = tabres(ji,jj) * r1_ hv_n(ji,jj) * r1_e1v(ji,jj)633 tabres(ji,jj) = tabres(ji,jj) * r1_e1v(ji,jj) 399 634 ! 400 635 ! Update "now" 3d velocities: … … 403 638 spgv(ji,jj) = spgv(ji,jj) + e3v_n(ji,jj,jk) * vn(ji,jj,jk) 404 639 END DO 405 spgv(ji,jj) = spgv(ji,jj) * r1_hv_n(ji,jj)406 640 ! 407 zcorr = tabres(ji,jj) - spgv(ji,jj)641 zcorr = (tabres(ji,jj) - spgv(ji,jj)) * r1_hv_n(ji,jj) 408 642 DO jk=1,jpkm1 409 643 vn(ji,jj,jk) = vn(ji,jj,jk) + zcorr * vmask(ji,jj,jk) … … 413 647 IF ( .NOT.ln_dynspg_ts .OR. (ln_dynspg_ts.AND.(.NOT.ln_bt_fw)) ) THEN 414 648 IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part 415 zcorr = tabres(ji,jj) - vn_b(ji,jj)649 zcorr = (tabres(ji,jj) - vn_b(ji,jj) * hv_a(ji,jj)) * r1_hv_b(ji,jj) 416 650 vb_b(ji,jj) = vb_b(ji,jj) + atfp * zcorr * vmask(ji,jj,1) 417 651 END IF 418 652 ENDIF 419 vn_b(ji,jj) = tabres(ji,jj) * vmask(ji,jj,1)653 vn_b(ji,jj) = tabres(ji,jj) * r1_hv_n(ji,jj) * vmask(ji,jj,1) 420 654 ! 421 655 ! Correct "before" velocities to hold correct bt component: … … 424 658 spgv(ji,jj) = spgv(ji,jj) + e3v_b(ji,jj,jk) * vb(ji,jj,jk) 425 659 END DO 426 spgv(ji,jj) = spgv(ji,jj) * r1_hv_b(ji,jj)427 660 ! 428 zcorr = vb_b(ji,jj) - spgv(ji,jj) 661 zcorr = vb_b(ji,jj) - spgv(ji,jj) * r1_hv_b(ji,jj) 429 662 DO jk=1,jpkm1 430 663 vb(ji,jj,jk) = vb(ji,jj,jk) + zcorr * vmask(ji,jj,jk) … … 433 666 END DO 434 667 END DO 668 ! 669 IF ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN 670 vb_b(i1:i2,j1:j2) = vn_b(i1:i2,j1:j2) 671 ENDIF 672 ! 435 673 ENDIF 436 674 ! … … 438 676 439 677 440 SUBROUTINE updateSSH( tabres, i1, i2, j1, j2, before )678 SUBROUTINE updateSSH( tabres, i1, i2, j1, j2, before, nb, ndir ) 441 679 !!--------------------------------------------- 442 680 !! *** ROUTINE updateSSH *** … … 445 683 REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: tabres 446 684 LOGICAL, INTENT(in) :: before 685 INTEGER, INTENT(in) :: nb, ndir 447 686 !! 687 LOGICAL :: western_side, eastern_side, southern_side, northern_side 448 688 INTEGER :: ji, jj 449 689 !!--------------------------------------------- … … 472 712 END DO 473 713 END DO 714 ! 715 IF ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN 716 sshb(i1:i2,j1:j2) = sshn(i1:i2,j1:j2) 717 ENDIF 718 ! 719 # if defined DECAL_FEEDBACK 720 ! western_side = (nb == 1).AND.(ndir == 1) 721 ! eastern_side = (nb == 1).AND.(ndir == 2) 722 ! southern_side = (nb == 2).AND.(ndir == 1) 723 ! northern_side = (nb == 2).AND.(ndir == 2) 724 ! ! 725 ! ! Asselin correction 726 ! IF ( ln_dynspg_ts.AND.ln_bt_fw ) THEN 727 ! IF (southern_side) THEN 728 ! DO ji=i1,i2 729 ! sshn(ji,j1-1) = sshn(ji,j1-1) - rdt * r1_e2t(ji,j1-1) * (vb2_b_s(ji,j1-1)-vb2_b(ji,j1-1)) 730 ! END DO 731 ! ENDIF 732 ! IF (northern_side) THEN 733 ! DO ji=i1,i2 734 ! sshn(ji,j1+1) = sshn(ji,j1+1) + rdt * r1_e2t(ji,j1+1) * (vb2_b_s(ji,j1)-vb2_b(ji,j1)) 735 ! END DO 736 ! ENDIF 737 ! IF (western_side) THEN 738 ! DO jj=j1,j2 739 ! sshn(i1-1,jj) = sshn(i1-1,jj) - rdt * r1_e2t(i1-1,jj) * (ub2_b_s(i1-1,jj)-ub2_b(i1-1,jj)) 740 ! END DO 741 ! ENDIF 742 ! IF (eastern_side) THEN 743 ! DO jj=j1,j2 744 ! sshn(i1+1,jj) = sshn(i1+1,jj) + rdt * r1_e2t(i1+1,jj) * (ub2_b_s(i1,jj)-ub2_b(i1,jj)) 745 ! END DO 746 ! ENDIF 747 ! ! 748 ! ENDIF 749 #endif 474 750 ENDIF 475 751 ! … … 486 762 !! 487 763 INTEGER :: ji, jj 488 REAL(wp) :: zrhoy 764 REAL(wp) :: zrhoy, za1 489 765 !!--------------------------------------------- 490 766 ! … … 498 774 tabres = zrhoy * tabres 499 775 ELSE 776 za1 = 1._wp / REAL(Agrif_rhot(), wp) 777 tabres(i1:i2,j1:j2) = tabres(i1:i2,j1:j2) * r1_e2u(i1:i2,j1:j2) 500 778 DO jj=j1,j2 501 779 DO ji=i1,i2 502 ub2_b(ji,jj) = tabres(ji,jj) / e2u(ji,jj) 780 ub2_i_b(ji,jj) = ub2_i_b(ji,jj) & 781 & + za1 * (tabres(ji,jj) - ub2_b(ji,jj)) 782 ! ub2_b_s(ji,jj) = ub2_b(ji,jj) 783 ub2_b(ji,jj) = tabres(ji,jj) 503 784 END DO 504 785 END DO … … 517 798 !! 518 799 INTEGER :: ji, jj 519 REAL(wp) :: zrhox 800 REAL(wp) :: zrhox, za1 520 801 !!--------------------------------------------- 521 802 ! … … 529 810 tabres = zrhox * tabres 530 811 ELSE 812 za1 = 1._wp / REAL(Agrif_rhot(), wp) 813 tabres(i1:i2,j1:j2) = tabres(i1:i2,j1:j2) * r1_e1v(i1:i2,j1:j2) 531 814 DO jj=j1,j2 532 815 DO ji=i1,i2 533 vb2_b(ji,jj) = tabres(ji,jj) / e1v(ji,jj) 816 vb2_i_b(ji,jj) = vb2_i_b(ji,jj) & 817 & + za1 * (tabres(ji,jj) - vb2_b(ji,jj)) 818 ! vb2_b_s(ji,jj) = vb2_b(ji,jj) 819 vb2_b(ji,jj) = tabres(ji,jj) 534 820 END DO 535 821 END DO … … 644 930 # endif /* key_zdftke */ 645 931 932 SUBROUTINE updatee3t( ptab, i1, i2, j1, j2, k1, k2, before ) 933 !!--------------------------------------------- 934 !! *** ROUTINE updatee3t *** 935 !!--------------------------------------------- 936 INTEGER, INTENT(in) :: i1, i2, j1, j2, k1, k2 937 REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: ptab 938 LOGICAL, INTENT(in) :: before 939 INTEGER :: ji,jj,jk 940 REAL(wp) :: zcoef 941 !!--------------------------------------------- 942 ! 943 IF (before) THEN 944 !> jc tmp: 945 ! ptab(i1:i2,j1:j2,k1:k2) = e3t_n(i1:i2,j1:j2,k1:k2) 946 ptab(i1:i2,j1:j2,k1:k2) = e3t_n(i1:i2,j1:j2,k1:k2) / e3t_0(i1:i2,j1:j2,k1:k2) * tmask(i1:i2,j1:j2,k1:k2) 947 !< jc tmp: 948 ELSE 949 ! 950 ! 1) Updates at BEFORE time step: 951 ! ------------------------------- 952 ! 953 !> jc tmp: 954 ! DO jk = 1, jpkm1 955 ! DO jj=j1,j2 956 ! DO ji=i1,i2 957 ! IF (tmask(ji,jj,jk)==1) THEN 958 ! ptab(ji,jj,jk) = ptab(ji,jj,jk) * e3t_0(ji,jj,jk) 959 ! ELSE 960 ! ptab(ji,jj,jk) = e3t_0(ji,jj,jk) 961 ! ENDIF 962 ! END DO 963 ! END DO 964 ! END DO 965 ptab(i1:i2,j1:j2,k1:k2) = ptab(i1:i2,j1:j2,k1:k2) * e3t_0(i1:i2,j1:j2,k1:k2) 966 !< jc tmp: 967 968 ! Save "old" scale factor (prior update) for subsequent asselin correction 969 ! of prognostic variables (needed to update initial state only) 970 e3t_a(i1:i2,j1:j2,k1:k2) = e3t_n(i1:i2,j1:j2,k1:k2) 971 ! hdivb(i1:i2,j1:j2,k1:k2) = e3t_b(i1:i2,j1:j2,k1:k2) 972 973 IF ( (.NOT.(lk_agrif_fstep.AND.(neuler==0)).AND.(ln_dynspg_exp)) & 974 & .OR.(.NOT.(lk_agrif_fstep.AND.(neuler==0)).AND.(ln_dynspg_ts & 975 & .AND.(.NOT.ln_bt_fw)))) THEN 976 977 DO jk = 1, jpkm1 978 DO jj=j1,j2 979 DO ji=i1,i2 980 e3t_b(ji,jj,jk) = e3t_b(ji,jj,jk) & 981 & + atfp * ( ptab(ji,jj,jk) - e3t_n(ji,jj,jk) ) 982 END DO 983 END DO 984 END DO 985 ! 986 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) 987 gdepw_b(i1:i2,j1:j2,1) = 0.0_wp 988 gdept_b(i1:i2,j1:j2,1) = 0.5_wp * e3w_b(i1:i2,j1:j2,1) 989 ! 990 DO jk = 2, jpk 991 DO jj = j1,j2 992 DO ji = i1,i2 993 zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk)) 994 e3w_b(ji,jj,jk) = e3w_0(ji,jj,jk) + ( 1.0_wp - 0.5_wp * tmask(ji,jj,jk) ) * & 995 & ( e3t_b(ji,jj,jk-1) - e3t_0(ji,jj,jk-1) ) & 996 & + 0.5_wp * tmask(ji,jj,jk) * & 997 & ( e3t_b(ji,jj,jk ) - e3t_0(ji,jj,jk ) ) 998 gdepw_b(ji,jj,jk) = gdepw_b(ji,jj,jk-1) + e3t_b(ji,jj,jk-1) 999 gdept_b(ji,jj,jk) = zcoef * ( gdepw_b(ji,jj,jk ) + 0.5 * e3w_b(ji,jj,jk)) & 1000 & + (1-zcoef) * ( gdept_b(ji,jj,jk-1) + e3w_b(ji,jj,jk)) 1001 END DO 1002 END DO 1003 END DO 1004 ! 1005 ENDIF 1006 ! 1007 ! 2) Updates at NOW time step: 1008 ! ---------------------------- 1009 ! 1010 ! Update vertical scale factor at T-points: 1011 e3t_n(i1:i2,j1:j2,k1:k2) = ptab(i1:i2,j1:j2,k1:k2) 1012 ! 1013 ! Update total depth: 1014 ht_n(i1:i2,j1:j2) = 0._wp 1015 DO jk = 1, jpkm1 1016 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) 1017 END DO 1018 ! 1019 ! Update vertical scale factor at W-points and depths: 1020 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) 1021 gdept_n(i1:i2,j1:j2,1) = 0.5_wp * e3w_n(i1:i2,j1:j2,1) 1022 gdepw_n(i1:i2,j1:j2,1) = 0.0_wp 1023 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 ssh 1024 ! 1025 DO jk = 2, jpk 1026 DO jj = j1,j2 1027 DO ji = i1,i2 1028 zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk)) 1029 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) ) & 1030 & + 0.5_wp * tmask(ji,jj,jk) * ( e3t_n(ji,jj,jk ) - e3t_0(ji,jj,jk ) ) 1031 gdepw_n(ji,jj,jk) = gdepw_n(ji,jj,jk-1) + e3t_n(ji,jj,jk-1) 1032 gdept_n(ji,jj,jk) = zcoef * ( gdepw_n(ji,jj,jk ) + 0.5 * e3w_n(ji,jj,jk)) & 1033 & + (1-zcoef) * ( gdept_n(ji,jj,jk-1) + e3w_n(ji,jj,jk)) 1034 gde3w_n(ji,jj,jk) = gdept_n(ji,jj,jk) - (ht_n(ji,jj)-ht_0(ji,jj)) ! Last term in the rhs is ssh 1035 END DO 1036 END DO 1037 END DO 1038 ! 1039 IF ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN 1040 e3t_b (i1:i2,j1:j2,1:jpk) = e3t_n (i1:i2,j1:j2,1:jpk) 1041 e3w_b (i1:i2,j1:j2,1:jpk) = e3w_n (i1:i2,j1:j2,1:jpk) 1042 gdepw_b(i1:i2,j1:j2,1:jpk) = gdepw_n(i1:i2,j1:j2,1:jpk) 1043 gdept_b(i1:i2,j1:j2,1:jpk) = gdept_n(i1:i2,j1:j2,1:jpk) 1044 ENDIF 1045 ! 1046 ENDIF 1047 ! 1048 END SUBROUTINE updatee3t 1049 646 1050 #else 647 1051 CONTAINS
Note: See TracChangeset
for help on using the changeset viewer.
