Changeset 4385


Ignore:
Timestamp:
2014-01-31T14:43:21+01:00 (7 years ago)
Author:
hliu
Message:

1) Modified Wetting and Drying flux limiter part
2) I will apply all the modifications to v3.6beta as soon as everything runs smoothly at this version.

Location:
branches/2013/dev_r4050_NOC_WaD/NEMOGCM/NEMO/OPA_SRC
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • branches/2013/dev_r4050_NOC_WaD/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90

    r4380 r4385  
    398398            ij_bump = jpjdta / 2 + 1                       ! j-index of the basin center 
    399399            r_bump  = 430620._wp                           ! basin radius (meters)        
    400             h_bump  =     50._wp                           ! basin depth  (meters) 
     400            h_bump  =     50._wp                           ! basin centre depth  (meters) 
    401401            h_oce   = gdepw_0(jpk)                         ! background ocean depth (meters) 
    402402            IF(lwp) WRITE(numout,*) '           basin characteristics: ' 
  • branches/2013/dev_r4050_NOC_WaD/NEMOGCM/NEMO/OPA_SRC/DYN/wadlmt.F90

    r4375 r4385  
    4343      !! ** Method  : - Prevent negative depth occurring (Not ready for Agrif)  
    4444      !! 
    45       !! ** Action  : - update: uwdlmt(:,:), vwdlmt(:,:)  
     45      !! ** Action  : - calculate flux limiter and W/D flag 
    4646      !!---------------------------------------------------------------------- 
    4747      INTEGER, INTENT(in) ::   kt       ! ocean time-step index 
    4848      ! 
    4949      INTEGER  ::   ji, jj, jk, jk1     ! dummy loop indices 
    50       INTEGER  ::   zflag, z2dt         ! local scalar 
     50      INTEGER  ::   z2dt, zflag         ! local scalar 
    5151      REAL(wp) ::   zcoef, zdep1, zdep2 ! local scalars 
     52      REAL(wp) ::   zzflxp, zzflxn      ! local scalars 
    5253      REAL(wp) ::   ztmp                ! local scalars 
    53       REAL(wp), POINTER,  DIMENSION(:,:) ::   zflxp,  zflxn   ! specific 2D workspace 
    54       REAL(wp), POINTER,  DIMENSION(:,:) ::   zflxu,  zflxv   ! specific 2D workspace 
    55       REAL(wp), POINTER,  DIMENSION(:,:) ::   zflxu1, zflxv1  ! specific 2D workspace 
    56       REAL(wp), POINTER,  DIMENSION(:,:) ::   wdlmt           !: W/D flux limiters 
     54      INTEGER,  POINTER,  DIMENSION(:,:) ::   iwdflg                   !: W/D cell flag 
     55      REAL(wp), POINTER,  DIMENSION(:,:) ::   zwdlmtu, zwdlmtv         !: W/D flux limiters 
     56      REAL(wp), POINTER,  DIMENSION(:,:) ::   zflxp,  zflxn            ! specific 2D workspace 
     57      REAL(wp), POINTER,  DIMENSION(:,:) ::   zflxu,  zflxv            ! specific 2D workspace 
     58      REAL(wp), POINTER,  DIMENSION(:,:) ::   zflxu1, zflxv1           ! specific 2D workspace 
    5759 
    5860      !!---------------------------------------------------------------------- 
     
    6365      IF(ln_wad) THEN 
    6466 
    65         CALL wrk_alloc( jpi, jpj, zflxp, zflxn, zflxu, zflxv, zflxu1, zflxv1, wdlmt ) 
     67        CALL wrk_alloc( jpi, jpj, zflxp, zflxn, zflxu, zflxv, zflxu1, zflxv1 ) 
     68        CALL wrk_alloc( jpi, jpj, zwdlmtu, zwdlmtv) 
     69        CALL wrk_alloc( jpi, jpj, iwdflg ) 
    6670        ! 
    6771        
     
    7276        IF( neuler == 0 .AND. kt == nit000 )  z2dt = rdt 
    7377        
    74         zflxp(:,:)  = 0._wp 
    75         zflxn(:,:)  = 0._wp 
    76         zflxu(:,:)  = 0._wp 
    77         zflxv(:,:)  = 0._wp 
    78         wdlmt(:,:)  = 1._wp 
     78        zflxp(:,:)   = 0._wp 
     79        zflxn(:,:)   = 0._wp 
     80        zflxu(:,:)   = 0._wp 
     81        zflxv(:,:)   = 0._wp 
     82 
     83        zwdlmtu(:,:)  = 1._wp 
     84        zwdlmtv(:,:)  = 1._wp 
    7985        
    8086        ! Horizontal Flux in u and v direction 
     
    110116 
    111117       
    112         !! start limiter iteration  
     118        !! start limiter iterations  
    113119        DO jk1 = 1, nn_waditr 
    114120        
    115121           zflag = 0     ! flag indicating if any further iteration is needed? 
    116122           
    117            zflxu1(:,:) = zflxu(:,:) * wdlmt(:,:) 
    118            zflxv1(:,:) = zflxv(:,:) * wdlmt(:,:) 
     123           zflxu1(:,:) = zflxu(:,:) * zwdlmtu(:,:) 
     124           zflxv1(:,:) = zflxv(:,:) * zwdlmtv(:,:) 
    119125           
    120126           DO jj = 2, jpjm1 
     
    123129                 IF(tmask(ji, jj, 1) < 0.5_wp) CYCLE  
    124130         
    125                  ztmp = e1t(ji,jj) * e2t(ji,jj) 
     131                 ztmp = e1t(ji,jj) * e2t(ji,jj)     !there must be an array ready for this 
    126132 
    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)  
     133                 zzflxp = max(zflxu1(ji,jj), 0._wp) - min(zflxu1(ji-1,jj),   0._wp) + & 
     134                        & max(zflxv1(ji,jj), 0._wp) - min(zflxv1(ji,  jj-1), 0._wp)  
     135                 zzflxn = min(zflxu1(ji,jj), 0._wp) - max(zflxu1(ji-1,jj),   0._wp) + & 
     136                        & min(zflxv1(ji,jj), 0._wp) - max(zflxv1(ji,  jj-1), 0._wp)  
    129137           
    130                  zdep1 = (zflxp(ji,jj) * wdlmt(ji,jj) + zflxn(ji,jj)) * z2dt / ztmp 
    131                  zdep2 = bathy(ji,jj) + sshb(ji,jj) - rn_wadmin  
     138                 zdep1 = (zzflxp + zzflxn) * z2dt / ztmp 
     139                 zdep2 = bathy(ji,jj) + sshb(ji,jj) - rn_wadmin   ! this one can be moved out of the loops 
    132140           
    133141                 IF(zdep1 > zdep2) THEN 
    134142                   zflag = 1 
    135                    zcoef = ( ( zdep2 - rn_wadmine ) * ztmp - zflxn(ji,jj) * z2dt ) / ( zflxp(ji,jj) * z2dt ) 
     143                   zcoef = ( ( zdep2 - rn_wadmine ) * ztmp - zzflxn * z2dt ) / ( zflxp(ji,jj) * z2dt ) 
    136144                   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 
     145                   iwdflg(ji, jj) = 1 
     146                   IF(zflxu1(ji,  jj) > 0._wp) zwdlmtu(ji  ,jj) = zcoef 
     147                   IF(zflxu1(ji-1,jj) < 0._wp) zwdlmtu(ji-1,jj) = zcoef 
     148                   IF(zflxv1(ji,  jj) > 0._wp) zwdlmtv(ji  ,jj) = zcoef 
     149                   IF(zflxv1(ji,jj-1) < 0._wp) zwdlmtv(ji-1,jj) = zcoef 
    141150                 END IF 
    142151              END DO ! ji loop 
    143152           END DO  ! jj loop 
    144153           
    145            CALL lbc_lnk( wdlmt, 'T', 1. ) 
     154           CALL lbc_lnk( zwdlmtu, 'U', 1. ) 
     155           CALL lbc_lnk( zwdlmtv, 'V', 1. ) 
    146156           
    147157           IF(zflag == 0) EXIT 
    148158           
    149159           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 
     160               
     161              IF(zflxu1(ji,  jj  ) >= 0._wp) zwdlmtu(ji,  jj  ) = 0._wp 
     162              IF(zflxu1(ji-1,jj  ) <  0._wp) zwdlmtu(ji-1,jj  ) = 0._wp 
     163              IF(zflxv1(ji,  jj  ) >= 0._wp) zwdlmtv(ji,  jj  ) = 0._wp 
     164              IF(zflxv1(ji,  jj-1) <  0._wp) zwdlmtv(ji,  jj-1) = 0._wp 
     165 
     166              CALL lbc_lnk( zwdlmtu, 'U', 1. ) 
     167              CALL lbc_lnk( zwdlmtv, 'V', 1. ) 
    154168           END IF 
    155169 
    156170        END DO  ! jk1 loop 
    157171        
    158         DO jk = 1, jpkm1   
    159            DO jj = 2, jpjm1 
    160               DO ji = fs_2, fs_jpim1   ! vector opt. 
    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)  
    165               END DO 
    166            END DO 
    167         END DO 
     172        un(:,:,:) = un(:,:,:) * zwdlmtu(ji, jj)  
     173        vn(:,:,:) = vn(:,:,:) * zwdlmtv(ji, jj)  
     174 
     175        CALL lbc_lnk( un, 'U', -1. ) 
     176        CALL lbc_lnk( vn, 'V', -1. ) 
    168177        
    169178        IF(zflag == 1 .AND. lwp) WRITE(numout,*) 'Need more iterations in wad_lmt!!!' 
     
    173182        ! 
    174183        ! 
    175         CALL wrk_dealloc( jpi, jpj, zflxp, zflxn, zflxu, zflxv, zflxu1, zflxv1, wdlmt ) 
     184        CALL wrk_dealloc( jpi, jpj, zflxp, zflxn, zflxu, zflxv, zflxu1, zflxv1 ) 
     185        CALL wrk_dealloc( jpi, jpj, zwdlmtu, zwdlmtv) 
     186        CALL wrk_dealloc( jpi, jpj, iwdflg ) 
    176187      ! 
    177188      END IF 
Note: See TracChangeset for help on using the changeset viewer.