Changeset 8741 for branches/2017/dev_r8624_AGRIF3_VVL/NEMOGCM/NEMO
- Timestamp:
- 2017-11-17T17:19:55+01:00 (7 years ago)
- Location:
- branches/2017/dev_r8624_AGRIF3_VVL/NEMOGCM/NEMO
- Files:
-
- 15 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2017/dev_r8624_AGRIF3_VVL/NEMOGCM/NEMO/NST_SRC/agrif_oce.F90
r5656 r8741 35 35 LOGICAL , PUBLIC :: spongedoneU = .FALSE. !: dynamics sponge layer indicator 36 36 LOGICAL , PUBLIC :: lk_agrif_fstep = .TRUE. !: if true: first step 37 LOGICAL , PUBLIC :: lk_agrif_doupd = .TRUE. !: if true: send update from current grid38 37 LOGICAL , PUBLIC :: lk_agrif_debug = .FALSE. !: if true: print debugging info 39 38 … … 65 64 INTEGER :: e3t_id, e1u_id, e2v_id, sshn_id 66 65 INTEGER :: scales_t_id 67 # if defined key_zdftke 66 # if defined key_zdftke || defined key_zdfgls 68 67 INTEGER :: avt_id, avm_id, en_id 69 68 # endif -
branches/2017/dev_r8624_AGRIF3_VVL/NEMOGCM/NEMO/NST_SRC/agrif_opa_interp.F90
r7646 r8741 24 24 USE agrif_oce 25 25 USE phycst 26 USE dynspg_ts, ONLY: un_adv, vn_adv 26 27 ! 27 28 USE in_out_manager … … 38 39 PUBLIC interpunb, interpvnb, interpub2b, interpvb2b 39 40 PUBLIC interpe3t, interpumsk, interpvmsk 40 # if defined key_zdftke 41 PUBLIC Agrif_ tke, interpavm41 # if defined key_zdftke || defined key_zdfgls 42 PUBLIC Agrif_avm, interpavm 42 43 # endif 43 44 … … 449 450 INTEGER :: ji, jj 450 451 LOGICAL :: ll_int_cons 451 REAL(wp) :: zrhot, zt452 452 !!---------------------------------------------------------------------- 453 453 ! … … 456 456 ll_int_cons = ln_bt_fw ! Assume conservative temporal integration in the forward case only 457 457 ! 458 zrhot = Agrif_rhot() 459 ! 460 ! "Central" time index for interpolation: 461 IF( ln_bt_fw ) THEN 462 zt = REAL( Agrif_NbStepint()+0.5_wp, wp ) / zrhot 463 ELSE 464 zt = REAL( Agrif_NbStepint() , wp ) / zrhot 465 ENDIF 466 ! 467 ! Linear interpolation of sea level 468 Agrif_SpecialValue = 0._wp 469 Agrif_UseSpecialValue = .TRUE. 470 CALL Agrif_Bc_variable( sshn_id, calledweight=zt, procname=interpsshn ) 471 Agrif_UseSpecialValue = .FALSE. 458 ! Enforce volume conservation if no time refinement: 459 IF ( Agrif_rhot()==1 ) ll_int_cons=.TRUE. 472 460 ! 473 461 ! Interpolate barotropic fluxes 474 Agrif_SpecialValue=0. 462 Agrif_SpecialValue=0._wp 475 463 Agrif_UseSpecialValue = ln_spc_dyn 476 464 ! … … 491 479 ubdy_n(:) = 0._wp ; vbdy_n(:) = 0._wp 492 480 ubdy_s(:) = 0._wp ; vbdy_s(:) = 0._wp 493 CALL Agrif_Bc_variable( unb_id, calledweight=zt,procname=interpunb )494 CALL Agrif_Bc_variable( vnb_id, calledweight=zt,procname=interpvnb )481 CALL Agrif_Bc_variable( unb_id, procname=interpunb ) 482 CALL Agrif_Bc_variable( vnb_id, procname=interpvnb ) 495 483 ENDIF 496 484 Agrif_UseSpecialValue = .FALSE. … … 501 489 SUBROUTINE Agrif_ssh( kt ) 502 490 !!---------------------------------------------------------------------- 503 !! *** ROUTINE Agrif_ DYN***491 !! *** ROUTINE Agrif_ssh *** 504 492 !!---------------------------------------------------------------------- 505 493 INTEGER, INTENT(in) :: kt 506 494 !! 495 INTEGER :: ji, jj 507 496 !!---------------------------------------------------------------------- 508 497 ! 509 498 IF( Agrif_Root() ) RETURN 499 ! 500 ! Linear interpolation in time of sea level 501 ! 502 Agrif_SpecialValue = 0._wp 503 Agrif_UseSpecialValue = .TRUE. 504 CALL Agrif_Bc_variable(sshn_id, procname=interpsshn ) 505 Agrif_UseSpecialValue = .FALSE. 510 506 ! 511 507 IF((nbondi == -1).OR.(nbondi == 2)) THEN 512 ssha(2,:)=ssha(3,:) 513 sshn(2,:)=sshn(3,:) 508 DO jj=1,jpj 509 ssha(2,jj) = hbdy_w(jj) 510 END DO 514 511 ENDIF 515 512 ! 516 513 IF((nbondi == 1).OR.(nbondi == 2)) THEN 517 ssha(nlci-1,:)=ssha(nlci-2,:) 518 sshn(nlci-1,:)=sshn(nlci-2,:) 514 DO jj=1,jpj 515 ssha(nlci-1,jj) = hbdy_e(jj) 516 END DO 519 517 ENDIF 520 518 ! 521 519 IF((nbondj == -1).OR.(nbondj == 2)) THEN 522 ssha(:,2)=ssha(:,3) 523 sshn(:,2)=sshn(:,3) 520 DO ji=1,jpi 521 ssha(ji,2) = hbdy_s(ji) 522 END DO 524 523 ENDIF 525 524 ! 526 525 IF((nbondj == 1).OR.(nbondj == 2)) THEN 527 ssha(:,nlcj-1)=ssha(:,nlcj-2) 528 sshn(:,nlcj-1)=sshn(:,nlcj-2) 526 DO ji=1,jpi 527 ssha(ji,nlcj-1) = hbdy_n(ji) 528 END DO 529 529 ENDIF 530 530 ! … … 541 541 !!---------------------------------------------------------------------- 542 542 ! 543 ! 544 IF( Agrif_Root() ) RETURN 545 ! 543 546 IF((nbondi == -1).OR.(nbondi == 2)) THEN 544 547 DO jj = 1, jpj … … 567 570 END SUBROUTINE Agrif_ssh_ts 568 571 569 # if defined key_zdftke 570 571 SUBROUTINE Agrif_ tke572 !!---------------------------------------------------------------------- 573 !! *** ROUTINE Agrif_ tke***572 # if defined key_zdftke || defined key_zdfgls 573 574 SUBROUTINE Agrif_avm 575 !!---------------------------------------------------------------------- 576 !! *** ROUTINE Agrif_avm *** 574 577 !!---------------------------------------------------------------------- 575 578 REAL(wp) :: zalpha 576 579 !!---------------------------------------------------------------------- 577 580 ! 578 zalpha = REAL( Agrif_NbStepint() + Agrif_IRhot() - 1, wp ) / REAL( Agrif_IRhot(), wp ) 579 IF( zalpha > 1. ) zalpha = 1. 581 IF( Agrif_Root() ) RETURN 582 ! 583 ! zalpha = REAL( Agrif_NbStepint() + Agrif_IRhot() - 1, wp ) / REAL( Agrif_IRhot(), wp ) 584 ! IF( zalpha > 1. ) zalpha = 1. 585 zalpha = 1._wp ! JC: proper time interpolation impossible 586 ! => use last available value from parent 580 587 ! 581 588 Agrif_SpecialValue = 0.e0 … … 586 593 Agrif_UseSpecialValue = .FALSE. 587 594 ! 588 END SUBROUTINE Agrif_ tke595 END SUBROUTINE Agrif_avm 589 596 590 597 # endif … … 781 788 ! 782 789 IF( before ) THEN 783 DO jk = k1, jpk790 DO jk = 1, jpkm1 784 791 ptab(i1:i2,j1:j2,jk) = e2u(i1:i2,j1:j2) * e3u_n(i1:i2,j1:j2,jk) * un(i1:i2,j1:j2,jk) 785 792 END DO … … 788 795 DO jk = 1, jpkm1 789 796 DO jj=j1,j2 790 ua(i1:i2,jj,jk) = ptab(i1:i2,jj,jk) / ( zrhoy * e2u(i1:i2,jj) * e3u_ n(i1:i2,jj,jk) )797 ua(i1:i2,jj,jk) = ptab(i1:i2,jj,jk) / ( zrhoy * e2u(i1:i2,jj) * e3u_a(i1:i2,jj,jk) ) 791 798 END DO 792 799 END DO … … 808 815 !!---------------------------------------------------------------------- 809 816 ! 810 IF( before ) THEN !interpv entre 1 et k2 et interpv2d en jpkp1811 DO jk = k1, jpk817 IF( before ) THEN 818 DO jk = 1, jpkm1 812 819 ptab(i1:i2,j1:j2,jk) = e1v(i1:i2,j1:j2) * e3v_n(i1:i2,j1:j2,jk) * vn(i1:i2,j1:j2,jk) 813 820 END DO … … 815 822 zrhox= Agrif_Rhox() 816 823 DO jk = 1, jpkm1 817 va(i1:i2,j1:j2,jk) = ptab(i1:i2,j1:j2,jk) / ( zrhox * e1v(i1:i2,j1:j2) * e3v_ n(i1:i2,j1:j2,jk) )824 va(i1:i2,j1:j2,jk) = ptab(i1:i2,j1:j2,jk) / ( zrhox * e1v(i1:i2,j1:j2) * e3v_a(i1:i2,j1:j2,jk) ) 818 825 END DO 819 826 ENDIF … … 978 985 !!---------------------------------------------------------------------- 979 986 IF( before ) THEN 980 ptab(i1:i2,j1:j2) = e2u(i1:i2,j1:j2) * ub2_b(i1:i2,j1:j2) 987 IF ( ln_bt_fw ) THEN 988 ptab(i1:i2,j1:j2) = e2u(i1:i2,j1:j2) * ub2_b(i1:i2,j1:j2) 989 ELSE 990 ptab(i1:i2,j1:j2) = e2u(i1:i2,j1:j2) * un_adv(i1:i2,j1:j2) 991 ENDIF 981 992 ELSE 982 993 western_side = (nb == 1).AND.(ndir == 1) … … 1016 1027 ! 1017 1028 IF( before ) THEN 1018 ptab(i1:i2,j1:j2) = e1v(i1:i2,j1:j2) * vb2_b(i1:i2,j1:j2) 1029 IF ( ln_bt_fw ) THEN 1030 ptab(i1:i2,j1:j2) = e1v(i1:i2,j1:j2) * vb2_b(i1:i2,j1:j2) 1031 ELSE 1032 ptab(i1:i2,j1:j2) = e1v(i1:i2,j1:j2) * vn_adv(i1:i2,j1:j2) 1033 ENDIF 1019 1034 ELSE 1020 1035 western_side = (nb == 1).AND.(ndir == 1) … … 1175 1190 END SUBROUTINE interpvmsk 1176 1191 1177 # if defined key_zdftke 1192 # if defined key_zdftke || defined key_zdfgls 1178 1193 1179 1194 SUBROUTINE interpavm( ptab, i1, i2, j1, j2, k1, k2, before ) … … 1189 1204 ptab (i1:i2,j1:j2,k1:k2) = avm_k(i1:i2,j1:j2,k1:k2) 1190 1205 ELSE 1191 avm _k(i1:i2,j1:j2,k1:k2) = ptab (i1:i2,j1:j2,k1:k2)1206 avm (i1:i2,j1:j2,k1:k2) = ptab (i1:i2,j1:j2,k1:k2) 1192 1207 ENDIF 1193 1208 ! 1194 1209 END SUBROUTINE interpavm 1195 1210 1196 # endif /* key_zdftke */1211 # endif /* key_zdftke || key_zdfgls */ 1197 1212 1198 1213 #else -
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 -
branches/2017/dev_r8624_AGRIF3_VVL/NEMOGCM/NEMO/NST_SRC/agrif_top_update.F90
r6140 r8741 33 33 CONTAINS 34 34 35 SUBROUTINE Agrif_Update_Trc( kt)35 SUBROUTINE Agrif_Update_Trc( ) 36 36 !!---------------------------------------------------------------------- 37 37 !! *** ROUTINE Agrif_Update_Trc *** 38 38 !!---------------------------------------------------------------------- 39 INTEGER, INTENT(in) :: kt40 !!----------------------------------------------------------------------41 39 ! 42 IF((Agrif_NbStepint() .NE. (Agrif_irhot()-1)).AND.(kt /= 0)) RETURN 40 IF (Agrif_Root()) RETURN 41 ! 43 42 #if defined TWO_WAY 44 43 Agrif_UseSpecialValueInUpdate = .TRUE. … … 66 65 67 66 68 SUBROUTINE updateTRC( ptab, i1, i2, j1, j2, k1, k2, n1, n2, before )67 SUBROUTINE updateTRC( ptab, i1, i2, j1, j2, k1, k2, n1, n2, before, nb, ndir ) 69 68 !!---------------------------------------------------------------------- 70 69 !! *** ROUTINE updateT *** … … 73 72 REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: ptab 74 73 LOGICAL , INTENT(in ) :: before 74 INTEGER, INTENT(in) :: nb, ndir 75 75 !! 76 INTEGER :: ji, jj, jk, jn 77 !!---------------------------------------------------------------------- 78 ! 79 IF( before ) THEN 80 ptab(i1:i2,j1:j2,k1:k2,n1:n2) = trn(i1:i2,j1:j2,k1:k2,n1:n2) 81 ELSE 82 IF( .NOT.(lk_agrif_fstep.AND.(neuler==0)) ) THEN 83 ! Add asselin part 84 DO jn = n1,n2 85 DO jk = k1, k2 86 DO jj = j1, j2 87 DO ji = i1, i2 88 IF( ptab(ji,jj,jk,jn) /= 0._wp ) THEN 89 trb(ji,jj,jk,jn) = trb(ji,jj,jk,jn) & 90 & + atfp * ( ptab(ji,jj,jk,jn) & 91 & - trn(ji,jj,jk,jn) ) * tmask(ji,jj,jk) 92 ENDIF 93 END DO 76 LOGICAL :: western_side, eastern_side, southern_side, northern_side 77 INTEGER :: ji,jj,jk,jn 78 REAL(wp) :: ztb, ztnu, ztno 79 !!---------------------------------------------------------------------- 80 ! 81 ! 82 IF (before) THEN 83 DO jn = n1,n2 84 DO jk=k1,k2 85 DO jj=j1,j2 86 DO ji=i1,i2 87 !> jc tmp 88 tabres(ji,jj,jk,jn) = trn(ji,jj,jk,jn) * e3t_n(ji,jj,jk) / e3t_0(ji,jj,jk) 89 ! tabres(ji,jj,jk,jn) = trn(ji,jj,jk,jn) * e3t_n(ji,jj,jk) 90 !< jc tmp 94 91 END DO 95 92 END DO 96 93 END DO 94 END DO 95 ELSE 96 !> jc tmp 97 DO jn = n1,n2 98 tabres(i1:i2,j1:j2,k1:k2,jn) = tabres(i1:i2,j1:j2,k1:k2,jn) * e3t_0(i1:i2,j1:j2,k1:k2) & 99 & * tmask(i1:i2,j1:j2,k1:k2) 100 ENDDO 101 !< jc tmp 102 IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN 103 ! Add asselin part 104 DO jn = n1,n2 105 DO jk=k1,k2 106 DO jj=j1,j2 107 DO ji=i1,i2 108 IF( tabres(ji,jj,jk,jn) .NE. 0. ) THEN 109 ztb = trb(ji,jj,jk,jn) * e3t_b(ji,jj,jk) ! fse3t_b prior update should be used 110 ztnu = tabres(ji,jj,jk,jn) 111 ztno = trn(ji,jj,jk,jn) * e3t_a(ji,jj,jk) 112 trb(ji,jj,jk,jn) = ( ztb + atfp * ( ztnu - ztno) ) & 113 & * tmask(ji,jj,jk) / e3t_b(ji,jj,jk) 114 ENDIF 115 ENDDO 116 ENDDO 117 ENDDO 118 ENDDO 97 119 ENDIF 98 DO jn = n1, 99 DO jk = k1,k2100 DO jj = j1,j2101 DO ji = i1,i2102 IF( ptab(ji,jj,jk,jn) /= 0._wp) THEN103 trn(ji,jj,jk,jn) = ptab(ji,jj,jk,jn) * tmask(ji,jj,jk)120 DO jn = n1,n2 121 DO jk=k1,k2 122 DO jj=j1,j2 123 DO ji=i1,i2 124 IF( tabres(ji,jj,jk,jn) .NE. 0. ) THEN 125 trn(ji,jj,jk,jn) = tabres(ji,jj,jk,jn) / e3t_n(ji,jj,jk) 104 126 END IF 105 127 END DO … … 107 129 END DO 108 130 END DO 131 ! 132 IF ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN 133 trb(i1:i2,j1:j2,k1:k2,n1:n2) = trn(i1:i2,j1:j2,k1:k2,n1:n2) 134 ENDIF 135 ! 136 ! 137 # if defined DECAL_FEEDBACK 138 IF (.NOT.ln_linssh) THEN 139 western_side = (nb == 1).AND.(ndir == 1) 140 eastern_side = (nb == 1).AND.(ndir == 2) 141 southern_side = (nb == 2).AND.(ndir == 1) 142 northern_side = (nb == 2).AND.(ndir == 2) 143 ! 144 ! Asselin correction 145 IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN 146 IF (southern_side) THEN 147 DO jn = n1,n2 148 DO jk=k1,k2 149 DO ji=i1,i2 150 ztb = trb(ji,j1-1,jk,jn) * e3t_b(ji,j1-1,jk) ! fse3t_b prior update should be used 151 ztnu = trn(ji,j1-1,jk,jn) * e3t_n(ji,j1-1,jk) 152 ztno = trn(ji,j1-1,jk,jn) * e3t_a(ji,j1-1,jk) 153 trb(ji,j1-1,jk,jn) = ( ztb + atfp * ( ztnu - ztno) ) & 154 & * tmask(ji,j1-1,jk) / e3t_b(ji,j1-1,jk) 155 END DO 156 ENDDO 157 ENDDO 158 ENDIF 159 IF (northern_side) THEN 160 DO jn = n1,n2 161 DO jk=k1,k2 162 DO ji=i1,i2 163 ztb = trb(ji,j2+1,jk,jn) * e3t_b(ji,j2+1,jk) ! fse3t_b prior update should be used 164 ztnu = trn(ji,j2+1,jk,jn) * e3t_n(ji,j2+1,jk) 165 ztno = trn(ji,j2+1,jk,jn) * e3t_a(ji,j2+1,jk) 166 trb(ji,j2+1,jk,jn) = ( ztb + atfp * ( ztnu - ztno) ) & 167 & * tmask(ji,j2+1,jk) / e3t_b(ji,j2+1,jk) 168 END DO 169 ENDDO 170 ENDDO 171 ENDIF 172 IF (western_side) THEN 173 DO jn = n1,n2 174 DO jk=k1,k2 175 DO jj=j1,j2 176 ztb = trb(i1-1,jj,jk,jn) * e3t_b(i1-1,jj,jk) ! fse3t_b prior update should be used 177 ztnu = trn(i1-1,jj,jk,jn) * e3t_n(i1-1,jj,jk) 178 ztno = trn(i1-1,jj,jk,jn) * e3t_a(i1-1,jj,jk) 179 trb(i1-1,jj,jk,jn) = ( ztb + atfp * ( ztnu - ztno) ) & 180 & * tmask(i1-1,jj,jk) / e3t_b(i1-1,jj,jk) 181 END DO 182 ENDDO 183 ENDDO 184 ENDIF 185 IF (eastern_side) THEN 186 DO jn = n1,n2 187 DO jk=k1,k2 188 DO jj=j1,j2 189 ztb = trb(i2+1,jj,jk,jn) * e3t_b(i2+1,jj,jk) ! fse3t_b prior update should be used 190 ztnu = trn(i2+1,jj,jk,jn) * e3t_n(i2+1,jj,jk) 191 ztno = trn(i2+1,jj,jk,jn) * e3t_a(i2+1,jj,jk) 192 trb(i2+1,jj,jk,jn) = ( ztb + atfp * ( ztnu - ztno) ) & 193 & * tmask(i2+1,jj,jk) / e3t_b(i2+1,jj,jk) 194 END DO 195 ENDDO 196 ENDDO 197 ENDIF 198 ENDIF ! Asselin correction 199 200 IF (southern_side) THEN 201 DO jn = n1,n2 202 DO jk=k1,k2 203 DO ji=i1,i2 204 trn(ji,j1-1,jk,jn) = trn(ji,j1-1,jk,jn) * e3t_a(ji,j1-1,jk) / e3t_n(ji,j1-1,jk) 205 END DO 206 ENDDO 207 ENDDO 208 ENDIF 209 IF (northern_side) THEN 210 DO jn = n1,n2 211 DO jk=k1,k2 212 DO ji=i1,i2 213 trn(ji,j2+1,jk,jn) = trn(ji,j2+1,jk,jn) * e3t_a(ji,j2+1,jk) / e3t_n(ji,j2+1,jk) 214 END DO 215 ENDDO 216 ENDDO 217 ENDIF 218 IF (western_side) THEN 219 DO jn = n1,n2 220 DO jk=k1,k2 221 DO jj=j1,j2 222 trn(i1-1,jj,jk,jn) = trn(i1-1,jj,jk,jn) * e3t_a(i1-1,jj,jk) / e3t_n(i1-1,jj,jk) 223 END DO 224 ENDDO 225 ENDDO 226 ENDIF 227 IF (eastern_side) THEN 228 DO jn = n1,n2 229 DO jk=k1,k2 230 DO jj=j1,j2 231 trn(i2+1,jj,jk,jn) = trn(i2+1,jj,jk,jn) * e3t_a(i2+1,jj,jk) / e3t_n(i2+1,jj,jk) 232 END DO 233 ENDDO 234 ENDDO 235 ENDIF 236 ENDIF 237 #endif 109 238 ENDIF 110 239 ! -
branches/2017/dev_r8624_AGRIF3_VVL/NEMOGCM/NEMO/NST_SRC/agrif_user.F90
r7761 r8741 1 #undef UPD_HIGH /* MIX HIGH UPDATE */ 1 2 #if defined key_agrif 2 3 !!---------------------------------------------------------------------- … … 88 89 # endif 89 90 ! 91 nbcline = 0 92 #if defined key_top 93 nbcline_trc = 0 94 #endif 95 ! 96 IF ( Agrif_Level().EQ.Agrif_MaxLevel() ) CALL agrif_Update_ini() 97 98 Agrif_UseSpecialValueInUpdate = .FALSE. 99 90 100 END SUBROUTINE Agrif_initvalues 91 101 … … 144 154 CALL Agrif_Set_bc(e2v_id,(/0,0/)) 145 155 146 ! 5. Update type156 ! 4. Update type 147 157 !--------------- 158 # if defined UPD_HIGH 159 CALL Agrif_Set_Updatetype(e1u_id,update1 = Agrif_Update_Average, update2=Agrif_Update_Full_Weighting) 160 CALL Agrif_Set_Updatetype(e2v_id,update1 = Agrif_Update_Full_Weighting, update2=Agrif_Update_Average) 161 #else 148 162 CALL Agrif_Set_Updatetype(e1u_id,update1 = Agrif_Update_Copy, update2=Agrif_Update_Average) 149 163 CALL Agrif_Set_Updatetype(e2v_id,update1 = Agrif_Update_Average, update2=Agrif_Update_Copy) 150 151 ! High order updates 152 ! CALL Agrif_Set_Updatetype(e1u_id,update1 = Agrif_Update_Average, update2=Agrif_Update_Full_Weighting) 153 ! CALL Agrif_Set_Updatetype(e2v_id,update1 = Agrif_Update_Full_Weighting, update2=Agrif_Update_Average) 154 ! 164 #endif 165 155 166 END SUBROUTINE agrif_declare_var_dom 156 167 … … 175 186 ! 176 187 LOGICAL :: check_namelist 177 CHARACTER(len=15) :: cl_check1, cl_check2, cl_check3 188 CHARACTER(len=15) :: cl_check1, cl_check2, cl_check3, cl_check4 178 189 !!---------------------------------------------------------------------- 179 190 … … 205 216 Agrif_UseSpecialValue = .TRUE. 206 217 CALL Agrif_Bc_variable(sshn_id,calledweight=1., procname=interpsshn ) 218 hbdy_w(:) = 0.e0 ; hbdy_e(:) = 0.e0 ; hbdy_n(:) = 0.e0 ; hbdy_s(:) = 0.e0 219 ssha(:,:) = 0.e0 207 220 208 221 IF ( ln_dynspg_ts ) THEN … … 212 225 CALL Agrif_Bc_variable(ub2b_interp_id,calledweight=1.,procname=interpub2b) 213 226 CALL Agrif_Bc_variable(vb2b_interp_id,calledweight=1.,procname=interpvb2b) 214 ubdy_w(:) = 0.e0 ; vbdy_w(:) = 0.e0 ; hbdy_w(:) =0.e0215 ubdy_e(:) = 0.e0 ; vbdy_e(:) = 0.e0 ; hbdy_e(:) =0.e0216 ubdy_n(:) = 0.e0 ; vbdy_n(:) = 0.e0 ; hbdy_n(:) =0.e0217 ubdy_s(:) = 0.e0 ; vbdy_s(:) = 0.e0 ; hbdy_s(:) =0.e0227 ubdy_w(:) = 0.e0 ; vbdy_w(:) = 0.e0 228 ubdy_e(:) = 0.e0 ; vbdy_e(:) = 0.e0 229 ubdy_n(:) = 0.e0 ; vbdy_n(:) = 0.e0 230 ubdy_s(:) = 0.e0 ; vbdy_s(:) = 0.e0 218 231 ENDIF 219 232 … … 234 247 WRITE(cl_check2,*) NINT(rdt) 235 248 WRITE(cl_check3,*) NINT(Agrif_Parent(rdt)/Agrif_Rhot()) 236 CALL ctl_stop( ' incompatible time step between ocean grids', &249 CALL ctl_stop( 'Incompatible time step between ocean grids', & 237 250 & 'parent grid value : '//cl_check1 , & 238 251 & 'child grid value : '//cl_check2 , & … … 245 258 WRITE(cl_check1,*) (Agrif_Parent(nit000)-1)*Agrif_IRhot() + 1 246 259 WRITE(cl_check2,*) Agrif_Parent(nitend) *Agrif_IRhot() 247 CALL ctl_warn( ' incompatible run length between grids', &248 & ' nit000 on fine grid will be changeto : '//cl_check1, &249 & ' nitend on fine grid will be changeto : '//cl_check2 )260 CALL ctl_warn( 'Incompatible run length between grids' , & 261 & 'nit000 on fine grid will be changed to : '//cl_check1, & 262 & 'nitend on fine grid will be changed to : '//cl_check2 ) 250 263 nit000 = (Agrif_Parent(nit000)-1)*Agrif_IRhot() + 1 251 264 nitend = Agrif_Parent(nitend) *Agrif_IRhot() 252 265 ENDIF 253 266 254 ! Check coordinates255 !SF IF( ln_zps ) THEN256 !SF ! check parameters for partial steps257 !SF IF( Agrif_Parent(e3zps_min) .NE. e3zps_min ) THEN258 !SF WRITE(*,*) 'incompatible e3zps_min between grids'259 !SF WRITE(*,*) 'parent grid :',Agrif_Parent(e3zps_min)260 !SF WRITE(*,*) 'child grid :',e3zps_min261 !SF WRITE(*,*) 'those values should be identical'262 !SF STOP263 !SF ENDIF264 !SF IF( Agrif_Parent(e3zps_rat) /= e3zps_rat ) THEN265 !SF WRITE(*,*) 'incompatible e3zps_rat between grids'266 !SF WRITE(*,*) 'parent grid :',Agrif_Parent(e3zps_rat)267 !SF WRITE(*,*) 'child grid :',e3zps_rat268 !SF WRITE(*,*) 'those values should be identical'269 !SF STOP270 !SF ENDIF271 !SF ENDIF272 273 267 ! Check free surface scheme 274 268 IF ( ( Agrif_Parent(ln_dynspg_ts ).AND.ln_dynspg_exp ).OR.& 275 269 & ( Agrif_Parent(ln_dynspg_exp).AND.ln_dynspg_ts ) ) THEN 276 WRITE(*,*) 'incompatible free surface scheme between grids' 277 WRITE(*,*) 'parent grid ln_dynspg_ts :', Agrif_Parent(ln_dynspg_ts ) 278 WRITE(*,*) 'parent grid ln_dynspg_exp :', Agrif_Parent(ln_dynspg_exp) 279 WRITE(*,*) 'child grid ln_dynspg_ts :', ln_dynspg_ts 280 WRITE(*,*) 'child grid ln_dynspg_exp :', ln_dynspg_exp 281 WRITE(*,*) 'those logicals should be identical' 270 WRITE(cl_check1,*) Agrif_Parent( ln_dynspg_ts ) 271 WRITE(cl_check2,*) ln_dynspg_ts 272 WRITE(cl_check3,*) Agrif_Parent( ln_dynspg_exp ) 273 WRITE(cl_check4,*) ln_dynspg_exp 274 CALL ctl_stop( 'Incompatible free surface scheme between grids' , & 275 & 'parent grid ln_dynspg_ts :'//cl_check1 , & 276 & 'child grid ln_dynspg_ts :'//cl_check2 , & 277 & 'parent grid ln_dynspg_exp :'//cl_check3 , & 278 & 'child grid ln_dynspg_exp :'//cl_check4 , & 279 & 'those logicals should be identical' ) 280 STOP 281 ENDIF 282 283 ! Check if identical linear free surface option 284 IF ( ( Agrif_Parent(ln_linssh ).AND.(.NOT.ln_linssh )).OR.& 285 & ( (.NOT.Agrif_Parent(ln_linssh)).AND.ln_linssh ) ) THEN 286 WRITE(cl_check1,*) Agrif_Parent(ln_linssh ) 287 WRITE(cl_check2,*) ln_linssh 288 CALL ctl_stop( 'Incompatible linearized fs option between grids', & 289 & 'parent grid ln_linssh :'//cl_check1 , & 290 & 'child grid ln_linssh :'//cl_check2 , & 291 & 'those logicals should be identical' ) 282 292 STOP 283 293 ENDIF … … 306 316 ENDIF 307 317 ! 308 ! Do update at initialisation because not done before writing restarts 309 ! This would indeed change boundary conditions values at initial time 310 ! hence produce restartability issues. 311 ! Note that update below is recursive (with lk_agrif_doupd=T): 312 ! 313 ! JC: I am not sure if Agrif_MaxLevel() is the "relative" 314 ! or the absolute maximum nesting level...TBC 315 IF ( Agrif_Level().EQ.Agrif_MaxLevel() ) THEN 316 ! NB: Do tracers first, dynamics after because nbcline incremented in dynamics 317 CALL Agrif_Update_tra() 318 CALL Agrif_Update_dyn() 319 ENDIF 320 ! 318 END SUBROUTINE Agrif_InitValues_cont 319 320 RECURSIVE SUBROUTINE Agrif_Update_ini( ) 321 !!---------------------------------------------------------------------- 322 !! *** ROUTINE agrif_Update_ini *** 323 !! 324 !! ** Purpose :: Recursive update done at initialization 325 !!---------------------------------------------------------------------- 326 USE dom_oce 327 USE agrif_opa_update 328 #if defined key_top 329 USE agrif_top_update 330 #endif 331 ! 332 IMPLICIT NONE 333 !!---------------------------------------------------------------------- 334 ! 335 IF (Agrif_Root()) RETURN 336 ! 337 IF (.NOT.ln_linssh) CALL Agrif_Update_vvl() 338 CALL Agrif_Update_tra() 339 #if defined key_top 340 CALL Agrif_Update_Trc() 341 #endif 342 CALL Agrif_Update_dyn() 321 343 # if defined key_zdftke 322 CALL Agrif_Update_tke(0) 323 # endif 324 ! 325 Agrif_UseSpecialValueInUpdate = .FALSE. 326 nbcline = 0327 lk_agrif_doupd = .FALSE.328 !329 END SUBROUTINE Agrif_InitValues_cont 330 344 ! JC remove update because this precludes from perfect restartability 345 !! CALL Agrif_Update_tke() 346 # endif 347 348 CALL Agrif_ChildGrid_To_ParentGrid() 349 CALL Agrif_Update_ini() 350 CALL Agrif_ParentGrid_To_ChildGrid() 351 352 END SUBROUTINE agrif_update_ini 331 353 332 354 SUBROUTINE agrif_declare_var … … 371 393 CALL agrif_declare_variable((/2,2/),(/3,3/),(/'x','y'/),(/1,1/),(/nlci,nlcj/),sshn_id) 372 394 373 # if defined key_zdftke 374 CALL agrif_declare_variable((/2,2,0/),(/3,3,0/),(/'x','y','N'/),(/1,1,1/),(/ jpi,jpj,jpk/), en_id)375 CALL agrif_declare_variable((/2,2,0/),(/3,3,0/),(/'x','y','N'/),(/1,1,1/),(/ jpi,jpj,jpk/),avt_id)376 CALL agrif_declare_variable((/2,2,0/),(/3,3,0/),(/'x','y','N'/),(/1,1,1/),(/ jpi,jpj,jpk/),avm_id)395 # if defined key_zdftke || defined key_zdfgls 396 CALL agrif_declare_variable((/2,2,0/),(/3,3,0/),(/'x','y','N'/),(/1,1,1/),(/nlci,nlcj,jpk/), en_id) 397 CALL agrif_declare_variable((/2,2,0/),(/3,3,0/),(/'x','y','N'/),(/1,1,1/),(/nlci,nlcj,jpk/),avt_id) 398 CALL agrif_declare_variable((/2,2,0/),(/3,3,0/),(/'x','y','N'/),(/1,1,1/),(/nlci,nlcj,jpk/),avm_id) 377 399 # endif 378 400 … … 400 422 CALL Agrif_Set_bcinterp(vmsk_id,interp=AGRIF_constant) 401 423 402 # if defined key_zdftke 424 # if defined key_zdftke || defined key_zdfgls 403 425 CALL Agrif_Set_bcinterp(avm_id ,interp=AGRIF_linear) 404 426 # endif … … 411 433 CALL Agrif_Set_bc(vn_interp_id,(/0,1/)) 412 434 413 ! CALL Agrif_Set_bc(tsn_sponge_id,(/-3*Agrif_irhox(),0/))414 ! CALL Agrif_Set_bc(un_sponge_id,(/-2*Agrif_irhox()-1,0/))415 ! CALL Agrif_Set_bc(vn_sponge_id,(/-2*Agrif_irhox()-1,0/))416 435 CALL Agrif_Set_bc(tsn_sponge_id,(/-nn_sponge_len*Agrif_irhox()-1,0/)) 417 436 CALL Agrif_Set_bc(un_sponge_id ,(/-nn_sponge_len*Agrif_irhox()-1,0/)) … … 428 447 CALL Agrif_Set_bc(vmsk_id,(/0,0/)) 429 448 449 # if defined key_zdftke || defined key_zdfgls 450 CALL Agrif_Set_bc(avm_id ,(/0,1/)) 451 # endif 452 453 ! 4. Update type 454 !--------------- 455 CALL Agrif_Set_Updatetype(scales_t_id, update = AGRIF_Update_Average) 456 457 # if defined UPD_HIGH 458 CALL Agrif_Set_Updatetype(tsn_id, update = Agrif_Update_Full_Weighting) 459 CALL Agrif_Set_Updatetype(un_update_id,update1 = Agrif_Update_Average, update2 = Agrif_Update_Full_Weighting) 460 CALL Agrif_Set_Updatetype(vn_update_id,update1 = Agrif_Update_Full_Weighting, update2 = Agrif_Update_Average) 461 462 CALL Agrif_Set_Updatetype(ub2b_update_id,update1 = Agrif_Update_Average, update2 = Agrif_Update_Full_Weighting) 463 CALL Agrif_Set_Updatetype(vb2b_update_id,update1 = Agrif_Update_Full_Weighting, update2 = Agrif_Update_Average) 464 CALL Agrif_Set_Updatetype(sshn_id, update = Agrif_Update_Full_Weighting) 465 CALL Agrif_Set_Updatetype(e3t_id, update = Agrif_Update_Full_Weighting) 466 430 467 # if defined key_zdftke 431 CALL Agrif_Set_bc(avm_id ,(/0,1/)) 432 # endif 433 434 ! 5. Update type 435 !--------------- 468 CALL Agrif_Set_Updatetype( en_id, update = AGRIF_Update_Full_Weighting) 469 CALL Agrif_Set_Updatetype(avt_id, update = AGRIF_Update_Full_Weighting) 470 CALL Agrif_Set_Updatetype(avm_id, update = AGRIF_Update_Full_Weighting) 471 # endif 472 473 #else 436 474 CALL Agrif_Set_Updatetype(tsn_id, update = AGRIF_Update_Average) 437 438 CALL Agrif_Set_Updatetype(scales_t_id, update = AGRIF_Update_Average)439 440 475 CALL Agrif_Set_Updatetype(un_update_id,update1 = Agrif_Update_Copy, update2 = Agrif_Update_Average) 441 476 CALL Agrif_Set_Updatetype(vn_update_id,update1 = Agrif_Update_Average, update2 = Agrif_Update_Copy) 442 477 443 CALL Agrif_Set_Updatetype(sshn_id, update = AGRIF_Update_Average)444 445 478 CALL Agrif_Set_Updatetype(ub2b_update_id,update1 = Agrif_Update_Copy, update2 = Agrif_Update_Average) 446 479 CALL Agrif_Set_Updatetype(vb2b_update_id,update1 = Agrif_Update_Average, update2 = Agrif_Update_Copy) 480 CALL Agrif_Set_Updatetype(sshn_id, update = AGRIF_Update_Average) 481 CALL Agrif_Set_Updatetype(e3t_id, update = AGRIF_Update_Average) 447 482 448 483 # if defined key_zdftke … … 452 487 # endif 453 488 454 ! High order updates 455 ! CALL Agrif_Set_Updatetype(tsn_id, update = Agrif_Update_Full_Weighting) 456 ! CALL Agrif_Set_Updatetype(un_update_id,update1 = Agrif_Update_Average, update2 = Agrif_Update_Full_Weighting) 457 ! CALL Agrif_Set_Updatetype(vn_update_id,update1 = Agrif_Update_Full_Weighting, update2 = Agrif_Update_Average) 458 ! 459 ! CALL Agrif_Set_Updatetype(ub2b_update_id,update1 = Agrif_Update_Average, update2 = Agrif_Update_Full_Weighting) 460 ! CALL Agrif_Set_Updatetype(vb2b_update_id,update1 = Agrif_Update_Full_Weighting, update2 = Agrif_Update_Average) 461 ! CALL Agrif_Set_Updatetype(sshn_id, update = Agrif_Update_Full_Weighting) 462 489 #endif 463 490 ! 464 491 END SUBROUTINE agrif_declare_var … … 733 760 ENDIF 734 761 735 ! Check coordinates736 IF( ln_zps ) THEN737 ! check parameters for partial steps738 IF( Agrif_Parent(e3zps_min) .NE. e3zps_min ) THEN739 WRITE(*,*) 'incompatible e3zps_min between grids'740 WRITE(*,*) 'parent grid :',Agrif_Parent(e3zps_min)741 WRITE(*,*) 'child grid :',e3zps_min742 WRITE(*,*) 'those values should be identical'743 STOP744 ENDIF745 IF( Agrif_Parent(e3zps_rat) .NE. e3zps_rat ) THEN746 WRITE(*,*) 'incompatible e3zps_rat between grids'747 WRITE(*,*) 'parent grid :',Agrif_Parent(e3zps_rat)748 WRITE(*,*) 'child grid :',e3zps_rat749 WRITE(*,*) 'those values should be identical'750 STOP751 ENDIF752 762 ENDIF 753 763 ! Check passive tracer cell … … 756 766 ENDIF 757 767 ENDIF 758 759 CALL Agrif_Update_trc(0)760 !761 Agrif_UseSpecialValueInUpdate = .FALSE.762 nbcline_trc = 0763 768 ! 764 769 END SUBROUTINE Agrif_InitValues_cont_top … … 792 797 !----------------------------- 793 798 CALL Agrif_Set_bc(trn_id,(/0,1/)) 794 ! CALL Agrif_Set_bc(trn_sponge_id,(/-3*Agrif_irhox(),0/))795 799 CALL Agrif_Set_bc(trn_sponge_id,(/-nn_sponge_len*Agrif_irhox()-1,0/)) 796 800 797 ! 5. Update type801 ! 4. Update type 798 802 !--------------- 803 # if defined UPD_HIGH 804 CALL Agrif_Set_Updatetype(trn_id, update = Agrif_Update_Full_Weighting) 805 #else 799 806 CALL Agrif_Set_Updatetype(trn_id, update = AGRIF_Update_Average) 800 801 ! Higher order update 802 ! CALL Agrif_Set_Updatetype(tsn_id, update = Agrif_Update_Full_Weighting) 803 807 #endif 804 808 ! 805 809 END SUBROUTINE agrif_declare_var_top … … 866 870 IF (lk_mpp) iminspon = MIN(iminspon,FLOOR(REAL(jpi-2)/REAL(Agrif_irhox())), FLOOR(REAL(jpj-2)/REAL(Agrif_irhox())) ) 867 871 IF (nn_sponge_len > iminspon) CALL ctl_stop('agrif sponge length is too large') 872 ! Check update frequency 873 IF (MOD((nitend-nit000+1), nbclineupdate).NE.0 ) CALL ctl_stop('number of time steps should be a multiple of nn_cln_update') 868 874 ! 869 875 IF( agrif_oce_alloc() > 0 ) CALL ctl_warn('agrif agrif_oce_alloc: allocation of arrays failed') … … 878 884 SUBROUTINE Agrif_InvLoc( indloc, nprocloc, i, indglob ) 879 885 !!---------------------------------------------------------------------- 880 !! *** ROUTINE Agrif_ detect***886 !! *** ROUTINE Agrif_InvLoc *** 881 887 !!---------------------------------------------------------------------- 882 888 USE dom_oce -
branches/2017/dev_r8624_AGRIF3_VVL/NEMOGCM/NEMO/OPA_SRC/DOM/dommsk.F90
r7753 r8741 276 276 ENDIF 277 277 END DO 278 #if defined key_agrif 279 IF( .NOT. AGRIF_Root() ) THEN 280 IF ((nbondi == 1).OR.(nbondi == 2)) fmask(nlci-1 , : ,jk) = 0.e0 ! east 281 IF ((nbondi == -1).OR.(nbondi == 2)) fmask(1 , : ,jk) = 0.e0 ! west 282 IF ((nbondj == 1).OR.(nbondj == 2)) fmask(: ,nlcj-1 ,jk) = 0.e0 ! north 283 IF ((nbondj == -1).OR.(nbondj == 2)) fmask(: ,1 ,jk) = 0.e0 ! south 284 ENDIF 285 #endif 278 286 END DO 279 287 ! -
branches/2017/dev_r8624_AGRIF3_VVL/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90
r7753 r8741 1029 1029 ! 1030 1030 #if defined key_agrif 1031 IF( .NOT.Agrif_Root() ) CALL ctl_stop( 'AGRIF not implemented with non-linear free surface' )1031 IF( (.NOT.Agrif_Root()).AND.(.NOT.ln_vvl_zstar) )CALL ctl_stop( 'AGRIF is implemented with zstar coordinate only' ) 1032 1032 #endif 1033 1033 ! -
branches/2017/dev_r8624_AGRIF3_VVL/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90
r7753 r8741 132 132 ! so that asselin contribution is removed at the same time 133 133 DO jk = 1, jpkm1 134 un(:,:,jk) = ( un(:,:,jk) - un_adv(:,:) + un_b(:,:) )*umask(:,:,jk)135 vn(:,:,jk) = ( vn(:,:,jk) - vn_adv(:,:) + vn_b(:,:) )*vmask(:,:,jk)134 un(:,:,jk) = ( un(:,:,jk) - un_adv(:,:)*r1_hu_n(:,:) + un_b(:,:) )*umask(:,:,jk) 135 vn(:,:,jk) = ( vn(:,:,jk) - vn_adv(:,:)*r1_hv_n(:,:) + vn_b(:,:) )*vmask(:,:,jk) 136 136 END DO 137 137 ENDIF -
branches/2017/dev_r8624_AGRIF3_VVL/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90
r7831 r8741 52 52 #if defined key_agrif 53 53 USE agrif_opa_interp ! agrif 54 USE agrif_oce 54 55 #endif 55 56 #if defined key_asminc … … 76 77 77 78 !! Time filtered arrays at baroclinic time step: 78 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: un_adv , vn_adv !: Advection vel.at "now" barocl. step79 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: un_adv , vn_adv !: Advection fluxes at "now" barocl. step 79 80 80 81 !! * Substitutions … … 127 128 !! -Update the filtered free surface at step "n+1" : ssha 128 129 !! -Update filtered barotropic velocities at step "n+1" : ua_b, va_b 129 !! -Compute barotropic advective velocities at step "n": un_adv, vn_adv130 !! -Compute barotropic advective fluxes at step "n" : un_adv, vn_adv 130 131 !! These are used to advect tracers and are compliant with discrete 131 132 !! continuity equation taken at the baroclinic time steps. This … … 1019 1020 ! 1020 1021 ! Set advection velocity correction: 1021 zwx(:,:) = un_adv(:,:) 1022 zwy(:,:) = vn_adv(:,:) 1023 IF( ( kt == nit000 .AND. neuler==0 ) .OR. .NOT.ln_bt_fw ) THEN 1024 un_adv(:,:) = zwx(:,:) * r1_hu_n(:,:) 1025 vn_adv(:,:) = zwy(:,:) * r1_hv_n(:,:) 1026 ELSE 1027 un_adv(:,:) = z1_2 * ( ub2_b(:,:) + zwx(:,:) ) * r1_hu_n(:,:) 1028 vn_adv(:,:) = z1_2 * ( vb2_b(:,:) + zwy(:,:) ) * r1_hv_n(:,:) 1029 END IF 1030 1031 IF( ln_bt_fw ) THEN ! Save integrated transport for next computation 1022 IF (ln_bt_fw) THEN 1023 zwx(:,:) = un_adv(:,:) 1024 zwy(:,:) = vn_adv(:,:) 1025 IF( .NOT.( kt == nit000 .AND. neuler==0 ) ) THEN 1026 un_adv(:,:) = z1_2 * ( ub2_b(:,:) + zwx(:,:) ) 1027 vn_adv(:,:) = z1_2 * ( vb2_b(:,:) + zwy(:,:) ) 1028 END IF 1029 ! Save integrated transport for next computation 1032 1030 ub2_b(:,:) = zwx(:,:) 1033 1031 vb2_b(:,:) = zwy(:,:) … … 1065 1063 DO jk = 1, jpkm1 1066 1064 ! Correct velocities: 1067 un(:,:,jk) = ( un(:,:,jk) + un_adv(:,:) - un_b(:,:) ) * umask(:,:,jk)1068 vn(:,:,jk) = ( vn(:,:,jk) + vn_adv(:,:) - vn_b(:,:) ) * vmask(:,:,jk)1065 un(:,:,jk) = ( un(:,:,jk) + un_adv(:,:)*r1_hu_n(:,:) - un_b(:,:) ) * umask(:,:,jk) 1066 vn(:,:,jk) = ( vn(:,:,jk) + vn_adv(:,:)*r1_hv_n(:,:) - vn_b(:,:) ) * vmask(:,:,jk) 1069 1067 ! 1070 1068 END DO 1071 1069 ! 1072 CALL iom_put( "ubar", un_adv(:,:) 1073 CALL iom_put( "vbar", vn_adv(:,:) 1070 CALL iom_put( "ubar", un_adv(:,:)*r1_hu_n(:,:) ) ! barotropic i-current 1071 CALL iom_put( "vbar", vn_adv(:,:)*r1_hv_n(:,:) ) ! barotropic i-current 1074 1072 ! 1075 1073 #if defined key_agrif … … 1297 1295 #if defined key_agrif 1298 1296 ! Restrict the use of Agrif to the forward case only 1299 IF( .NOT.ln_bt_fw .AND. .NOT.Agrif_Root() ) CALL ctl_stop( 'AGRIF not implemented if ln_bt_fw=.FALSE.' )1297 !!! IF( .NOT.ln_bt_fw .AND. .NOT.Agrif_Root() ) CALL ctl_stop( 'AGRIF not implemented if ln_bt_fw=.FALSE.' ) 1300 1298 #endif 1301 1299 ! -
branches/2017/dev_r8624_AGRIF3_VVL/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90
r7753 r8741 109 109 ! 110 110 ssha(:,:) = ( sshb(:,:) - z2dt * ( zcoef * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) ) ) * ssmask(:,:) 111 111 ! 112 #if defined key_agrif 113 CALL agrif_ssh( kt ) 114 #endif 115 ! 112 116 IF ( .NOT.ln_dynspg_ts ) THEN 113 ! These lines are not necessary with time splitting since114 ! boundary condition on sea level is set during ts loop115 # if defined key_agrif116 CALL agrif_ssh( kt )117 # endif118 117 IF( ln_bdy ) THEN 119 118 CALL lbc_lnk( ssha, 'T', 1. ) ! Not sure that's necessary … … 214 213 END DO 215 214 ENDIF 215 ! 216 #if defined key_agrif 217 IF( .NOT. AGRIF_Root() ) THEN 218 IF ((nbondi == 1).OR.(nbondi == 2)) wn(nlci-1 , : ,:) = 0.e0 ! east 219 IF ((nbondi == -1).OR.(nbondi == 2)) wn(2 , : ,:) = 0.e0 ! west 220 IF ((nbondj == 1).OR.(nbondj == 2)) wn(: ,nlcj-1 ,:) = 0.e0 ! north 221 IF ((nbondj == -1).OR.(nbondj == 2)) wn(: ,2 ,:) = 0.e0 ! south 222 ENDIF 223 #endif 216 224 ! 217 225 IF( nn_timing == 1 ) CALL timing_stop('wzv') -
branches/2017/dev_r8624_AGRIF3_VVL/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfgls.F90
r7646 r8741 33 33 USE iom ! I/O manager library 34 34 USE timing ! Timing 35 USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 35 USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 36 #if defined key_agrif 37 USE agrif_opa_interp ! Set bc on avm 38 #endif 36 39 37 40 IMPLICIT NONE … … 204 207 DO jj = 2, jpjm1 205 208 DO ji = fs_2, fs_jpim1 ! vector opt. 206 avmu(ji,jj,jk) = avmu(ji,jj,jk) * ( un(ji,jj,jk-1) - un(ji,jj,jk) ) & 207 & * ( ub(ji,jj,jk-1) - ub(ji,jj,jk) ) & 208 & / ( e3uw_n(ji,jj,jk) * e3uw_b(ji,jj,jk) ) 209 avmv(ji,jj,jk) = avmv(ji,jj,jk) * ( vn(ji,jj,jk-1) - vn(ji,jj,jk) ) & 210 & * ( vb(ji,jj,jk-1) - vb(ji,jj,jk) ) & 211 & / ( e3vw_n(ji,jj,jk) * e3vw_b(ji,jj,jk) ) 209 avmu(ji,jj,jk) = 0.5_wp * ( avm(ji,jj,jk) + avm(ji+1,jj,jk) ) & 210 & * ( un(ji,jj,jk-1) - un(ji,jj,jk) ) & 211 & * ( ub(ji,jj,jk-1) - ub(ji,jj,jk) ) & 212 & / ( e3uw_n(ji,jj,jk) * e3uw_b(ji,jj,jk) ) 213 avmv(ji,jj,jk) = 0.5_wp * ( avm(ji,jj,jk) + avm(ji,jj+1,jk) ) & 214 & * ( vn(ji,jj,jk-1) - vn(ji,jj,jk) ) & 215 & * ( vb(ji,jj,jk-1) - vb(ji,jj,jk) ) & 216 & / ( e3vw_n(ji,jj,jk) * e3vw_b(ji,jj,jk) ) 212 217 eps(ji,jj,jk) = rc03 * en(ji,jj,jk) * SQRT(en(ji,jj,jk)) / mxln(ji,jj,jk) 213 218 END DO … … 800 805 ! Lateral boundary conditions (sign unchanged) 801 806 avt(:,:,1) = 0._wp 807 ! 808 #if defined key_agrif 809 CALL Agrif_avm 810 #endif 811 ! 802 812 CALL lbc_lnk( avm, 'W', 1. ) ; CALL lbc_lnk( avt, 'W', 1. ) 803 813 -
branches/2017/dev_r8624_AGRIF3_VVL/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90
r7813 r8741 173 173 !!---------------------------------------------------------------------- 174 174 ! 175 #if defined key_agrif176 ! interpolation parent grid => child grid for avm_k ( ex : at west border: update column 1 and 2)177 IF( .NOT.Agrif_Root() ) CALL Agrif_Tke178 #endif179 !180 175 IF( kt /= nit000 ) THEN ! restore before value to compute tke 181 176 avt (:,:,:) = avt_k (:,:,:) … … 196 191 #if defined key_agrif 197 192 ! Update child grid f => parent grid 198 IF( .NOT.Agrif_Root() ) CALL Agrif_Update_Tke( kt ) ! children only 193 !!! JC: suppress update since restartability is not possible in that case 194 !!! IF( .NOT.Agrif_Root() ) CALL Agrif_Update_Tke( ) ! children only 199 195 #endif 200 196 ! … … 666 662 END DO 667 663 END DO 664 ! 665 #if defined key_agrif 666 CALL Agrif_avm 667 #endif 668 ! 668 669 CALL lbc_lnk( avm, 'W', 1. ) ! Lateral boundary conditions (sign unchanged) 669 670 ! -
branches/2017/dev_r8624_AGRIF3_VVL/NEMOGCM/NEMO/OPA_SRC/step.F90
r7753 r8741 315 315 IF( Agrif_NbStepint() == 0 ) THEN ! AGRIF Update 316 316 !!jc in fact update is useless at last time step, but do it for global diagnostics 317 CALL Agrif_Update_Tra() ! Update active tracers 318 CALL Agrif_Update_Dyn() ! Update momentum 317 IF(.NOT.ln_linssh) CALL Agrif_Update_vvl() ! Update vertical scale factors 318 CALL Agrif_Update_Tra() ! Update active tracers 319 CALL Agrif_Update_Dyn() ! Update momentum 320 #if defined key_top 321 CALL Agrif_Update_Trc() ! Update passive tracers 322 #endif 319 323 ENDIF 320 324 #endif -
branches/2017/dev_r8624_AGRIF3_VVL/NEMOGCM/NEMO/OPA_SRC/step_oce.F90
r7646 r8741 113 113 USE agrif_opa_sponge ! Momemtum and tracers sponges 114 114 USE agrif_opa_update ! Update (2-way nesting) 115 #if defined key_top 116 USE agrif_top_update ! passive tracers update (2-way nesting) 117 #endif 115 118 #endif 116 119 #if defined key_top -
branches/2017/dev_r8624_AGRIF3_VVL/NEMOGCM/NEMO/TOP_SRC/TRP/trctrp.F90
r7646 r8741 30 30 #if defined key_agrif 31 31 USE agrif_top_sponge ! tracers sponges 32 USE agrif_top_update ! tracers updates33 32 #endif 34 33 … … 83 82 IF( ln_trcdmp_clo ) CALL trc_dmp_clo( kt ) ! internal damping trends on closed seas only 84 83 85 #if defined key_agrif86 IF( .NOT.Agrif_Root()) CALL Agrif_Update_Trc( kt ) ! Update tracer at AGRIF zoom boundaries : children only87 #endif88 84 ! 89 85 ELSE ! 1D vertical configuration
Note: See TracChangeset
for help on using the changeset viewer.