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 5727 for branches/UKMO/2014_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90 – NEMO

Ignore:
Timestamp:
2015-09-10T19:05:13+02:00 (9 years ago)
Author:
rfurner
Message:

some bug fixes for wetting and drying elements...still not working though

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/UKMO/2014_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90

    r5066 r5727  
    371371      !                                   !* Right-Hand-Side of the barotropic momentum equation 
    372372      !                                   ! ---------------------------------------------------- 
     373 
    373374      IF( lk_vvl ) THEN                         ! Variable volume : remove surface pressure gradient 
    374375        IF(ln_wd) THEN                    ! Calculating and applying W/D gravity filters 
    375376          DO jj = 2, jpjm1 
    376377             DO ji = 2, jpim1 
    377                 ll_tmp1 = MIN(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) 
    378                 ll_tmp2 = MAX(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) + & 
     378                IF ( tmask(ji+1,jj,1) == 0._wp ) THEN 
     379                   zcpx = 1._wp 
     380                ELSE 
     381                   ll_tmp1 = MIN(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj)*tmask(ji,jj,1), -bathy(ji+1,jj)*tmask(ji+1,jj,1)) & 
     382                           & .and. MAX(sshn(ji,jj) + bathy(ji,jj)*tmask(ji,jj,1), sshn(ji+1,jj) + bathy(ji+1,jj)*tmask(ji+1,jj,1)) & 
     383                           &  > rn_wdmin1 + rn_wdmin2 
     384                   ll_tmp2 = MAX(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj)*tmask(ji,jj,1), -bathy(ji+1,jj)*tmask(ji+1,jj,1)) + & 
    379385                                                          & rn_wdmin1 + rn_wdmin2 
    380                 IF(ll_tmp1) THEN 
    381                   zcpx(ji,jj) = 1.0_wp 
    382                 ELSE IF(ll_tmp2) THEN 
    383                    ! no worries about sshn(ji+1,jj)-sshn(ji,jj) = 0, it won't happen ! here 
    384                   zcpx(ji,jj) = ABS((sshn(ji+1,jj) + bathy(ji+1,jj) - sshn(ji,jj) - bathy(ji,jj)) /& 
    385                               &     (sshn(ji+1,jj) - sshn(ji,jj))) 
     386                   IF(ll_tmp1) THEN 
     387                      zcpx(ji,jj) = 1.0_wp 
     388                   ELSE IF(ll_tmp2) THEN 
     389                      ! no worries about sshn(ji+1,jj)-sshn(ji,jj) = 0, it won't happen ! here 
     390                      zcpx(ji,jj) = ABS((sshn(ji+1,jj) + bathy(ji+1,jj)*tmask(ji+1,jj,1) - sshn(ji,jj) - bathy(ji,jj)*tmask(ji,jj,1)) /& 
     391                                  &     (sshn(ji+1,jj) - sshn(ji,jj))) 
     392                   ELSE 
     393                      zcpx(ji,jj) = 0._wp 
     394                   END IF 
     395                ENDIF 
     396 
     397                IF ( tmask(ji,jj+1,1) == 0._wp ) THEN 
     398                   zcpy = 1._wp 
    386399                ELSE 
    387                    zcpx(ji,jj) = 0._wp 
    388                 END IF 
    389  
    390                 ll_tmp1 = MIN(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) 
    391                 ll_tmp2 = MAX(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) + & 
    392                                                           & rn_wdmin1 + rn_wdmin2 
    393                 IF(ll_tmp1) THEN 
    394                    zcpy(ji,jj) = 1.0_wp 
    395                 ELSE IF(ll_tmp2) THEN 
    396                    ! no worries about sshn(ji,jj+1)-sshn(ji,jj) = 0, it won't happen ! here 
    397                   zcpy(ji,jj) = ABS((sshn(ji,jj+1) + bathy(ji,jj+1) - sshn(ji,jj) + bathy(ji,jj)) /& 
    398                               &     (sshn(ji,jj+1) - sshn(ji,jj))) 
    399                 ELSE 
    400                   zcpy(ji,jj) = 0._wp 
    401                 END IF 
     400                   ll_tmp1 = MIN(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj)*tmask(ji,jj,1), -bathy(ji,jj+1)*tmask(ji,jj+1,1)) & 
     401                           & .and. MAX(sshn(ji,jj) + bathy(ji,jj)*tmask(ji,jj,1), sshn(ji,jj+1) + bathy(ji,jj+1)*tmask(ji,jj+1,1)) & 
     402                           &  > rn_wdmin1 + rn_wdmin2 
     403                   ll_tmp2 = MAX(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj)*tmask(ji,jj,1), -bathy(ji,jj+1)*tmask(ji,jj+1,1)) + & 
     404                                                             & rn_wdmin1 + rn_wdmin2 
     405                   IF(ll_tmp1) THEN 
     406                      zcpy(ji,jj) = 1.0_wp 
     407                   ELSE IF(ll_tmp2) THEN 
     408                      ! no worries about sshn(ji,jj+1)-sshn(ji,jj) = 0, it won't happen ! here 
     409                     zcpy(ji,jj) = ABS((sshn(ji,jj+1) + bathy(ji,jj+1)*tmask(ji,jj+1,1) - sshn(ji,jj) - bathy(ji,jj)*tmask(ji,jj,1)) /& 
     410                                 &     (sshn(ji,jj+1) - sshn(ji,jj))) 
     411                   ELSE 
     412                     zcpy(ji,jj) = 0._wp 
     413                   END IF 
     414                ENDIF 
    402415             END DO 
    403416           END DO 
     
    453466      ! 
    454467      ! Note that the "unclipped" bottom friction parameter is used even with explicit drag 
    455       zu_frc(:,:) = zu_frc(:,:) + hur(:,:) * bfrua(:,:) * zwx(:,:) 
    456       zv_frc(:,:) = zv_frc(:,:) + hvr(:,:) * bfrva(:,:) * zwy(:,:) 
     468      IF(ln_wd) THEN 
     469        zu_frc(:,:) = zu_frc(:,:) + MAX(hur(:,:) * bfrua(:,:),-1._wp / rdtbt) * zwx(:,:) 
     470        zv_frc(:,:) = zv_frc(:,:) + MAX(hvr(:,:) * bfrva(:,:),-1._wp / rdtbt) * zwy(:,:) 
     471      ELSE 
     472        zu_frc(:,:) = zu_frc(:,:) + hur(:,:) * bfrua(:,:) * zwx(:,:) 
     473        zv_frc(:,:) = zv_frc(:,:) + hvr(:,:) * bfrva(:,:) * zwy(:,:) 
     474      END IF 
    457475      !        
    458476      IF (ln_bt_fw) THEN                        ! Add wind forcing 
     
    518536      IF(ln_wd) THEN      !preserve the positivity of water depth 
    519537                          !ssh[b,n,a] should have already been processed for this 
    520          sshbb_e(:,:) = MAX(sshbb_e(:,:), rn_wdmin1 - bathy(:,:)) 
    521          sshb_e(:,:)  = MAX(sshb_e(:,:) , rn_wdmin1 - bathy(:,:)) 
     538         sshbb_e(:,:jj) = MAX( sshbb_e(:,:), (rn_wdmin1 - bathy(:,:)) ) *tmask(:,:,1) 
     539         sshb_e(:,:jj)  = MAX( sshb_e(:,:) , (rn_wdmin1 - bathy(:,:)) ) *tmask(:,:,1) 
    522540      ENDIF 
    523541 
     
    637655         ENDIF 
    638656#endif 
     657 
    639658         IF(ln_wd) CALL wad_lmt_bt(zwx, zwy, sshn_e, zssh_frc, rdtbt) 
    640659         ! 
     
    652671         END DO 
    653672         ssha_e(:,:) = (  sshn_e(:,:) - rdtbt * ( zssh_frc(:,:) + zhdiv(:,:) )  ) * tmask(:,:,1) 
    654          IF(ln_wd) ssha_e(:,:) = MAX(ssha_e(:,:), rn_wdmin1 - bathy(:,:))  
     673         IF(ln_wd) ssha_e(:,:) = MAX(ssha_e(:,:), (rn_wdmin1 - bathy(:,:)) ) * tmask(:,:,1) 
    655674         CALL lbc_lnk( ssha_e, 'T',  1._wp ) 
    656675 
     
    785804         ! 
    786805         ! Add bottom stresses: 
    787          zu_trd(:,:) = zu_trd(:,:) + bfrua(:,:) * zun_e(:,:) * hur_e(:,:) 
    788          zv_trd(:,:) = zv_trd(:,:) + bfrva(:,:) * zvn_e(:,:) * hvr_e(:,:) 
    789          ! 
     806         IF(ln_wd) THEN 
     807           zu_trd(:,:) = zu_trd(:,:) + MAX(bfrua(:,:) * hur_e(:,:), -1._wp / rdtbt) * zun_e(:,:) 
     808           zv_trd(:,:) = zv_trd(:,:) + MAX(bfrva(:,:) * hvr_e(:,:), -1._wp / rdtbt) * zvn_e(:,:) 
     809         ELSE 
     810           zu_trd(:,:) = zu_trd(:,:) + bfrua(:,:) * zun_e(:,:) * hur_e(:,:) 
     811           zv_trd(:,:) = zv_trd(:,:) + bfrva(:,:) * zvn_e(:,:) * hvr_e(:,:) 
     812         END IF 
     813          
    790814         IF(ln_wd) THEN                    ! Calculating and applying W/D gravity filters 
    791815           DO jj = 2, jpjm1 
    792816              DO ji = 2, jpim1 
    793                  ll_tmp1 = MIN(zsshp2_e(ji,jj), zsshp2_e(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) 
    794                  ll_tmp2 = MAX(zsshp2_e(ji,jj), zsshp2_e(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) + & 
     817                 IF ( tmask(ji+1,jj,1) == 0._wp ) THEN 
     818                    zcpx = 1._wp 
     819                 ELSE 
     820                    ll_tmp1 = MIN(zsshp2_e(ji,jj), zsshp2_e(ji+1,jj)) > MAX(-bathy(ji,jj)*tmask(ji,jj,1), -bathy(ji+1,jj)*tmask(ji+1,jj,1)) & 
     821                           & .and. MAX(zsshp2_e(ji,jj) + bathy(ji,jj)*tmask(ji,jj,1), zsshp2_e(ji+1,jj) + bathy(ji+1,jj)*tmask(ji+1,jj,1))  & 
     822                           &  > rn_wdmin1 + rn_wdmin2 
     823                    ll_tmp2 = MAX(zsshp2_e(ji,jj), zsshp2_e(ji+1,jj)) > MAX(-bathy(ji,jj)*tmask(ji,jj,1), -bathy(ji+1,jj)*tmask(ji+1,jj,1)) + & 
    795824                                                           & rn_wdmin1 + rn_wdmin2 
    796                  IF(ll_tmp1) THEN 
    797                    zcpx(ji,jj) = 1.0_wp 
    798                  ELSE IF(ll_tmp2) THEN 
    799                     ! no worries about zsshp2_e(ji+1,jj)-zsshp2_e(ji,jj) = 0, it won't happen ! here 
    800                    zcpx(ji,jj) = ABS((zsshp2_e(ji+1,jj) + bathy(ji+1,jj) - zsshp2_e(ji,jj) - bathy(ji,jj)) /& 
    801                                & (zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj))) 
     825                    IF(ll_tmp1) THEN 
     826                      zcpx(ji,jj) = 1.0_wp 
     827                    ELSE IF(ll_tmp2) THEN 
     828                       ! no worries about zsshp2_e(ji+1,jj)-zsshp2_e(ji,jj) = 0, it won't happen ! here 
     829                      zcpx(ji,jj) = ABS((zsshp2_e(ji+1,jj) + bathy(ji+1,jj)*tmask(ji+1,jj,1) - zsshp2_e(ji,jj) - bathy(ji,jj)*tmask(ji,jj,1)) /& 
     830                                 & (zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj))) 
     831                    ELSE 
     832                       zcpx(ji,jj) = 0._wp 
     833                    END IF 
     834                 ENDIF 
     835 
     836                 IF ( tmask(ji,jj+1,1) == 0._wp ) THEN 
     837                    zcpy = 1._wp 
    802838                 ELSE 
    803                     zcpx(ji,jj) = 0._wp 
    804                  END IF 
    805  
    806                  ll_tmp1 = MIN(zsshp2_e(ji,jj), zsshp2_e(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) 
    807                  ll_tmp2 = MAX(zsshp2_e(ji,jj), zsshp2_e(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) + & 
     839                    ll_tmp1 = MIN(zsshp2_e(ji,jj), zsshp2_e(ji,jj+1)) > MAX(-bathy(ji,jj)*tmask(ji,jj,1), -bathy(ji,jj+1)*tmask(ji,jj+1,1))& 
     840                           & .and. MAX(zsshp2_e(ji,jj) + bathy(ji,jj)*tmask(ji,jj,1), zsshp2_e(ji,jj+1) + bathy(ji,jj+1)*tmask(ji,jj+1,1)) & 
     841                           &  > rn_wdmin1 + rn_wdmin2 
     842                    ll_tmp2 = MAX(zsshp2_e(ji,jj), zsshp2_e(ji,jj+1)) > MAX(-bathy(ji,jj)*tmask(ji,jj,1), -bathy(ji,jj+1)*tmask(ji,jj+1,1)) + & 
    808843                                                           & rn_wdmin1 + rn_wdmin2 
    809                  IF(ll_tmp1) THEN 
    810                     zcpy(ji,jj) = 1.0_wp 
    811                  ELSE IF(ll_tmp2) THEN 
    812                     ! no worries about zsshp2_e(ji,jj+1)-zsshp2_e(ji,jj) = 0, it won't happen ! here 
    813                    zcpy(ji,jj) = ABS((zsshp2_e(ji,jj+1) + bathy(ji,jj+1) - zsshp2_e(ji,jj) - bathy(ji,jj)) /& 
    814                                & (zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj))) 
    815                  ELSE 
    816                    zcpy(ji,jj) = 0._wp 
    817                  END IF 
     844                    IF(ll_tmp1) THEN 
     845                       zcpy(ji,jj) = 1.0_wp 
     846                    ELSE IF(ll_tmp2) THEN 
     847                       ! no worries about zsshp2_e(ji,jj+1)-zsshp2_e(ji,jj) = 0, it won't happen ! here 
     848                      zcpy(ji,jj) = ABS((zsshp2_e(ji,jj+1) + bathy(ji,jj+1)*tmask(ji,jj+1,1) - zsshp2_e(ji,jj) - bathy(ji,jj)*tmask(ji,jj,1)) /& 
     849                                  & (zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj))) 
     850                    ELSE 
     851                      zcpy(ji,jj) = 0._wp 
     852                    END IF 
     853                 ENDIF 
    818854              END DO 
    819855            END DO 
     
    894930            !                                          !  ----------------------------------------------         
    895931            IF(ln_wd) THEN 
    896               hu_e (:,:) = MAX(hu_0(:,:) + zsshu_a(:,:), rn_wdmin1) 
    897               hv_e (:,:) = MAX(hv_0(:,:) + zsshv_a(:,:), rn_wdmin1) 
     932              hu_e (:,:) = MAX(hu_0(:,:) + zsshu_a(:,:), rn_wdmin1 * umask(:,:,1) ) 
     933              hv_e (:,:) = MAX(hv_0(:,:) + zsshv_a(:,:), rn_wdmin1 * vmask(:,:,1) ) 
    898934            ELSE 
    899935              hu_e (:,:) = hu_0(:,:) + zsshu_a(:,:) 
     
    10321068      CALL wrk_dealloc( jpi, jpj, zsshu_a, zsshv_a                                   ) 
    10331069      CALL wrk_dealloc( jpi, jpj, zhf ) 
    1034       IF(ln_wd) CALL wrk_alloc( jpi, jpj, zcpx, zcpy ) 
     1070      IF(ln_wd) CALL wrk_dealloc( jpi, jpj, zcpx, zcpy ) 
    10351071      ! 
    10361072      IF( nn_timing == 1 )  CALL timing_stop('dyn_spg_ts') 
Note: See TracChangeset for help on using the changeset viewer.