IMPORTANT: Wetting and drying revisited

[arango]

There have been several issues with the wetting and drying (WET_DRY) algorithms and regular and perfect restart. In some applications, the model blows-up after the restart. Hopefully, this update fixes all these problems. A lot of files were modified.

What it is new:

• The wetting and drying algorithm in step2d_LF_AM3.h was updated and it is the same as that in COAWST. Many thanks to John Warner for making the algorithm more robust.

• The file wetdry.F was split into four subroutines: wetdry_tile, wetdry_ini_tile, wetdry_mask_tile, and wetdry_avg_mask_tile. This allows to have access to different parts of the wetting and drying algorithm from various places in ROMS without repeating code.

• A new routine get_wetdry.F is introduced to read wet/dry mask arrays (pmask_wet, rmask_wet, umask_wet, and vmask_wet) from restart NetCDF file. In initial.F, we now have:

  #ifdef WET_DRY
  !
  !-----------------------------------------------------------------------
  !  Process initial wet/dry masks.
  !-----------------------------------------------------------------------
  !
        DO ng=1,Ngrids
  !
  !  If restart, read in wet/dry masks.
  !
          IF (nrrec(ng).ne.0) THEN
  !$OMP MASTER
            CALL get_wetdry (ng, iNLM, IniRec(ng))
  !$OMP END MASTER
  # ifdef DISTRIBUTE
            CALL mp_bcasti (ng, iNLM, exit_flag)
  # endif
  !$OMP BARRIER
            IF (exit_flag.ne.NoError) RETURN
          ELSE
            DO tile=first_tile(ng),last_tile(ng),+1
              CALL wetdry (ng, tile, Tindex(ng), .TRUE.)
            END DO
  !$OMP BARRIER
          END IF
        END DO
  #endif
  

• The wetting and drying mask have the following range values:

     pmask_wet, wetdry_mask_psi:    0=land  1=water  2=no-slip boundary
     rmask_wet, wetdry_mask_rho:    0=land  1=water
     umask_wet, wetdry_mask_u:      0=land  1=water
     vmask-wet, wetdry_mask_v:      0=land  1=water
  

• In wetting and drying, we need to activate both MASKING and WET_DRY since there is permanent land/sea mask and time dependent wetting and drying mask. Both directives need to be separate because the wetting and drying are adjointable. For example, we need to have:

        DO j=Jstr,Jend
          DO i=Istr,Iend
            ...
  #ifdef MASKING
            Rarray(i,j)=Rarray(i,j)*rmask(i,j)
  #endif
  #ifdef WET_DRY
            Rarray(i,j)=Rarray(i,j)*rmask_wet(i,j)
  #endif
            ...
          END DO
        END DO
  

• The internal arrays pmask_io, rmask_io, umask_io, and vmask_io are renamed to a more appropriate names pmask_full, rmask_full, umask_full, and vmask_full, respectively. They are computed in wetdry.F as:

          DO j=JstrR,JendR
            DO i=IstrR,IendR
              rmask_full(i,j)=rmask_wet(i,j)*rmask(i,j)
            END DO
          END DO
          DO j=Jstr,JendR
            DO i=Istr,IendR
              pmask_full(i,j)=pmask_wet(i,j)*pmask(i,j)
            END DO
          END DO
          DO j=JstrR,JendR
            DO i=Istr,IendR
              umask_full(i,j)=umask_wet(i,j)*umask(i,j)
            END DO
          END DO
          DO j=Jstr,JendR
            DO i=IstrR,IendR
              vmask_full(i,j)=vmask_wet(i,j)*vmask(i,j)
            END DO
          END DO
  

• The wetting and drying masks are now written to the restart file as: wetdry_mask_psi, wetdry_mask_rho, wetdry_mask_u, and wetdry_mask_v.

• Several calls to NetCDF writing routines in wrt_rst.F use the optional argument SetFillVal = .FALSE. to avoid setting large special value (1E+37) on land points. This is very important in perfect restart. Many thanks to Kate Hedstrom for her help debugging this. The tracer values cannot be masked when wetting and drying is activated in the numerical kernel or in the manipulation of output fields.

• If writing only water points (WRITE_WATER option), we need to use only the permanent land/sea masking arrays as an argument to the NetCDF output calls. The full mask arrays (pmask_full, rmask_full, umask_full, and vmask_full) cannot be used to compact the data into one spatial vector because the time dependency of the wetting and drying masks. Therefore, the output routines were changed accordingly.

---

In addition, I corrected few typos and small bugs that have been reported in the ROMS forum. Thank you to everybody for bringing these issues to my attention.

---

Many thanks to Lyon Lanerolle for reporting these problems with perfect restart and for providing his Cook Inlet application for testing:

The grid is huge (1044x724x30) and it was blowing-up after restarting the model. It was used to test the changes to the wetting and drying algorithms. I having running this huge test case in my 8-CPUs desktop in case that I need to examine problem in the TotalView debugger. I have been able to restart the model several times and it is running successfully:

   STEP   Day HH:MM:SS  KINETIC_ENRG   POTEN_ENRG    TOTAL_ENRG    NET_VOLUME
          C => (i,j,k)       Cu            Cv            Cw         Max Speed

      0     0 00:00:00  0.000000E+00  7.600860E+02  7.600860E+02  6.396795E+12
         (000,0000,00)  0.000000E+00  0.000000E+00  0.000000E+00  0.000000E+00
      DEF_HIS   - creating history file, Grid 01: cook_inlet_his_0001.nc
      WRT_HIS   - wrote history  fields (Index=1,1) into time record = 0000001
      DEF_STATION - creating stations file, Grid 01: cook_inlet_sta.nc
      1     0 00:00:05  8.664375E-16  7.600860E+02  7.600860E+02  6.396795E+12
         (386,0113,01)  1.100680E-09  7.030073E-09  0.000000E+00  2.549995E-06
      2     0 00:00:10  2.402784E-10  7.600860E+02  7.600860E+02  6.396795E+12
         (181,0640,01)  5.215735E-04  5.364590E-05  2.222223E-11  2.209626E-02
      3     0 00:00:15  8.536992E-12  7.600860E+02  7.600860E+02  6.396795E+12
         (181,0640,30)  0.000000E+00  0.000000E+00  2.634209E+01  1.285078E-02
....

 NLM: GET_STATE - Read state initial conditions,             t =     0 12:00:00
                   (Grid 01, File: cook_inlet_rst.nc, Rec=0002, Index=1)

   STEP   Day HH:MM:SS  KINETIC_ENRG   POTEN_ENRG    TOTAL_ENRG    NET_VOLUME
          C => (i,j,k)       Cu            Cv            Cw         Max Speed

   8640     0 12:00:00  6.990734E-05  7.599744E+02  7.599745E+02  6.392954E+12
         (177,0823,30)  1.022442E-02  6.039011E-03  0.000000E+00  2.274181E-01
      DEF_STATION - inquiring stations file, Grid 01: cook_inlet_sta.nc
   8641     0 12:00:05  6.995225E-05  7.599746E+02  7.599746E+02  6.392955E+12
         (178,0637,30)  0.000000E+00  0.000000E+00  9.772995E+00  2.270133E-01
   8642     0 12:00:10  6.999734E-05  7.599747E+02  7.599748E+02  6.392957E+12
         (228,0667,30)  2.714938E-06  0.000000E+00  1.091613E+01  2.266488E-01
   8643     0 12:00:15  7.004253E-05  7.599748E+02  7.599749E+02  6.392958E+12
...

 NLM: GET_STATE - Read state initial conditions,             t =     1 00:00:00
                   (Grid 01, File: cook_inlet_rst.nc, Rec=0002, Index=1)

   STEP   Day HH:MM:SS  KINETIC_ENRG   POTEN_ENRG    TOTAL_ENRG    NET_VOLUME
          C => (i,j,k)       Cu            Cv            Cw         Max Speed

  17280     1 00:00:00  1.628109E-03  7.595062E+02  7.595078E+02  6.385508E+12
         (172,0592,30)  3.321712E-02  1.269262E-02  0.000000E+00  7.566996E-01
      DEF_STATION - inquiring stations file, Grid 01: cook_inlet_sta.nc
  17281     1 00:00:05  1.628419E-03  7.595064E+02  7.595080E+02  6.385506E+12
         (177,0636,30)  0.000000E+00  0.000000E+00  3.088772E+01  7.568594E-01
  17282     1 00:00:10  1.628729E-03  7.595065E+02  7.595081E+02  6.385503E+12
         (170,0619,30)  0.000000E+00  0.000000E+00  1.418715E+01  7.570190E-01
...

 NLM: GET_STATE - Read state initial conditions,             t =     1 12:00:00
                   (Grid 01, File: cook_inlet_rst.nc, Rec=0002, Index=1)

   STEP   Day HH:MM:SS  KINETIC_ENRG   POTEN_ENRG    TOTAL_ENRG    NET_VOLUME
          C => (i,j,k)       Cu            Cv            Cw         Max Speed

  25920     1 12:00:00  7.244958E-04  7.594712E+02  7.594719E+02  6.389436E+12
         (616,0511,30)  2.413032E-02  1.226487E-01  0.000000E+00  2.061889E+00
      DEF_STATION - inquiring stations file, Grid 01: cook_inlet_sta.nc
  25921     1 12:00:05  7.252824E-04  7.594714E+02  7.594722E+02  6.389437E+12
         (179,0638,30)  1.161984E-03  3.170779E-06  8.484577E+00  2.060438E+00
  25922     1 12:00:10  7.260429E-04  7.594717E+02  7.594724E+02  6.389439E+12
         (207,0968,30)  6.414001E-03  0.000000E+00  9.312273E+00  2.059241E+00
...

 NLM: GET_STATE - Read state initial conditions,             t =     2 00:00:00
                   (Grid 01, File: cook_inlet_rst.nc, Rec=0002, Index=1)

   STEP   Day HH:MM:SS  KINETIC_ENRG   POTEN_ENRG    TOTAL_ENRG    NET_VOLUME
          C => (i,j,k)       Cu            Cv            Cw         Max Speed

  34560     2 00:00:00  5.806858E-03  7.586865E+02  7.586923E+02  6.375351E+12
         (617,0507,30)  6.565633E-03  1.592644E-01  0.000000E+00  2.073913E+00
      DEF_STATION - inquiring stations file, Grid 01: cook_inlet_sta.nc
  34561     2 00:00:05  5.808851E-03  7.586864E+02  7.586922E+02  6.375339E+12
         (212,0974,30)  1.397617E-02  4.390969E-02  1.038334E+01  2.073022E+00
  34562     2 00:00:10  5.810832E-03  7.586862E+02  7.586920E+02  6.375327E+12
         (211,0975,30)  2.637443E-02  4.564990E-02  2.213804E+01  2.072133E+00
...
and so on

As you can see, the model is going stable restarting with wetting and drying and perfect restart.