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

Changeset 8865


Ignore:
Timestamp:
2017-12-01T05:41:32+01:00 (7 years ago)
Author:
deazer
Message:

Moving Changes from CS15mini config into NEMO main src
Updating TEST configs to run wit this version of the code
all sette tests pass at this revision other than AGRIF
Includes updates to dynnxt and tranxt required for 3D rives in WAD case to be conservative.

Next commit will update naming conventions and tidy the code.

Location:
branches/UKMO/ROMS_WAD_7832/NEMOGCM
Files:
1 added
8 deleted
19 edited

Legend:

Unmodified
Added
Removed

  • branches/UKMO/ROMS_WAD_7832/NEMOGCM/CONFIG/TEST_CASES/WAD/EXP00/namelist_cfg

    r8841 r8865  
    88   rn_dx = 1000.0 
    99   rn_dz = 1.0 
    10    nn_wad_test = 7 
     10   nn_wad_test = 1 
    1111/ 
    1212!----------------------------------------------------------------------- 
    … …  
    142142&nambdy        !  unstructured open boundaries 
    143143!----------------------------------------------------------------------- 
    144     ln_bdy         = .true.             
    145     nb_bdy         = 1                    !  number of open boundary sets 
     144    ln_bdy         = .false.             
     145    nb_bdy         = 0                    !  number of open boundary sets 
    146146    ln_coords_file = .false.              !  =T : read bdy coordinates from file 
    147147    cn_coords_file = 'coordinates.bdy.nc' !  bdy coordinates files 
    … …  
    394394!----------------------------------------------------------------------- 
    395395   ln_glo_trd  = .false.   ! (T) global domain averaged diag for T, T^2, KE, and PE 
    396    ln_dyn_trd  = .true.   ! (T) 3D momentum trend output 
     396   ln_dyn_trd  = .false.   ! (T) 3D momentum trend output 
    397397   ln_dyn_mxl  = .FALSE.   ! (T) 2D momentum trends averaged over the mixed layer (not coded yet) 
    398398   ln_vor_trd  = .FALSE.   ! (T) 2D barotropic vorticity trends (not coded yet) 
    399399   ln_KE_trd   = .false.   ! (T) 3D Kinetic   Energy     trends 
    400400   ln_PE_trd   = .false.   ! (T) 3D Potential Energy     trends 
    401    ln_tra_trd  = .true.    ! (T) 3D tracer trend output 
     401   ln_tra_trd  = .false.    ! (T) 3D tracer trend output 
    402402   ln_tra_mxl  = .false.   ! (T) 2D tracer trends averaged over the mixed layer (not coded yet) 
    403403   nn_trd      = 365       !  print frequency (ln_glo_trd=T) (unit=time step) 
    … …  
    444444&namwad  !   Wetting and drying 
    445445!----------------------------------------------------------------------- 
    446    ln_wd             = .false.  ! T/F activation of NOC  wetting and drying scheme 
     446   ln_wd             = .false  ! T/F activation of NOC  wetting and drying scheme 
    447447   ln_rwd            = .true.   ! T/F activation of ROMS wetting and drying scheme 
     448   ln_rwd_bc         = .true. 
    448449   ln_rwd_rmp        = .true.   ! Turn on the limiter 
    449    ln_rwd_bc         = .true.   ! ROMS Baroclinic option 
    450450   ln_wd_diag        = .false.   ! T/F activation of diagnostics for ROMS wd scheme 
    451    rn_wdmin0         =  0.50     ! Rmp value for NOCL option 
    452    rn_wdmin1         =  0.150    ! Minimum wet depth on dried cells 
    453    rn_wdmin2         =  0.001    ! Tolerance of min wet depth on dried cells 
    454    rn_ssh_ref        =  3.0     ! reference level 
     451   rn_wdmin0         =  0.30     ! Minimum wet depth on dried cells 
     452   rn_wdmin1         =  0.2     ! Minimum wet depth on dried cells 
     453   rn_wdmin2         =  0.0001  ! Tolerance of min wet depth on dried cells 
    455454   rn_wdld           =  2.5     ! Land elevation below which wetting/drying is allowed 
    456455   nn_wdit           =   20     ! Max iterations for W/D limiter 
    457456   jn_wd_i           =  22      ! i index of diagnostics 
    458457   jn_wd_j           =   3      ! j index of diagnostics 
    459    jn_wd_k           =   3      ! k index of diagnostics 
    460 / 
     458  
     459/ 

  • branches/UKMO/ROMS_WAD_7832/NEMOGCM/CONFIG/TEST_CASES/WAD/MY_SRC/domain.F90

    r8403 r8865  
    674674      ENDIF 
    675675      ! 
    676       IF( ln_wd  .or. ln_rwd ) THEN              ! wetting and drying domain 
     676      !IF( ln_wd  .or. ln_rwd ) THEN              ! wetting and drying domain 
     677      IF( ln_wd   ) THEN              ! wetting and drying domain 
    677678         CALL iom_rstput( 0, 0, inum, 'ht_0'   , ht_0   , ktype = jp_r8 ) 
    678          CALL iom_rstput( 0, 0, inum, 'ht_wd'  , ht_wd  , ktype = jp_r8 ) 
    679679      ENDIF 
    680680      ! 

  • branches/UKMO/ROMS_WAD_7832/NEMOGCM/CONFIG/TEST_CASES/WAD/MY_SRC/usrdef_istate.F90

    r8403 r8865  
    171171 
    172172! subtract the height of z=0 above the geoid (this allows z = 0 to be higher than all points that may become wet)    
    173       pssh(:,:) = pssh(:,:) - rn_ssh_ref 
     173      pssh(:,:) = pssh(:,:)  - rn_ssh_ref 
    174174 
    175175      ! 

  • branches/UKMO/ROMS_WAD_7832/NEMOGCM/CONFIG/TEST_CASES/WAD/MY_SRC/usrdef_zgr.F90

    r8403 r8865  
    1616   USE oce            ! ocean variables 
    1717   USE dom_oce ,  ONLY: ln_zco, ln_zps, ln_sco   ! ocean space and time domain 
    18    USE dom_oce ,  ONLY: mi0, mi1, nimpp, njmpp,  & 
    19                       & mj0, mj1, glamt, gphit   ! ocean space and time domain 
     18   USE dom_oce ,  ONLY: ht_0,mi0, mi1, nimpp, njmpp,  & 
     19                      & mj0, mj1, glamt, gphit, ht_0   ! ocean space and time domain 
    2020   USE usrdef_nam     ! User defined : namelist variables 
    21    USE wet_dry ,  ONLY: rn_wdmin1, rn_wdmin2, rn_wdld, ht_wd, rn_ht_0,rn_ssh_ref 
     21   USE wet_dry ,  ONLY: rn_wdmin1, rn_wdmin2, rn_wdld, rn_ht_0,rn_ssh_ref 
    2222   ! 
    2323   USE in_out_manager ! I/O manager 
    … …  
    234234               
    235235! increase the depth of the bathymetry by rn_ssh_ref and rn_ht_0    
    236       zht(:,:) = zht(:,:) + rn_ssh_ref + rn_ht_0   
     236      !zht(:,:) = zht(:,:) + rn_ssh_ref + rn_ht_0   
     237      !zht(:,:) = zht(:,:) + rn_ssh_ref + rn_ht_0   
    237238 
    238239      ! at u-point: averaging zht 
    … …  
    285286      IF ( ln_sco ) THEN      !==  s-coordinate  ==!   (terrain-following coordinate) 
    286287         ! 
    287          ht_wd = zht 
     288         ht_0 = zht 
    288289         k_bot(:,:) = jpkm1 * k_top(:,:)  !* bottom ocean = jpk-1 (here use k_top as a land mask) 
    289290         DO jj = 1, jpj 

  • branches/UKMO/ROMS_WAD_7832/NEMOGCM/CONFIG/TEST_CASES/cfg.txt

    r8544 r8865  
    1010WAD7 OPA_SRC 
    1111WAD7LONG OPA_SRC 
     12WAD2 OPA_SRC 

  • branches/UKMO/ROMS_WAD_7832/NEMOGCM/NEMO/OPA_SRC/BDY/bdydyn2d.F90

    r7646 r8865  
    2121   USE phycst          ! physical constants 
    2222   USE lbclnk          ! ocean lateral boundary conditions (or mpp link) 
     23   USE wet_dry         ! Use wet dry to get reference ssh level 
    2324   USE in_out_manager  ! 
    2425 
    … …  
    169170         ii = idx%nbi(jb,igrd) 
    170171         ij = idx%nbj(jb,igrd) 
    171          spgu(ii, ij) = dta%ssh(jb) 
     172         IF( ln_wd .OR. ln_rwd ) THEN 
     173            spgu(ii, ij) = dta%ssh(jb)  - rn_ssh_ref  
     174         ELSE 
     175            spgu(ii, ij) = dta%ssh(jb) 
     176         ENDIF 
    172177      END DO 
    173178 

  • branches/UKMO/ROMS_WAD_7832/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90

    r8841 r8865  
    154154         CALL iom_put( "e3tdef"  , ( ( e3t_n(:,:,:) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100 * tmask(:,:,:) ) ** 2 ) 
    155155 
    156       CALL iom_put( "ssh" , sshn )                 ! sea surface height 
     156      IF( (ln_wd .OR. ln_rwd)) THEN 
     157         CALL iom_put( "ssh" , sshn+rn_ssh_ref )                 ! sea surface height !bring it back to the reference need wad if here 
     158      ELSE 
     159         CALL iom_put( "ssh" , sshn )                 ! sea surface height !bring it back to the reference need wad if here 
     160      ENDIF 
     161 
    157162      IF( iom_use("wetdep") )   &                  ! wet depth 
    158          CALL iom_put( "wetdep" , ht_wd(:,:) + sshn(:,:) ) 
     163         CALL iom_put( "wetdep" , ht_0(:,:) + sshn(:,:) ) 
    159164       
    160165      CALL iom_put( "toce", tsn(:,:,:,jp_tem) )    ! 3D temperature 

  • branches/UKMO/ROMS_WAD_7832/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90

    r7822 r8865  
    674674      ENDIF 
    675675      ! 
    676       IF( ln_wd ) THEN              ! wetting and drying domain 
     676      IF( ln_wd  .or. ln_rwd ) THEN              ! wetting and drying domain 
    677677         CALL iom_rstput( 0, 0, inum, 'ht_0'   , ht_0   , ktype = jp_r8 ) 
    678          CALL iom_rstput( 0, 0, inum, 'ht_wd'  , ht_wd  , ktype = jp_r8 ) 
    679678      ENDIF 
    680679      ! 

  • branches/UKMO/ROMS_WAD_7832/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90

    r7753 r8865  
    876876         ELSE                                   !* Initialize at "rest" 
    877877            ! 
    878             IF( ln_wd .AND. ( cn_cfg == 'wad' ) ) THEN 
     878 
     879! MJB ln_rwd edits start here - these are essential  
     880 
     881            IF( (ln_wd .OR. ln_rwd)) THEN  
     882 
    879883              ! Wetting and drying test case 
    880               CALL usr_def_istate( gdept_b, tmask, tsb, ub, vb, sshb  ) 
    881                        tsn  (:,:,:,:) = tsb (:,:,:,:)       ! set now values from to before ones 
    882                        sshn (:,:)     = sshb(:,:) 
    883                        un   (:,:,:)   = ub  (:,:,:) 
    884                        vn   (:,:,:)   = vb  (:,:,:) 
    885                                                 ! uniform T-S fields and initial ssh slope 
    886                ! needs to be called here and in istate which is called later. 
    887                ! Adjust vertical metrics 
     884               IF( cn_cfg == 'wad' ) THEN 
     885                  CALL usr_def_istate( gdept_b, tmask, tsb, ub, vb, sshb  ) 
     886                  tsn  (:,:,:,:) = tsb (:,:,:,:)       ! set now values from to before ones 
     887                  sshn (:,:)     = sshb(:,:) 
     888                  un   (:,:,:)   = ub  (:,:,:) 
     889                  vn   (:,:,:)   = vb  (:,:,:) 
     890               ELSEIF( ln_wd .or. ln_rwd   ) THEN ! if not test case 
     891                  sshn(:,:) = -rn_ssh_ref 
     892                  sshb(:,:) = -rn_ssh_ref 
     893 
     894                  DO jj = 1, jpj 
     895                     DO ji = 1, jpi 
     896                        IF( ht_0(ji,jj)-rn_ssh_ref <  rn_wdmin1 ) THEN ! if total depth is less than min depth 
     897 
     898                           sshb(ji,jj) = rn_wdmin1 - (ht_0(ji,jj) ) 
     899                           sshn(ji,jj) = rn_wdmin1 - (ht_0(ji,jj) ) 
     900                           ssha(ji,jj) = rn_wdmin1 - (ht_0(ji,jj) ) 
     901                        ENDIF 
     902                     ENDDO 
     903                  ENDDO 
     904               ENDIF !If wad elseif ln_wd or ln_rwd 
     905 
     906               ! Adjust vertical metrics for all wad 
    888907               DO jk = 1, jpk 
    889                   e3t_n(:,:,jk) =  e3t_0(:,:,jk) * ( ht_0(:,:) + sshn(:,:) ) & 
    890                     &                            / ( ht_0(:,:) + 1._wp - ssmask(:,:) ) * tmask(:,:,jk)   & 
    891                     &            + e3t_0(:,:,jk)                               * (1._wp -tmask(:,:,jk)) 
     908                    e3t_n(:,:,jk) =  e3t_0(:,:,jk) * ( ht_0(:,:) - rn_ht_0 + sshn(:,:) ) & 
     909                      &                            / ( ht_0(:,:) + 1._wp - ssmask(:,:) ) * tmask(:,:,jk)   & 
     910                      &            + e3t_0(:,:,jk)                               * (1._wp -tmask(:,:,jk)) 
    892911               END DO 
    893912               e3t_b(:,:,:) = e3t_n(:,:,:) 
    894                ! 
    895             ELSEIF( ln_wd ) THEN 
    896                ! 
    897               DO jj = 1, jpj 
    898                 DO ji = 1, jpi 
    899                   IF( e3t_0(ji,jj,1) <= 0.5_wp * rn_wdmin1 ) THEN 
    900                      ! potential bug 
    901                      ! Warning this assumes 2 layers only over wetting locations. needs investigating 
    902                      e3t_b(ji,jj,:) = 0.5_wp * rn_wdmin1 
    903                      e3t_n(ji,jj,:) = 0.5_wp * rn_wdmin1 
    904                      e3t_a(ji,jj,:) = 0.5_wp * rn_wdmin1 
    905                      sshb(ji,jj) = rn_wdmin1 - ht_wd(ji,jj)           !!gm I don't understand that ! 
    906                      sshn(ji,jj) = rn_wdmin1 - ht_wd(ji,jj) 
    907                      ssha(ji,jj) = rn_wdmin1 - ht_wd(ji,jj) 
    908                   ENDIF 
    909                 ENDDO 
    910               ENDDO 
     913 
     914               DO ji = 1, jpi 
     915                  DO jj = 1, jpj 
     916                     IF ( ht_0(ji,jj) .LE. 0.0 .AND. NINT(ssmask(ji,jj)) .EQ. 1) THEN 
     917                       CALL ctl_stop( 'dom_vvl_rst: ht_0 must be positive at potentially wet points' ) 
     918                     ENDIF 
     919                  END DO  
     920               END DO  
     921  
    911922               ! 
    912923            ELSE 

  • branches/UKMO/ROMS_WAD_7832/NEMOGCM/NEMO/OPA_SRC/DOM/domwri.F90

    r7646 r8865  
    1818   USE dom_oce         ! ocean space and time domain 
    1919   USE phycst ,   ONLY :   rsmall 
    20    USE wet_dry,   ONLY :   ln_wd, ht_wd 
     20   USE wet_dry,   ONLY :   ln_wd, ln_rwd 
    2121   ! 
    2222   USE in_out_manager  ! I/O manager 
    … …  
    198198      ENDIF 
    199199      ! 
    200       IF( ln_wd ) THEN                                          ! wetting and drying domain 
    201          CALL iom_rstput( 0, 0, inum, 'ht_0'   , ht_0   , ktype = jp_r8 ) 
    202          CALL iom_rstput( 0, 0, inum, 'ht_wd'  , ht_wd  , ktype = jp_r8 ) 
    203       ENDIF 
     200      IF( ln_wd .OR. ln_rwd ) CALL iom_rstput( 0, 0, inum, 'ht_0'   , ht_0   , ktype = jp_r8 ) 
     201 
    204202      !                                     ! ============================ 
    205203      CALL iom_close( inum )                !        close the files  

  • branches/UKMO/ROMS_WAD_7832/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90

    r7753 r8865  
    3030   USE usrdef_zgr     ! user defined vertical coordinate system 
    3131   USE depth_e3       ! depth <=> e3 
    32    USE wet_dry, ONLY: ln_wd, ht_wd 
     32   USE wet_dry, ONLY: ln_wd, ln_rwd, rn_ssh_ref 
    3333   ! 
    3434   USE in_out_manager ! I/O manager 
    … …  
    258258      k_bot(:,:) = INT( z2d(:,:) ) 
    259259      ! 
    260       ! bathymetry with orography (wetting and drying only) 
    261       IF( ln_wd )  CALL iom_get( inum, jpdom_data, 'ht_wd' , ht_wd  , lrowattr=ln_use_jattr ) 
     260      ! reference depth for negative bathy (wetting and drying only) 
     261      IF( ln_wd .OR. ln_rwd )  CALL iom_get( inum,  'rn_wd_ref_depth' , rn_ssh_ref  ) 
    262262      ! 
    263263      CALL iom_close( inum ) 

  • branches/UKMO/ROMS_WAD_7832/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90

    r7761 r8865  
    455455           DO ji = 2, jpim1  
    456456             ll_tmp1 = MIN(   sshn(ji,jj)               ,   sshn(ji+1,jj) ) >                & 
    457                   &    MAX( -ht_wd(ji,jj)               , -ht_wd(ji+1,jj) ) .AND.            & 
    458                   &    MAX(   sshn(ji,jj) + ht_wd(ji,jj),   sshn(ji+1,jj) + ht_wd(ji+1,jj) ) & 
     457                  &    MAX( -ht_0(ji,jj)               , -ht_0(ji+1,jj) ) .AND.            & 
     458                  &    MAX(   sshn(ji,jj) + ht_0(ji,jj),   sshn(ji+1,jj) + ht_0(ji+1,jj) ) & 
    459459                  &                                                         > rn_wdmin1 + rn_wdmin2 
    460460             ll_tmp2 = ( ABS( sshn(ji,jj)               -   sshn(ji+1,jj) ) > 1.E-12 ) .AND. (             & 
    461461                  &    MAX(   sshn(ji,jj)               ,   sshn(ji+1,jj) ) >                & 
    462                   &    MAX( -ht_wd(ji,jj)               , -ht_wd(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 ) 
     462                  &    MAX( -ht_0(ji,jj)               , -ht_0(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 ) 
    463463 
    464464             IF(ll_tmp1) THEN 
    … …  
    466466             ELSE IF(ll_tmp2) THEN 
    467467               ! no worries about  sshn(ji+1,jj) -  sshn(ji  ,jj) = 0, it won't happen ! here 
    468                zcpx(ji,jj) = ABS( (sshn(ji+1,jj) + ht_wd(ji+1,jj) - sshn(ji,jj) - ht_wd(ji,jj)) & 
     468               zcpx(ji,jj) = ABS( (sshn(ji+1,jj) + ht_0(ji+1,jj) - sshn(ji,jj) - ht_0(ji,jj)) & 
    469469                           &    / (sshn(ji+1,jj) -  sshn(ji  ,jj)) ) 
    470470             ELSE 
    … …  
    473473       
    474474             ll_tmp1 = MIN(   sshn(ji,jj)               ,   sshn(ji,jj+1) ) >                & 
    475                   &    MAX( -ht_wd(ji,jj)               , -ht_wd(ji,jj+1) ) .AND.            & 
    476                   &    MAX(   sshn(ji,jj) + ht_wd(ji,jj),   sshn(ji,jj+1) + ht_wd(ji,jj+1) ) & 
     475                  &    MAX( -ht_0(ji,jj)               , -ht_0(ji,jj+1) ) .AND.            & 
     476                  &    MAX(   sshn(ji,jj) + ht_0(ji,jj),   sshn(ji,jj+1) + ht_0(ji,jj+1) ) & 
    477477                  &                                                         > rn_wdmin1 + rn_wdmin2 
    478478             ll_tmp2 = ( ABS( sshn(ji,jj)               -   sshn(ji,jj+1) ) > 1.E-12 ) .AND. (             & 
    479479                  &    MAX(   sshn(ji,jj)               ,   sshn(ji,jj+1) ) >                & 
    480                   &    MAX( -ht_wd(ji,jj)               , -ht_wd(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 ) 
     480                  &    MAX( -ht_0(ji,jj)               , -ht_0(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 ) 
    481481 
    482482             IF(ll_tmp1) THEN 
    … …  
    484484             ELSE IF(ll_tmp2) THEN 
    485485               ! no worries about  sshn(ji,jj+1) -  sshn(ji,jj  ) = 0, it won't happen ! here 
    486                zcpy(ji,jj) = ABS( (sshn(ji,jj+1) + ht_wd(ji,jj+1) - sshn(ji,jj) - ht_wd(ji,jj)) & 
     486               zcpy(ji,jj) = ABS( (sshn(ji,jj+1) + ht_0(ji,jj+1) - sshn(ji,jj) - ht_0(ji,jj)) & 
    487487                           &    / (sshn(ji,jj+1) -  sshn(ji,jj  )) ) 
    488488             ELSE 
    … …  
    707707           DO ji = 2, jpim1  
    708708             ll_tmp1 = MIN(   sshn(ji,jj)               ,   sshn(ji+1,jj) ) >                & 
    709                   &    MAX( -ht_wd(ji,jj)               , -ht_wd(ji+1,jj) ) .AND.            & 
    710                   &    MAX(   sshn(ji,jj) + ht_wd(ji,jj),   sshn(ji+1,jj) + ht_wd(ji+1,jj) ) & 
     709                  &    MAX( -ht_0(ji,jj)               , -ht_0(ji+1,jj) ) .AND.            & 
     710                  &    MAX(   sshn(ji,jj) + ht_0(ji,jj),   sshn(ji+1,jj) + ht_0(ji+1,jj) ) & 
    711711                  &                                                         > rn_wdmin1 + rn_wdmin2 
    712712             ll_tmp2 = ( ABS( sshn(ji,jj)               -   sshn(ji+1,jj) ) > 1.E-12 ) .AND. (             & 
    713713                  &    MAX(   sshn(ji,jj)               ,   sshn(ji+1,jj) ) >                & 
    714                   &    MAX( -ht_wd(ji,jj)               , -ht_wd(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 ) 
     714                  &    MAX( -ht_0(ji,jj)               , -ht_0(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 ) 
    715715 
    716716             IF(ll_tmp1) THEN 
    … …  
    718718             ELSE IF(ll_tmp2) THEN 
    719719               ! no worries about  sshn(ji+1,jj) -  sshn(ji  ,jj) = 0, it won't happen ! here 
    720                zcpx(ji,jj) = ABS( (sshn(ji+1,jj) + ht_wd(ji+1,jj) - sshn(ji,jj) - ht_wd(ji,jj)) & 
     720               zcpx(ji,jj) = ABS( (sshn(ji+1,jj) + ht_0(ji+1,jj) - sshn(ji,jj) - ht_0(ji,jj)) & 
    721721                           &    / (sshn(ji+1,jj) -  sshn(ji  ,jj)) ) 
    722722             ELSE 
    … …  
    725725       
    726726             ll_tmp1 = MIN(   sshn(ji,jj)               ,   sshn(ji,jj+1) ) >                & 
    727                   &    MAX( -ht_wd(ji,jj)               , -ht_wd(ji,jj+1) ) .AND.            & 
    728                   &    MAX(   sshn(ji,jj) + ht_wd(ji,jj),   sshn(ji,jj+1) + ht_wd(ji,jj+1) ) & 
     727                  &    MAX( -ht_0(ji,jj)               , -ht_0(ji,jj+1) ) .AND.            & 
     728                  &    MAX(   sshn(ji,jj) + ht_0(ji,jj),   sshn(ji,jj+1) + ht_0(ji,jj+1) ) & 
    729729                  &                                                         > rn_wdmin1 + rn_wdmin2 
    730730             ll_tmp2 = ( ABS( sshn(ji,jj)               -   sshn(ji,jj+1) ) > 1.E-12 ) .AND. (             & 
    731731                  &    MAX(   sshn(ji,jj)               ,   sshn(ji,jj+1) ) >                & 
    732                   &    MAX( -ht_wd(ji,jj)               , -ht_wd(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 ) 
     732                  &    MAX( -ht_0(ji,jj)               , -ht_0(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 ) 
    733733 
    734734             IF(ll_tmp1) THEN 
    … …  
    736736             ELSE IF(ll_tmp2) THEN 
    737737               ! no worries about  sshn(ji,jj+1) -  sshn(ji,jj  ) = 0, it won't happen ! here 
    738                zcpy(ji,jj) = ABS( (sshn(ji,jj+1) + ht_wd(ji,jj+1) - sshn(ji,jj) - ht_wd(ji,jj)) & 
     738               zcpy(ji,jj) = ABS( (sshn(ji,jj+1) + ht_0(ji,jj+1) - sshn(ji,jj) - ht_0(ji,jj)) & 
    739739                           &    / (sshn(ji,jj+1) -  sshn(ji,jj  )) ) 
    740740             ELSE 
    … …  
    10061006           DO ji = 2, jpim1  
    10071007             ll_tmp1 = MIN(   sshn(ji,jj)               ,   sshn(ji+1,jj) ) >                & 
    1008                   &    MAX( -ht_wd(ji,jj)               , -ht_wd(ji+1,jj) ) .AND.            & 
    1009                   &    MAX(   sshn(ji,jj) + ht_wd(ji,jj),   sshn(ji+1,jj) + ht_wd(ji+1,jj) ) & 
     1008                  &    MAX( -ht_0(ji,jj)               , -ht_0(ji+1,jj) ) .AND.            & 
     1009                  &    MAX(   sshn(ji,jj) + ht_0(ji,jj),   sshn(ji+1,jj) + ht_0(ji+1,jj) ) & 
    10101010                  &                                                         > rn_wdmin1 + rn_wdmin2 
    10111011             ll_tmp2 = ( ABS( sshn(ji,jj)               -   sshn(ji+1,jj) ) > 1.E-12 ) .AND. (             & 
    10121012                  &    MAX(   sshn(ji,jj)               ,   sshn(ji+1,jj) ) >                & 
    1013                   &    MAX( -ht_wd(ji,jj)               , -ht_wd(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 ) 
     1013                  &    MAX( -ht_0(ji,jj)               , -ht_0(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 ) 
    10141014 
    10151015             IF(ll_tmp1) THEN 
    … …  
    10171017             ELSE IF(ll_tmp2) THEN 
    10181018               ! no worries about  sshn(ji+1,jj) -  sshn(ji  ,jj) = 0, it won't happen ! here 
    1019                zcpx(ji,jj) = ABS( (sshn(ji+1,jj) + ht_wd(ji+1,jj) - sshn(ji,jj) - ht_wd(ji,jj)) & 
     1019               zcpx(ji,jj) = ABS( (sshn(ji+1,jj) + ht_0(ji+1,jj) - sshn(ji,jj) - ht_0(ji,jj)) & 
    10201020                           &    / (sshn(ji+1,jj) -  sshn(ji  ,jj)) ) 
     1021               
     1022                zcpx(ji,jj) = max(min( zcpx(ji,jj) , 1.0_wp),0.0_wp) 
    10211023             ELSE 
    10221024               zcpx(ji,jj) = 0._wp 
    … …  
    10241026       
    10251027             ll_tmp1 = MIN(   sshn(ji,jj)               ,   sshn(ji,jj+1) ) >                & 
    1026                   &    MAX( -ht_wd(ji,jj)               , -ht_wd(ji,jj+1) ) .AND.            & 
    1027                   &    MAX(   sshn(ji,jj) + ht_wd(ji,jj),   sshn(ji,jj+1) + ht_wd(ji,jj+1) ) & 
     1028                  &    MAX( -ht_0(ji,jj)               , -ht_0(ji,jj+1) ) .AND.            & 
     1029                  &    MAX(   sshn(ji,jj) + ht_0(ji,jj),   sshn(ji,jj+1) + ht_0(ji,jj+1) ) & 
    10281030                  &                                                         > rn_wdmin1 + rn_wdmin2 
    10291031             ll_tmp2 = ( ABS( sshn(ji,jj)               -   sshn(ji,jj+1) ) > 1.E-12 ) .AND. (             & 
    10301032                  &    MAX(   sshn(ji,jj)               ,   sshn(ji,jj+1) ) >                & 
    1031                   &    MAX( -ht_wd(ji,jj)               , -ht_wd(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 ) 
     1033                  &    MAX( -ht_0(ji,jj)               , -ht_0(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 ) 
    10321034 
    10331035             IF(ll_tmp1) THEN 
    … …  
    10351037             ELSE IF(ll_tmp2) THEN 
    10361038               ! no worries about  sshn(ji,jj+1) -  sshn(ji,jj  ) = 0, it won't happen ! here 
    1037                zcpy(ji,jj) = ABS( (sshn(ji,jj+1) + ht_wd(ji,jj+1) - sshn(ji,jj) - ht_wd(ji,jj)) & 
     1039               zcpy(ji,jj) = ABS( (sshn(ji,jj+1) + ht_0(ji,jj+1) - sshn(ji,jj) - ht_0(ji,jj)) & 
    10381040                           &    / (sshn(ji,jj+1) -  sshn(ji,jj  )) ) 
     1041                zcpy(ji,jj) = max(min( zcpy(ji,jj) , 1.0_wp),0.0_wp) 
     1042 
    10391043             ELSE 
    10401044               zcpy(ji,jj) = 0._wp 
    … …  
    12281232               ENDIF 
    12291233               IF( ln_wd ) THEN 
    1230                   zdpdx1 = zdpdx1 * zcpx(ji,jj) 
    1231                   zdpdx2 = zdpdx2 * zcpx(ji,jj) 
     1234                  zdpdx1 = zdpdx1 * zcpx(ji,jj) * wdrampu(ji,jj) 
     1235                  zdpdx2 = zdpdx2 * zcpx(ji,jj) * wdrampu(ji,jj) 
    12321236               ENDIF 
    12331237               ua(ji,jj,jk) = ua(ji,jj,jk) + (zdpdx1 + zdpdx2) * umask(ji,jj,jk)  
    … …  
    12871291               ENDIF 
    12881292               IF( ln_wd ) THEN 
    1289                   zdpdy1 = zdpdy1 * zcpy(ji,jj) 
    1290                   zdpdy2 = zdpdy2 * zcpy(ji,jj) 
     1293                  zdpdy1 = zdpdy1 * zcpy(ji,jj) * wdrampv(ji,jj)  
     1294                  zdpdy2 = zdpdy2 * zcpy(ji,jj) * wdrampv(ji,jj)  
    12911295               ENDIF 
    12921296 

  • branches/UKMO/ROMS_WAD_7832/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90

    r7753 r8865  
    2828   USE dom_oce        ! ocean space and time domain 
    2929   USE sbc_oce        ! Surface boundary condition: ocean fields 
     30   USE sbcrnf         ! river runoffs 
    3031   USE phycst         ! physical constants 
    3132   USE dynadv         ! dynamics: vector invariant versus flux form 
    … …  
    220221               IF ( .NOT. ln_isf ) THEN   ! if no ice shelf melting 
    221222                  e3t_b(:,:,1) = e3t_b(:,:,1) - zcoef * ( emp_b(:,:) - emp(:,:) & 
    222                                  &                      - rnf_b(:,:) + rnf(:,:) ) * tmask(:,:,1) 
     223                                 &                      ) * tmask(:,:,1) 
     224                  IF( ln_rnf_depth ) THEN  
     225                     DO jk = 1, jpkm1 ! Deal with Rivers separetely, as can be through depth too, not sure for ice shelf case yet 
     226                        DO jj = 1, jpj 
     227                           DO ji = 1, jpi 
     228                              IF( mikt(ji,jj) <= jk .and. jk <=  nk_rnf(ji,jj)  ) THEN 
     229                                 e3t_b(ji,jj,jk) = e3t_b(ji,jj,jk) - zcoef *  (rnf_b(ji,jj) - rnf(ji,jj))*(e3t_n(ji,jj,jk)/h_rnf(ji,jj)  )*tmask(ji,jj,jk) 
     230                              ENDIF 
     231                           ENDDO 
     232                        ENDDO 
     233                     ENDDO 
     234                  ELSE 
     235                      e3t_b(:,:,1) = e3t_b(:,:,1) - zcoef *  (rnf_b(ji,jj) - rnf(ji,jj))*tmask(ji,jj,1) 
     236                  ENDIF 
    223237               ELSE                     ! if ice shelf melting 
    224238                  DO jj = 1, jpj 

  • branches/UKMO/ROMS_WAD_7832/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90

    r7831 r8865  
    11MODULE dynspg_ts 
     2 
     3   !! Includes ROMS wd scheme with diagnostic outputs ; un and ua updates are commented out !  
     4 
    25   !!====================================================================== 
    36   !!                   ***  MODULE  dynspg_ts  *** 
    … …  
    150153      REAL(wp) ::   zhura, zhvra          !   -      - 
    151154      REAL(wp) ::   za0, za1, za2, za3    !   -      - 
     155      REAL(wp) ::   zwdramp                     ! local scalar - only used if ln_rwd = .True.  
     156 
     157      INTEGER  :: iwdg, jwdg, kwdg   ! short-hand values for the indices of the output point 
     158 
    152159      ! 
    153160      REAL(wp), POINTER, DIMENSION(:,:) :: zsshp2_e 
    … …  
    158165      REAL(wp), POINTER, DIMENSION(:,:) :: zhf 
    159166      REAL(wp), POINTER, DIMENSION(:,:) :: zcpx, zcpy                 ! Wetting/Dying gravity filter coef. 
     167      REAL(wp), POINTER, DIMENSION(:,:) :: ztwdmask, zuwdmask, zvwdmask ! ROMS wetting and drying masks at t,u,v points 
     168      REAL(wp), POINTER, DIMENSION(:,:) :: zuwdav2, zvwdav2    ! averages over the sub-steps of zuwdmask and zvwdmask 
    160169      !!---------------------------------------------------------------------- 
    161170      ! 
    … …  
    170179      CALL wrk_alloc( jpi,jpj,   zhf ) 
    171180      IF( ln_wd ) CALL wrk_alloc( jpi, jpj, zcpx, zcpy ) 
     181      IF( ln_rwd ) CALL wrk_alloc( jpi, jpj, ztwdmask, zuwdmask, zvwdmask, zuwdav2, zvwdav2) 
     182 
     183      IF ( ln_wd_diag ) THEN  
     184         iwdg = jn_wd_i ; jwdg = jn_wd_j ; kwdg = jn_wd_k 
     185         WRITE(numout,*) 'kt, iwdg, jwdg, kwdg = ', kt, iwdg, jwdg, kwdg  
     186      END IF  
     187 
    172188      ! 
    173189      zmdi=1.e+20                               !  missing data indicator for masking 
    … …  
    178194      z1_2  = 0.5_wp      
    179195      zraur = 1._wp / rau0 
     196      zwdramp = 1._wp / rn_wdmin1                 ! simplest ramp  
     197!     zwdramp = 1._wp / (rn_wdmin2 - rn_wdmin1)   ! more general ramp 
    180198      !                                            ! reciprocal of baroclinic time step  
    181199      IF( kt == nit000 .AND. neuler == 0 ) THEN   ;   z2dt_bf =          rdt 
    … …  
    403421              DO ji = 2, jpim1  
    404422                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) ) & 
     423                     &    MAX( -ht_0(ji,jj)               , -ht_0(ji+1,jj) ) .AND.            & 
     424                     &    MAX(   sshn(ji,jj) + ht_0(ji,jj),   sshn(ji+1,jj) + ht_0(ji+1,jj) ) & 
    407425                     &                                                         > rn_wdmin1 + rn_wdmin2 
    408426                ll_tmp2 = ( ABS( sshn(ji+1,jj)             -   sshn(ji  ,jj))  > 1.E-12 ).AND.( & 
    409427                     &    MAX(   sshn(ji,jj)               ,   sshn(ji+1,jj) ) >                & 
    410                      &    MAX( -ht_wd(ji,jj)               , -ht_wd(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 ) 
     428                     &    MAX( -ht_0(ji,jj)               , -ht_0(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 ) 
    411429    
    412430                IF(ll_tmp1) THEN 
    … …  
    414432                ELSE IF(ll_tmp2) THEN 
    415433                  ! 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)) & 
     434                  zcpx(ji,jj) = ABS( (sshn(ji+1,jj) + ht_0(ji+1,jj) - sshn(ji,jj) - ht_0(ji,jj)) & 
    417435                              &    / (sshn(ji+1,jj) -  sshn(ji  ,jj)) ) 
     436                  zcpx(ji,jj) = max(min( zcpx(ji,jj) , 1.0_wp),0.0_wp) 
     437 
    418438                ELSE 
    419439                  zcpx(ji,jj) = 0._wp 
    … …  
    421441          
    422442                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) ) & 
     443                     &    MAX( -ht_0(ji,jj)               , -ht_0(ji,jj+1) ) .AND.            & 
     444                     &    MAX(   sshn(ji,jj) + ht_0(ji,jj),   sshn(ji,jj+1) + ht_0(ji,jj+1) ) & 
    425445                     &                                                         > rn_wdmin1 + rn_wdmin2 
    426446                ll_tmp2 = ( ABS( sshn(ji,jj)               -   sshn(ji,jj+1))  > 1.E-12 ).AND.( & 
    427447                     &    MAX(   sshn(ji,jj)               ,   sshn(ji,jj+1) ) >                & 
    428                      &    MAX( -ht_wd(ji,jj)               , -ht_wd(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 ) 
     448                     &    MAX( -ht_0(ji,jj)               , -ht_0(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 ) 
    429449    
    430450                IF(ll_tmp1) THEN 
    … …  
    432452                ELSE IF(ll_tmp2) THEN 
    433453                  ! 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)) & 
     454                  zcpy(ji,jj) = ABS( (sshn(ji,jj+1) + ht_0(ji,jj+1) - sshn(ji,jj) - ht_0(ji,jj)) & 
    435455                              &    / (sshn(ji,jj+1) -  sshn(ji,jj  )) ) 
     456                  zcpy(ji,jj) = max(min( zcpy(ji,jj) , 1.0_wp),0.0_wp) 
     457 
    436458                ELSE 
    437459                  zcpy(ji,jj) = 0._wp 
    … …  
    443465              DO ji = 2, jpim1 
    444466                 zu_trd(ji,jj) = zu_trd(ji,jj) - grav * ( sshn(ji+1,jj  ) - sshn(ji  ,jj ) )   & 
    445                         &                        * r1_e1u(ji,jj) * zcpx(ji,jj) 
     467                        &                        * r1_e1u(ji,jj) * zcpx(ji,jj)  * wdrampu(ji,jj)  !jth 
    446468                 zv_trd(ji,jj) = zv_trd(ji,jj) - grav * ( sshn(ji  ,jj+1) - sshn(ji  ,jj ) )   & 
    447                         &                        * r1_e2v(ji,jj) * zcpy(ji,jj) 
     469                        &                        * r1_e2v(ji,jj) * zcpy(ji,jj)  * wdrampv(ji,jj)  !jth 
     470 
    448471              END DO 
    449472           END DO 
    … …  
    491514      ! Note that the "unclipped" bottom friction parameter is used even with explicit drag 
    492515      IF( ln_wd ) THEN 
    493         zu_frc(:,:) = zu_frc(:,:) + MAX(r1_hu_n(:,:) * bfrua(:,:),-1._wp / rdtbt) * zwx(:,:) 
    494         zv_frc(:,:) = zv_frc(:,:) + MAX(r1_hv_n(:,:) * bfrva(:,:),-1._wp / rdtbt) * zwy(:,:) 
     516        zu_frc(:,:) = zu_frc(:,:) + MAX(r1_hu_n(:,:) * bfrua(:,:),-1._wp / rdtbt) * zwx(:,:) *  wdrampu(ji,jj) 
     517        zv_frc(:,:) = zv_frc(:,:) + MAX(r1_hv_n(:,:) * bfrva(:,:),-1._wp / rdtbt) * zwy(:,:) *  wdrampv(ji,jj) 
    495518      ELSE 
    496519        zu_frc(:,:) = zu_frc(:,:) + r1_hu_n(:,:) * bfrua(:,:) * zwx(:,:) 
    … …  
    627650      vn_adv(:,:) = 0._wp 
    628651      !                                             ! ==================== ! 
     652 
     653      IF (ln_rwd) THEN 
     654         zuwdmask(:,:) = 0._wp  ! set to zero for definiteness (not sure this is necessary)  
     655         zvwdmask(:,:) = 0._wp  !  
     656         zuwdav2(:,:) =  0._wp  
     657         zvwdav2(:,:) =  0._wp    
     658      END IF  
     659 
     660 
    629661      DO jn = 1, icycle                             !  sub-time-step loop  ! 
    630662         !                                          ! ==================== ! 
    … …  
    654686            ! Extrapolate Sea Level at step jit+0.5: 
    655687            zsshp2_e(:,:) = za1 * sshn_e(:,:)  + za2 * sshb_e(:,:) + za3 * sshbb_e(:,:) 
    656             ! 
     688             
     689            ! set wetting & drying mask at tracer points for this barotropic sub-step  
     690            IF ( ln_rwd ) THEN  
     691 
     692               IF ( ln_rwd_rmp ) THEN  
     693                  DO jj = 1, jpj                                  
     694                     DO ji = 1, jpi   ! vector opt.   
     695                        IF ( zsshp2_e(ji,jj) + ht_0(ji,jj) >  2._wp * rn_wdmin1 ) THEN  
     696!                        IF ( zsshp2_e(ji,jj) + ht_0(ji,jj) >  rn_wdmin2 ) THEN  
     697                           ztwdmask(ji,jj) = 1._wp 
     698                        ELSE IF ( zsshp2_e(ji,jj) + ht_0(ji,jj) >  rn_wdmin1 ) THEN 
     699                           ztwdmask(ji,jj) = (tanh(50._wp*( ( zsshp2_e(ji,jj) + ht_0(ji,jj) -  rn_wdmin1 )/rn_wdmin1)) )  
     700                        ELSE  
     701                           ztwdmask(ji,jj) = 0._wp 
     702                        END IF 
     703                     END DO 
     704                  END DO 
     705               ELSE 
     706                  DO jj = 1, jpj                                  
     707                     DO ji = 1, jpi   ! vector opt.   
     708                        IF ( zsshp2_e(ji,jj) + ht_0(ji,jj) >  rn_wdmin1 ) THEN  
     709                           ztwdmask(ji,jj) = 1._wp 
     710                        ELSE  
     711                           ztwdmask(ji,jj) = 0._wp 
     712                        END IF 
     713                     END DO 
     714                  END DO 
     715               END IF  
     716 
     717               IF ( ln_wd_diag ) WRITE(numout,*) 'kt, jn = ', kt, jn  
     718               IF ( ln_wd_diag ) WRITE(numout, *) 'zsshp2_e: (i,j), (i+1,j), (i,j+1) = ', zsshp2_e(iwdg,jwdg), zsshp2_e(iwdg+1,jwdg), zsshp2_e(iwdg,jwdg+1) 
     719               IF ( ln_wd_diag ) WRITE(numout, *) 'ht_0:     (i,j), (i+1,j), (i,j+1) = ', ht_0(iwdg,jwdg), ht_0(iwdg+1,jwdg), (iwdg,jwdg+1) 
     720               IF ( ln_wd_diag ) WRITE(numout, *) 'ztwdmask: (i,j), (i+1,j), (i,j+1) = ', ztwdmask(iwdg,jwdg), ztwdmask(iwdg+1,jwdg), ztwdmask(iwdg,jwdg+1)  
     721            END IF  
     722            
     723 
    657724            DO jj = 2, jpjm1                                    ! Sea Surface Height at u- & v-points 
    658725               DO ji = 2, fs_jpim1   ! Vector opt. 
    … …  
    707774#endif 
    708775         IF( ln_wd ) CALL wad_lmt_bt(zwx, zwy, sshn_e, zssh_frc, rdtbt) 
    709          ! 
     776 
     777         IF ( ln_rwd ) THEN  
     778 
     779            IF ( ln_wd_diag ) THEN  
     780               WRITE(numout, *) 'zwx: (i,j), (i+1,j) = ', zwx(iwdg,jwdg), zwx(iwdg+1,jwdg) 
     781               WRITE(numout, *) 'zwy: (i,j), (i,j+1) = ', zwy(iwdg,jwdg), zwx(iwdg,jwdg+1) 
     782            END IF  
     783 
     784! un_e and vn_e are set to zero at faces where the direction of the flow is from dry cells  
     785 
     786            DO jj = 1, jpjm1                                  
     787               DO ji = 1, jpim1    
     788                  IF ( zwx(ji,jj) > 0.0 ) THEN  
     789                     zuwdmask(ji, jj) = ztwdmask(ji  ,jj)  
     790                  ELSE  
     791                     zuwdmask(ji, jj) = ztwdmask(ji+1,jj)   
     792                  END IF  
     793                  zwx(ji, jj) = zuwdmask(ji,jj)*zwx(ji, jj) 
     794                  un_e(ji,jj) = zuwdmask(ji,jj)*un_e(ji,jj) 
     795 
     796                  IF ( zwy(ji,jj) > 0.0 ) THEN 
     797                     zvwdmask(ji, jj) = ztwdmask(ji, jj  ) 
     798                  ELSE  
     799                     zvwdmask(ji, jj) = ztwdmask(ji, jj+1)   
     800                  END IF  
     801                  zwy(ji, jj) = zvwdmask(ji,jj)*zwy(ji,jj)  
     802                  vn_e(ji,jj) = zvwdmask(ji,jj)*vn_e(ji,jj) 
     803               END DO 
     804            END DO 
     805 
     806            IF ( ln_wd_diag ) THEN  
     807               WRITE(numout, *) 'zuwdmask: (i,j)     = ', zuwdmask(iwdg,jwdg) 
     808               WRITE(numout, *) 'zwx: (i,j)          = ', zwx(iwdg,jwdg) 
     809               WRITE(numout, *) 'e2u: (i,j)          = ', e2u(iwdg,jwdg) 
     810               WRITE(numout, *) 'ua_e: (i,j)         = ', ua_e(iwdg,jwdg) 
     811               WRITE(numout, *) 'un_e: (i,j)         = ', un_e(iwdg,jwdg) 
     812               WRITE(numout, *) 'zhup2_e: (i,j)      = ', zhup2_e(iwdg,jwdg) 
     813               WRITE(numout, *) 'zvwdmask: (i,j)     = ', zvwdmask(iwdg,jwdg) 
     814               WRITE(numout, *) 'zwy: (i,j)          = ', zwy(iwdg,jwdg)  
     815            END IF  
     816 
     817         END IF     
     818          
    710819         ! Sum over sub-time-steps to compute advective velocities 
    711820         za2 = wgtbtp2(jn) 
    712821         un_adv(:,:) = un_adv(:,:) + za2 * zwx(:,:) * r1_e2u(:,:) 
    713822         vn_adv(:,:) = vn_adv(:,:) + za2 * zwy(:,:) * r1_e1v(:,:) 
    714          ! 
     823          
     824         ! sum over sub-time-steps to decide which baroclinic velocities to set to zero (zuwdav2 is only used when ln_rwd_bc = True)  
     825         IF ( ln_rwd_bc ) THEN 
     826            zuwdav2(:,:) = zuwdav2(:,:) + za2 * zuwdmask(:,:) 
     827            zvwdav2(:,:) = zvwdav2(:,:) + za2 * zvwdmask(:,:) 
     828 
     829            IF ( ln_wd_diag ) THEN  
     830               WRITE(numout, *) 'za2, r1_e2u(i,j)     = ', za2, r1_e2u(iwdg,jwdg)  
     831               WRITE(numout, *) 'un_adv: (i,j)        = ', un_adv(iwdg,jwdg) 
     832               WRITE(numout, *) 'zuwdav2: (i,j)       = ', zuwdav2(iwdg,jwdg) 
     833               WRITE(numout, *) 'zvwdav2: (i,j)       = ', zvwdav2(iwdg,jwdg) 
     834            END IF  
     835 
     836         END IF  
     837 
    715838         ! Set next sea level: 
    716839         DO jj = 2, jpjm1                                  
    … …  
    770893              DO ji = 2, jpim1  
    771894                ll_tmp1 = MIN( zsshp2_e(ji,jj)               , zsshp2_e(ji+1,jj) ) >                & 
    772                      &    MAX(   -ht_wd(ji,jj)               ,   -ht_wd(ji+1,jj) ) .AND.            & 
    773                      &    MAX( zsshp2_e(ji,jj) + ht_wd(ji,jj), zsshp2_e(ji+1,jj) + ht_wd(ji+1,jj) ) & 
     895                     &    MAX(   -ht_0(ji,jj)               ,   -ht_0(ji+1,jj) ) .AND.            & 
     896                     &    MAX( zsshp2_e(ji,jj) + ht_0(ji,jj), zsshp2_e(ji+1,jj) + ht_0(ji+1,jj) ) & 
    774897                     &                                                             > rn_wdmin1 + rn_wdmin2 
    775898                ll_tmp2 = (ABS(zsshp2_e(ji,jj)               - zsshp2_e(ji+1,jj))  > 1.E-12 ).AND.( & 
    776899                     &    MAX( zsshp2_e(ji,jj)               , zsshp2_e(ji+1,jj) ) >                & 
    777                      &    MAX(   -ht_wd(ji,jj)               ,   -ht_wd(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 ) 
     900                     &    MAX(   -ht_0(ji,jj)               ,   -ht_0(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 ) 
    778901    
    779902                IF(ll_tmp1) THEN 
    … …  
    781904                ELSE IF(ll_tmp2) THEN 
    782905                  ! no worries about  zsshp2_e(ji+1,jj) - zsshp2_e(ji  ,jj) = 0, it won't happen ! here 
    783                   zcpx(ji,jj) = ABS( (zsshp2_e(ji+1,jj) +    ht_wd(ji+1,jj) - zsshp2_e(ji,jj) - ht_wd(ji,jj)) & 
     906                  zcpx(ji,jj) = ABS( (zsshp2_e(ji+1,jj) +    ht_0(ji+1,jj) - zsshp2_e(ji,jj) - ht_0(ji,jj)) & 
    784907                              &    / (zsshp2_e(ji+1,jj) - zsshp2_e(ji  ,jj)) ) 
    785908                ELSE 
    … …  
    788911          
    789912                ll_tmp1 = MIN( zsshp2_e(ji,jj)               , zsshp2_e(ji,jj+1) ) >                & 
    790                      &    MAX(   -ht_wd(ji,jj)               ,   -ht_wd(ji,jj+1) ) .AND.            & 
    791                      &    MAX( zsshp2_e(ji,jj) + ht_wd(ji,jj), zsshp2_e(ji,jj+1) + ht_wd(ji,jj+1) ) & 
     913                     &    MAX(   -ht_0(ji,jj)               ,   -ht_0(ji,jj+1) ) .AND.            & 
     914                     &    MAX( zsshp2_e(ji,jj) + ht_0(ji,jj), zsshp2_e(ji,jj+1) + ht_0(ji,jj+1) ) & 
    792915                     &                                                             > rn_wdmin1 + rn_wdmin2 
    793916                ll_tmp2 = (ABS(zsshp2_e(ji,jj)               - zsshp2_e(ji,jj+1))  > 1.E-12 ).AND.( & 
    794917                     &    MAX( zsshp2_e(ji,jj)               , zsshp2_e(ji,jj+1) ) >                & 
    795                      &    MAX(   -ht_wd(ji,jj)               ,   -ht_wd(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 ) 
     918                     &    MAX(   -ht_0(ji,jj)               ,   -ht_0(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 ) 
    796919    
    797920                IF(ll_tmp1) THEN 
    … …  
    799922                ELSE IF(ll_tmp2) THEN 
    800923                  ! no worries about  zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj  ) = 0, it won't happen ! here 
    801                   zcpy(ji,jj) = ABS( (zsshp2_e(ji,jj+1) +    ht_wd(ji,jj+1) - zsshp2_e(ji,jj) - ht_wd(ji,jj)) & 
     924                  zcpy(ji,jj) = ABS( (zsshp2_e(ji,jj+1) +    ht_0(ji,jj+1) - zsshp2_e(ji,jj) - ht_0(ji,jj)) & 
    802925                              &    / (zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj  )) ) 
    803926                ELSE 
    … …  
    8861009         ! 
    8871010         ! Add bottom stresses: 
    888          zu_trd(:,:) = zu_trd(:,:) + bfrua(:,:) * un_e(:,:) * hur_e(:,:) 
    889          zv_trd(:,:) = zv_trd(:,:) + bfrva(:,:) * vn_e(:,:) * hvr_e(:,:) 
     1011!jth do implicitly instead 
     1012!         zu_trd(:,:) = zu_trd(:,:) + bfrua(:,:) * un_e(:,:) * hur_e(:,:) 
     1013!         zv_trd(:,:) = zv_trd(:,:) + bfrva(:,:) * vn_e(:,:) * hvr_e(:,:) 
    8901014         ! 
    8911015         ! Add top stresses: 
    … …  
    9331057                            &                                 + zv_frc(ji,jj) ) & 
    9341058                            &   ) * ssvmask(ji,jj) 
     1059  
     1060!jth implicit bottom friction: 
     1061               ua_e(ji,jj) =  ua_e(ji,jj) /(1.0 -   rdtbt * bfrua(ji,jj) * hur_e(ji,jj)) 
     1062               va_e(ji,jj) =  va_e(ji,jj) /(1.0 -   rdtbt * bfrva(ji,jj) * hvr_e(ji,jj)) 
     1063 
    9351064               END DO 
    9361065            END DO 
    … …  
    9411070 
    9421071                  IF( ln_wd ) THEN 
    943                     zhura = MAX(hu_0(ji,jj) + zsshu_a(ji,jj), rn_wdmin1) 
    944                     zhvra = MAX(hv_0(ji,jj) + zsshv_a(ji,jj), rn_wdmin1) 
     1072                    zhura = hu_0(ji,jj) + zsshu_a(ji,jj) 
     1073                    zhvra = hv_0(ji,jj) + zsshv_a(ji,jj) 
    9451074                  ELSE 
    9461075                    zhura = hu_0(ji,jj) + zsshu_a(ji,jj) 
    … …  
    9641093            END DO 
    9651094         ENDIF 
    966          ! 
     1095 
     1096! if ln_rwd: ua_e and va_e should not be masked ; they are used to determine the direction of flow into all cells 
     1097 
     1098!         IF ( ln_rwd) THEN  
     1099!            IF ( ln_wd_diag ) THEN  
     1100!               WRITE(numout, *) 'ua_e: (i,j)         = ', ua_e(iwdg,jwdg) 
     1101!               WRITE(numout, *) 'va_e: (i,j)         = ', va_e(iwdg,jwdg) 
     1102!            END IF  
     1103!            ua_e(:,:) = ua_e(:,:) * zuwdmask(:,:)  
     1104!            va_e(:,:) = va_e(:,:) * zvwdmask(:,:)  
     1105!            IF ( ln_wd_diag ) THEN  
     1106!               WRITE(numout, *) 'ua_e: (i,j)         = ', ua_e(iwdg,jwdg) 
     1107!               WRITE(numout, *) 'va_e: (i,j)         = ', va_e(iwdg,jwdg) 
     1108!            END IF  
     1109!         END IF  
     1110 
     1111          
    9671112         IF( .NOT.ln_linssh ) THEN                     !* Update ocean depth (variable volume case only) 
    9681113            IF( ln_wd ) THEN 
    969               hu_e (:,:) = MAX(hu_0(:,:) + zsshu_a(:,:), rn_wdmin1) 
    970               hv_e (:,:) = MAX(hv_0(:,:) + zsshv_a(:,:), rn_wdmin1) 
     1114              hu_e (:,:) = hu_0(:,:) + zsshu_a(:,:) 
     1115              hv_e (:,:) = hv_0(:,:) + zsshv_a(:,:) 
    9711116            ELSE 
    9721117              hu_e (:,:) = hu_0(:,:) + zsshu_a(:,:) 
    … …  
    10061151            va_b  (:,:) = va_b  (:,:) + za1 * va_e  (:,:)  
    10071152         ELSE                                              ! Sum transports 
    1008             ua_b  (:,:) = ua_b  (:,:) + za1 * ua_e  (:,:) * hu_e (:,:) 
    1009             va_b  (:,:) = va_b  (:,:) + za1 * va_e  (:,:) * hv_e (:,:) 
     1153            IF (.NOT.ln_rwd) THEN   
     1154               ua_b  (:,:) = ua_b  (:,:) + za1 * ua_e  (:,:) * hu_e (:,:) 
     1155               va_b  (:,:) = va_b  (:,:) + za1 * va_e  (:,:) * hv_e (:,:) 
     1156            ELSE  
     1157               ua_b  (:,:) = ua_b  (:,:) + za1 * ua_e  (:,:) * hu_e (:,:) * zuwdmask(:,:) 
     1158               va_b  (:,:) = va_b  (:,:) + za1 * va_e  (:,:) * hv_e (:,:) * zvwdmask(:,:) 
     1159            END IF  
    10101160         ENDIF 
    10111161         !                                   ! Sum sea level 
    10121162         ssha(:,:) = ssha(:,:) + za1 * ssha_e(:,:) 
     1163 
    10131164         !                                                 ! ==================== ! 
    10141165      END DO                                               !        end loop      ! 
    … …  
    10331184         vb2_b(:,:) = zwy(:,:) 
    10341185      ENDIF 
     1186 
     1187      IF ( ln_wd_diag ) THEN  
     1188         WRITE(numout, *) 'ub2_b: (i,j)        = ', ub2_b(iwdg,jwdg) 
     1189         WRITE(numout, *) 'r1_hu_n: (i,j)      = ', r1_hu_n(iwdg,jwdg) 
     1190         WRITE(numout, *) 'zwx: (i,j)          = ', zwx(iwdg,jwdg) 
     1191         WRITE(numout, *) 'un_adv: (i,j)        = ', un_adv(iwdg,jwdg) 
     1192      END IF  
     1193 
    10351194      ! 
    10361195      ! Update barotropic trend: 
    … …  
    10621221         va_b(:,:) =  va_b(:,:) / ( hv_0(:,:) + zsshv_a(:,:) + 1._wp - ssvmask(:,:) ) 
    10631222      ENDIF 
    1064       ! 
     1223 
     1224      IF ( ln_wd_diag ) THEN  
     1225         WRITE(numout, *) 'ua_b: (i,j) A        = ', ua_b(iwdg,jwdg) 
     1226         WRITE(numout, *) 'va_b: (i,j) B        = ', va_b(iwdg,jwdg) 
     1227      END IF  
     1228 
     1229! temporary debugging code  
     1230      IF ( ln_wd_diag ) THEN  
     1231         WRITE(numout, *) 'ua: (i,j,k)  B       = ', ua(iwdg,jwdg,kwdg) 
     1232         WRITE(numout, *) 'ua_b: (i,j)  B       = ', ua_b(iwdg,jwdg) 
     1233         WRITE(numout, *) 'un: (i,j,k)          = ', un(iwdg,jwdg,kwdg) 
     1234         WRITE(numout, *) 'un_b: (i,j)          = ', un_b(iwdg,jwdg) 
     1235         WRITE(numout, *) 'un_adv: (i,j)        = ', un_adv(iwdg,jwdg) 
     1236         WRITE(numout, *) 'va: (i,j,k)          = ', va(iwdg,jwdg,kwdg) 
     1237         WRITE(numout, *) 'va_b: (i,j,k)        = ', va_b(iwdg,jwdg) 
     1238         WRITE(numout, *) 'vn: (i,j,k)          = ', vn(iwdg,jwdg,kwdg) 
     1239         WRITE(numout, *) 'vn_b: (i,j)          = ', vn_b(iwdg,jwdg) 
     1240         WRITE(numout, *) 'vn_adv: (i,j)        = ', vn_adv(iwdg,jwdg) 
     1241      END IF  
     1242 
     1243      ! Correct velocities so that the barotropic velocity equals (un_adv, vn_adv) (in all cases)   
    10651244      DO jk = 1, jpkm1 
    1066          ! Correct velocities: 
    10671245         un(:,:,jk) = ( un(:,:,jk) + un_adv(:,:) - un_b(:,:) ) * umask(:,:,jk) 
    10681246         vn(:,:,jk) = ( vn(:,:,jk) + vn_adv(:,:) - vn_b(:,:) ) * vmask(:,:,jk) 
    1069          ! 
    10701247      END DO 
    1071       ! 
     1248 
     1249      IF ( ln_rwd .and. ln_rwd_bc) THEN  
     1250         DO jk = 1, jpkm1 
     1251            un(:,:,jk) = ( un_adv(:,:) + zuwdav2(:,:)*(un(:,:,jk) - un_adv(:,:)) ) * umask(:,:,jk)  
     1252            vn(:,:,jk) = ( vn_adv(:,:) + zvwdav2(:,:)*(vn(:,:,jk) - vn_adv(:,:)) ) * vmask(:,:,jk)   
     1253         END DO 
     1254      END IF  
     1255 
     1256      IF ( ln_wd_diag ) THEN  
     1257         WRITE(numout, *) 'ua: (i,j,k)          = ', ua(iwdg,jwdg,kwdg) 
     1258         WRITE(numout, *) 'ua_b: (i,j,k)        = ', ua_b(iwdg,jwdg) 
     1259         WRITE(numout, *) 'un: (i,j,k)          = ', un(iwdg,jwdg,kwdg) 
     1260         WRITE(numout, *) 'va: (i,j,k)          = ', va(iwdg,jwdg,kwdg) 
     1261         WRITE(numout, *) 'va_b: (i,j,k)        = ', va_b(iwdg,jwdg) 
     1262         WRITE(numout, *) 'vn: (i,j,k)          = ', vn(iwdg,jwdg,kwdg) 
     1263      END IF  
     1264       
    10721265      CALL iom_put(  "ubar", un_adv(:,:)      )    ! barotropic i-current 
    10731266      CALL iom_put(  "vbar", vn_adv(:,:)      )    ! barotropic i-current 
    … …  
    10981291      CALL wrk_dealloc( jpi,jpj,   zhf ) 
    10991292      IF( ln_wd ) CALL wrk_dealloc( jpi, jpj, zcpx, zcpy ) 
     1293      IF( ln_rwd ) CALL wrk_dealloc( jpi, jpj, ztwdmask, zuwdmask, zvwdmask, zuwdav2, zvwdav2 ) 
    11001294      ! 
    11011295      IF ( ln_diatmb ) THEN 

  • branches/UKMO/ROMS_WAD_7832/NEMOGCM/NEMO/OPA_SRC/DYN/wet_dry.F90

    r7646 r8865  
    11MODULE wet_dry 
     2 
     3   !! includes updates to namelist namwad for diagnostic outputs of ROMS wetting and drying 
     4 
    25   !!============================================================================== 
    36   !!                       ***  MODULE  wet_dry  *** 
    … …  
    3235 
    3336   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) ::   wdmask         !: u- and v- limiter  
    34    REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) ::   ht_wd          !: wetting and drying t-pt depths 
    3537                                                                     !  (can include negative depths) 
     38   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) ::   wdramp, wdrampu, wdrampv ! for hpg limiting 
    3639 
    3740   LOGICAL,  PUBLIC  ::   ln_wd       !: Wetting/drying activation switch (T:on,F:off) 
     41   REAL(wp), PUBLIC  ::   rn_wdmin0   !: depth at which wetting/drying starts 
     42   LOGICAL,  PUBLIC  ::   ln_rwd      !: ROMS Wetting/drying activation switch (T:on,F:off) 
    3843   REAL(wp), PUBLIC  ::   rn_wdmin1   !: minimum water depth on dried cells 
    39    REAL(wp), PUBLIC  ::   rn_wdmin2   !: tolerrance of minimum water depth on dried cells 
     44   REAL(wp), PUBLIC  ::   rn_wdmin2   !: tolerance of minimum water depth on dried cells 
    4045   REAL(wp), PUBLIC  ::   rn_wdld     !: land elevation below which wetting/drying  
    4146                                           !: will be considered 
    4247   INTEGER , PUBLIC  ::   nn_wdit     !: maximum number of iteration for W/D limiter 
     48   LOGICAL,  PUBLIC  ::   ln_rwd_bc   !: ROMS scheme: True implies 3D velocities are set to the barotropic values at points  
     49                                      !: where the flow is from wet points on less than half the barotropic sub-steps   
     50   LOGICAL,  PUBLIC  ::   ln_rwd_rmp  !: use a ramp for the rwd flux limiter between 2 rn_wdmin1 and rn_wdmin1 (rather than a cut-off at rn_wdmin1)       
     51   REAL(wp), PUBLIC  ::   rn_ssh_ref  !: height of z=0 with respect to the geoid;  
     52   REAL(wp), PUBLIC  ::   rn_ht_0     !: the height at which ht_0 = 0    
     53 
     54   LOGICAL,  PUBLIC  ::   ln_wd_diag  ! True implies wad diagnostic at chosen point are printed out 
     55   INTEGER , PUBLIC  ::   jn_wd_i, jn_wd_j, jn_wd_k   ! indices at which diagnostic outputs are generated  
    4356 
    4457   PUBLIC   wad_init                  ! initialisation routine called by step.F90 
    … …  
    5871      !! ** input   : - namwad namelist 
    5972      !!---------------------------------------------------------------------- 
    60       NAMELIST/namwad/ ln_wd, rn_wdmin1, rn_wdmin2, rn_wdld, nn_wdit 
     73      NAMELIST/namwad/ ln_wd, ln_rwd, rn_wdmin0, ln_rwd, rn_wdmin1, rn_wdmin2, rn_wdld, nn_wdit, ln_wd_diag, ln_rwd_bc, & 
     74                     &  rn_ht_0, jn_wd_i, jn_wd_j, jn_wd_k,ln_rwd_rmp 
     75 
    6176      INTEGER  ::   ios                 ! Local integer output status for namelist read 
    6277      INTEGER  ::   ierr                ! Local integer status array allocation  
    … …  
    7893         WRITE(numout,*) '~~~~~~~~' 
    7994         WRITE(numout,*) '   Namelist namwad' 
    80          WRITE(numout,*) '      Logical activation                 ln_wd      = ', ln_wd 
     95         WRITE(numout,*) '      Logical for NOC  wd scheme         ln_wd      = ', ln_wd 
     96         WRITE(numout,*) '      Logical for ROMS wd scheme         ln_rwd     = ', ln_rwd 
     97         WRITE(numout,*) '      Depth at which wet/drying starts rn_wdmin0    = ', rn_wdmin0 
    8198         WRITE(numout,*) '      Minimum wet depth on dried cells rn_wdmin1    = ', rn_wdmin1 
    8299         WRITE(numout,*) '      Tolerance of min wet depth     rn_wdmin2      = ', rn_wdmin2 
    83100         WRITE(numout,*) '      land elevation threshold         rn_wdld      = ', rn_wdld 
    84101         WRITE(numout,*) '      Max iteration for W/D limiter    nn_wdit      = ', nn_wdit 
     102         WRITE(numout,*) '      Logical for WAD diagnostics      ln_wd_diag   = ', ln_wd_diag 
     103         WRITE(numout,*) '      T => baroclinic u,v = 0 at dry pts: ln_rwd_bc = ', ln_rwd_bc      
     104         WRITE(numout,*) '      use a ramp for rwd limiter:      ln_rwd_rmp   = ', ln_rwd_rmp 
     105         WRITE(numout,*) '      the height (z) at which ht_0 = 0:rn_ht_0      = ', rn_ht_0   
     106         WRITE(numout,*) '      i-index for diagnostic point     jn_wd_i      = ', jn_wd_i 
     107         WRITE(numout,*) '      j-index for diagnostic point     jn_wd_j      = ', jn_wd_j 
     108         WRITE(numout,*) '      k-index for diagnostic point     jn_wd_k      = ', jn_wd_k 
    85109      ENDIF 
    86110      ! 
    87       IF(ln_wd) THEN 
    88          ALLOCATE( wdmask(jpi,jpj), ht_wd(jpi,jpj),  STAT=ierr ) 
     111      IF(ln_wd .OR. ln_rwd) THEN 
     112         ALLOCATE( wdmask(jpi,jpj),   STAT=ierr ) 
     113         ALLOCATE( wdramp(jpi,jpj), wdrampu(jpi,jpj), wdrampv(jpi,jpj), STAT=ierr )  
    89114         IF( ierr /= 0 ) CALL ctl_stop('STOP', 'wad_init : Array allocation error') 
    90115      ENDIF 
    … …  
    161186 
    162187             IF( tmask(ji, jj, 1) < 0.5_wp ) CYCLE   ! we don't care about land cells 
    163              IF( ht_wd(ji,jj) > zdepwd )      CYCLE   ! and cells which are unlikely to dry 
     188             IF( ht_0(ji,jj)-rn_ssh_ref > zdepwd )      CYCLE   ! and cells which are unlikely to dry 
    164189 
    165190              zflxp(ji,jj) = max(zflxu(ji,jj), 0._wp) - min(zflxu(ji-1,jj),   0._wp) + & 
    … …  
    168193                           & min(zflxv(ji,jj), 0._wp) - max(zflxv(ji,  jj-1), 0._wp)  
    169194 
    170               zdep2 = ht_wd(ji,jj) + sshb1(ji,jj) - rn_wdmin1 
     195              zdep2 = ht_0(ji,jj) + sshb1(ji,jj) - rn_wdmin1 
    171196              IF(zdep2 .le. 0._wp) THEN  !add more safty, but not necessary 
    172                 sshb1(ji,jj) = rn_wdmin1 - ht_wd(ji,jj) 
     197                sshb1(ji,jj) = rn_wdmin1 - ht_0(ji,jj) 
    173198                IF(zflxu(ji,  jj) > 0._wp) zwdlmtu(ji  ,jj) = 0._wp 
    174199                IF(zflxu(ji-1,jj) < 0._wp) zwdlmtu(ji-1,jj) = 0._wp 
    … …  
    180205        END DO 
    181206 
     207 
     208! slwa 
     209! HPG limiter from jholt 
     210      wdramp(:,:) = min((ht_0(:,:) + sshb1(:,:) - rn_wdmin1)/(rn_wdmin0 - rn_wdmin1),1.0_wp) 
     211!      write(6,*)'wdramp ',wdramp(10,10),wdramp(10,11) 
     212!jth assume don't need a lbc_lnk here 
     213       DO jj = 1, jpjm1 
     214         DO ji = 1, jpim1 
     215           wdrampu(ji,jj) = min(wdramp(ji,jj),wdramp(ji+1,jj)) 
     216           wdrampv(ji,jj) = min(wdramp(ji,jj),wdramp(ji,jj+1)) 
     217         ENDDO 
     218       ENDDO 
     219        !wdrampu(:,:)=1.0_wp 
     220        !wdrampv(:,:)=1.0_wp 
     221! end HPG limiter 
     222 
     223 
    182224       
    183225        !! start limiter iterations  
    … …  
    193235         
    194236                 IF( tmask(ji, jj, 1) < 0.5_wp ) CYCLE  
    195                  IF( ht_wd(ji,jj) > zdepwd )      CYCLE 
     237                 IF( ht_0(ji,jj) > zdepwd )      CYCLE 
    196238         
    197239                 ztmp = e1e2t(ji,jj) 
    … …  
    203245           
    204246                 zdep1 = (zzflxp + zzflxn) * z2dt / ztmp 
    205                  zdep2 = ht_wd(ji,jj) + sshb1(ji,jj) - rn_wdmin1 - z2dt * sshemp(ji,jj) 
     247                 zdep2 = ht_0(ji,jj) + sshb1(ji,jj) - rn_wdmin1 - z2dt * sshemp(ji,jj) 
    206248           
    207249                 IF( zdep1 > zdep2 ) THEN 
    … …  
    317359 
    318360             IF( tmask(ji, jj, 1 ) < 0.5_wp) CYCLE   ! we don't care about land cells 
    319              IF( ht_wd(ji,jj) > zdepwd )      CYCLE   ! and cells which are unlikely to dry 
     361             IF( ht_0(ji,jj) > zdepwd )      CYCLE   ! and cells which are unlikely to dry 
    320362 
    321363              zflxp(ji,jj) = max(zflxu(ji,jj), 0._wp) - min(zflxu(ji-1,jj),   0._wp) + & 
    … …  
    324366                           & min(zflxv(ji,jj), 0._wp) - max(zflxv(ji,  jj-1), 0._wp)  
    325367 
    326               zdep2 = ht_wd(ji,jj) + sshn_e(ji,jj) - rn_wdmin1 
     368              zdep2 = ht_0(ji,jj) + sshn_e(ji,jj) - rn_wdmin1 
    327369              IF(zdep2 .le. 0._wp) THEN  !add more safety, but not necessary 
    328                 sshn_e(ji,jj) = rn_wdmin1 - ht_wd(ji,jj) 
     370                sshn_e(ji,jj) = rn_wdmin1 - ht_0(ji,jj) 
    329371                IF(zflxu(ji,  jj) > 0._wp) zwdlmtu(ji  ,jj) = 0._wp 
    330372                IF(zflxu(ji-1,jj) < 0._wp) zwdlmtu(ji-1,jj) = 0._wp 
    … …  
    348390         
    349391                 IF( tmask(ji, jj, 1 ) < 0.5_wp) CYCLE  
    350                  IF( ht_wd(ji,jj) > zdepwd )      CYCLE 
     392                 IF( ht_0(ji,jj) > zdepwd )      CYCLE 
    351393         
    352394                 ztmp = e1e2t(ji,jj) 
    … …  
    358400           
    359401                 zdep1 = (zzflxp + zzflxn) * z2dt / ztmp 
    360                  zdep2 = ht_wd(ji,jj) + sshn_e(ji,jj) - rn_wdmin1 - z2dt * zssh_frc(ji,jj) 
     402                 zdep2 = ht_0(ji,jj) + sshn_e(ji,jj) - rn_wdmin1 - z2dt * zssh_frc(ji,jj) 
    361403           
    362404                 IF(zdep1 > zdep2) THEN 

  • branches/UKMO/ROMS_WAD_7832/NEMOGCM/NEMO/OPA_SRC/SBC/sbcrnf.F90

    r8841 r8865  
    3939   !                                                !!* namsbc_rnf namelist * 
    4040   CHARACTER(len=100)         ::   cn_dir            !: Root directory for location of rnf files 
    41    LOGICAL                   ::   ln_rnf_depth      !: depth       river runoffs attribute specified in a file 
     41   LOGICAL           , PUBLIC ::   ln_rnf_depth      !: depth       river runoffs attribute specified in a file 
    4242   LOGICAL                    ::      ln_rnf_depth_ini  !: depth       river runoffs  computed at the initialisation 
    4343   REAL(wp)                   ::      rn_rnf_max        !: maximum value of the runoff climatologie (ln_rnf_depth_ini =T) 
    … …  
    133133            END WHERE 
    134134         ELSE                                                        ! use SST as runoffs temperature 
    135             rnf_tsc(:,:,jp_tem) = sst_m(:,:) * rnf(:,:) * r1_rau0 
     135         rnf_tsc(:,:,jp_tem) = MAX(sst_m(:,:),0.0_wp) * rnf(:,:) * r1_rau0 !CEOD River is fresh water so must at least be 0 unless we consider ice 
    136136         ENDIF 
    137137         !                                                           ! use runoffs salinity data 

  • branches/UKMO/ROMS_WAD_7832/NEMOGCM/NEMO/OPA_SRC/TRA/tranxt.F90

    r8841 r8865  
    317317                  IF( jk == mikt(ji,jj) ) THEN           ! first level  
    318318                     ze3t_f = ze3t_f - zfact2 * ( (emp_b(ji,jj)    - emp(ji,jj)   )  & 
    319                             &                   - (rnf_b(ji,jj)    - rnf(ji,jj)   )  & 
    320319                            &                   + (fwfisf_b(ji,jj) - fwfisf(ji,jj))  ) 
    321320                     ztc_f  = ztc_f  - zfact1 * ( psbc_tc(ji,jj,jn) - psbc_tc_b(ji,jj,jn) ) 
    322321                  ENDIF 
     322                  ! Rivers can be not just at the surface must go down to nk_rnd(ji,jj) 
     323                  IF( ln_rnf_depth ) THEN 
     324                     IF( mikt(ji,jj) <=jk .and. jk <= nk_rnf(ji,jj)  ) THEN 
     325                      ze3t_f = ze3t_f - zfact2 * ( - (rnf_b(ji,jj)    - rnf(ji,jj)   )  )*(e3t_n(ji,jj,jk)/h_rnf(ji,jj) ) ! as we have sigma can do that here change later 
     326                     ENDIF 
     327                  ELSE 
     328                      IF( jk == mikt(ji,jj) ) THEN           ! first level  
     329                         ze3t_f = ze3t_f - zfact2 * ( - (rnf_b(ji,jj)    - rnf(ji,jj)   ) )  
     330                      ENDIF 
     331                  ENDIF 
     332 
    323333                  ! 
    324334                  ! solar penetration (temperature only) 

  • branches/UKMO/ROMS_WAD_7832/NEMOGCM/NEMO/OPA_SRC/TRA/trasbc.F90

    r7788 r8865  
    3434   USE wrk_nemo       ! Memory Allocation 
    3535   USE timing         ! Timing 
     36   USE wet_dry 
    3637 
    3738   IMPLICIT NONE 
    … …  
    122123      DO jj = 2, jpj 
    123124         DO ji = fs_2, fs_jpim1   ! vector opt. 
    124             sbc_tsc(ji,jj,jp_tem) = r1_rau0_rcp * qns(ji,jj)   ! non solar heat flux 
     125            IF ( ln_rwd_rmp ) THEN    ! If near WAD point limite the flux for now 
     126               IF ( sshn(ji,jj) + ht_0(ji,jj) >  2._wp * rn_wdmin1 ) THEN 
     127                  sbc_tsc(ji,jj,jp_tem) = r1_rau0_rcp * qns(ji,jj)   ! non solar heat flux 
     128               ELSE IF ( sshn(ji,jj) + ht_0(ji,jj) >  rn_wdmin1 ) THEN 
     129                  sbc_tsc(ji,jj,jp_tem) = r1_rau0_rcp * qns(ji,jj) * (tanh(5._wp*( ( sshn(ji,jj) + ht_0(ji,jj) -  rn_wdmin1 )/rn_wdmin1)) ) 
     130               ELSE 
     131                  sbc_tsc(ji,jj,jp_tem) = 0._wp 
     132               ENDIF 
     133            ELSE  
     134               sbc_tsc(ji,jj,jp_tem) = r1_rau0_rcp * qns(ji,jj)   ! non solar heat flux 
     135            ENDIF 
     136 
    125137            sbc_tsc(ji,jj,jp_sal) = r1_rau0     * sfx(ji,jj)   ! salt flux due to freezing/melting 
    126138         END DO 

  • branches/UKMO/ROMS_WAD_7832/NEMOGCM/NEMO/OPA_SRC/stpctl.F90

    r8841 r8865  
    2222   USE lib_mpp         ! distributed memory computing 
    2323   USE lib_fortran     ! Fortran routines library  
     24   USE wet_dry,  ONLY: ln_wd, ln_rwd, rn_ssh_ref    ! reference depth for negative bathy 
    2425 
    2526   IMPLICIT NONE 
    … …  
    150151      DO jj = 1, jpj 
    151152         DO ji = 1, jpi 
    152             IF( tmask(ji,jj,1) == 1) zsshmax = MAX( zsshmax, ABS(sshn(ji,jj)) ) 
     153            IF( (ln_wd .OR. ln_rwd)) THEN 
     154               IF( tmask(ji,jj,1) == 1) zsshmax = MAX( zsshmax, ABS(sshn(ji,jj)+rn_ssh_ref) ) 
     155            ELSE 
     156               IF( tmask(ji,jj,1) == 1) zsshmax = MAX( zsshmax, ABS(sshn(ji,jj)) ) 
     157            ENDIF 
    153158         END DO 
    154159      END DO 
Note: See TracChangeset for help on using the changeset viewer.