metno / emep-ctm

Open Source EMEP/MSC-W model
GNU General Public License v3.0
30 stars 20 forks source link

PODs calculation with ad-hoc SGS and EGS #114

Open mvieno opened 1 year ago

mvieno commented 1 year ago

Hi, I am using the latest rv4.45 to calculate some PODs over the UK. I have an issue when I want to use a later SGS (let say for WHEAT) SGS=141 and earlier end of growing season EGS= 180. These options do trigger the 20 days check (the comment in the cods says 30 days) for the minimum growing season at specific altitude and latitude combination.

EMEP error and stop: LUtopo:iTopo problem too short GS: 33.9057426452637 187 183 67.8665699207478 46.4848414259221

My workaround (not fully tested) is as follow:

        ! Just a crude check that the growing season remains > 30 days:
        ! 30 days would work for W. Europe, but v. big montains in Asian
        ! part complicate things. Shouldn't use latitude function here
        if ( LandCover(i,j)%EGS(ilu)-LandCover(i,j)%SGS(ilu) < 20 ) then
            print *, dtxt//"iTopo problem too short GS: ", model_surf_elevation(i,j), &
            LandCover(i,j)%SGS(ilu), LandCover(i,j)%EGS(ilu), &
             glat(i,j), glon(i,j)

            !!! MV 13092023 call StopAll( dtxt//'iTopo SGS EGS problem')
            !!! workaround for GRID squares where GS may be less than 20 days
            !!! add + and minus 10 days centered to the SGS 

            LandCover(i,j)%EGS(ilu) = LandCover(i,j)%SGS(ilu)

            LandCover(i,j)%SGS(ilu) = LandCover(i,j)%SGS(ilu) -10
            LandCover(i,j)%EGS(ilu) = LandCover(i,j)%EGS(ilu) +10

            print *, dtxt//"Assumed GS of 20 days centered on the new SGS ", model_surf_elevation(i,j), &
            LandCover(i,j)%SGS(ilu), LandCover(i,j)%EGS(ilu), &
             glat(i,j), glon(i,j)   

        end if

Workaround: new SGS and EGS : LUtopo: Assumed GS of 20 days centred on the new SGS 33.9057426452637 177 197 67.8665699207478 46.4848414259221

Massimo

mifads commented 1 year ago

Hi Massimo, In the first printout it is saying that SGS = 187 and EGS = 183, ie before SGS. That just suggests a problem in the input tables. Are you sure you set SGS, DSGS etc properly?

Also, the printout says lat = 67 degN. That is far north for a UK site and wheat growing, and is I suspect the source of your problem ;-)

/Dave

mvieno commented 1 year ago

Hi Dave, thanks for your reply.

The issue in EMEPrv4.45 is that if the SGS50 and EGS50 are not distant enough the calculated SGS and EGS that are adjusted for the latitude and elevation are less than 20 days apart. If not a long list of locations (across Europe) violate the 20 days.

SGS=187 and EGS=183 have been after the calculation.

As an example, I used the follow setup for POTATO: PFT hveg Alb eNH4 SGS50 DSGS EGS50 DEGS LAImin LAImax ... POTATO POTATO ECR NOLPJ 1 20 1 146 0 216 0 0 4.2 35 65 ...

With this (or similar) in some places the calculated SGS and EGS stops the model.

mifads commented 1 year ago

With those potato settings the GS should always be 71 days (216-146+1), since the DSGS and DEGS are zero. For wheat the GS should also stay constant, since we be default set DSGS = DEGS = 2.57. Thus SGS gets later as one progresses away from 50 deg N, but so does EGS. (This is different for e.g. forests, where the GS gets shorter with latitude.)

mvieno commented 1 year ago

.....mmmmm...

Do the Input_LandDefs and Input_DO3SE file need to have the same structure? maybe that's my problem...the location of POTATO is different. I'll do more tests before using more of your time :)

mifads commented 1 year ago

The land-cover need to be in the same order.

mifads commented 1 year ago

Hi Massimo, is this sorted out now? If so, we can close the issue.

mvieno commented 1 year ago

Hi Massimo, is this sorted out now? If so, we can close the issue.

Hi Dave, in short yes. I managed to output the PODs in rv4.45 if I append my PODs to the one already included in the config. This was superseded by the release of rv5.0 where I managed to output all the PODs following your suggestion to re-name my additional land types with the prefix "IAM_".

mvieno commented 1 year ago

I guess the "IAM_" prefix is checked here: Line 639 of Landuse_mod.f90 (rv5.0)

... else if ( index(OutputVegO3(i)%txtLC, 'IAM_') >0 ) then j = find_index(OutputVegO3(i)%txtLC, EXTRA_FLUX_VEGS) end if ...

this is how my additional entry in the Input_LandDefs.csv looks like: Name,code,type,PFT,hveg,Alb,eNH4,SGS50,DSGS,EGS50,DEGS,LAImin,.... IAM_WHEAT,IAM_WHEAT,ECR,NOLPJ,1,20,1,1,2.57,197,0,0,....

I suppose only the name need to have the prefix "IAM_"?