Changeset 4375


Ignore:
Timestamp:
2014-01-28T14:55:35+01:00 (7 years ago)
Author:
hliu
Message:

updated gravity filters

Location:
branches/2013/dev_r4050_NOC_WaD/NEMOGCM
Files:
8 edited

Legend:

Unmodified
Added
Removed
  • branches/2013/dev_r4050_NOC_WaD/NEMOGCM/CONFIG/cfg.txt

    r3905 r4375  
    99ORCA2_OFF_PISCES OPA_SRC OFF_SRC TOP_SRC 
    1010ORCA2_LIM_PISCES OPA_SRC LIM_SRC_2 NST_SRC TOP_SRC 
     11PBASIN OPA_SRC 
  • branches/2013/dev_r4050_NOC_WaD/NEMOGCM/NEMO/OPA_SRC/DOM/dom_oce.F90

    r4355 r4375  
    6060   REAL(wp), PUBLIC ::   rn_wadmine = 0.01_wp   !: tolerrance of minimum water depth on dried cells 
    6161   REAL(wp), PUBLIC ::   rn_landele = 20.0_wp   !: land elevation below which wetting/drying will be considered 
     62   INTEGER , PUBLIC ::   nn_waditr  =   10      !: Maximum number of iteration for W/D limiter 
    6263 
    6364   !                                                  !!! associated variables 
  • branches/2013/dev_r4050_NOC_WaD/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90

    r4355 r4375  
    280280      !!---------------------------------------------------------------------- 
    281281      USE ioipsl 
    282       NAMELIST/namwad/ ln_wad , rn_wadmin, rn_wadmine, rn_landele 
     282      NAMELIST/namwad/ ln_wad , rn_wadmin, rn_wadmine, rn_landele, nn_waditr 
    283283      !!---------------------------------------------------------------------- 
    284284 
     
    295295         WRITE(numout,*) '      tollerance of minimual depth    rn_wadmine    = ', rn_wadmine 
    296296         WRITE(numout,*) ' threshold elevation for land masking rn_landele    = ', rn_landele 
     297         WRITE(numout,*) ' Maximum number of iteration for W/D limiter        = ', nn_waditr 
    297298      ENDIF 
    298299 
  • branches/2013/dev_r4050_NOC_WaD/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90

    r4006 r4375  
    14111411!! 
    14121412      ! HYBRID :  
    1413       DO jj = 1, jpj 
    1414          DO ji = 1, jpi 
    1415             DO jk = 1, jpkm1 
    1416                IF( scobot(ji,jj) >= fsdept(ji,jj,jk) )   mbathy(ji,jj) = MAX( 2, jk ) 
    1417                IF( scobot(ji,jj) == 0._wp            )   mbathy(ji,jj) = 0 
    1418             END DO 
    1419          END DO 
    1420       END DO 
     1413      IF(ln_wad) THEN 
     1414      !! Wetting/Drying case 
     1415         DO jj = 1, jpj 
     1416            DO ji = 1, jpi 
     1417               DO jk = 1, jpkm1 
     1418                  IF( scobot(ji,jj) >= fsdept(ji,jj,jk) ) THEN 
     1419                    mbathy(ji,jj) = MAX( 2, jk ) 
     1420                  ELSE IF(scobot(ji,jj) <= rn_landele) THEN 
     1421                    mbathy(ji,jj) = 0 
     1422                  ELSE 
     1423                    mbathy(ji,jj) = 2 
     1424                  END IF 
     1425               END DO 
     1426            END DO 
     1427         END DO 
     1428      ELSE 
     1429         DO jj = 1, jpj 
     1430            DO ji = 1, jpi 
     1431               DO jk = 1, jpkm1 
     1432                  IF( scobot(ji,jj) >= fsdept(ji,jj,jk) )   mbathy(ji,jj) = MAX( 2, jk ) 
     1433                  IF( scobot(ji,jj) == 0._wp            )   mbathy(ji,jj) = 0 
     1434               END DO 
     1435            END DO 
     1436         END DO 
     1437      END IF 
    14211438      IF( nprint == 1 .AND. lwp ) WRITE(numout,*) ' MIN val mbathy h90 ', MINVAL( mbathy(:,:) ),   & 
    14221439         &                                                       ' MAX ', MAXVAL( mbathy(:,:) ) 
  • branches/2013/dev_r4050_NOC_WaD/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90

    r3764 r4375  
    7676      !!             - Save the trend (l_trddyn=T) 
    7777      !!---------------------------------------------------------------------- 
    78       INTEGER, INTENT(in) ::   kt   ! ocean time-step index 
     78      INTEGER, INTENT(in) ::   kt                               ! ocean time-step index 
     79      INTEGER             ::   ji, jj, jk                       ! dummy loop indices 
     80      REAL                ::   zcpx, zcpy                       ! gravity filter for W/D 
     81      REAL                ::   zij, zim1j, zijm1, zim1jm1   ! local variables 
     82      REAL                ::   dij, dim1j, dijm1, dim1jm1   ! local variables 
    7983      REAL(wp), POINTER, DIMENSION(:,:,:) ::  ztrdu, ztrdv 
    8084      !!---------------------------------------------------------------------- 
     
    8286      IF( nn_timing == 1 )  CALL timing_start('dyn_hpg') 
    8387      ! 
    84       IF( l_trddyn ) THEN                    ! Temporary saving of ua and va trends (l_trddyn) 
     88 
     89      IF( l_trddyn .OR. ln_wad ) THEN               ! Temporary saving of ua and va trends (l_trddyn) 
    8590         CALL wrk_alloc( jpi,jpj,jpk, ztrdu, ztrdv ) 
    8691         ztrdu(:,:,:) = ua(:,:,:) 
    8792         ztrdv(:,:,:) = va(:,:,:) 
    8893      ENDIF 
     94 
    8995      ! 
    9096      SELECT CASE ( nhpg )      ! Hydrostatic pressure gradient computation 
     
    95101      CASE (  4 )   ;   CALL hpg_prj    ( kt )      ! s-coordinate (Pressure Jacobian scheme) 
    96102      END SELECT 
     103 
     104      !! Gravity filter for W/D 
     105      IF(ln_wad) THEN 
     106         DO jk = 2, jpkm1 
     107            DO jj = 2, jpjm1 
     108               DO ji = fs_2, fs_jpim1   ! vector opt. 
     109                  zij     =  sshn(ji   , jj   ) 
     110                  zim1j   =  sshn(ji- 1, jj   ) 
     111                  zijm1   =  sshn(ji   , jj- 1) 
     112                  zim1jm1 =  sshn(ji-1 , jj- 1) 
     113                  dij     = bathy(ji   , jj   ) 
     114                  dim1j   = bathy(ji- 1, jj   ) 
     115                  dijm1   = bathy(ji   , jj- 1) 
     116                  dim1jm1 = bathy(ji-1 , jj- 1) 
     117 
     118                  zcpx = 0.5_wp + sign(0.5_wp, min(zij,zim1j) - max(-dij,-dim1j)) 
     119                  zcpy = 0.5_wp + sign(0.5_wp, min(zij,zijm1) - max(-dij,-dijm1)) 
     120                   
     121                  ua(ji,jj,jk) = ztrdu(ji,jj,jk) + zcpx * ( ua(ji,jj,jk) - ztrdu(ji,jj,jk)) 
     122                  va(ji,jj,jk) = ztrdv(ji,jj,jk) + zcpy * ( va(ji,jj,jk) - ztrdv(ji,jj,jk)) 
     123               END DO 
     124            END DO 
     125         END DO 
     126      END IF 
    97127      ! 
    98128      IF( l_trddyn ) THEN      ! save the hydrostatic pressure gradient trends for momentum trend diagnostics 
     
    100130         ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 
    101131         CALL trd_mod( ztrdu, ztrdv, jpdyn_trd_hpg, 'DYN', kt ) 
     132      ENDIF 
     133 
     134      IF( l_trddyn .OR. ln_wad ) THEN  
    102135         CALL wrk_dealloc( jpi,jpj,jpk, ztrdu, ztrdv ) 
    103136      ENDIF 
     137 
    104138      ! 
    105139      IF(ln_ctl)   CALL prt_ctl( tab3d_1=ua, clinfo1=' hpg  - Ua: ', mask1=umask,   & 
  • branches/2013/dev_r4050_NOC_WaD/NEMOGCM/NEMO/OPA_SRC/DYN/wadlmt.F90

    r4356 r4375  
    4949      INTEGER  ::   ji, jj, jk, jk1     ! dummy loop indices 
    5050      INTEGER  ::   zflag, z2dt         ! local scalar 
    51       REAL(wp) ::   zflxp, zflxn        ! local scalars 
    5251      REAL(wp) ::   zcoef, zdep1, zdep2 ! local scalars 
    5352      REAL(wp) ::   ztmp                ! local scalars 
     53      REAL(wp), POINTER,  DIMENSION(:,:) ::   zflxp,  zflxn   ! specific 2D workspace 
    5454      REAL(wp), POINTER,  DIMENSION(:,:) ::   zflxu,  zflxv   ! specific 2D workspace 
    5555      REAL(wp), POINTER,  DIMENSION(:,:) ::   zflxu1, zflxv1  ! specific 2D workspace 
    56       REAL(wp), POINTER,  DIMENSION(:,:) ::   uwdlmt, vwdlmt  !: W/D flux limiters 
     56      REAL(wp), POINTER,  DIMENSION(:,:) ::   wdlmt           !: W/D flux limiters 
     57 
    5758      !!---------------------------------------------------------------------- 
    5859      ! 
     
    6263      IF(ln_wad) THEN 
    6364 
    64         CALL wrk_alloc( jpi, jpj, zflxu, zflxv, zflxu1, zflxv1, uwdlmt, vwdlmt) 
     65        CALL wrk_alloc( jpi, jpj, zflxp, zflxn, zflxu, zflxv, zflxu1, zflxv1, wdlmt ) 
    6566        ! 
    6667        
     
    7172        IF( neuler == 0 .AND. kt == nit000 )  z2dt = rdt 
    7273        
    73         
     74        zflxp(:,:)  = 0._wp 
     75        zflxn(:,:)  = 0._wp 
    7476        zflxu(:,:)  = 0._wp 
    7577        zflxv(:,:)  = 0._wp 
    76         uwdlmt(:,:) = 1._wp 
    77         vwdlmt(:,:) = 1._wp 
    78         
     78        wdlmt(:,:)  = 1._wp 
    7979        
    8080        ! Horizontal Flux in u and v direction 
     
    9191        zflxv(:,:) = zflxv(:,:) * e1v(:,:) 
    9292        
    93         DO jk1 = 1, 10 
    94         
    95         zflag = 0     ! flag marking that if we need further iteration of W/D limiters? 
    96         
    97         zflxu1(:,:) = zflxu(:,:) * uwdlmt(:,:) 
    98         zflxv1(:,:) = zflxv(:,:) * vwdlmt(:,:) 
    99         
    100         
    10193        DO jj = 2, jpjm1 
    10294           DO ji = fs_2, fs_jpim1   ! vector opt. 
    103               zflxp = max(zflxu1(ji,jj), 0._wp) - min(zflxu1(ji-1,jj), 0._wp) + & 
    104                     & max(zflxv1(ji,jj), 0._wp) - min(zflxv1(ji,jj-1), 0._wp)  
    105               zflxn = min(zflxu1(ji,jj), 0._wp) - max(zflxu1(ji-1,jj), 0._wp) + & 
    106                     & min(zflxv1(ji,jj), 0._wp) - max(zflxv1(ji,jj-1), 0._wp)  
    107         
     95 
    10896              IF(tmask(ji, jj, 1) < 0.5_wp) CYCLE  
    109         
    110               !IF( .NOT. AGRIF_Root() ) THEN    !Wetting and Drying but not for AGRIF 
    111                  IF (((nbondi ==  1).OR.(nbondi == 2)).AND.(ji == nlci-1)) CYCLE ! east 
    112                  IF (((nbondi == -1).OR.(nbondi == 2)).AND.(ji == 2     )) CYCLE ! west  
    113                  IF (((nbondj ==  1).OR.(nbondj == 2)).AND.(jj == nlcj-1)) CYCLE ! north 
    114                  IF (((nbondj == -1).OR.(nbondj == 2)).AND.(jj == 2     )) CYCLE ! south 
    115               !END IF 
    116         
    117                
    118               ztmp = e1t(ji,jj) * e2t(ji,jj) 
    119         
    120               zdep1 = (zflxp + zflxn) * z2dt / ztmp 
     97 
     98              zflxp(ji,jj) = max(zflxu(ji,jj), 0._wp) - min(zflxu(ji-1,jj),   0._wp) + & 
     99                           & max(zflxv(ji,jj), 0._wp) - min(zflxv(ji,  jj-1), 0._wp)  
     100              zflxn(ji,jj) = min(zflxu(ji,jj), 0._wp) - max(zflxu(ji-1,jj),   0._wp) + & 
     101                           & min(zflxv(ji,jj), 0._wp) - max(zflxv(ji,  jj-1), 0._wp)  
     102 
    121103              zdep2 = bathy(ji,jj) + sshb(ji,jj) - rn_wadmin  
    122         
    123104              IF(zdep2 < 0._wp) THEN 
    124                 !WRITE(numout,*) 'depth less than minimum depth, cell(ji,jj):', ji,jj 
    125105                zdep2 = 0._wp 
    126106                sshb(ji,jj) = rn_wadmin - bathy(ji,jj) 
    127107              END IF 
     108           ENDDO 
     109        END DO 
     110 
     111       
     112        !! start limiter iteration  
     113        DO jk1 = 1, nn_waditr 
    128114        
    129               IF(zdep1 > zdep2) THEN 
    130                 zflag = 1 
    131                 zcoef = ( ( zdep2 + rn_wadmine ) * ztmp + zflxn * z2dt ) / ( zflxp * z2dt ) 
    132                 IF(zflxu1(ji,  jj) >= 0._wp) uwdlmt(ji,jj  ) = zcoef 
    133                 IF(zflxu1(ji-1,jj) <  0._wp) uwdlmt(ji-1,jj) = zcoef 
    134                 IF(zflxv1(ji,  jj) >= 0._wp) vwdlmt(ji,jj  ) = zcoef 
    135                 IF(zflxv1(ji,jj-1) <  0._wp) vwdlmt(ji,jj-1) = zcoef 
    136               END IF 
    137            END DO ! ji loop 
    138         END DO  ! jj loop 
    139         
    140         CALL lbc_lnk( uwdlmt, 'U', 1. )   ;   CALL lbc_lnk( vwdlmt , 'V', 1. )  
    141         
    142         IF(zflag == 0) EXIT 
    143         
     115           zflag = 0     ! flag indicating if any further iteration is needed? 
     116           
     117           zflxu1(:,:) = zflxu(:,:) * wdlmt(:,:) 
     118           zflxv1(:,:) = zflxv(:,:) * wdlmt(:,:) 
     119           
     120           DO jj = 2, jpjm1 
     121              DO ji = fs_2, fs_jpim1   ! vector opt. 
     122         
     123                 IF(tmask(ji, jj, 1) < 0.5_wp) CYCLE  
     124         
     125                 ztmp = e1t(ji,jj) * e2t(ji,jj) 
     126 
     127                 zflxn(ji,jj) = min(zflxu1(ji,jj), 0._wp) - max(zflxu1(ji-1,jj),   0._wp) + & 
     128                              & min(zflxv1(ji,jj), 0._wp) - max(zflxv1(ji,  jj-1), 0._wp)  
     129           
     130                 zdep1 = (zflxp(ji,jj) * wdlmt(ji,jj) + zflxn(ji,jj)) * z2dt / ztmp 
     131                 zdep2 = bathy(ji,jj) + sshb(ji,jj) - rn_wadmin  
     132           
     133                 IF(zdep1 > zdep2) THEN 
     134                   zflag = 1 
     135                   zcoef = ( ( zdep2 - rn_wadmine ) * ztmp - zflxn(ji,jj) * z2dt ) / ( zflxp(ji,jj) * z2dt ) 
     136                   zcoef = max(zcoef, 0._wp) 
     137                   IF(zflxu1(ji,  jj  ) >= 0._wp) wdlmt(ji,  jj  ) = zcoef 
     138                   IF(zflxu1(ji-1,jj  ) <  0._wp) wdlmt(ji-1,jj  ) = zcoef 
     139                   IF(zflxv1(ji,  jj  ) >= 0._wp) wdlmt(ji,  jj  ) = zcoef 
     140                   IF(zflxv1(ji,  jj-1) <  0._wp) wdlmt(ji,  jj-1) = zcoef 
     141                 END IF 
     142              END DO ! ji loop 
     143           END DO  ! jj loop 
     144           
     145           CALL lbc_lnk( wdlmt, 'T', 1. ) 
     146           
     147           IF(zflag == 0) EXIT 
     148           
     149           IF(jk1 == nn_waditr) THEN 
     150              IF(zflxu1(ji,  jj  ) >= 0._wp) wdlmt(ji,  jj  ) = 0._wp 
     151              IF(zflxu1(ji-1,jj  ) <  0._wp) wdlmt(ji-1,jj  ) = 0._wp 
     152              IF(zflxv1(ji,  jj  ) >= 0._wp) wdlmt(ji,  jj  ) = 0._wp 
     153              IF(zflxv1(ji,  jj-1) <  0._wp) wdlmt(ji,  jj-1) = 0._wp 
     154           END IF 
     155 
    144156        END DO  ! jk1 loop 
    145157        
     
    147159           DO jj = 2, jpjm1 
    148160              DO ji = fs_2, fs_jpim1   ! vector opt. 
    149                un(ji,jj,jk)   = 0.5_wp * ( sign(1._wp, un(ji,jj,jk))   + 1._wp ) * uwdlmt(ji,jj)  
    150                vn(ji,jj,jk)   = 0.5_wp * ( sign(1._wp, vn(ji,jj,jk))   + 1._wp ) * vwdlmt(ji,jj)  
    151                un(ji-1,jj,jk) = 0.5_wp * ( 1._wp - sign(1._wp, un(ji-1,jj,jk) )) * uwdlmt(ji-1,jj)  
    152                vn(ji,jj-1,jk) = 0.5_wp * ( 1._wp - sign(1._wp, vn(ji,jj-1,jk) )) * vwdlmt(ji,jj-1)  
     161               un(ji,  jj,  jk) = ( sign(0.5_wp, un(ji,  jj,  jk))  + 0.5_wp ) * wdlmt(ji  ,jj)  
     162               vn(ji,  jj,  jk) = ( sign(0.5_wp, vn(ji,  jj,  jk))  + 0.5_wp ) * wdlmt(ji  ,jj)  
     163               un(ji-1,jj,  jk) = (-sign(0.5_wp, un(ji-1,jj,  jk))  + 0.5_wp ) * wdlmt(ji-1,jj)  
     164               vn(ji,  jj-1,jk) = (-sign(0.5_wp, vn(ji,  jj-1,jk))  + 0.5_wp ) * wdlmt(ji,  jj-1)  
    153165              END DO 
    154166           END DO 
     
    161173        ! 
    162174        ! 
    163         CALL wrk_dealloc( jpi, jpj, zflxu, zflxv, zflxu1, zflxv1, uwdlmt, vwdlmt) 
     175        CALL wrk_dealloc( jpi, jpj, zflxp, zflxn, zflxu, zflxv, zflxu1, zflxv1, wdlmt ) 
    164176      ! 
    165177      END IF 
  • branches/2013/dev_r4050_NOC_WaD/NEMOGCM/NEMO/OPA_SRC/par_oce.F90

    r3294 r4375  
    8686   !!--------------------------------------------------------------------- 
    8787#             include "par_AMM_12km.h90" 
     88#elif defined key_pbasin 
     89   !!--------------------------------------------------------------------- 
     90   !!   'key_pbasin':                             Parabolic Basin : PBASIN 
     91   !!--------------------------------------------------------------------- 
     92#             include "par_PBASIN.h90" 
    8893#else 
    8994   !!--------------------------------------------------------------------- 
  • branches/2013/dev_r4050_NOC_WaD/NEMOGCM/NEMO/OPA_SRC/step.F90

    r4355 r4375  
    8787 
    8888      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    89       !  Ocean dynamics : Wetting/Drying Limiter 
    90       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    91                          CALL wad_lmt( kstp )         ! wetting/drying limiter 
    92  
    93       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    9489      ! Update data, open boundaries, surface boundary condition (including sea-ice) 
    9590      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
     
    10095      IF( lk_obc     )   CALL obc_rad( kstp )         ! compute phase velocities at open boundaries 
    10196      IF( lk_bdy     )   CALL bdy_dta( kstp, time_offset=+1 ) ! update dynamic and tracer data at open boundaries 
     97 
     98      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
     99      !  Ocean dynamics : Wetting/Drying Limiter 
     100      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
     101                         CALL wad_lmt( kstp )         ! wetting/drying limiter 
    102102 
    103103      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
Note: See TracChangeset for help on using the changeset viewer.