Skip to content
Snippets Groups Projects
Commit f11c8749 authored by Dethinne Thomas's avatar Dethinne Thomas
Browse files

working on MODIS per sector

parent 776a8318
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
MAR/code_nestor/src/MODlai.f 0 → 100644
+ 241
− 0
View file @ f11c8749
C +-------------------------------------------------------------------+
C | Subroutine MODlai Apr 2023 NESTING |
C +-------------------------------------------------------------------+
SUBROUTINE MODlai
IMPLICIT none
C +---General variables
C + -----------------
INCLUDE 'NSTdim.inc'
INCLUDE 'NSTvar.inc'
INCLUDE 'LOCfil.inc'
INCLUDE 'NetCDF.inc'
INCLUDE 'NESTOR.inc'
real ,parameter :: reso_y=0.0042 !ny/cy
real ,parameter :: reso_x=0.0042 !nx/cx
real ,parameter :: MOD_lon=-0.4996
real ,parameter :: MOD_lat=46.9046
integer,parameter :: cx = 2750 ! I of MODIS
integer,parameter :: cy = 1700 ! J of MODIS
integer ii,jj,i,j,k,l, i1,i2,j1,j2,notin,n
integer NET_ID_LAI,NETcid,Rcode,start(4),count(4)
integer nbr_day,i_cent,j_cent,G_nx,G_ny, nerror
real AUXlon,AUXlat,debug
real AUXlo1,AUXlo2,AUXla1,AUXla2
real MODIS_lai(cx,cy),nsamp(3), laisum(3),MODIS_lu(cx,cy)
integer DATiyr,DATmma,DATjda,DATjhu
CALL DATcnv (DATiyr,DATmma,DATjda,DATjhu,DATtim,.false.)
call CF_read2D("input/VEGE/MODIS_lu.nc","dominant_value",
. 1,cx,cy,1,MODIS_lu)
nbr_day=0
do i=1,DATmma-1
if(i==1.or.i==3.or.i==5.or.i==7.or.i==8.or.i==10.or.i==12)
. nbr_day=nbr_day+31
if(i==4.or.i==6.or.i==9.or.i==11) nbr_day=nbr_day+30
if(i==2) nbr_day=nbr_day+28
enddo
nbr_day=nbr_day+DATjda
!-----------------------------------------------------------------------
NETcid = NCOPN("input/VEGE/Climato_non_leap_year.nc"
. ,NCNOWRIT,Rcode)
NET_ID_LAI = NCVID(NETcid,'LAI',Rcode)
start(1)=1
start(2)=1
start(3)=nbr_day ! time step
count(1)=cx
count(2)=cy
count(3)=1
Rcode = nf_get_vara_real(NETcid,NET_ID_LAI,start,count,MODIS_lai)
CALL NCCLOS(NETcid, RCODE)
!-----------------------------------------------------------------------
write(6,*) 'MERRA2 LAI-GLF data set'
write(6,*) '~~~~~~~~~~~~~~~~~~~~~~~'
write(6,*) ' '
nerror = 0
DO j=1,my !Loop on each Y of MAR
DO i=1,mx !Loop on each x of MAR
C + *****
IF(NSTsol(i,j)>=3.and.NST__y(i,j)>-60.)THEN !if not on ice
C + *****
if(NSTsvt(i,j,1)<=0.or.NSTsvt(i,j,1)==13) then
DO l=1,nvx
NSTlai(i,j,l) = 0.0
NSTglf(i,j,l) = 0.0
enddo
else
AUXlon = NST__x(i,j)
AUXlat = NST__y(i,j)
C +---Search for the closest point in data file
C + -----------------------------------------
i_cent=NINT((AUXlon-REAL(MOD_lon))/reso_x)+1
j_cent=NINT((AUXlat-REAL(MOD_lat))/reso_y)+1
C +---Compute the resolution of the considered NST cell
C + -------------------------------------------------
ii = MAX(i,2)
jj = MAX(j,2)
AUXlo1 = NST__x(ii ,jj )
AUXla1 = NST__y(ii ,jj )
AUXlo2 = NST__x(ii-1,jj-1)
AUXla2 = NST__y(ii-1,jj-1)
C +---Define the data points to be read around (i_cent,j_cent)
C + --------------------------------------------------------
G_nx=MAX(NINT(ABS(AUXlo1-AUXlo2)/reso_x),0)
G_ny=MAX(NINT(ABS(AUXla1-AUXla2)/reso_y),0)
i1=i_cent-G_nx
i2=i_cent+G_nx
j1=j_cent-G_ny
j2=j_cent+G_ny
! not to go out of domain
i1=MAX(i1,1)
i2=MIN(i2,cx)
j1=MAX(j1,1)
j2=MIN(j2,cy)
C +---Read subset of data
C + -------------------
nsamp =0
laisum=0.
DO l=j1,j2 ! Loop on data grid points
DO k=i1,i2 ! contained in the (i,j) NST cell
notin = 0
IF (MODIS_lai(k,cy-l+1).ge.0
. .and. MODIS_lai(k,cy-l+1).le.20) THEN
DO n=1,nvx
IF (MODIS_lu(k,cy-l+1) .eq. NSTsvt(k,cy-l+1,n)) Then
laisum(n) = laisum(n)+MAX(0.,MODIS_lai(k,cy-l+1))
nsamp(n) = nsamp(n)+1
else
notin = notin +1
endif
enddo
if (notin .eq. 3) nerror = nerror +1
ENDIF
ENDDO
ENDDO
DO l=1,nvx
NSTlai(i,j,l)=min(10.,laisum(l)/nsamp(l)) ! interpolate to NST grid
if (l.eq.3) NSTlai(i,j,l) = 0.75
if ((NSTsvt(i,j,l)>0.and.NSTsvt(i,j,l)<13)
. .and. NSTlai(i,j,l).eq.0) then
write(6,*) "No modis_lu NST !!!", i , j, l
endif
ENDDO
! -----------------------------------------------
DO l=1,nvx !For each vegetation type, we define a LAI max --> is
!also in GLOveg.f
IF (NSTsvt(i,j,l).eq. 0) NSTlmx(i,j,l) = 0.0
IF (NSTsvt(i,j,l).eq. 1) NSTlmx(i,j,l) = 0.6
IF (NSTsvt(i,j,l).eq. 2) NSTlmx(i,j,l) = 0.9
IF (NSTsvt(i,j,l).eq. 3) NSTlmx(i,j,l) = 1.2
IF (NSTsvt(i,j,l).eq. 4) NSTlmx(i,j,l) = 0.7
IF (NSTsvt(i,j,l).eq. 5) NSTlmx(i,j,l) = 1.4
IF (NSTsvt(i,j,l).eq. 6) NSTlmx(i,j,l) = 2.0
IF (NSTsvt(i,j,l).eq. 7.or.NSTsvt(i,j,l).eq.10)
. NSTlmx(i,j,l) = 3.0
IF (NSTsvt(i,j,l).eq. 8.or.NSTsvt(i,j,l).eq.11)
. NSTlmx(i,j,l) = 4.5
IF (NSTsvt(i,j,l).eq. 9.or.NSTsvt(i,j,l).eq.12)
. NSTlmx(i,j,l) = 6.0
ENDDO
! -----------------------------------------------
DO l=1,nvx
! NSTlai(i,j,l) = NSTlai(i,j,l) *
!. max(1.,min(2.,(1.+(NSTlmx(i,j,l)-3.)/12.)))
! MERRA lai = mean lai over 50 x 50 km2
! it is a bit corrected here in fct of vegetation.
!NSTlai(i,j,l) =max(0.,min(1.25*NSTlmx(i,j,l),NSTlai(i,j,l)))
! maximum values are a bit too low in respect to literature
if(NSTsvt(i,j,l)<=0.or.NSTsvt(i,j,l)==13) then
NSTlai(i,j,l) = 0.0
NSTglf(i,j,l) = 0.0
endif ! city or bare soil or ice
ENDDO
endif
! -----------------------------------------------
C + ****
ELSE ! Ocean, ice, snow
C + ****
DO l=1,nvx
NSTlai(i,j,l) = 0.0
NSTglf(i,j,l) = 0.0
ENDDO
C + *****
ENDIF ! Continental areas
C + *****
ENDDO
ENDDO
write(6,*)"Number of errors", nerror
END SUBROUTINE
!--------------------------------------------------------------------------------------------------------------------------
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment