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 7527 for branches/2016/dev_merge_2016/NEMOGCM/CONFIG/WAD_TEST_CASES/MY_SRC/usrdef_zgr.F90 – NEMO

Ignore:
Timestamp:
2017-01-05T18:26:14+01:00 (7 years ago)
Author:
acc
Message:

Branch 2016/dev_merge_2016. Changes to WAD_TEST_CASES to cope with wetting and drying above sea level

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/2016/dev_merge_2016/NEMOGCM/CONFIG/WAD_TEST_CASES/MY_SRC/usrdef_zgr.F90

    r7514 r7527  
    6767      REAL(wp) ::   ze3min            ! local scalar 
    6868      REAL(wp) ::   zi, zj, zbathy    ! local scalar 
    69       REAL(wp) ::   ztmpu, ztmpv, ztmpf, ztmpu1, ztmpv1, ztmpf1 
     69      REAL(wp) ::   ztmpu, ztmpv, ztmpf, ztmpu1, ztmpv1, ztmpf1, zwet 
    7070      REAL(wp), DIMENSION(jpi,jpj) ::   zht, zhu, zhv, z2d   ! 2D workspace 
    7171      !!---------------------------------------------------------------------- 
     
    105105               DO ji = 1, jpi 
    106106                 zi = MIN((glamt(ji,1) - 10.0)/40.0, 1.0 ) 
    107                  zht(ji,:) = MAX(zbathy*zi, 1.5)  
     107                 zht(ji,:) = MAX(zbathy*zi, -2.0)  
    108108               END DO 
    109109               zht(mi0(1):mi1(1),:) = -4._wp 
     
    120120               ! 
    121121               DO ji = 1, jpi 
    122                  zi = MAX(1.0-((glamt(ji,1)-25.0)**2)/484.0, 0.0 ) 
    123                  zht(ji,:) = MAX(zbathy*zi, 1.5) 
     122                 zi = MAX(1.0-((glamt(ji,1)-25.0)**2)/484.0, -0.3 ) 
     123                 zht(ji,:) = MAX(zbathy*zi, -2.0) 
    124124               END DO 
    125125               zht(mi0(1):mi1(1),:) = -4._wp 
     
    139139               DO jj = 1, jpj 
    140140                 zj = MAX(1.0-((gphit(1,jj)-17.0)**2)/196.0, -2.0 ) 
    141                  zht(ji,jj) = MAX(zbathy*zi*zj, 1.5) 
     141                 zht(ji,jj) = MAX(zbathy*zi*zj, -2.0) 
    142142               END DO 
    143143               END DO 
     
    156156               DO ji = 1, jpi 
    157157                 zi = MIN(glamt(ji,1)/45.0, 1.0 ) 
    158                  zht(ji,:) = MAX(zbathy*zi, 1.5) 
     158                 zht(ji,:) = MAX(zbathy*zi, -2.0) 
    159159                 IF( glamt(ji,1) >= 46.0 ) THEN 
    160160                   zht(ji,:) = 10.0 
     
    165165                   zht(ji,:) = 2.5 
    166166                 ELSE IF( glamt(ji,1) >= 4.0 .AND. glamt(ji,1) < 15.0 ) THEN 
    167                    zi = 2.0/11.0 
    168                    zht(ji,:) = MAX(2.5 - zi*(16.0-glamt(ji,1)), 1.5) 
    169                  ELSE IF( glamt(ji,1) >= 1.0 .AND. glamt(ji,1) < 4.0 ) THEN 
    170                    zht(ji,:) = 1.5 
     167                   zi = 4.5/11.0 
     168                   zht(ji,:) = MAX(2.5 - zi*(16.0-glamt(ji,1)), -2.0) 
     169                 ELSE IF( glamt(ji,1) >= 0.0 .AND. glamt(ji,1) < 4.0 ) THEN 
     170                   zht(ji,:) = -2.0 
    171171                 ENDIF 
    172172               END DO 
     
    186186               DO ji = 1, jpi 
    187187                 zi = MAX(1.0-((glamt(ji,1)-25.0)**2)/484.0, -2.0 ) 
    188                  zj = 0.95*MAX(EXP(-1.0*((glamt(ji,1)-25.0)**2)/32.0) , 0.0 ) 
    189                  zht(ji,:) = MAX(zbathy*(zi-zj), 1.0) 
     188                 zj = 1.075*MAX(EXP(-1.0*((glamt(ji,1)-25.0)**2)/32.0) , 0.0 ) 
     189                 zht(ji,:) = MAX(zbathy*(zi-zj), -2.0) 
    190190               END DO 
    191191               zht(mi0(1):mi1(1),:) = -4._wp 
     
    222222      END DO 
    223223      CALL lbc_lnk( zhv, 'V', 1. )     ! boundary condition: this mask the surrounding grid-points 
    224       !avoid the zero depth on T- (U-,V-,F-) points 
    225         DO jj = 1, jpj 
    226           DO ji = 1, jpi 
    227             IF(ABS(zht(ji,jj)) < rn_wdmin1) & 
    228               & zht(ji,jj) = SIGN(1._wp, zht(ji,jj)) * rn_wdmin1 
    229             IF(ABS(zhu(ji,jj)) < rn_wdmin1) & 
    230               & zhu(ji,jj) = SIGN(1._wp, zhu(ji,jj)) * rn_wdmin1 
    231             IF(ABS(zhv(ji,jj)) < rn_wdmin1) & 
    232               & zhv(ji,jj) = SIGN(1._wp, zhv(ji,jj)) * rn_wdmin1 
    233           END DO 
    234         END DO 
     224      DO jj = mj0(1), mj1(1)   ! first  row of global domain only 
     225         zhv(:,jj) = zht(:,jj) 
     226      END DO 
     227      DO jj = mj0(jpjglo), mj1(jpjglo)   ! last  row of global domain only 
     228         zhv(:,jj) = zht(:,jj) 
     229      END DO 
    235230      !      
    236231      CALL zgr_z1d( pdept_1d, pdepw_1d, pe3t_1d , pe3w_1d )   ! Reference z-coordinate system 
     
    264259              IF( zht(ji,jj) <= -(rn_wdld - rn_wdmin2)) THEN 
    265260                k_bot(ji,jj) = 0 
    266                 zht(ji,jj) = rn_wdmin1 * REAL( jpkm1 , wp) * 0.5_wp 
    267               ELSEIF(zht(ji,jj) <= rn_wdmin1) THEN 
    268                 k_bot(ji,jj) = 2 
    269                 zht(ji,jj) = rn_wdmin1 * REAL( jpkm1 , wp) * 0.5_wp 
     261                k_top(ji,jj) = 0 
    270262              ENDIF 
    271263           END DO 
    272264         END DO 
    273       DO ji = 1, jpim1 
    274          zhu(ji,:) = 0.5_wp * ( zht(ji,:) + zht(ji+1,:) ) 
    275       END DO 
    276       !                                ! ==>>>  set by hand non-zero value on 
    277       !                                first/last columns & rows  
    278       CALL lbc_lnk( zhu, 'U', 1. )     ! boundary condition: this mask the surrounding grid-points 
    279       DO ji = mi0(1), mi1(1)              ! first row of global domain only 
    280          zhu(ji,:) = zht(ji,:) 
    281       END DO 
    282       DO ji = mi0(jpiglo), mi1(jpiglo)   ! last  row of global domain only 
    283          zhu(ji,:) = zht(ji,:) 
    284       END DO 
    285       ! at v-point: averaging zht 
    286       zhv = 0._wp 
    287       DO jj = 1, jpjm1 
    288          zhv(:,jj) = 0.5_wp * ( zht(:,jj) + zht(:,jj+1) ) 
    289       END DO 
    290       CALL lbc_lnk( zhv, 'V', 1. ) 
    291       DO jj = mj0(1), mj1(1)   ! first  row of global domain only 
    292          zhv(:,jj) = zht(:,jj) 
    293       END DO 
    294       DO jj = mj0(jpjglo), mj1(jpjglo)   ! last  row of global domain only 
    295          zhv(:,jj) = zht(:,jj) 
    296       END DO 
    297265         ! 
    298266         !                                !* terrain-following coordinate with e3.(k)=cst) 
    299267         !                                !  OVERFLOW case : identical with j-index (T=V, U=F) 
    300          z1_jpkm1 = 1._wp / REAL( jpkm1 , wp) 
    301          DO jk = 1, jpk 
    302                   pdept(:,:,jk) = zht(:,:) * z1_jpkm1 * ( REAL( jk   , wp ) - 0.5_wp ) 
    303                   pdepw(:,:,jk) = zht(:,:) * z1_jpkm1 * ( REAL( jk-1 , wp )          ) 
    304                   pe3t (:,:,jk) = zht(:,:) * z1_jpkm1 
    305                   pe3u (:,:,jk) = zhu(:,:) * z1_jpkm1 
    306                   pe3v (:,:,jk) = zhv(:,:) * z1_jpkm1 
    307                   pe3f (:,:,jk) = zhu(:,:) * z1_jpkm1 
    308                   pe3w (:,:,jk) = zht(:,:) * z1_jpkm1 
    309                   pe3uw(:,:,jk) = zhu(:,:) * z1_jpkm1 
    310                   pe3vw(:,:,jk) = zhv(:,:) * z1_jpkm1 
    311          END DO       
    312268         DO jj = 1, jpjm1 
    313269            DO ji = 1, jpim1 
    314                ztmpu  = zht(ji,jj)+zht(ji+1,jj) 
    315                ztmpv  = zht(ji,jj)+zht(ji,jj+1) 
    316                ztmpf  = zht(ji,jj)+zht(ji+1,jj)+zht(ji,jj+1)+zht(ji+1,jj+1) 
    317                ztmpu1 = zht(ji,jj)*zht(ji+1,jj) 
    318                ztmpv1 = zht(ji,jj)*zht(ji,jj+1) 
    319                ztmpf1 = MIN(zht(ji,jj), zht(ji+1,jj), zht(ji,jj+1), zht(ji+1,jj+1)) * & 
    320                       & MAX(zht(ji,jj), zht(ji+1,jj), zht(ji,jj+1), zht(ji+1,jj+1)) 
    321                IF(  (ztmpu1 < 0._wp.OR.ABS(ztmpu) < rn_wdmin1) ) THEN 
    322                  pe3u (ji,jj,:) = 0.5_wp* ( pe3t(ji,jj,:) + pe3t(ji+1,jj,:) ) 
    323                  pe3uw(ji,jj,:) = 0.5_wp* ( pe3w(ji,jj,:) + pe3w(ji+1,jj,:) ) 
    324                ENDIF 
    325                IF(  (ztmpv1 < 0._wp.OR.ABS(ztmpv) < rn_wdmin1) ) THEN 
    326                  pe3v (ji,jj,:) = 0.5_wp* ( pe3t(ji,jj,:) + pe3t(ji,jj+1,:) ) 
    327                  pe3vw(ji,jj,:) = 0.5_wp* ( pe3w(ji,jj,:) + pe3w(ji,jj+1,:) ) 
    328                ENDIF 
    329                IF(  (ztmpf1 < 0._wp.OR.ABS(ztmpf) < rn_wdmin1) ) THEN 
    330                  pe3f (ji,jj,:) = 0.25_wp* ( pe3t(ji,jj,:)   + pe3t(ji+1,jj,:) & 
    331                                 &          + pe3t(ji,jj+1,:) + pe3t(ji+1,jj+1,:) ) 
    332                ENDIF 
    333                  pe3f (ji,jj,:) = 0.25_wp* ( pe3t(ji,jj,:)   + pe3t(ji+1,jj,:) & 
    334                                 &          + pe3t(ji,jj+1,:) + pe3t(ji+1,jj+1,:) ) 
    335            END DO 
    336          END DO 
     270              z1_jpkm1 = 1._wp / REAL( k_bot(ji,jj) - k_top(ji,jj) + 1 , wp) 
     271              DO jk = 1, jpk 
     272                  zwet = MAX( zht(ji,jj), rn_wdmin1 ) 
     273                  pdept(ji,jj,jk) = zwet * z1_jpkm1 * ( REAL( jk   , wp ) - 0.5_wp ) 
     274                  pdepw(ji,jj,jk) = zwet * z1_jpkm1 * ( REAL( jk-1 , wp )          ) 
     275                  pe3t (ji,jj,jk) = zwet * z1_jpkm1 
     276                  pe3w (ji,jj,jk) = zwet * z1_jpkm1 
     277                  zwet = MAX( zhu(ji,jj), rn_wdmin1 ) 
     278                  pe3u (ji,jj,jk) = zwet * z1_jpkm1 
     279                  pe3f (ji,jj,jk) = zwet * z1_jpkm1 
     280                  pe3uw(ji,jj,jk) = zwet * z1_jpkm1 
     281                  zwet = MAX( zhv(ji,jj), rn_wdmin1 ) 
     282                  pe3vw(ji,jj,jk) = zwet * z1_jpkm1 
     283                  pe3v (ji,jj,jk) = zwet * z1_jpkm1 
     284              END DO       
     285           END DO       
     286         END DO       
    337287         CALL lbc_lnk( pe3u , 'U', 1. ) 
    338288         CALL lbc_lnk( pe3uw, 'U', 1. ) 
Note: See TracChangeset for help on using the changeset viewer.