Evaporation + precipitation mass conservation

[ashbre]

Hi all,

I'm currently running a lake model where E and P are important in terms of removing and adding water to the domain. However I found that ROMS does not support this function.

I wondered if anybody added this feature to their model code and if so, could they share their experience.

As far as I can tell, I can follow 2 recipes. One in the routine "step2d_LF_AM3.h"

Code: Select all

     IF (LuvSrc(ng)) THEN
        DO is=1,Nsrc(ng)
          i=SOURCES(ng)%Isrc(is)
          j=SOURCES(ng)%Jsrc(is)
          IF (((IstrR.le.i).and.(i.le.IendR)).and.                      &
     &        ((JstrR.le.j).and.(j.le.JendR))) THEN
            IF (INT(SOURCES(ng)%Dsrc(is)).eq.0) THEN
              cff=1.0_r8/(on_u(i,j)*                                    &
     &                    0.5_r8*(zeta(i-1,j,knew)+h(i-1,j)+            &
     &                            zeta(i  ,j,knew)+h(i  ,j)))
              ubar(i,j,knew)=SOURCES(ng)%Qbar(is)*cff
            ELSE
              cff=1.0_r8/(om_v(i,j)*                                    &
     &                    0.5_r8*(zeta(i,j-1,knew)+h(i,j-1)+            &
     &                            zeta(i,j  ,knew)+h(i,j  )))
              vbar(i,j,knew)=SOURCES(ng)%Qbar(is)*cff

            END IF
          END IF
        END DO
      END IF

Replacing SOURCES(ng)%Qbar(is) with my equivalent evaporation/precipitation calculation - though I imagine this will introduce an unwanted horizontal velocity.

Or I can follow the recipe of omega.F

! Apply mass point sources (volume vertical influx), if any.
!

Code: Select all

       IF (LwSrc(ng)) THEN
          DO is=1,Nsrc(ng)
            ii=SOURCES(ng)%Isrc(is)
            jj=SOURCES(ng)%Jsrc(is)
            IF (((IstrR.le.ii).and.(ii.le.IendR)).and.                  &
     &          ((JstrR.le.jj).and.(jj.le.JendR)).and.                  &
     &          (j.eq.jj)) THEN



              DO k=1,N(ng)
                W(ii,jj,k)=W(ii,jj,k)+SOURCES(ng)%Qsrc(is,k)
              END DO
            END IF
          END DO
        END IF

Again with an appropriate Evap/precip value replacing the SOURCES(ng)%Qsrc(is,k) value.

However, I might be thinking about this problem in the wrong way, so if anybody has any ideas, it would be greatly appreciated!

Many thanks,
Ash

[jcwarner]

does this help
viewtopic.php?f=14&t=5556

[arango]

I think that there is a confusion here with the forcing terminology. The surface and bottom heat (Q) and freshwater fluxes (E-P) are the boundary conditions to the vertical tracer (temperature and salinity, respectively) diffusion terms in the governing equations. We need to provide boundary conditions to both surface and bottom fluxes. Although, it will be tuff to have freshwater precipitation at the bottom! Such fluxes need to be specified at every horizontal grid point! By default, the bottom tracer flux is set to zero in ROMS with the analytical routine ana_btflux.h. Also, If we have sea ice at the top, we will need a salt flux, but that will be taken care of in the sea-ice coupling interface. By the way, the 3D momentum equations also have a surface (wind stress) and bottom boundary conditions to the vertical viscosity terms.

I can foresee applications that may have a bottom flux to the diffusion tracer equation (heat, freshwater pipe, nutrient, etc.) that are different than zero. The user can specify such fluxes with analytical functions or reading from a NetCDF file. ROMS provides comprehensive metadata for such forcing.

Now, the code sections shown above are for additional point sources in the transport or concentration of tracers. It is a different term in the governing equations. We usually use point sources to model river runoff forcing. Such point sources can occur anywhere in the water column. It is up to the user.

[ashbre]

Hi there,

Thank you very much for your quick responses!

Apologies I think I didn't explain my problem very clearly.

My problem is that when I add evaporation and precipitation, the amount of water that I have in the system doesn't change. Actually in my case I'm not solving for salinity, as everything is fresh in my lake. So I wondered, if it rains, what's the best way of adding the water that should be coming into the system.

Best wishes,
Ash

[wilkin]

I suggest you look at the code for handling river point sources in the LwSrc option. This adds the volume flux as a cell centered divergence. You obviously don't want to create Lm times Mm sources, but a few modifications to omega.F and step3d_t.F might be sufficient.

If you look in omega.F at the block inside ...

Code: Select all

 IF (LwSrc(ng)) THEN

you'll see how for cells where there is a source the volume flux divergence is augmented with the river discharge split across the cells in the vertical. There is a recursion in the vertical so the code is written to accumulate flux in the vertical.

You just want to impose a volume flux in/out of the surface cell only, so in the code block above LwSrc which reads ...

Code: Select all

        DO k=1,N(ng)
          DO i=Istr,Iend
            W(i,j,k)=W(i,j,k-1)-                                        &
     &               (Huon(i+1,j,k)-Huon(i,j,k)+                        &
     &                Hvom(i,j+1,k)-Hvom(i,j,k))
          END DO
        END DO

I think you can simply modify the block to something like this ...

Code: Select all

       DO k=1,N(ng)-1
          DO i=Istr,Iend
            W(i,j,k)=W(i,j,k-1)-                                        &
     &               (Huon(i+1,j,k)-Huon(i,j,k)+                        &
     &                Hvom(i,j+1,k)-Hvom(i,j,k))
          END DO
          W(i,j,N(ng))=W(i,j,N(ng)-1)-                                  &
     &             (Huon(i+1,j,k)-Huon(i,j,k)+                          &
     &              Hvom(i,j+1,k)-Hvom(i,j,k))+                         &
     &             [EminusP as cubic meters per second]
        END DO

Check my logic here on the units of EminusP. And you'll have to figure out how to make EminusP available to omega.F.

It's important you also locate the LwSrc code in step2d_LF_AM3.h and make a complimentary adjustment there or else the 2-D zeta equation will be inconsistent with omega.

Even though you have all fresh water, you should still check the LwSrc code in step3d_t.F because this is there is make sure that the tracers (temp and salt) do the right thing. Since the source is changing the volume of water, keeping the concentration constant would imply a change in mass of tracer. This might be irrelevant to salinity in your case, but you don't want to create an implicit heat flux because you are "diluting" temperature. You'll have to think it through and do some careful tests on conservation in idealized cases.

[ashbre]

Thank you for your help on this.

I made the changes as you recommended and it works perfectly!

Thanks again,
Ash