New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Changeset 7016 for branches/2016/dev_r6393_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90 – NEMO

Ignore:
Timestamp:
2016-10-10T18:24:04+02:00 (8 years ago)
Author:
acc
Message:

Branch 2016/dev_r6393_NOC_WAD. Reinstated and fixed limiter in barotropic loop to resolve conservation issues with wetting and drying. Now maintains constant salinity in WAD_TEST_CASES 1,2 and 3.

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/2016/dev_r6393_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90

    r6986 r7016  
    156156      REAL(wp), POINTER, DIMENSION(:,:) :: zhf 
    157157      REAL(wp), POINTER, DIMENSION(:,:) :: zcpx, zcpy                 ! Wetting/Dying gravity filter coef. 
    158       REAL(wp), POINTER, DIMENSION(:,:) :: wduflt1, wdvflt1           ! Wetting/Dying velocity filter coef. 
    159       REAL(wp), POINTER, DIMENSION(:,:) :: sshai                      ! Wetting/Dying velocity filter coef. 
    160158      !!---------------------------------------------------------------------- 
    161159      ! 
     
    169167      CALL wrk_alloc( jpi,jpj,   zsshu_a, zsshv_a                  ) 
    170168      CALL wrk_alloc( jpi,jpj,   zhf ) 
    171       IF( ln_wd ) CALL wrk_alloc( jpi, jpj, zcpx, zcpy, wduflt1, wdvflt1,sshai ) 
     169      IF( ln_wd ) CALL wrk_alloc( jpi, jpj, zcpx, zcpy ) 
    172170      ! 
    173171      zmdi=1.e+20                               !  missing data indicator for masking 
     
    375373      IF( .NOT.ln_linssh ) THEN                 ! Variable volume : remove surface pressure gradient 
    376374        IF( ln_wd ) THEN                        ! Calculating and applying W/D gravity filters 
    377           wduflt1(:,:) = 1.0_wp 
    378           wdvflt1(:,:) = 1.0_wp 
    379375          DO jj = 2, jpjm1 
    380376             DO ji = 2, jpim1 
     
    392388                ELSE 
    393389                  zcpx(ji,jj)    = 0._wp 
    394                   wduflt1(ji,jj) = 0.0_wp 
    395390                END IF 
    396391 
     
    408403                ELSE 
    409404                  zcpy(ji,jj)    = 0._wp 
    410                   wdvflt1(ji,jj) = 0.0_wp 
    411405                ENDIF 
    412406 
     
    414408           END DO 
    415409 
    416            CALL lbc_lnk( zcpx, 'U', 1._wp )    ;   CALL lbc_lnk( zcpy, 'V', 1._wp ) 
    417410 
    418411           DO jj = 2, jpjm1 
    419412              DO ji = 2, jpim1 
    420413                 zu_trd(ji,jj) = ( zu_trd(ji,jj) - grav * ( sshn(ji+1,jj  ) - sshn(ji  ,jj ) )   & 
    421                         &                        * r1_e1u(ji,jj) ) * zcpx(ji,jj) * wduflt1(ji,jj) 
     414                        &                        * r1_e1u(ji,jj) ) * zcpx(ji,jj) 
    422415                 zv_trd(ji,jj) = ( zv_trd(ji,jj) - grav * ( sshn(ji  ,jj+1) - sshn(ji  ,jj ) )   & 
    423                         &                        * r1_e2v(ji,jj) ) * zcpy(ji,jj) * wdvflt1(ji,jj) 
     416                        &                        * r1_e2v(ji,jj) ) * zcpy(ji,jj) 
    424417              END DO 
    425418           END DO 
     
    568561      ENDIF 
    569562 
    570       IF( ln_wd ) THEN      !preserve the positivity of water depth 
    571                           !ssh[b,n,a] should have already been processed for this 
    572          sshbb_e(:,:) = MAX(sshbb_e(:,:), rn_wdmin1 - bathy(:,:)) 
    573          sshb_e(:,:)  = MAX(sshb_e(:,:) , rn_wdmin1 - bathy(:,:)) 
    574       ENDIF 
    575563      ! 
    576564      IF (ln_bt_fw) THEN                  ! FORWARD integration: start from NOW fields                     
     
    597585      ! 
    598586      ! Initialize sums: 
    599       IF(ln_wd) sshai(:,:) = ssha(:,:) 
    600587      ua_b  (:,:) = 0._wp       ! After barotropic velocities (or transport if flux form)           
    601588      va_b  (:,:) = 0._wp 
     
    648635            zhup2_e (:,:) = hu_0(:,:) + zwx(:,:)                ! Ocean depth at U- and V-points 
    649636            zhvp2_e (:,:) = hv_0(:,:) + zwy(:,:) 
    650             IF( ln_wd ) THEN 
    651               zhup2_e(:,:) = MAX(zhup2_e (:,:), rn_wdmin1) 
    652               zhvp2_e(:,:) = MAX(zhvp2_e (:,:), rn_wdmin1) 
    653             END IF 
    654637         ELSE 
    655638            zhup2_e (:,:) = hu_n(:,:) 
     
    689672         ENDIF 
    690673#endif 
    691          !IF( ln_wd ) CALL wad_lmt_bt(zwx, zwy, sshn_e, zssh_frc, rdtbt) 
     674         IF( ln_wd ) CALL wad_lmt_bt(zwx, zwy, sshn_e, zssh_frc, rdtbt) 
    692675         ! 
    693676         ! Sum over sub-time-steps to compute advective velocities 
     
    704687         END DO 
    705688 
    706          IF(ln_wd) THEN 
    707            ssha_e(:,:) = (  sshn_e(:,:) - rdtbt * ( zssh_frc(:,:) + zhdiv(:,:) ) * wdmask(:,:)    & 
    708            &             + ( 1._wp - wdmask(:,:) ) * ( sshai(:,:) - sshn_e(:,:) ) ) * ssmask(:,:) 
    709          ELSE 
    710            ssha_e(:,:) = (  sshn_e(:,:) - rdtbt * ( zssh_frc(:,:) + zhdiv(:,:) )  ) * ssmask(:,:) 
    711          END IF 
     689         ssha_e(:,:) = (  sshn_e(:,:) - rdtbt * ( zssh_frc(:,:) + zhdiv(:,:) )  ) * ssmask(:,:) 
    712690          
    713691         CALL lbc_lnk( ssha_e, 'T',  1._wp ) 
     
    758736          &            + za2 *  sshb_e(:,:) + za3 *  sshbb_e(:,:) 
    759737         IF( ln_wd ) THEN                   ! Calculating and applying W/D gravity filters 
    760            wduflt1(:,:) = 1._wp 
    761            wdvflt1(:,:) = 1._wp 
    762738           DO jj = 2, jpjm1 
    763739              DO ji = 2, jpim1 
     
    775751                 ELSE 
    776752                    zcpx(ji,jj)    = 0._wp 
    777                     wduflt1(ji,jj) = 0._wp 
    778753                 END IF 
    779754 
     
    791766                 ELSE 
    792767                    zcpy(ji,jj)    = 0._wp 
    793                     wdvflt1(ji,jj) = 0._wp 
    794768                 END IF 
    795769              END DO 
    796770            END DO 
    797             CALL lbc_lnk( zcpx, 'U', 1._wp )    ;   CALL lbc_lnk( zcpy, 'V', 1._wp ) 
    798771         ENDIF 
    799772         ! 
     
    813786               END DO 
    814787            END DO 
    815  
    816             IF( ln_wd ) THEN 
    817               zhust_e(:,:) = MAX(zhust_e (:,:), rn_wdmin1 ) 
    818               zhvst_e(:,:) = MAX(zhvst_e (:,:), rn_wdmin1 ) 
    819             END IF 
    820788 
    821789         ENDIF 
     
    935903               DO ji = fs_2, fs_jpim1   ! vector opt. 
    936904 
    937                   IF( ln_wd ) THEN 
    938                     zhura = MAX(hu_0(ji,jj) + zsshu_a(ji,jj), rn_wdmin1) 
    939                     zhvra = MAX(hv_0(ji,jj) + zsshv_a(ji,jj), rn_wdmin1) 
    940                   ELSE 
    941                     zhura = hu_0(ji,jj) + zsshu_a(ji,jj) 
    942                     zhvra = hv_0(ji,jj) + zsshv_a(ji,jj) 
    943                   END IF 
     905                  zhura = hu_0(ji,jj) + zsshu_a(ji,jj) 
     906                  zhvra = hv_0(ji,jj) + zsshv_a(ji,jj) 
    944907                  zhura = ssumask(ji,jj)/(zhura + 1._wp - ssumask(ji,jj)) 
    945908                  zhvra = ssvmask(ji,jj)/(zhvra + 1._wp - ssvmask(ji,jj)) 
     
    961924         ! 
    962925         IF( .NOT.ln_linssh ) THEN                     !* Update ocean depth (variable volume case only) 
    963             IF( ln_wd ) THEN 
    964               hu_e (:,:) = MAX(hu_0(:,:) + zsshu_a(:,:), rn_wdmin1) 
    965               hv_e (:,:) = MAX(hv_0(:,:) + zsshv_a(:,:), rn_wdmin1) 
    966             ELSE 
    967               hu_e (:,:) = hu_0(:,:) + zsshu_a(:,:) 
    968               hv_e (:,:) = hv_0(:,:) + zsshv_a(:,:) 
    969             END IF 
     926            hu_e (:,:) = hu_0(:,:) + zsshu_a(:,:) 
     927            hv_e (:,:) = hv_0(:,:) + zsshv_a(:,:) 
    970928            hur_e(:,:) = ssumask(:,:) / ( hu_e(:,:) + 1._wp - ssumask(:,:) ) 
    971929            hvr_e(:,:) = ssvmask(:,:) / ( hv_e(:,:) + 1._wp - ssvmask(:,:) ) 
     
    1008966         !                                   ! Sum sea level 
    1009967         ssha(:,:) = ssha(:,:) + za1 * ssha_e(:,:) 
    1010          IF(ln_wd) ssha(:,:) = ssha(:,:) * wdmask(:,:) + (1._wp - wdmask(:,:)) * sshai(:,:) 
    1011968         !                                                 ! ==================== ! 
    1012969      END DO                                               !        end loop      ! 
     
    11261083      CALL wrk_dealloc( jpi,jpj,   zsshu_a, zsshv_a                                   ) 
    11271084      CALL wrk_dealloc( jpi,jpj,   zhf ) 
    1128       IF( ln_wd ) CALL wrk_dealloc( jpi, jpj, zcpx, zcpy, wduflt1, wdvflt1 ) 
     1085      IF( ln_wd ) CALL wrk_dealloc( jpi, jpj, zcpx, zcpy ) 
    11291086      ! 
    11301087      IF ( ln_diatmb ) THEN 
Note: See TracChangeset for help on using the changeset viewer.