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 7433 – NEMO

Changeset 7433


Ignore:
Timestamp:
2016-12-02T11:46:27+01:00 (6 years ago)
Author:
timgraham
Message:

#1811 Fix for wetting and drying

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/2016/dev_merge_2016/NEMOGCM/NEMO/OPA_SRC/DYN/wet_dry.F90

    r7421 r7433  
    3131   !! --------------------------------------------------------------------- 
    3232 
    33    REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) ::   wduflt, wdvflt !: u- and v- filter 
    3433   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) ::   wdmask         !: u- and v- limiter  
    3534 
     
    4443   PUBLIC   wad_lmt                   ! routine called by sshwzv.F90 
    4544   PUBLIC   wad_lmt_bt                ! routine called by dynspg_ts.F90 
     45   PUBLIC   wad_istate                ! routine called by istate.F90 and domvvl.F90 
    4646 
    4747   !! * Substitutions 
     
    8585      ! 
    8686      IF(ln_wd) THEN 
    87          ALLOCATE( wduflt(jpi,jpj), wdvflt(jpi,jpj), wdmask(jpi,jpj), STAT=ierr ) 
     87         ALLOCATE( wdmask(jpi,jpj), STAT=ierr ) 
    8888         IF( ierr /= 0 ) CALL ctl_stop('STOP', 'wad_init : Array allocation error') 
    8989      ENDIF 
     
    144144        ! Horizontal Flux in u and v direction 
    145145        DO jk = 1, jpkm1   
    146            DO jj = 1, jpjm1 
    147               DO ji = 1, jpim1 
     146           DO jj = 1, jpj 
     147              DO ji = 1, jpi 
    148148                 zflxu(ji,jj) = zflxu(ji,jj) + e3u_n(ji,jj,jk) * un(ji,jj,jk) * umask(ji,jj,jk) 
    149149                 zflxv(ji,jj) = zflxv(ji,jj) + e3v_n(ji,jj,jk) * vn(ji,jj,jk) * vmask(ji,jj,jk) 
     
    155155        zflxv(:,:) = zflxv(:,:) * e1v(:,:) 
    156156        
    157          DO jj = 2, jpjm1 
    158            DO ji = 2, jpim1  
     157        wdmask(:,:) = 1 
     158        DO jj = 2, jpj 
     159           DO ji = 2, jpi  
    159160 
    160161             IF( tmask(ji,jj,1) == 0._wp  )   CYCLE   ! we don't care about land cells 
     
    182183           zflxv1(:,:) = zflxv(:,:) * zwdlmtv(:,:) 
    183184           
    184            DO jj = 2, jpjm1 
    185               DO ji = 2, jpim1  
     185           DO jj = 2, jpj 
     186              DO ji = 2, jpi  
    186187         
    187188                 wdmask(ji,jj) = 0 
     
    201202                 IF(zdep1 > zdep2) THEN 
    202203                   zflag = 1 
    203                    wdmask(ji, jj) = 1 
     204                   wdmask(ji, jj) = 0 
    204205                   zcoef = ( ( zdep2 - rn_wdmin2 ) * ztmp - zzflxn * z2dt ) / ( zflxp(ji,jj) * z2dt ) 
    205206                   zcoef = max(zcoef, 0._wp) 
     
    208209                   IF(zflxu1(ji-1,jj) < 0._wp) zwdlmtu(ji-1,jj) = zcoef 
    209210                   IF(zflxv1(ji,  jj) > 0._wp) zwdlmtv(ji  ,jj) = zcoef 
    210                    IF(zflxv1(ji,jj-1) < 0._wp) zwdlmtv(ji-1,jj) = zcoef 
     211                   IF(zflxv1(ji,jj-1) < 0._wp) zwdlmtv(ji,jj-1) = zcoef 
    211212                 END IF 
    212213              END DO ! ji loop 
     
    230231        CALL lbc_lnk( un, 'U', -1. ) 
    231232        CALL lbc_lnk( vn, 'V', -1. ) 
     233      ! 
     234        un_b(:,:) = un_b(:,:) * zwdlmtu(:, :) 
     235        vn_b(:,:) = vn_b(:,:) * zwdlmtv(:, :) 
     236        CALL lbc_lnk( un_b, 'U', -1. ) 
     237        CALL lbc_lnk( vn_b, 'V', -1. ) 
    232238        
    233239        IF(zflag == 1 .AND. lwp) WRITE(numout,*) 'Need more iterations in wad_lmt!!!' 
     
    239245        CALL wrk_dealloc( jpi, jpj, zflxp, zflxn, zflxu, zflxv, zflxu1, zflxv1 ) 
    240246        CALL wrk_dealloc( jpi, jpj, zwdlmtu, zwdlmtv) 
    241          ! 
     247        ! 
    242248      ENDIF 
    243249      ! 
     
    298304        ! Horizontal Flux in u and v direction 
    299305        
    300         !zflxu(:,:) = zflxu(:,:) * e2u(:,:) 
    301         !zflxv(:,:) = zflxv(:,:) * e1v(:,:) 
    302         
    303         DO jj = 2, jpjm1 
    304            DO ji = 2, jpim1  
     306        DO jj = 2, jpj 
     307           DO ji = 2, jpi  
    305308 
    306309             IF(tmask(ji,jj,1) < 0.5_wp) CYCLE   ! we don't care about land cells 
     
    313316 
    314317              zdep2 = ht_0(ji,jj) + sshn_e(ji,jj) - rn_wdmin1 
    315               IF(zdep2 < 0._wp) THEN  !add more safty, but not necessary 
    316                 !zdep2 = 0._wp 
    317                sshn_e(ji,jj) = rn_wdmin1 - ht_0(ji,jj) 
    318               END IF 
    319318           ENDDO 
    320319        END DO 
     
    328327           zflxv1(:,:) = zflxv(:,:) * zwdlmtv(:,:) 
    329328           
    330            DO jj = 2, jpjm1 
    331               DO ji = 2, jpim1  
     329           DO jj = 2, jpj 
     330              DO ji = 2, jpi  
    332331         
    333                  wdmask(ji,jj) = 0 
    334332                 IF(tmask(ji,jj,1) < 0.5_wp) CYCLE  
    335333                 IF(ht_0 (ji,jj)   > zdepwd) CYCLE 
     
    355353                   IF(zflxu1(ji-1,jj) < 0._wp) zwdlmtu(ji-1,jj) = zcoef 
    356354                   IF(zflxv1(ji,  jj) > 0._wp) zwdlmtv(ji  ,jj) = zcoef 
    357                    IF(zflxv1(ji,jj-1) < 0._wp) zwdlmtv(ji-1,jj) = zcoef 
     355                   IF(zflxv1(ji,jj-1) < 0._wp) zwdlmtv(ji,jj-1) = zcoef 
    358356                 END IF 
    359357              END DO ! ji loop 
     
    384382        CALL wrk_dealloc( jpi, jpj, zflxp, zflxn, zflxu1, zflxv1 ) 
    385383        CALL wrk_dealloc( jpi, jpj, zwdlmtu, zwdlmtv) 
    386          ! 
     384        ! 
    387385      END IF 
    388       ! 
     386 
    389387      IF( nn_timing == 1 )  CALL timing_stop('wad_lmt') 
    390       ! 
    391388   END SUBROUTINE wad_lmt_bt 
     389 
     390   SUBROUTINE wad_istate 
     391      !!---------------------------------------------------------------------- 
     392      !!                   ***  ROUTINE wad_istate  *** 
     393      !!  
     394      !! ** Purpose :   Initialization of the dynamics and tracers for WAD test 
     395      !!      configurations (channels or bowls with initial ssh gradients) 
     396      !! 
     397      !! ** Method  : - set temperature field 
     398      !!              - set salinity field 
     399      !!              - set ssh slope (needs to be repeated in domvvl_rst_init to 
     400      !!                               set vertical metrics ) 
     401      !!---------------------------------------------------------------------- 
     402      ! 
     403      INTEGER  ::   ji, jj            ! dummy loop indices 
     404      REAL(wp) ::   zi, zj 
     405      !!---------------------------------------------------------------------- 
     406      ! 
     407      ! Uniform T & S in all test cases 
     408      tsn(:,:,:,jp_tem) = 10._wp 
     409      tsb(:,:,:,jp_tem) = 10._wp 
     410      tsn(:,:,:,jp_sal) = 35._wp 
     411      tsb(:,:,:,jp_sal) = 35._wp 
     412      SELECT CASE ( jp_cfg )  
     413         !                                        ! ==================== 
     414         CASE ( 1 )                               ! WAD 1 configuration 
     415            !                                     ! ==================== 
     416            ! 
     417            IF(lwp) WRITE(numout,*) 
     418            IF(lwp) WRITE(numout,*) 'istate_wad : Closed box with EW linear bottom slope' 
     419            IF(lwp) WRITE(numout,*) '~~~~~~~~~~' 
     420            ! 
     421            do ji = 1,jpi 
     422             sshn(ji,:) = ( -5.5_wp + 5.5_wp*FLOAT(mig(ji))/FLOAT(jpidta-1))*tmask(ji,:,1) 
     423            end do 
     424            !                                     ! ==================== 
     425         CASE ( 2 )                               ! WAD 2 configuration 
     426            !                                     ! ==================== 
     427            ! 
     428            IF(lwp) WRITE(numout,*) 
     429            IF(lwp) WRITE(numout,*) 'istate_wad : Parobolic EW channel, mid-range initial ssh slope' 
     430            IF(lwp) WRITE(numout,*) '~~~~~~~~~~' 
     431            ! 
     432            do ji = 1,jpi 
     433             sshn(ji,:) = ( -5.5_wp + 3.9_wp*FLOAT(jpidta - mig(ji))/FLOAT(jpidta-1))*tmask(ji,:,1) 
     434            end do 
     435            !                                     ! ==================== 
     436         CASE ( 3 )                               ! WAD 3 configuration 
     437            !                                     ! ==================== 
     438            ! 
     439            IF(lwp) WRITE(numout,*) 
     440            IF(lwp) WRITE(numout,*) 'istate_wad : Parobolic EW channel, extreme initial ssh slope'  
     441            IF(lwp) WRITE(numout,*) '~~~~~~~~~~' 
     442            ! 
     443            do ji = 1,jpi 
     444             sshn(ji,:) = ( -7.5_wp + 6.9_wp*FLOAT(jpidta - mig(ji))/FLOAT(jpidta-1))*tmask(ji,:,1) 
     445            end do 
     446 
     447            ! 
     448            !                                     ! ==================== 
     449         CASE ( 4 )                               ! WAD 4 configuration 
     450            !                                     ! ==================== 
     451            ! 
     452            IF(lwp) WRITE(numout,*) 
     453            IF(lwp) WRITE(numout,*) 'istate_wad : Parobolic bowl, mid-range initial ssh slope'  
     454            IF(lwp) WRITE(numout,*) '~~~~~~~~~~' 
     455            ! 
     456            DO ji = 1, jpi 
     457               zi = MAX(1.0-FLOAT((mig(ji)-25)**2)/400.0, 0.0 ) 
     458               DO jj = 1, jpj 
     459                  zj = MAX(1.0-FLOAT((mjg(jj)-17)**2)/144.0, 0.0 ) 
     460                  sshn(ji,jj) = -8.5_wp + 8.5_wp*zi*zj 
     461               END DO 
     462            END DO 
     463 
     464            ! 
     465            !                                    ! =========================== 
     466         CASE ( 5 )                              ! WAD 5 configuration 
     467            !                                    ! ==================== 
     468            ! 
     469            IF(lwp) WRITE(numout,*) 
     470            IF(lwp) WRITE(numout,*) 'istate_wad : Double slope with shelf' 
     471            IF(lwp) WRITE(numout,*) '~~~~~~~~~~' 
     472            ! 
     473            ! Needed rn_wdmin2 increased to 0.01 for this case? 
     474            do ji = 1,jpi 
     475             sshn(ji,:) = ( -5.5_wp + 9.0_wp*FLOAT(mig(ji))/FLOAT(jpidta-1))*tmask(ji,:,1) 
     476            end do 
     477 
     478            ! 
     479            !                                     ! =========================== 
     480         CASE ( 6 )                               ! WAD 6 configuration 
     481            !                                     ! ==================== 
     482            ! 
     483            IF(lwp) WRITE(numout,*) 
     484            IF(lwp) WRITE(numout,*) 'istate_wad : Parobolic EW channel with gaussian ridge'  
     485            IF(lwp) WRITE(numout,*) '~~~~~~~~~~' 
     486            ! 
     487            do ji = 1,jpi 
     488             !6a 
     489             sshn(ji,:) = ( -5.5_wp + 9.0_wp*FLOAT(jpidta - mig(ji))/FLOAT(jpidta-1))*tmask(ji,:,1) 
     490             !Some variations in initial slope that have been tested 
     491             !6b 
     492             !sshn(ji,:) = ( -5.5_wp + 6.5_wp*FLOAT(jpidta - mig(ji))/FLOAT(jpidta-1))*tmask(ji,:,1) 
     493             !6c 
     494             !sshn(ji,:) = ( -5.5_wp + 7.5_wp*FLOAT(jpidta - mig(ji))/FLOAT(jpidta-1))*tmask(ji,:,1) 
     495             !6d 
     496             !sshn(ji,:) = ( -4.5_wp + 8.0_wp*FLOAT(jpidta - mig(ji))/FLOAT(jpidta-1))*tmask(ji,:,1) 
     497            end do 
     498 
     499            ! 
     500            !                                    ! =========================== 
     501         CASE DEFAULT                            ! NONE existing configuration 
     502            !                                    ! =========================== 
     503            WRITE(ctmp1,*) 'WAD test with a ', jp_cfg,' option is not coded' 
     504            ! 
     505            CALL ctl_stop( ctmp1 ) 
     506            ! 
     507      END SELECT 
     508      ! 
     509      ! Apply minimum wetdepth criterion 
     510      ! 
     511      do jj = 1,jpj 
     512         do ji = 1,jpi 
     513            IF( bathy(ji,jj) + sshn(ji,jj) < rn_wdmin1 ) THEN 
     514               sshn(ji,jj) = tmask(ji,jj,1)*( rn_wdmin1 - bathy(ji,jj) ) 
     515            ENDIF 
     516         end do 
     517      end do 
     518      sshb = sshn 
     519      ssha = sshn 
     520      ! 
     521   END SUBROUTINE wad_istate 
    392522 
    393523   !!============================================================================== 
Note: See TracChangeset for help on using the changeset viewer.