- Timestamp:
- 2017-12-13T18:08:50+01:00 (6 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90
r9019 r9023 1 1 MODULE dynspg_ts 2 3 !! Includes ROMS wd scheme with diagnostic outputs ; un and ua updates are commented out ! 4 2 5 !!====================================================================== 3 6 !! *** MODULE dynspg_ts *** … … 57 60 USE restart ! only for lrst_oce 58 61 USE timing ! Timing 62 USE diatmb ! Top,middle,bottom output 63 #if defined key_agrif 64 USE agrif_opa_interp ! agrif 65 USE agrif_oce 66 #endif 67 #if defined key_asminc 68 USE asminc ! Assimilation increment 69 #endif 59 70 60 71 IMPLICIT NONE … … 68 79 !! Time filtered arrays at baroclinic time step: 69 80 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: un_adv , vn_adv !: Advection vel. at "now" barocl. step 70 71 INTEGER , SAVE :: icycle ! Number of barotropic sub-steps for each internal step nn_baro <= 2.5 nn_baro 72 REAL(wp), SAVE :: rdtbt ! Barotropic time step 81 INTEGER, SAVE :: icycle ! Number of barotropic sub-steps for each internal step nn_baro <= 2.5 nn_baro 82 REAL(wp),SAVE :: rdtbt ! Barotropic time step 73 83 ! 74 84 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:) :: wgtbtp1, wgtbtp2 ! 1st & 2nd weights used in time filtering of barotropic fields … … 131 141 !! -Update the filtered free surface at step "n+1" : ssha 132 142 !! -Update filtered barotropic velocities at step "n+1" : ua_b, va_b 133 !! -Compute barotropic advective velocities at step "n": un_adv, vn_adv143 !! -Compute barotropic advective fluxes at step "n" : un_adv, vn_adv 134 144 !! These are used to advect tracers and are compliant with discrete 135 145 !! continuity equation taken at the baroclinic time steps. This … … 159 169 REAL(wp), DIMENSION(jpi,jpj) :: zCdU_u, zCdU_v ! top/bottom stress at u- & v-points 160 170 ! 171 REAL(wp) :: zwdramp ! local scalar - only used if ln_wd_dl = .True. 172 173 INTEGER :: iwdg, jwdg, kwdg ! short-hand values for the indices of the output point 174 175 REAL(wp) :: zepsilon, zgamma ! - - 161 176 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zcpx, zcpy ! Wetting/Dying gravity filter coef. 177 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: ztwdmask, zuwdmask, zvwdmask ! ROMS wetting and drying masks at t,u,v points 178 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zuwdav2, zvwdav2 ! averages over the sub-steps of zuwdmask and zvwdmask 162 179 !!---------------------------------------------------------------------- 163 180 ! 164 181 IF( ln_timing ) CALL timing_start('dyn_spg_ts') 165 182 ! 166 IF( ln_wd ) ALLOCATE( zcpx(jpi,jpj), zcpy(jpi,jpj) ) 183 IF( ln_wd_il ) ALLOCATE( zcpx(jpi,jpj), zcpy(jpi,jpj) ) 184 ! !* Allocate temporary arrays 185 IF( ln_wd_dl ) ALLOCATE( ztwdmask(jpi,jpj), zuwdmask(jpi,jpj), zvwdmask(jpi,jpj), zuwdav2(jpi,jpj), zvwdav2(jpi,jpj)) 167 186 ! 168 187 zmdi=1.e+20 ! missing data indicator for masking 169 188 ! 170 ! ! reciprocal of baroclinic time step 189 zwdramp = r_rn_wdmin1 ! simplest ramp 190 ! zwdramp = 1._wp / (rn_wdmin2 - rn_wdmin1) ! more general ramp 191 ! ! reciprocal of baroclinic time step 171 192 IF( kt == nit000 .AND. neuler == 0 ) THEN ; z2dt_bf = rdt 172 193 ELSE ; z2dt_bf = 2.0_wp * rdt … … 174 195 z1_2dt_b = 1.0_wp / z2dt_bf 175 196 ! 176 ll_init = ln_bt_av 197 ll_init = ln_bt_av ! if no time averaging, then no specific restart 177 198 ll_fw_start = .FALSE. 178 ! 199 ! ! time offset in steps for bdy data update 179 200 IF( .NOT.ln_bt_fw ) THEN ; noffset = - nn_baro 180 201 ELSE ; noffset = 0 181 202 ENDIF 182 203 ! 183 IF( kt == nit000 ) THEN 204 IF( kt == nit000 ) THEN !* initialisation 184 205 ! 185 206 IF(lwp) WRITE(numout,*) … … 405 426 ! ! ---------------------------------------------------- 406 427 IF( .NOT.ln_linssh ) THEN ! Variable volume : remove surface pressure gradient 407 IF( ln_wd ) THEN ! Calculating and applying W/D gravity filters 408 DO jj = 2, jpjm1 409 DO ji = 2, jpim1 410 ll_tmp1 = MIN( sshn(ji,jj) , sshn(ji+1,jj) ) > & 411 & MAX( -ht_wd(ji,jj) , -ht_wd(ji+1,jj) ) .AND. & 412 & MAX( sshn(ji,jj) + ht_wd(ji,jj), sshn(ji+1,jj) + ht_wd(ji+1,jj) ) & 413 & > rn_wdmin1 + rn_wdmin2 414 ll_tmp2 = ( ABS( sshn(ji+1,jj) - sshn(ji ,jj)) > 1.E-12 ).AND.( & 415 & MAX( sshn(ji,jj) , sshn(ji+1,jj) ) > & 416 & MAX( -ht_wd(ji,jj) , -ht_wd(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 ) 417 ! 418 IF(ll_tmp1) THEN 419 zcpx(ji,jj) = 1.0_wp 420 ELSE IF(ll_tmp2) THEN ! no worries about sshn(ji+1,jj) - sshn(ji ,jj) = 0, it won't happen ! here 421 zcpx(ji,jj) = ABS( (sshn(ji+1,jj) + ht_wd(ji+1,jj) - sshn(ji,jj) - ht_wd(ji,jj)) & 422 & / (sshn(ji+1,jj) - sshn(ji ,jj)) ) 423 ELSE 424 zcpx(ji,jj) = 0._wp 425 ENDIF 426 ! 427 ll_tmp1 = MIN( sshn(ji,jj) , sshn(ji,jj+1) ) > & 428 & MAX( -ht_wd(ji,jj) , -ht_wd(ji,jj+1) ) .AND. & 429 & MAX( sshn(ji,jj) + ht_wd(ji,jj), sshn(ji,jj+1) + ht_wd(ji,jj+1) ) & 430 & > rn_wdmin1 + rn_wdmin2 431 ll_tmp2 = ( ABS( sshn(ji,jj) - sshn(ji,jj+1)) > 1.E-12 ).AND.( & 432 & MAX( sshn(ji,jj) , sshn(ji,jj+1) ) > & 433 & MAX( -ht_wd(ji,jj) , -ht_wd(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 ) 434 ! 435 IF(ll_tmp1) THEN 436 zcpy(ji,jj) = 1.0_wp 437 ELSE IF(ll_tmp2) THEN 438 ! no worries about sshn(ji,jj+1) - sshn(ji,jj ) = 0, it won't happen ! here 439 zcpy(ji,jj) = ABS( (sshn(ji,jj+1) + ht_wd(ji,jj+1) - sshn(ji,jj) - ht_wd(ji,jj)) & 440 & / (sshn(ji,jj+1) - sshn(ji,jj )) ) 441 ELSE 442 zcpy(ji,jj) = 0._wp 443 ENDIF 444 END DO 445 END DO 446 ! 447 DO jj = 2, jpjm1 448 DO ji = 2, jpim1 449 zu_trd(ji,jj) = zu_trd(ji,jj) - grav * ( sshn(ji+1,jj ) - sshn(ji ,jj ) ) & 450 & * r1_e1u(ji,jj) * zcpx(ji,jj) 451 zv_trd(ji,jj) = zv_trd(ji,jj) - grav * ( sshn(ji ,jj+1) - sshn(ji ,jj ) ) & 452 & * r1_e2v(ji,jj) * zcpy(ji,jj) 453 END DO 454 END DO 455 ! 428 IF( ln_wd_il ) THEN ! Calculating and applying W/D gravity filters 429 DO jj = 2, jpjm1 430 DO ji = 2, jpim1 431 ll_tmp1 = MIN( sshn(ji,jj) , sshn(ji+1,jj) ) > & 432 & MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) .AND. & 433 & MAX( sshn(ji,jj) + ht_0(ji,jj) , sshn(ji+1,jj) + ht_0(ji+1,jj) ) & 434 & > rn_wdmin1 + rn_wdmin2 435 ll_tmp2 = ( ABS( sshn(ji+1,jj) - sshn(ji ,jj)) > 1.E-12 ).AND.( & 436 & MAX( sshn(ji,jj) , sshn(ji+1,jj) ) > & 437 & MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 ) 438 439 IF(ll_tmp1) THEN 440 zcpx(ji,jj) = 1.0_wp 441 ELSE IF(ll_tmp2) THEN 442 ! no worries about sshn(ji+1,jj) - sshn(ji ,jj) = 0, it won't happen ! here 443 zcpx(ji,jj) = ABS( (sshn(ji+1,jj) + ht_0(ji+1,jj) - sshn(ji,jj) - ht_0(ji,jj)) & 444 & / (sshn(ji+1,jj) - sshn(ji ,jj)) ) 445 zcpx(ji,jj) = max(min( zcpx(ji,jj) , 1.0_wp),0.0_wp) 446 447 ELSE 448 zcpx(ji,jj) = 0._wp 449 END IF 450 451 ll_tmp1 = MIN( sshn(ji,jj) , sshn(ji,jj+1) ) > & 452 & MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) .AND. & 453 & MAX( sshn(ji,jj) + ht_0(ji,jj) , sshn(ji,jj+1) + ht_0(ji,jj+1) ) & 454 & > rn_wdmin1 + rn_wdmin2 455 ll_tmp2 = ( ABS( sshn(ji,jj) - sshn(ji,jj+1)) > 1.E-12 ).AND.( & 456 & MAX( sshn(ji,jj) , sshn(ji,jj+1) ) > & 457 & MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 ) 458 459 IF(ll_tmp1) THEN 460 zcpy(ji,jj) = 1.0_wp 461 ELSE IF(ll_tmp2) THEN 462 ! no worries about sshn(ji,jj+1) - sshn(ji,jj ) = 0, it won't happen ! here 463 zcpy(ji,jj) = ABS( (sshn(ji,jj+1) + ht_0(ji,jj+1) - sshn(ji,jj) - ht_0(ji,jj)) & 464 & / (sshn(ji,jj+1) - sshn(ji,jj )) ) 465 zcpy(ji,jj) = max(min( zcpy(ji,jj) , 1.0_wp),0.0_wp) 466 467 ELSE 468 zcpy(ji,jj) = 0._wp 469 END IF 470 END DO 471 END DO 472 473 DO jj = 2, jpjm1 474 DO ji = 2, jpim1 475 zu_trd(ji,jj) = zu_trd(ji,jj) - grav * ( sshn(ji+1,jj ) - sshn(ji ,jj ) ) & 476 & * r1_e1u(ji,jj) * zcpx(ji,jj) * wdrampu(ji,jj) !jth 477 zv_trd(ji,jj) = zv_trd(ji,jj) - grav * ( sshn(ji ,jj+1) - sshn(ji ,jj ) ) & 478 & * r1_e2v(ji,jj) * zcpy(ji,jj) * wdrampv(ji,jj) !jth 479 480 END DO 481 END DO 482 456 483 ELSE 457 484 ! … … 473 500 END DO 474 501 ! 475 ! ! Add BOTTOMstress contribution from baroclinic velocities:476 IF ( ln_bt_fw) THEN502 ! ! Add bottom stress contribution from baroclinic velocities: 503 IF (ln_bt_fw) THEN 477 504 DO jj = 2, jpjm1 478 505 DO ji = fs_2, fs_jpim1 ! vector opt. … … 495 522 ! 496 523 ! Note that the "unclipped" bottom friction parameter is used even with explicit drag 497 IF( ln_wd ) THEN 498 zztmp = - 1._wp / rdtbt 499 DO jj = 2, jpjm1 524 IF( ln_wd_il ) THEN 525 zu_frc(:,:) = zu_frc(:,:) + MAX(r1_hu_n(:,:) * bfrua(:,:),-1._wp / rdtbt) * zwx(:,:) * wdrampu(ji,jj) 526 zv_frc(:,:) = zv_frc(:,:) + MAX(r1_hv_n(:,:) * bfrva(:,:),-1._wp / rdtbt) * zwy(:,:) * wdrampv(ji,jj) 527 ELSE 528 zu_frc(:,:) = zu_frc(:,:) + r1_hu_n(:,:) * bfrua(:,:) * zwx(:,:) 529 zv_frc(:,:) = zv_frc(:,:) + r1_hv_n(:,:) * bfrva(:,:) * zwy(:,:) 530 END IF 531 ! 532 ! ! Add top stress contribution from baroclinic velocities: 533 IF( ln_bt_fw ) THEN 534 DO jj = 2, jpjm1 500 535 DO ji = fs_2, fs_jpim1 ! vector opt. 501 536 zu_frc(ji,jj) = zu_frc(ji,jj) + MAX( r1_hu_n(ji,jj) * 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) , zztmp ) * zwx(ji,jj) … … 657 692 vn_adv(:,:) = 0._wp 658 693 ! ! ==================== ! 694 695 IF (ln_wd_dl) THEN 696 zuwdmask(:,:) = 0._wp ! set to zero for definiteness (not sure this is necessary) 697 zvwdmask(:,:) = 0._wp ! 698 zuwdav2(:,:) = 0._wp 699 zvwdav2(:,:) = 0._wp 700 END IF 701 702 659 703 DO jn = 1, icycle ! sub-time-step loop ! 660 704 ! ! ==================== ! … … 684 728 ! Extrapolate Sea Level at step jit+0.5: 685 729 zsshp2_e(:,:) = za1 * sshn_e(:,:) + za2 * sshb_e(:,:) + za3 * sshbb_e(:,:) 686 ! 730 731 ! set wetting & drying mask at tracer points for this barotropic sub-step 732 IF ( ln_wd_dl ) THEN 733 734 IF ( ln_wd_dl_rmp ) THEN 735 DO jj = 1, jpj 736 DO ji = 1, jpi ! vector opt. 737 IF ( zsshp2_e(ji,jj) + ht_0(ji,jj) > 2._wp * rn_wdmin1 ) THEN 738 ! IF ( zsshp2_e(ji,jj) + ht_0(ji,jj) > rn_wdmin2 ) THEN 739 ztwdmask(ji,jj) = 1._wp 740 ELSE IF ( zsshp2_e(ji,jj) + ht_0(ji,jj) > rn_wdmin1 ) THEN 741 ztwdmask(ji,jj) = (tanh(50._wp*( ( zsshp2_e(ji,jj) + ht_0(ji,jj) - rn_wdmin1 )*r_rn_wdmin1)) ) 742 ELSE 743 ztwdmask(ji,jj) = 0._wp 744 END IF 745 END DO 746 END DO 747 ELSE 748 DO jj = 1, jpj 749 DO ji = 1, jpi ! vector opt. 750 IF ( zsshp2_e(ji,jj) + ht_0(ji,jj) > rn_wdmin1 ) THEN 751 ztwdmask(ji,jj) = 1._wp 752 ELSE 753 ztwdmask(ji,jj) = 0._wp 754 END IF 755 END DO 756 END DO 757 END IF 758 759 END IF 760 761 687 762 DO jj = 2, jpjm1 ! Sea Surface Height at u- & v-points 688 763 DO ji = 2, fs_jpim1 ! Vector opt. … … 736 811 ENDIF 737 812 #endif 738 IF( ln_wd ) CALL wad_lmt_bt(zwx, zwy, sshn_e, zssh_frc, rdtbt) 739 ! 813 IF( ln_wd_il ) CALL wad_lmt_bt(zwx, zwy, sshn_e, zssh_frc, rdtbt) 814 815 IF ( ln_wd_dl ) THEN 816 817 818 ! un_e and vn_e are set to zero at faces where the direction of the flow is from dry cells 819 820 DO jj = 1, jpjm1 821 DO ji = 1, jpim1 822 IF ( zwx(ji,jj) > 0.0 ) THEN 823 zuwdmask(ji, jj) = ztwdmask(ji ,jj) 824 ELSE 825 zuwdmask(ji, jj) = ztwdmask(ji+1,jj) 826 END IF 827 zwx(ji, jj) = zuwdmask(ji,jj)*zwx(ji, jj) 828 un_e(ji,jj) = zuwdmask(ji,jj)*un_e(ji,jj) 829 830 IF ( zwy(ji,jj) > 0.0 ) THEN 831 zvwdmask(ji, jj) = ztwdmask(ji, jj ) 832 ELSE 833 zvwdmask(ji, jj) = ztwdmask(ji, jj+1) 834 END IF 835 zwy(ji, jj) = zvwdmask(ji,jj)*zwy(ji,jj) 836 vn_e(ji,jj) = zvwdmask(ji,jj)*vn_e(ji,jj) 837 END DO 838 END DO 839 840 841 END IF 842 740 843 ! Sum over sub-time-steps to compute advective velocities 741 844 za2 = wgtbtp2(jn) 742 845 un_adv(:,:) = un_adv(:,:) + za2 * zwx(:,:) * r1_e2u(:,:) 743 846 vn_adv(:,:) = vn_adv(:,:) + za2 * zwy(:,:) * r1_e1v(:,:) 744 ! 847 848 ! sum over sub-time-steps to decide which baroclinic velocities to set to zero (zuwdav2 is only used when ln_wd_dl_bc = True) 849 IF ( ln_wd_dl_bc ) THEN 850 zuwdav2(:,:) = zuwdav2(:,:) + za2 * zuwdmask(:,:) 851 zvwdav2(:,:) = zvwdav2(:,:) + za2 * zvwdmask(:,:) 852 END IF 853 745 854 ! Set next sea level: 746 855 DO jj = 2, jpjm1 … … 788 897 za3= 0._wp 789 898 ELSE ! AB3-AM4 Coefficients; bet=0.281105 ; eps=0.013 ; gam=0.0880 790 za0=0.614_wp ! za0 = 1/2 + gam + 2*eps 791 za1=0.285_wp ! za1 = 1/2 - 2*gam - 3*eps 792 za2=0.088_wp ! za2 = gam 793 za3=0.013_wp ! za3 = eps 899 IF (rn_bt_alpha==0._wp) THEN 900 za0=0.614_wp ! za0 = 1/2 + gam + 2*eps 901 za1=0.285_wp ! za1 = 1/2 - 2*gam - 3*eps 902 za2=0.088_wp ! za2 = gam 903 za3=0.013_wp ! za3 = eps 904 ELSE 905 zepsilon = 0.00976186_wp - 0.13451357_wp * rn_bt_alpha 906 zgamma = 0.08344500_wp - 0.51358400_wp * rn_bt_alpha 907 za0 = 0.5_wp + zgamma + 2._wp * rn_bt_alpha + 2._wp * zepsilon 908 za1 = 1._wp - za0 - zgamma - zepsilon 909 za2 = zgamma 910 za3 = zepsilon 911 ENDIF 794 912 ENDIF 795 913 ! 796 914 zsshp2_e(:,:) = za0 * ssha_e(:,:) + za1 * sshn_e (:,:) & 797 915 & + za2 * sshb_e(:,:) + za3 * sshbb_e(:,:) 798 IF( ln_wd ) THEN ! Calculating and applying W/D gravity filters916 IF( ln_wd_il ) THEN ! Calculating and applying W/D gravity filters 799 917 DO jj = 2, jpjm1 800 918 DO ji = 2, jpim1 801 919 ll_tmp1 = MIN( zsshp2_e(ji,jj) , zsshp2_e(ji+1,jj) ) > & 802 & MAX( -ht_wd(ji,jj) , -ht_wd(ji+1,jj) ) .AND. &803 & MAX( zsshp2_e(ji,jj) + ht_ wd(ji,jj), zsshp2_e(ji+1,jj) + ht_wd(ji+1,jj) ) &920 & MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) .AND. & 921 & MAX( zsshp2_e(ji,jj) + ht_0(ji,jj) , zsshp2_e(ji+1,jj) + ht_0(ji+1,jj) ) & 804 922 & > rn_wdmin1 + rn_wdmin2 805 923 ll_tmp2 = (ABS(zsshp2_e(ji,jj) - zsshp2_e(ji+1,jj)) > 1.E-12 ).AND.( & 806 924 & MAX( zsshp2_e(ji,jj) , zsshp2_e(ji+1,jj) ) > & 807 & MAX( -ht_wd(ji,jj) , -ht_wd(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 )925 & MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 ) 808 926 809 927 IF(ll_tmp1) THEN … … 811 929 ELSE IF(ll_tmp2) THEN 812 930 ! no worries about zsshp2_e(ji+1,jj) - zsshp2_e(ji ,jj) = 0, it won't happen ! here 813 zcpx(ji,jj) = ABS( (zsshp2_e(ji+1,jj) + ht_wd(ji+1,jj) - zsshp2_e(ji,jj) - ht_wd(ji,jj)) &931 zcpx(ji,jj) = ABS( (zsshp2_e(ji+1,jj) + ht_0(ji+1,jj) - zsshp2_e(ji,jj) - ht_0(ji,jj)) & 814 932 & / (zsshp2_e(ji+1,jj) - zsshp2_e(ji ,jj)) ) 815 933 ELSE … … 818 936 819 937 ll_tmp1 = MIN( zsshp2_e(ji,jj) , zsshp2_e(ji,jj+1) ) > & 820 & MAX( -ht_wd(ji,jj) , -ht_wd(ji,jj+1) ) .AND. &821 & MAX( zsshp2_e(ji,jj) + ht_ wd(ji,jj), zsshp2_e(ji,jj+1) + ht_wd(ji,jj+1) ) &938 & MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) .AND. & 939 & MAX( zsshp2_e(ji,jj) + ht_0(ji,jj) , zsshp2_e(ji,jj+1) + ht_0(ji,jj+1) ) & 822 940 & > rn_wdmin1 + rn_wdmin2 823 941 ll_tmp2 = (ABS(zsshp2_e(ji,jj) - zsshp2_e(ji,jj+1)) > 1.E-12 ).AND.( & 824 942 & MAX( zsshp2_e(ji,jj) , zsshp2_e(ji,jj+1) ) > & 825 & MAX( -ht_wd(ji,jj) , -ht_wd(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 )943 & MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 ) 826 944 827 945 IF(ll_tmp1) THEN … … 829 947 ELSE IF(ll_tmp2) THEN 830 948 ! no worries about zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj ) = 0, it won't happen ! here 831 zcpy(ji,jj) = ABS( (zsshp2_e(ji,jj+1) + ht_wd(ji,jj+1) - zsshp2_e(ji,jj) - ht_wd(ji,jj)) &949 zcpy(ji,jj) = ABS( (zsshp2_e(ji,jj+1) + ht_0(ji,jj+1) - zsshp2_e(ji,jj) - ht_0(ji,jj)) & 832 950 & / (zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj )) ) 833 951 ELSE … … 839 957 ! 840 958 ! Compute associated depths at U and V points: 841 IF( .NOT.ln_linssh .AND. .NOT.ln_dynadv_vec ) THEN 959 IF( .NOT.ln_linssh .AND. .NOT.ln_dynadv_vec ) THEN !* Vector form 842 960 ! 843 961 DO jj = 2, jpjm1 … … 915 1033 ENDIF 916 1034 ! 917 DO jj = 2, jpjm1 ! Add top/bottom stresses:918 DO ji = fs_2, fs_jpim1 ! vector opt. 919 zu_trd(ji,jj) = zu_trd(ji,jj) + zCdU_u(ji,jj) * un_e(ji,jj) * hur_e(ji,jj)920 zv_trd(ji,jj) = zv_trd(ji,jj) + zCdU_v(ji,jj) * vn_e(ji,jj) * hvr_e(ji,jj)921 END DO922 END DO1035 ! Add bottom stresses: 1036 !jth do implicitly instead 1037 IF ( .NOT. ll_wd ) THEN ! Revert to explicit for bit comparison tests in non wad runs 1038 zu_trd(:,:) = zu_trd(:,:) + bfrua(:,:) * un_e(:,:) * hur_e(:,:) 1039 zv_trd(:,:) = zv_trd(:,:) + bfrva(:,:) * vn_e(:,:) * hvr_e(:,:) 1040 ENDIF 923 1041 ! 924 1042 ! Surface pressure trend: 925 IF( ln_wd ) THEN1043 IF( ln_wd_il ) THEN 926 1044 DO jj = 2, jpjm1 927 1045 DO ji = 2, jpim1 … … 929 1047 zu_spg = - grav * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) * r1_e1u(ji,jj) 930 1048 zv_spg = - grav * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) * r1_e2v(ji,jj) 931 zwx(ji,jj) = zu_spg * zcpx(ji,jj)932 zwy(ji,jj) = zv_spg * zcpy(ji,jj)1049 zwx(ji,jj) = (1._wp - rn_scal_load) * zu_spg * zcpx(ji,jj) 1050 zwy(ji,jj) = (1._wp - rn_scal_load) * zv_spg * zcpy(ji,jj) 933 1051 END DO 934 1052 END DO … … 939 1057 zu_spg = - grav * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) * r1_e1u(ji,jj) 940 1058 zv_spg = - grav * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) * r1_e2v(ji,jj) 941 zwx(ji,jj) = zu_spg942 zwy(ji,jj) = zv_spg1059 zwx(ji,jj) = (1._wp - rn_scal_load) * zu_spg 1060 zwy(ji,jj) = (1._wp - rn_scal_load) * zv_spg 943 1061 END DO 944 1062 END DO … … 947 1065 ! 948 1066 ! Set next velocities: 949 IF( ln_dynadv_vec .OR. ln_linssh ) THEN 1067 IF( ln_dynadv_vec .OR. ln_linssh ) THEN !* Vector form 950 1068 DO jj = 2, jpjm1 951 1069 DO ji = fs_2, fs_jpim1 ! vector opt. … … 961 1079 & + zv_frc(ji,jj) ) & 962 1080 & ) * ssvmask(ji,jj) 1081 1082 !jth implicit bottom friction: 1083 IF ( ll_wd ) THEN ! revert to explicit for bit comparison tests in non wad runs 1084 ua_e(ji,jj) = ua_e(ji,jj) /(1.0 - rdtbt * bfrua(ji,jj) * hur_e(ji,jj)) 1085 va_e(ji,jj) = va_e(ji,jj) /(1.0 - rdtbt * bfrva(ji,jj) * hvr_e(ji,jj)) 1086 ENDIF 1087 963 1088 END DO 964 1089 END DO 965 1090 ! 966 ELSE 1091 ELSE !* Flux form 967 1092 DO jj = 2, jpjm1 968 1093 DO ji = fs_2, fs_jpim1 ! vector opt. 969 1094 970 IF( ln_wd ) THEN 971 zhura = MAX(hu_0(ji,jj) + zsshu_a(ji,jj), rn_wdmin1) 972 zhvra = MAX(hv_0(ji,jj) + zsshv_a(ji,jj), rn_wdmin1) 973 ELSE 974 zhura = hu_0(ji,jj) + zsshu_a(ji,jj) 975 zhvra = hv_0(ji,jj) + zsshv_a(ji,jj) 976 END IF 1095 zhura = hu_0(ji,jj) + zsshu_a(ji,jj) 1096 zhvra = hv_0(ji,jj) + zsshv_a(ji,jj) 1097 977 1098 zhura = ssumask(ji,jj)/(zhura + 1._wp - ssumask(ji,jj)) 978 1099 zhvra = ssvmask(ji,jj)/(zhvra + 1._wp - ssvmask(ji,jj)) … … 992 1113 END DO 993 1114 ENDIF 994 ! 1115 1116 995 1117 IF( .NOT.ln_linssh ) THEN !* Update ocean depth (variable volume case only) 996 IF( ln_wd ) THEN 997 hu_e (:,:) = MAX(hu_0(:,:) + zsshu_a(:,:), rn_wdmin1) 998 hv_e (:,:) = MAX(hv_0(:,:) + zsshv_a(:,:), rn_wdmin1) 999 ELSE 1000 hu_e (:,:) = hu_0(:,:) + zsshu_a(:,:) 1001 hv_e (:,:) = hv_0(:,:) + zsshv_a(:,:) 1002 END IF 1118 hu_e (:,:) = hu_0(:,:) + zsshu_a(:,:) 1119 hv_e (:,:) = hv_0(:,:) + zsshv_a(:,:) 1003 1120 hur_e(:,:) = ssumask(:,:) / ( hu_e(:,:) + 1._wp - ssumask(:,:) ) 1004 1121 hvr_e(:,:) = ssvmask(:,:) / ( hv_e(:,:) + 1._wp - ssvmask(:,:) ) … … 1033 1150 ua_b (:,:) = ua_b (:,:) + za1 * ua_e (:,:) 1034 1151 va_b (:,:) = va_b (:,:) + za1 * va_e (:,:) 1035 ELSE ! Sum transports 1036 ua_b (:,:) = ua_b (:,:) + za1 * ua_e (:,:) * hu_e (:,:) 1037 va_b (:,:) = va_b (:,:) + za1 * va_e (:,:) * hv_e (:,:) 1038 ENDIF 1039 ! ! Sum sea level 1152 ELSE ! Sum transports 1153 IF ( .NOT.ln_wd_dl ) THEN 1154 ua_b (:,:) = ua_b (:,:) + za1 * ua_e (:,:) * hu_e (:,:) 1155 va_b (:,:) = va_b (:,:) + za1 * va_e (:,:) * hv_e (:,:) 1156 ELSE 1157 ua_b (:,:) = ua_b (:,:) + za1 * ua_e (:,:) * hu_e (:,:) * zuwdmask(:,:) 1158 va_b (:,:) = va_b (:,:) + za1 * va_e (:,:) * hv_e (:,:) * zvwdmask(:,:) 1159 END IF 1160 ENDIF 1161 ! ! Sum sea level 1040 1162 ssha(:,:) = ssha(:,:) + za1 * ssha_e(:,:) 1163 1041 1164 ! ! ==================== ! 1042 1165 END DO ! end loop ! … … 1047 1170 ! 1048 1171 ! Set advection velocity correction: 1049 zwx(:,:) = un_adv(:,:) 1050 zwy(:,:) = vn_adv(:,:) 1051 IF( ( kt == nit000 .AND. neuler==0 ) .OR. .NOT.ln_bt_fw ) THEN 1052 un_adv(:,:) = zwx(:,:) * r1_hu_n(:,:) 1053 vn_adv(:,:) = zwy(:,:) * r1_hv_n(:,:) 1054 ELSE 1055 un_adv(:,:) = r1_2 * ( ub2_b(:,:) + zwx(:,:) ) * r1_hu_n(:,:) 1056 vn_adv(:,:) = r1_2 * ( vb2_b(:,:) + zwy(:,:) ) * r1_hv_n(:,:) 1057 END IF 1058 1059 IF( ln_bt_fw ) THEN ! Save integrated transport for next computation 1172 IF (ln_bt_fw) THEN 1173 zwx(:,:) = un_adv(:,:) 1174 zwy(:,:) = vn_adv(:,:) 1175 IF( .NOT.( kt == nit000 .AND. neuler==0 ) ) THEN 1176 un_adv(:,:) = z1_2 * ( ub2_b(:,:) + zwx(:,:) - atfp * un_bf(:,:) ) 1177 vn_adv(:,:) = z1_2 * ( vb2_b(:,:) + zwy(:,:) - atfp * vn_bf(:,:) ) 1178 ! 1179 ! Update corrective fluxes for next time step: 1180 un_bf(:,:) = atfp * un_bf(:,:) + (zwx(:,:) - ub2_b(:,:)) 1181 vn_bf(:,:) = atfp * vn_bf(:,:) + (zwy(:,:) - vb2_b(:,:)) 1182 ELSE 1183 un_bf(:,:) = 0._wp 1184 vn_bf(:,:) = 0._wp 1185 END IF 1186 ! Save integrated transport for next computation 1060 1187 ub2_b(:,:) = zwx(:,:) 1061 1188 vb2_b(:,:) = zwy(:,:) 1062 1189 ENDIF 1190 1191 1063 1192 ! 1064 1193 ! Update barotropic trend: … … 1090 1219 va_b(:,:) = va_b(:,:) / ( hv_0(:,:) + zsshv_a(:,:) + 1._wp - ssvmask(:,:) ) 1091 1220 ENDIF 1092 ! 1221 1222 1223 ! Correct velocities so that the barotropic velocity equals (un_adv, vn_adv) (in all cases) 1093 1224 DO jk = 1, jpkm1 1094 ! Correct velocities: 1095 un(:,:,jk) = ( un(:,:,jk) + un_adv(:,:) - un_b(:,:) ) * umask(:,:,jk) 1096 vn(:,:,jk) = ( vn(:,:,jk) + vn_adv(:,:) - vn_b(:,:) ) * vmask(:,:,jk) 1097 ! 1225 un(:,:,jk) = ( un(:,:,jk) + un_adv(:,:)*r1_hu_n(:,:) - un_b(:,:) ) * umask(:,:,jk) 1226 vn(:,:,jk) = ( vn(:,:,jk) + vn_adv(:,:)*r1_hv_n(:,:) - vn_b(:,:) ) * vmask(:,:,jk) 1098 1227 END DO 1099 ! 1100 CALL iom_put( "ubar", un_adv(:,:) ) ! barotropic i-current 1101 CALL iom_put( "vbar", vn_adv(:,:) ) ! barotropic i-current 1228 1229 IF ( ln_wd_dl .and. ln_wd_dl_bc) THEN 1230 DO jk = 1, jpkm1 1231 un(:,:,jk) = ( un_adv(:,:) + zuwdav2(:,:)*(un(:,:,jk) - un_adv(:,:)) ) * umask(:,:,jk) 1232 vn(:,:,jk) = ( vn_adv(:,:) + zvwdav2(:,:)*(vn(:,:,jk) - vn_adv(:,:)) ) * vmask(:,:,jk) 1233 END DO 1234 END IF 1235 1236 1237 CALL iom_put( "ubar", un_adv(:,:)*r1_hu_n(:,:) ) ! barotropic i-current 1238 CALL iom_put( "vbar", vn_adv(:,:)*r1_hv_n(:,:) ) ! barotropic i-current 1102 1239 ! 1103 1240 #if defined key_agrif … … 1119 1256 IF( lrst_oce .AND.ln_bt_fw ) CALL ts_rst( kt, 'WRITE' ) 1120 1257 ! 1121 IF( ln_wd ) DEALLOCATE( zcpx, zcpy ) 1258 IF( ln_wd_il ) DEALLOCATE( zcpx, zcpy ) 1259 IF( ln_wd_dl ) DEALLOCATE( ztwdmask, zuwdmask, zvwdmask, zuwdav2, zvwdav2 ) 1122 1260 ! 1123 1261 IF( ln_diatmb ) THEN … … 1222 1360 CALL iom_get( numror, jpdom_autoglo, 'ub2_b' , ub2_b (:,:) ) 1223 1361 CALL iom_get( numror, jpdom_autoglo, 'vb2_b' , vb2_b (:,:) ) 1362 CALL iom_get( numror, jpdom_autoglo, 'un_bf' , un_bf (:,:) ) 1363 CALL iom_get( numror, jpdom_autoglo, 'vn_bf' , vn_bf (:,:) ) 1224 1364 IF( .NOT.ln_bt_av ) THEN 1225 1365 CALL iom_get( numror, jpdom_autoglo, 'sshbb_e' , sshbb_e(:,:) ) … … 1241 1381 CALL iom_rstput( kt, nitrst, numrow, 'ub2_b' , ub2_b (:,:) ) 1242 1382 CALL iom_rstput( kt, nitrst, numrow, 'vb2_b' , vb2_b (:,:) ) 1383 CALL iom_rstput( kt, nitrst, numrow, 'un_bf' , un_bf (:,:) ) 1384 CALL iom_rstput( kt, nitrst, numrow, 'vn_bf' , vn_bf (:,:) ) 1243 1385 ! 1244 1386 IF (.NOT.ln_bt_av) THEN … … 1291 1433 rdtbt = rdt / REAL( nn_baro , wp ) 1292 1434 zcmax = zcmax * rdtbt 1293 1435 ! Print results 1294 1436 IF(lwp) WRITE(numout,*) 1295 1437 IF(lwp) WRITE(numout,*) 'dyn_spg_ts : split-explicit free surface' … … 1317 1459 #if defined key_agrif 1318 1460 ! Restrict the use of Agrif to the forward case only 1319 IF( .NOT.ln_bt_fw .AND. .NOT.Agrif_Root() ) CALL ctl_stop( 'AGRIF not implemented if ln_bt_fw=.FALSE.' )1461 !!! IF( .NOT.ln_bt_fw .AND. .NOT.Agrif_Root() ) CALL ctl_stop( 'AGRIF not implemented if ln_bt_fw=.FALSE.' ) 1320 1462 #endif 1321 1463 ! … … 1333 1475 IF(lwp) WRITE(numout,*) ' Maximum Courant number is :', zcmax 1334 1476 ! 1477 IF(lwp) WRITE(numout,*) ' Time diffusion parameter rn_bt_alpha: ', rn_bt_alpha 1478 IF ((ln_bt_av.AND.nn_bt_flt/=0).AND.(rn_bt_alpha>0._wp)) THEN 1479 CALL ctl_stop( 'dynspg_ts ERROR: if rn_bt_alpha > 0, remove temporal averaging' ) 1480 ENDIF 1481 ! 1335 1482 IF( .NOT.ln_bt_av .AND. .NOT.ln_bt_fw ) THEN 1336 1483 CALL ctl_stop( 'dynspg_ts ERROR: No time averaging => only forward integration is possible' )
Note: See TracChangeset
for help on using the changeset viewer.
