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

Update file make_eco_txt copy.f90

parent 8237c3bb
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
couplage/MAR/make_eco_txt copy.f90 deleted 100644 → 0
+ 0
− 312
View file @ 8237c3bb
! +------------------------------------------------------------------+
! | Subroutine Make_eco 12/04/2024 INTERPOL |
! +------------------------------------------------------------------+
! | |
! | This routine purpose us to create the CARIB domain from the MAR |
! | one. It create the ecotxt file that is read by CARAIB. The info |
! | that are copied are the LON, LAT, FAO number, %of argile, silt, |
! | sand as well as SH and two other undifined vars. |
! | |
! | Input : mar_lon(mx, my) : Input MAR points position x(i,j) |
! | ^^^^^^^ mar_lat(mx, my) : " " " " y(i,j) |
! | FAO_nbr(mx, my) : FAO number of soil type |
! | clay(mx, my) : % of clay in soil |
! | sand(mx, my) : % of sand in soil |
! | silt(mx, my) : % of silt in soil |
! | sh(mx, my) : Elevation above see level |
! | avg_col(mx, my) : average color of soil |
! | fland(mx, my) : ? |
! | |
! | |
! | Output: ecotxt.dat |
! | ^^^^^^^ |
! | |
! +------------------------------------------------------------------+
!SUBROUTINE
program mar2car!(mar_lon,mar_lat,FAO_nbr,clay,sand,silt,sh,avg_col,fland)
use mardim
IMPLICIT NONE
! +---General and local variables
! + ---------------------------
INTEGER i,j,nbr_pixel,iostat
REAL mar_lon(mx,my), mar_lat(mx,my), FAO_nbr(mx,my), &
clay(mx,my), sand(mx,my), silt(mx,my), &
sh(mx,my), avg_col(mx,my), fland(mx,my)
character(len=50) :: filename = '../CARAIB_couplage/clim/Belgium/ecotxt.dat'
data ngr_lk /4320,576/
data cla / 9.,30.,66.,20.,38.,48.,35.,48.,30./
data sil / 8.,33.,17.,20.,12.,25.,19.,25.,33./
data san /83.,37.,17.,60.,50.,27.,46.,27.,37./
data hwsd_lat /-60.,-40.,-30.,-20.,-10.,0.,10.,20.,30.,40.,50.
& ,60.,70.,90./
data hwsd_npix / 1935182, 5754888,11868114,11329216,12087991
& ,11723098,12801602,17929486,21392072,25909962
& ,28717217,33829821,10052296/
pi = 3.1415926535
! + ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
! +---mock fill
! + -----------------
nbr_pixel = 1
inquire(file=filename, exist=iostat)
if (iostat == 0) then
! File does not exist, create it
open(unit=20, file=filename, status='unknown', action='write', iostat=iostat)
if (iostat /= 0) then
print *, 'Error creating file ', trim(filename)
stop
else
ngr_1=15347 ! number of line in global data
ngr= cx*cy ! number of pixel in grid
& ,nlat=1080,nlong=2160
& ,nlat1=180,nlong1=360
& ,nlat_hwsd=20*120,nlong_hwsd=360*120
& ,nlake=2,ntile=13)
write(*,*) 'Initialization of variables'
do ilt=1,nlat
do ilg=1,nlong
elev_land(ilg,ilt)=-999.
elev_pix(ilg,ilt)=-999.
f_land(ilg,ilt)=-999.
end do
end do
do ilt1=0,nlat1+1
do ilg1=0,nlong1+1
col(ilg1,ilt1)=-999.
sand(ilg1,ilt1)=-999.
clay(ilg1,ilt1)=-999.
silt(ilg1,ilt1)=-999.
isu(ilg1,ilt1)=-99
end do
end do
print *, "creating the file"
open(9,file=trim(dirhome)//'../data/ecotx.1x1')
write(*,*)'Reads texture file at 1x1 deg and assigns to array'
reso_lg1=360./float(nlong1)
reso_lt1=180./float(nlat1)
xlg10=-180.-reso_lg1
xlt10=90.+reso_lt1
do ilg1=0,nlong1+1
xlg_1(ilg1) = -180.-reso_lg1/2.+float(ilg1)*reso_lg1
end do
do ilt1=0,nlat1+1
xlt_1(ilt1) = 90.+reso_lt1/2.-float(ilt1)*reso_lt1
end do
do igr = 1, ngr_1
read(9,*)ig,xlg1,xlt1,isun,itxt,elv,icol
ilg1 = int((xlg1-xlg10)/reso_lg1+0.001)
ilt1 = int((xlt10-xlt1)/reso_lt1+0.001)
! ------ colour: soil colour (0:dark - 1:light) ----------------
colour = 0.5
if((icol.gt.10).and.(icol.lt.34))then
if((icol.le.16).or.(icol.eq.29))then
colour=1.
elseif ((icol.le.22).or.(icol.eq.30))then
colour=0.5
elseif ((icol.le.28).or.(icol.eq.31))then
colour=0.
endif
endif
! -------- soil texture -------------
! ------ glace permanente: itxt=9 (medium texture:2) -----
if (itxt.eq.9) itxt=2
if ((itxt.eq.999).or.(itxt.eq.250)) itxt=9
sandy = san(itxt)
silty = sil(itxt)
clayy = cla(itxt)
! -------------------------------------------
if ((ilg1.gt.nlong1).or.(ilg1.lt.1))write(*,*)'ilg1:',ilg1
if ((ilt1.gt.nlat1).or.(ilt1.lt.1))write(*,*)'ilt1:',ilt1
col(ilg1,ilt1)=colour
sand(ilg1,ilt1)=sandy
clay(ilg1,ilt1)=clayy
silt(ilg1,ilt1)=silty
isu(ilg1,ilt1)=isun
end do
call fillborder(col,nlong1,nlat1)
call fillborder(sand,nlong1,nlat1)
call fillborder(clay,nlong1,nlat1)
call fillborder(silt,nlong1,nlat1)
DO i=1,mx
DO j=1,my
mar_lon(mx,my) = 2
mar_lat(mx,my) = 3
sh(mx,my) = 8
xlg = xlg_WC(mx,my)
xlt = xlt_WC(mx,my)
ilg = int((xlg-xlg0)/reso_lg+0.001)
ilt = int((xlt0-xlt)/reso_lt+0.001)
elv_p=elev_pix(ilg,ilt)
elv_l=elev_land(ilg,ilt)
f_l=1.
if(elv_l.lt.-1.)then
write(21,120)xlg,xlt,f_l,elv_p
elv_l=0.
endif
ilg1 = int((xlg-xlg10)/reso_lg1+0.001)
ilt1 = int((xlt10-xlt)/reso_lt1+0.001)
isuy=isu(ilg1,ilt1)
if (isuy.lt.-1.) isuy=10
if (xlg.lt.xlg_1(ilg1)) then
ilg1_bef=ilg1-1
ilg1_aft=ilg1
else
ilg1_bef=ilg1
ilg1_aft=ilg1+1
endif
if (xlt.ge.xlt_1(ilt1)) then
ilt1_bef=ilt1-1
ilt1_aft=ilt1
else
ilt1_bef=ilt1
ilt1_aft=ilt1+1
endif
av_w = 0.
av_col = 0.
av_san = 0.
av_cla = 0.
do jlg1 = ilg1_bef, ilg1_aft
do jlt1 = ilt1_bef, ilt1_aft
if (col(jlg1,jlt1).ge.-0.001) then
dist=sqrt( (xlt-xlt_1(jlt1))**2.
& +(cos(xlt*pi/180.)*(xlg-xlg_1(jlg1)))**2.)
if (dist.le.1.e-5) then
av_col=col(jlg1,jlt1)
av_san=sand(jlg1,jlt1)
av_cla=clay(jlg1,jlt1)
go to 2000
else
av_w=av_w+1./dist
av_col=av_col+col(jlg1,jlt1)/dist
av_san=av_san+sand(jlg1,jlt1)/dist
av_cla=av_cla+clay(jlg1,jlt1)/dist
endif
endif
end do
end do
if (av_w.lt.1.e-10) then
av_col=0.5
av_san=san(2)
av_cla=cla(2)
write(22,120)xlg,xlt,f_l,elv_l
else
av_col=av_col/av_w
av_san=av_san/av_w
av_cla=av_cla/av_w
endif
2000 continue
if (av_col.lt.0.) av_col=0.
if (av_col.gt.1.) av_col=1.
if (av_san.lt.0.) av_san=0.
if (av_cla.lt.0.) av_cla=0.
av_sil = 100.-av_san-av_cla
if (av_sil.lt.0.) then
av_tot = av_san+av_cla
av_san = av_san/av_tot
av_cla = av_cla/av_tot
av_sil = 0.
endif
if (hwsd_clay(mx,my) .lt.0.) then
write(20,123)igr,xlg,xlt,isuy,av_cla,av_sil
& ,av_san,elv_WC(mx,my)
& ,av_col,f_l
else
write(20,123)igr,xlg,xlt,isuy,hwsd_clay(mx,my)
& ,hwsd_silt(mx,my) ,hwsd_sand(mx,my)
& ,elv_WC(mx,my)
& ,av_col,f_l
endif
FAO_nbr(mx,my) = 4
clay(mx,my) = 5
sand(mx,my) = 6
silt(mx,my) = 7
avg_col(mx,my) = 9
fland(mx,my) = 10
! write(20,123)nbr_pixel ,mar_lon(mx,my),mar_lat(mx,my),FAO_nbr(mx,my), &
! clay(mx,my),sand(mx,my),silt(mx,my), &
! sh(mx,my),avg_col(mx,my),fland(mx,my)
123 format(i7,2(1x,f9.4),1x,i3,3(1x,f6.2),1x,f5.0,1x,f5.2,1x,f6.3)
nbr_pixel = nbr_pixel + 1
enddo
enddo
endif
else
print *, 'File ', trim(filename), ' already exists.'
endif
end program
!end subroutine
!=======================================================================
subroutine fillborder(vari,nlong,nlat)
!=======================================================================
real*4 vari(0:nlong+1,0:nlat+1)
do ilong = 1, nlong
vari(ilong,0) = vari(ilong,1)
vari(ilong,nlat+1) = vari(ilong,nlat)
end do
do ilat = 0, nlat+1
vari(0,ilat) = vari(nlong,ilat)
vari(nlong+1,ilat) = vari(1,ilat)
end do
return
end
\ No newline at end of file
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