Opened 6 years ago

Closed 6 years ago

Last modified 6 years ago

#648 closed upgrade (Done)

IMPORTANT: Wetting and drying revisited

Reported by: arango Owned by: arango
Priority: major Milestone: Release ROMS/TOMS 3.7
Component: Nonlinear Version: 3.7
Keywords: Cc:

Description (last modified by 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:

https://www.myroms.org/trac/CookInlet.png

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.

Change History (2)

comment:1 Changed 6 years ago by arango

  • Description modified (diff)
  • Resolution set to Done
  • Status changed from new to closed

comment:2 Changed 6 years ago by arango

  • Description modified (diff)
Note: See TracTickets for help on using tickets.