Changeset 7016


Ignore:
Timestamp:
2016-10-10T18:24:04+02:00 (4 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.

Location:
branches/2016/dev_r6393_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/DYN
Files:
2 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 
  • branches/2016/dev_r6393_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/DYN/wet_dry.F90

    r6995 r7016  
    298298        zflxp(:,:)   = 0._wp 
    299299        zflxn(:,:)   = 0._wp 
    300         !zflxu(:,:)   = 0._wp 
    301         !zflxv(:,:)   = 0._wp 
    302300 
    303301        zwdlmtu(:,:)  = 1._wp 
     
    306304        ! Horizontal Flux in u and v direction 
    307305        
    308         !zflxu(:,:) = zflxu(:,:) * e2u(:,:) 
    309         !zflxv(:,:) = zflxv(:,:) * e1v(:,:) 
    310         
    311         DO jj = 2, jpjm1 
    312            DO ji = 2, jpim1  
     306        DO jj = 2, jpj 
     307           DO ji = 2, jpi  
    313308 
    314309             IF(tmask(ji, jj, 1) < 0.5_wp) CYCLE   ! we don't care about land cells 
     
    321316 
    322317              zdep2 = bathy(ji,jj) + sshn_e(ji,jj) - rn_wdmin1 
    323               IF(zdep2 < 0._wp) THEN  !add more safty, but not necessary 
    324                 !zdep2 = 0._wp 
    325                sshn_e(ji,jj) = rn_wdmin1 - bathy(ji,jj) 
    326               END IF 
    327318           ENDDO 
    328319        END DO 
     
    336327           zflxv1(:,:) = zflxv(:,:) * zwdlmtv(:,:) 
    337328           
    338            DO jj = 2, jpjm1 
    339               DO ji = 2, jpim1  
     329           DO jj = 2, jpj 
     330              DO ji = 2, jpi  
    340331         
    341                  !wdmask(ji,jj) = 0 
    342332                 IF(tmask(ji, jj, 1) < 0.5_wp) CYCLE  
    343333                 IF(bathy(ji,jj) > zdepwd) CYCLE 
     
    356346                 IF(zdep1 > zdep2) THEN 
    357347                   zflag = 1 
    358                    !wdmask(ji, jj) = 1 
    359348                   zcoef = ( ( zdep2 - rn_wdmin2 ) * ztmp - zzflxn * z2dt ) / ( zflxp(ji,jj) * z2dt ) 
    360349                   zcoef = max(zcoef, 0._wp) 
     
    386375        IF(zflag == 1 .AND. lwp) WRITE(numout,*) 'Need more iterations in wad_lmt_bt!!!' 
    387376        
    388         !IF( ln_rnf      )   CALL sbc_rnf_div( hdivn )          ! runoffs (update hdivn field) 
    389         !IF( nn_cla == 1 )   CALL cla_div    ( kt )             ! Cross Land Advection (update hdivn field) 
    390377        ! 
    391378        ! 
Note: See TracChangeset for help on using the changeset viewer.