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 7646 for trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90 – NEMO

Ignore:
Timestamp:
2017-02-06T10:25:03+01:00 (7 years ago)
Author:
timgraham
Message:

Merge of dev_merge_2016 into trunk. UPDATE TO ARCHFILES NEEDED for XIOS2.
LIM_SRC_s/limrhg.F90 to follow in next commit due to change of kind (I'm unable to do it in this commit).
Merged using the following steps:

1) svn merge --reintegrate ​svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/trunk .
2) Resolve minor conflicts in sette.sh and namelist_cfg for ORCA2LIM3 (due to a change in trunk after branch was created)
3) svn commit
4) svn switch ​svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/trunk
5) svn merge ​svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/branches/2016/dev_merge_2016 .
6) At this stage I checked out a clean copy of the branch to compare against what is about to be committed to the trunk.
6) svn commit #Commit code to the trunk

In this commit I have also reverted a change to Fcheck_archfile.sh which was causing problems on the Paris machine.

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90

    r6152 r7646  
    1515   !!             3.6  ! 2013-11  (A. Coward) Update for z-tilde compatibility 
    1616   !!             3.7  ! 2015-11  (J. Chanut) free surface simplification 
     17   !!              -   ! 2016-12  (G. Madec, E. Clementi) update for Stoke-Drift divergence 
    1718   !!--------------------------------------------------------------------- 
    1819 
    … …  
    3334   USE dynvor          ! vorticity term 
    3435   USE wet_dry         ! wetting/drying flux limter 
    35    USE bdy_par         ! for lk_bdy 
     36   USE bdy_oce         ! open boundary 
    3637   USE bdytides        ! open boundary condition data 
    3738   USE bdydyn2d        ! open boundary conditions on barotropic variables 
    3839   USE sbctide         ! tides 
    3940   USE updtide         ! tide potential 
     41   USE sbcwave         ! surface wave 
    4042   ! 
    4143   USE in_out_manager  ! I/O manager 
    … …  
    6971   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:) ::   wgtbtp1, wgtbtp2   !: 1st & 2nd weights used in time filtering of barotropic fields 
    7072 
    71    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::  zwz          !: ff/h at F points 
     73   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::  zwz          !: ff_f/h at F points 
    7274   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::  ftnw, ftne   !: triad of coriolis parameter 
    7375   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::  ftsw, ftse   !: (only used with een vorticity scheme) 
    … …  
    156158      REAL(wp), POINTER, DIMENSION(:,:) :: zhf 
    157159      REAL(wp), POINTER, DIMENSION(:,:) :: zcpx, zcpy                 ! Wetting/Dying gravity filter coef. 
    158       REAL(wp), POINTER, DIMENSION(:,:) :: wduflt1, wdvflt1           ! Wetting/Dying velocity filter coef. 
    159160      !!---------------------------------------------------------------------- 
    160161      ! 
    … …  
    168169      CALL wrk_alloc( jpi,jpj,   zsshu_a, zsshv_a                  ) 
    169170      CALL wrk_alloc( jpi,jpj,   zhf ) 
    170       IF( ln_wd ) CALL wrk_alloc( jpi, jpj, zcpx, zcpy, wduflt1, wdvflt1 ) 
     171      IF( ln_wd ) CALL wrk_alloc( jpi, jpj, zcpx, zcpy ) 
    171172      ! 
    172173      zmdi=1.e+20                               !  missing data indicator for masking 
    … …  
    220221      IF( kt == nit000 .OR. .NOT.ln_linssh ) THEN 
    221222         IF( ln_dynvor_een ) THEN               !==  EEN scheme  ==! 
    222             SELECT CASE( nn_een_e3f )              !* ff/e3 at F-point 
     223            SELECT CASE( nn_een_e3f )              !* ff_f/e3 at F-point 
    223224            CASE ( 0 )                                   ! original formulation  (masked averaging of e3t divided by 4) 
    224225               DO jj = 1, jpjm1 
    … …  
    226227                     zwz(ji,jj) =   ( ht_n(ji  ,jj+1) + ht_n(ji+1,jj+1) +                    & 
    227228                        &             ht_n(ji  ,jj  ) + ht_n(ji+1,jj  )   ) * 0.25_wp   
    228                      IF( zwz(ji,jj) /= 0._wp )   zwz(ji,jj) = ff(ji,jj) / zwz(ji,jj) 
     229                     IF( zwz(ji,jj) /= 0._wp )   zwz(ji,jj) = ff_f(ji,jj) / zwz(ji,jj) 
    229230                  END DO 
    230231               END DO 
    … …  
    236237                        &       / ( MAX( 1._wp, tmask(ji  ,jj+1, 1) + tmask(ji+1,jj+1, 1) +    & 
    237238                        &                       tmask(ji  ,jj  , 1) + tmask(ji+1,jj  , 1) ) ) 
    238                      IF( zwz(ji,jj) /= 0._wp )   zwz(ji,jj) = ff(ji,jj) / zwz(ji,jj) 
     239                     IF( zwz(ji,jj) /= 0._wp )   zwz(ji,jj) = ff_f(ji,jj) / zwz(ji,jj) 
    239240                  END DO 
    240241               END DO 
    … …  
    255256            zwz(:,:) = 0._wp 
    256257            zhf(:,:) = 0._wp 
    257             IF ( .not. ln_sco ) THEN 
    258  
    259 !!gm  agree the JC comment  : this should be done in a much clear way 
    260  
    261 ! JC: It not clear yet what should be the depth at f-points over land in z-coordinate case 
    262 !     Set it to zero for the time being  
    263 !              IF( rn_hmin < 0._wp ) THEN    ;   jk = - INT( rn_hmin )                                      ! from a nb of level 
    264 !              ELSE                          ;   jk = MINLOC( gdepw_0, mask = gdepw_0 > rn_hmin, dim = 1 )  ! from a depth 
    265 !              ENDIF 
    266 !              zhf(:,:) = gdepw_0(:,:,jk+1) 
    267             ELSE 
    268                zhf(:,:) = hbatf(:,:) 
    269             END IF 
    270  
    271             DO jj = 1, jpjm1 
    272                zhf(:,jj) = zhf(:,jj) * (1._wp- umask(:,jj,1) * umask(:,jj+1,1)) 
    273             END DO 
     258             
     259!!gm  assume 0 in both cases (xhich is almost surely WRONG ! ) as hvatf has been removed  
     260!!gm    A priori a better value should be something like : 
     261!!gm          zhf(i,j) = masked sum of  ht(i,j) , ht(i+1,j) , ht(i,j+1) , (i+1,j+1)  
     262!!gm                     divided by the sum of the corresponding mask  
     263!!gm  
     264!!             
     265              IF ( .not. ln_sco ) THEN 
     266   
     267   !!gm  agree the JC comment  : this should be done in a much clear way 
     268   
     269   ! JC: It not clear yet what should be the depth at f-points over land in z-coordinate case 
     270   !     Set it to zero for the time being  
     271   !              IF( rn_hmin < 0._wp ) THEN    ;   jk = - INT( rn_hmin )                                      ! from a nb of level 
     272   !              ELSE                          ;   jk = MINLOC( gdepw_0, mask = gdepw_0 > rn_hmin, dim = 1 )  ! from a depth 
     273   !              ENDIF 
     274   !              zhf(:,:) = gdepw_0(:,:,jk+1) 
     275               ELSE 
     276                 !zhf(:,:) = hbatf(:,:) 
     277                 DO jj = 1, jpjm1 
     278                   DO ji = 1, jpim1 
     279                     zhf(ji,jj) = MAX( 0._wp,                                & 
     280                                & ( ht_0(ji  ,jj  )*tmask(ji  ,jj  ,1) +     & 
     281                                &   ht_0(ji+1,jj  )*tmask(ji+1,jj  ,1) +     & 
     282                                &   ht_0(ji  ,jj+1)*tmask(ji  ,jj+1,1) +     & 
     283                                &   ht_0(ji+1,jj+1)*tmask(ji+1,jj+1,1) ) /   & 
     284                                & ( tmask(ji  ,jj  ,1) + tmask(ji+1,jj  ,1) +& 
     285                                &   tmask(ji  ,jj+1,1) + tmask(ji+1,jj+1,1) +& 
     286                                &   rsmall  ) ) 
     287                   END DO 
     288                 END DO 
     289              END IF 
     290   
     291              DO jj = 1, jpjm1 
     292                 zhf(:,jj) = zhf(:,jj) * (1._wp- umask(:,jj,1) * umask(:,jj+1,1)) 
     293              END DO 
     294!!gm end 
    274295 
    275296            DO jk = 1, jpkm1 
    … …  
    285306               END DO 
    286307            END DO 
    287             zwz(:,:) = ff(:,:) * zwz(:,:) 
     308            zwz(:,:) = ff_f(:,:) * zwz(:,:) 
    288309         ENDIF 
    289310      ENDIF 
    … …  
    324345         END DO 
    325346      END DO 
     347       
     348!!gm  Question here when removing the Vertically integrated trends, we remove the vertically integrated NL trends on momentum.... 
     349!!gm  Is it correct to do so ?   I think so... 
     350       
     351       
    326352      !                                   !* barotropic Coriolis trends (vorticity scheme dependent) 
    327353      !                                   ! -------------------------------------------------------- 
    … …  
    374400      IF( .NOT.ln_linssh ) THEN                 ! Variable volume : remove surface pressure gradient 
    375401        IF( ln_wd ) THEN                        ! Calculating and applying W/D gravity filters 
    376           wduflt1(:,:) = 1.0_wp 
    377           wdvflt1(:,:) = 1.0_wp 
    378           DO jj = 2, jpjm1 
    379              DO ji = 2, jpim1 
    380                 ll_tmp1 = MIN(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj))   & 
    381                         & .and. MAX(sshn(ji,jj) + bathy(ji,jj), sshn(ji+1,jj) + bathy(ji+1,jj))   & 
    382                         &  > rn_wdmin1 + rn_wdmin2 
    383                 ll_tmp2 = MAX(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj))   & 
    384                         &                                   + rn_wdmin1 + rn_wdmin2 
     402           DO jj = 2, jpjm1 
     403              DO ji = 2, jpim1  
     404                ll_tmp1 = MIN(   sshn(ji,jj)               ,   sshn(ji+1,jj) ) >                & 
     405                     &    MAX( -ht_wd(ji,jj)               , -ht_wd(ji+1,jj) ) .AND.            & 
     406                     &    MAX(   sshn(ji,jj) + ht_wd(ji,jj),   sshn(ji+1,jj) + ht_wd(ji+1,jj) ) & 
     407                     &                                                         > rn_wdmin1 + rn_wdmin2 
     408                ll_tmp2 = ( ABS( sshn(ji+1,jj)             -   sshn(ji  ,jj))  > 1.E-12 ).AND.( & 
     409                     &    MAX(   sshn(ji,jj)               ,   sshn(ji+1,jj) ) >                & 
     410                     &    MAX( -ht_wd(ji,jj)               , -ht_wd(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 ) 
     411    
    385412                IF(ll_tmp1) THEN 
    386                   zcpx(ji,jj)    = 1.0_wp 
    387                 ELSEIF(ll_tmp2) THEN 
    388                    ! no worries about sshn(ji+1,jj)-sshn(ji,jj) = 0, it won't happen here 
    389                   zcpx(ji,jj) = ABS((sshn(ji+1,jj) + bathy(ji+1,jj) - sshn(ji,jj) - bathy(ji,jj)) & 
    390                         &          /(sshn(ji+1,jj) - sshn(ji,jj))) 
     413                  zcpx(ji,jj) = 1.0_wp 
     414                ELSE IF(ll_tmp2) THEN 
     415                  ! no worries about  sshn(ji+1,jj) -  sshn(ji  ,jj) = 0, it won't happen ! here 
     416                  zcpx(ji,jj) = ABS( (sshn(ji+1,jj) + ht_wd(ji+1,jj) - sshn(ji,jj) - ht_wd(ji,jj)) & 
     417                              &    / (sshn(ji+1,jj) -  sshn(ji  ,jj)) ) 
    391418                ELSE 
    392                   zcpx(ji,jj)    = 0._wp 
    393                   wduflt1(ji,jj) = 0.0_wp 
     419                  zcpx(ji,jj) = 0._wp 
    394420                END IF 
    395  
    396                 ll_tmp1 = MIN(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1))   & 
    397                         & .and. MAX(sshn(ji,jj) + bathy(ji,jj), sshn(ji,jj+1) + bathy(ji,jj+1))   & 
    398                         &  > rn_wdmin1 + rn_wdmin2 
    399                 ll_tmp2 = MAX(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1))   & 
    400                         &                                   + rn_wdmin1 + rn_wdmin2 
     421          
     422                ll_tmp1 = MIN(   sshn(ji,jj)               ,   sshn(ji,jj+1) ) >                & 
     423                     &    MAX( -ht_wd(ji,jj)               , -ht_wd(ji,jj+1) ) .AND.            & 
     424                     &    MAX(   sshn(ji,jj) + ht_wd(ji,jj),   sshn(ji,jj+1) + ht_wd(ji,jj+1) ) & 
     425                     &                                                         > rn_wdmin1 + rn_wdmin2 
     426                ll_tmp2 = ( ABS( sshn(ji,jj)               -   sshn(ji,jj+1))  > 1.E-12 ).AND.( & 
     427                     &    MAX(   sshn(ji,jj)               ,   sshn(ji,jj+1) ) >                & 
     428                     &    MAX( -ht_wd(ji,jj)               , -ht_wd(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 ) 
     429    
    401430                IF(ll_tmp1) THEN 
    402                    zcpy(ji,jj)    = 1.0_wp 
    403                 ELSEIF(ll_tmp2) THEN 
    404                    ! no worries about sshn(ji,jj+1)-sshn(ji,jj) = 0, it won't happen here 
    405                   zcpy(ji,jj) = ABS((sshn(ji,jj+1) + bathy(ji,jj+1) - sshn(ji,jj) - bathy(ji,jj)) & 
    406                         &          /(sshn(ji,jj+1) - sshn(ji,jj))) 
     431                  zcpy(ji,jj) = 1.0_wp 
     432                ELSE IF(ll_tmp2) THEN 
     433                  ! no worries about  sshn(ji,jj+1) -  sshn(ji,jj  ) = 0, it won't happen ! here 
     434                  zcpy(ji,jj) = ABS( (sshn(ji,jj+1) + ht_wd(ji,jj+1) - sshn(ji,jj) - ht_wd(ji,jj)) & 
     435                              &    / (sshn(ji,jj+1) -  sshn(ji,jj  )) ) 
    407436                ELSE 
    408                   zcpy(ji,jj)    = 0._wp 
    409                   wdvflt1(ji,jj) = 0.0_wp 
    410                 ENDIF 
    411  
    412              END DO 
     437                  zcpy(ji,jj) = 0._wp 
     438                END IF 
     439              END DO 
    413440           END DO 
    414  
    415            CALL lbc_lnk( zcpx, 'U', 1._wp )    ;   CALL lbc_lnk( zcpy, 'V', 1._wp ) 
    416  
     441  
    417442           DO jj = 2, jpjm1 
    418443              DO ji = 2, jpim1 
    419                  zu_trd(ji,jj) = ( zu_trd(ji,jj) - grav * ( sshn(ji+1,jj  ) - sshn(ji  ,jj ) )   & 
    420                         &                        * r1_e1u(ji,jj) ) * zcpx(ji,jj) * wduflt1(ji,jj) 
    421                  zv_trd(ji,jj) = ( zv_trd(ji,jj) - grav * ( sshn(ji  ,jj+1) - sshn(ji  ,jj ) )   & 
    422                         &                        * r1_e2v(ji,jj) ) * zcpy(ji,jj) * wdvflt1(ji,jj) 
     444                 zu_trd(ji,jj) = zu_trd(ji,jj) - grav * ( sshn(ji+1,jj  ) - sshn(ji  ,jj ) )   & 
     445                        &                        * r1_e1u(ji,jj) * zcpx(ji,jj) 
     446                 zv_trd(ji,jj) = zv_trd(ji,jj) - grav * ( sshn(ji  ,jj+1) - sshn(ji  ,jj ) )   & 
     447                        &                        * r1_e2v(ji,jj) * zcpy(ji,jj) 
    423448              END DO 
    424449           END DO 
    … …  
    474499      ! 
    475500      !                                         ! Add top stress contribution from baroclinic velocities:       
    476       IF (ln_bt_fw) THEN 
     501      IF( ln_bt_fw ) THEN 
    477502         DO jj = 2, jpjm1 
    478503            DO ji = fs_2, fs_jpim1   ! vector opt. 
    … …  
    538563                &                        + fwfisf(:,:) + fwfisf_b(:,:)                     ) 
    539564      ENDIF 
     565      ! 
     566      IF( ln_sdw ) THEN                         ! Stokes drift divergence added if necessary 
     567         zssh_frc(:,:) = zssh_frc(:,:) + div_sd(:,:) 
     568      ENDIF 
     569      ! 
    540570#if defined key_asminc 
    541571      !                                         ! Include the IAU weighted SSH increment 
    … …  
    567597      ENDIF 
    568598 
    569       IF( ln_wd ) THEN      !preserve the positivity of water depth 
    570                           !ssh[b,n,a] should have already been processed for this 
    571          sshbb_e(:,:) = MAX(sshbb_e(:,:), rn_wdmin1 - bathy(:,:)) 
    572          sshb_e(:,:)  = MAX(sshb_e(:,:) , rn_wdmin1 - bathy(:,:)) 
    573       ENDIF 
    574599      ! 
    575600      IF (ln_bt_fw) THEN                  ! FORWARD integration: start from NOW fields                     
    … …  
    608633         ! Update only tidal forcing at open boundaries 
    609634#if defined key_tide 
    610          IF( lk_bdy      .AND. lk_tide )   CALL bdy_dta_tides( kt, kit=jn, time_offset= noffset+1 ) 
    611          IF( ln_tide_pot .AND. lk_tide )   CALL upd_tide     ( kt, kit=jn, time_offset= noffset   ) 
     635         IF( ln_bdy      .AND. ln_tide )   CALL bdy_dta_tides( kt, kit=jn, time_offset= noffset+1 ) 
     636         IF( ln_tide_pot .AND. ln_tide )   CALL upd_tide     ( kt, kit=jn, time_offset= noffset   ) 
    612637#endif 
    613638         ! 
    … …  
    646671            zhup2_e (:,:) = hu_0(:,:) + zwx(:,:)                ! Ocean depth at U- and V-points 
    647672            zhvp2_e (:,:) = hv_0(:,:) + zwy(:,:) 
    648             IF( ln_wd ) THEN 
    649               zhup2_e(:,:) = MAX(zhup2_e (:,:), rn_wdmin1) 
    650               zhvp2_e(:,:) = MAX(zhvp2_e (:,:), rn_wdmin1) 
    651             END IF 
    652673         ELSE 
    653674            zhup2_e (:,:) = hu_n(:,:) 
    … …  
    702723         END DO 
    703724         ssha_e(:,:) = (  sshn_e(:,:) - rdtbt * ( zssh_frc(:,:) + zhdiv(:,:) )  ) * ssmask(:,:) 
    704          IF( ln_wd ) ssha_e(:,:) = MAX(ssha_e(:,:), rn_wdmin1 - bathy(:,:))  
     725          
    705726         CALL lbc_lnk( ssha_e, 'T',  1._wp ) 
    706727 
    707 #if defined key_bdy 
    708728         ! Duplicate sea level across open boundaries (this is only cosmetic if linssh=T) 
    709          IF( lk_bdy )   CALL bdy_ssh( ssha_e ) 
    710 #endif 
     729         IF( ln_bdy )   CALL bdy_ssh( ssha_e ) 
    711730#if defined key_agrif 
    712731         IF( .NOT.Agrif_Root() )   CALL agrif_ssh_ts( jn ) 
    … …  
    750769          &            + za2 *  sshb_e(:,:) + za3 *  sshbb_e(:,:) 
    751770         IF( ln_wd ) THEN                   ! Calculating and applying W/D gravity filters 
    752            wduflt1(:,:) = 1._wp 
    753            wdvflt1(:,:) = 1._wp 
    754771           DO jj = 2, jpjm1 
    755               DO ji = 2, jpim1 
    756                  ll_tmp1 = MIN( zsshp2_e(ji,jj), zsshp2_e(ji+1,jj) ) > MAX( -bathy(ji,jj), -bathy(ji+1,jj) ) & 
    757                         & .AND. MAX( zsshp2_e(ji,jj) + bathy(ji,jj), zsshp2_e(ji+1,jj) + bathy(ji+1,jj) )    & 
    758                         &                                  > rn_wdmin1 + rn_wdmin2 
    759                  ll_tmp2 = MAX( zsshp2_e(ji,jj), zsshp2_e(ji+1,jj) ) > MAX( -bathy(ji,jj), -bathy(ji+1,jj) ) & 
    760                         &                                  + rn_wdmin1 + rn_wdmin2 
    761                  IF(ll_tmp1) THEN 
    762                     zcpx(ji,jj) = 1._wp 
    763                  ELSE IF(ll_tmp2) THEN 
    764                     ! no worries about zsshp2_e(ji+1,jj)-zsshp2_e(ji,jj) = 0, it won't happen here 
    765                     zcpx(ji,jj) = ABS( (zsshp2_e(ji+1,jj) + bathy(ji+1,jj) - zsshp2_e(ji,jj) - bathy(ji,jj)) & 
    766                         &             / (zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj)) ) 
    767                  ELSE 
    768                     zcpx(ji,jj)    = 0._wp 
    769                     wduflt1(ji,jj) = 0._wp 
    770                  END IF 
    771  
    772                  ll_tmp1 = MIN( zsshp2_e(ji,jj), zsshp2_e(ji,jj+1) ) > MAX( -bathy(ji,jj), -bathy(ji,jj+1) ) & 
    773                         & .AND. MAX( zsshp2_e(ji,jj) + bathy(ji,jj), zsshp2_e(ji,jj+1) + bathy(ji,jj+1) )    & 
    774                         &                                  > rn_wdmin1 + rn_wdmin2 
    775                  ll_tmp2 = MAX( zsshp2_e(ji,jj), zsshp2_e(ji,jj+1) ) > MAX( -bathy(ji,jj), -bathy(ji,jj+1) ) & 
    776                         &                                  + rn_wdmin1 + rn_wdmin2 
    777                  IF(ll_tmp1) THEN 
    778                     zcpy(ji,jj) = 1._wp 
    779                  ELSE IF(ll_tmp2) THEN 
    780                     ! no worries about zsshp2_e(ji,jj+1)-zsshp2_e(ji,jj) = 0, it won't happen here 
    781                     zcpy(ji,jj) = ABS( (zsshp2_e(ji,jj+1) + bathy(ji,jj+1) - zsshp2_e(ji,jj) - bathy(ji,jj)) & 
    782                         &             / (zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj)) ) 
    783                  ELSE 
    784                     zcpy(ji,jj)    = 0._wp 
    785                     wdvflt1(ji,jj) = 0._wp 
    786                  END IF 
     772              DO ji = 2, jpim1  
     773                ll_tmp1 = MIN( zsshp2_e(ji,jj)               , zsshp2_e(ji+1,jj) ) >                & 
     774                     &    MAX(   -ht_wd(ji,jj)               ,   -ht_wd(ji+1,jj) ) .AND.            & 
     775                     &    MAX( zsshp2_e(ji,jj) + ht_wd(ji,jj), zsshp2_e(ji+1,jj) + ht_wd(ji+1,jj) ) & 
     776                     &                                                             > rn_wdmin1 + rn_wdmin2 
     777                ll_tmp2 = (ABS(zsshp2_e(ji,jj)               - zsshp2_e(ji+1,jj))  > 1.E-12 ).AND.( & 
     778                     &    MAX( zsshp2_e(ji,jj)               , zsshp2_e(ji+1,jj) ) >                & 
     779                     &    MAX(   -ht_wd(ji,jj)               ,   -ht_wd(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 ) 
     780    
     781                IF(ll_tmp1) THEN 
     782                  zcpx(ji,jj) = 1.0_wp 
     783                ELSE IF(ll_tmp2) THEN 
     784                  ! no worries about  zsshp2_e(ji+1,jj) - zsshp2_e(ji  ,jj) = 0, it won't happen ! here 
     785                  zcpx(ji,jj) = ABS( (zsshp2_e(ji+1,jj) +    ht_wd(ji+1,jj) - zsshp2_e(ji,jj) - ht_wd(ji,jj)) & 
     786                              &    / (zsshp2_e(ji+1,jj) - zsshp2_e(ji  ,jj)) ) 
     787                ELSE 
     788                  zcpx(ji,jj) = 0._wp 
     789                END IF 
     790          
     791                ll_tmp1 = MIN( zsshp2_e(ji,jj)               , zsshp2_e(ji,jj+1) ) >                & 
     792                     &    MAX(   -ht_wd(ji,jj)               ,   -ht_wd(ji,jj+1) ) .AND.            & 
     793                     &    MAX( zsshp2_e(ji,jj) + ht_wd(ji,jj), zsshp2_e(ji,jj+1) + ht_wd(ji,jj+1) ) & 
     794                     &                                                             > rn_wdmin1 + rn_wdmin2 
     795                ll_tmp2 = (ABS(zsshp2_e(ji,jj)               - zsshp2_e(ji,jj+1))  > 1.E-12 ).AND.( & 
     796                     &    MAX( zsshp2_e(ji,jj)               , zsshp2_e(ji,jj+1) ) >                & 
     797                     &    MAX(   -ht_wd(ji,jj)               ,   -ht_wd(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 ) 
     798    
     799                IF(ll_tmp1) THEN 
     800                  zcpy(ji,jj) = 1.0_wp 
     801                ELSE IF(ll_tmp2) THEN 
     802                  ! no worries about  zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj  ) = 0, it won't happen ! here 
     803                  zcpy(ji,jj) = ABS( (zsshp2_e(ji,jj+1) +    ht_wd(ji,jj+1) - zsshp2_e(ji,jj) - ht_wd(ji,jj)) & 
     804                              &    / (zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj  )) ) 
     805                ELSE 
     806                  zcpy(ji,jj) = 0._wp 
     807                END IF 
    787808              END DO 
    788             END DO 
    789             CALL lbc_lnk( zcpx, 'U', 1._wp )    ;   CALL lbc_lnk( zcpy, 'V', 1._wp ) 
    790          ENDIF 
     809           END DO 
     810         END IF 
    791811         ! 
    792812         ! Compute associated depths at U and V points: 
    … …  
    806826            END DO 
    807827 
    808             IF( ln_wd ) THEN 
    809               zhust_e(:,:) = MAX(zhust_e (:,:), rn_wdmin1 ) 
    810               zhvst_e(:,:) = MAX(zhvst_e (:,:), rn_wdmin1 ) 
    811             END IF 
    812  
    813828         ENDIF 
    814829         ! 
    … …  
    861876         ! 
    862877         ! Add tidal astronomical forcing if defined 
    863          IF ( lk_tide.AND.ln_tide_pot ) THEN 
     878         IF ( ln_tide .AND. ln_tide_pot ) THEN 
    864879            DO jj = 2, jpjm1 
    865880               DO ji = fs_2, fs_jpim1   ! vector opt. 
    … …  
    967982         CALL lbc_lnk_multi( ua_e, 'U', -1._wp, va_e , 'V', -1._wp ) 
    968983         ! 
    969 #if defined key_bdy   
    970984         !                                                 ! open boundaries 
    971          IF( lk_bdy )   CALL bdy_dyn2d( jn, ua_e, va_e, un_e, vn_e, hur_e, hvr_e, ssha_e ) 
    972 #endif 
     985         IF( ln_bdy )   CALL bdy_dyn2d( jn, ua_e, va_e, un_e, vn_e, hur_e, hvr_e, ssha_e ) 
    973986#if defined key_agrif                                                            
    974987         IF( .NOT.Agrif_Root() )  CALL agrif_dyn_ts( jn )  ! Agrif 
    … …  
    10861099      CALL wrk_dealloc( jpi,jpj,   zsshu_a, zsshv_a                                   ) 
    10871100      CALL wrk_dealloc( jpi,jpj,   zhf ) 
    1088       IF( ln_wd ) CALL wrk_dealloc( jpi, jpj, zcpx, zcpy, wduflt1, wdvflt1 ) 
     1101      IF( ln_wd ) CALL wrk_dealloc( jpi, jpj, zcpx, zcpy ) 
    10891102      ! 
    10901103      IF ( ln_diatmb ) THEN 
    … …  
    12481261            zxr2 = r1_e1t(ji,jj) * r1_e1t(ji,jj) 
    12491262            zyr2 = r1_e2t(ji,jj) * r1_e2t(ji,jj) 
    1250             zcu(ji,jj) = SQRT( grav * ht_0(ji,jj) * (zxr2 + zyr2) ) 
     1263            zcu(ji,jj) = SQRT( grav * MAX(ht_0(ji,jj),0._wp) * (zxr2 + zyr2) ) 
    12511264         END DO 
    12521265      END DO 
Note: See TracChangeset for help on using the changeset viewer.