diff --git a/couplage/CARAIB/ver01_Iv_couplage/caraib_main_01_Iv_couplage.F b/couplage/CARAIB/ver01_Iv_couplage/caraib_main_01_Iv_couplage.F
index 3871ea7029e576a00e9b23e48dd5e239baf1cf7b..c6f8f5e035e13d4cde9d31cd94fecb1741eb588f 100644
--- a/couplage/CARAIB/ver01_Iv_couplage/caraib_main_01_Iv_couplage.F
+++ b/couplage/CARAIB/ver01_Iv_couplage/caraib_main_01_Iv_couplage.F
@@ -47,6 +47,10 @@ c ('parameter.common' file)
c
c=====================================================================72
USE MOD_NETCDFCARAIB
+ use marctr
+ use mar_ge
+ use marcar
+
IMPLICIT NONE
c// BEGIN
@@ -84,6 +88,7 @@ c-----------------------------------------------
include './com_18/ecopro.common'
include './com_18/envi.common'
include './com_18/files_car.common'
+ include './com_18/files_ext.common'
include './com_18/files_ibm.common'
include './com_18/fileunits.common'
include './com_18/flux_w.common'
@@ -125,18 +130,34 @@ c default value of nd (number of days in year)
c time step for hydrological module (IBM) integration (day)
stept = 1.
+ nyear = 1 !It will always be runed 1Y until changed
+
call open_input5(nyrmax,itmt,stept,iread)
write(*,*) 'open input file finish'
+
+ ! If this is not the first iteration, then you can read the spinup data
+ if (first_iter.eq.0) then
+ ifrac = 0
+ ifrac_rd = 1
+ ilai_rd = 1
+ iread = 1
- if (ifull.eq.1) open(24,file=filres)
+ endif
+ if (ifull.eq.1) open(99924,file=filres)
+ write(*,*) 'ifull file finish'
c=======================================================================
c constants
c=======================================================================
- if (iterun .eq.1) then
+ call open_file(iread) ! open_file opens all the files and
+ ! initialises the longitude (zlon)
+ ! and latitude (zlat) arrays
+
+ if (first_iter.eq.1) then
+
call ctgen
write(*,*) 'ctgen finish'
-
+
call monthparam
write(*,*) 'monthparam finish'
@@ -148,7 +169,7 @@ c=======================================================================
call crops_seas
write(*,*) 'crop_seas finish'
-
+
c=======================================================================
c calculates hourly (only) dependent parameters.
c=======================================================================
@@ -158,33 +179,31 @@ c=======================================================================
if(itmt.ne.0) tmt_it = nyear/20.
- call open_file(iread) ! open_file opens all the files and
- ! initialises the longitude (zlon)
- ! and latitude (zlat) arrays
-
+
+
+
if ((imig.eq.1).and.(igtyp.ne.0)) then
write(*,*)'A regular rectangular grid should be used '
write(*,*)'when the migration module is called : '
write(*,*)'parameters igtyp=',igtyp,' and imig=',imig
write(*,*)'are not compatible - Program stop'
- write(61,*)'A regular rectangular grid should be used '
- write(61,*)'when the migration module is called : '
- write(61,*)'parameters igtyp=',igtyp,' and imig=',imig
- write(61,*)'are not compatible - Program stop'
+ write(99961,*)'A regular rectangular grid should be used '
+ write(99961,*)'when the migration module is called : '
+ write(99961,*)'parameters igtyp=',igtyp,' and imig=',imig
+ write(99961,*)'are not compatible - Program stop'
stop
endif
DO ngt = 1, n_pix ! Read in soil and lon-lat data
call read_eco(ngt)
ENDDO
-
c=======================================================================
c Initializes the CO2 pressure of the previous years
c=======================================================================
if ((iread.ge.1)) then
- read(325,*)(co2_prev(iy),iy=1,20)
+ read(999325,*)(co2_prev(iy),iy=1,20)
endif
c=======================================================================
@@ -199,18 +218,17 @@ c=======================================================================
iprint=0
call generator(graine,iprint)
else
- read(23) rapportP
- read(23) rapportT
- read(23) rapportDT
- read(23) ioccu
- read(23) nombrejp
- read(23) nombrejpS
+ read(99923) rapportP
+ read(99923) rapportT
+ read(99923) rapportDT
+ read(99923) ioccu
+ read(99923) nombrejp
+ read(99923) nombrejpS
endif
endif
- endif !(iterun=1)
+ endif
-
c=======================================================================
c Calculates the orbital parameters for a given situation/year
c=======================================================================
@@ -225,7 +243,6 @@ c=======================================================================
ELSE
ntimedata = nd
ENDIF
-
ALLOCATE(tcel1year(ntimedata, n_pix))
ALLOCATE(temax1year(ntimedata, n_pix))
@@ -241,7 +258,7 @@ c=======================================================================
if (num_ncdf.ge.1) then
- IF (iterun == 1) THEN ! At the beginning of the first year
+ IF (first_iter.eq.1) THEN ! At the beginning of the first year
ALLOCATE(ngt4ilonjlat(nclen_lon,nclen_lat))
ALLOCATE(ilon4ngt(n_pix))
@@ -294,7 +311,7 @@ c=======================================================================
endif !(num_ncdf)
c=======================================================================
-c i dont know
+c Read or generate temperature ?
c=======================================================================
@@ -317,7 +334,7 @@ c=======================================================================
endif !(incdf_tem)
c=======================================================================
-c i dont know
+c Read or generate temperature max ?
c=======================================================================
if (incdf_dta.eq.1) then
@@ -340,7 +357,7 @@ c=======================================================================
endif !(incdf_dta)
c=======================================================================
-c i dont know
+c Read or generate temperature min ?
c=======================================================================
@@ -363,7 +380,7 @@ c=======================================================================
endif !(idtem)
c=======================================================================
-c i dont know
+c Read or generate precip ?
c=======================================================================
@@ -387,7 +404,7 @@ c=======================================================================
endif
c=======================================================================
-c i dont know
+c Read or generate sunshine hour ?
c=======================================================================
if (incdf_shi.eq.1) then
@@ -408,7 +425,7 @@ c=======================================================================
endif
c=======================================================================
-c i dont know
+c Read or generate Relative humidity ?
c=======================================================================
if (incdf_rhu.eq.1) then
@@ -428,7 +445,7 @@ c=======================================================================
endif
c=======================================================================
-c i dont know
+c Read or generate windspeed ?
c=======================================================================
if (incdf_win.eq.1) then
@@ -453,10 +470,10 @@ c=======================================================================
c#######################################################################
c=======================================================================
-c Loop over the pixels for ....
+c Loop over the pixels for preprocess
c=======================================================================
c#######################################################################
-
+ write(*,*) "first pixel loop"
DO ngt = 1, n_pix
c=======================================================================
@@ -464,20 +481,30 @@ c calculates or reads :
c - coordinates of pixel corners
c - pixel area
c=======================================================================
- if (iterun.eq. 1) call pixel_corners(ngt)
-
+ if (first_iter.eq.1) call pixel_corners(ngt)
+
+
c=======================================================================
c reads land use file and
c initialization of variables for all pixels
c=======================================================================
call read_lu(ngt,iread)
-
- if (iterun.eq. 1) then
- call init_frac_and_lai(ngt)
- call read_init(iread,ngt)
- endif
-
+
+ ! if ((first_iter.eq.1) .and. (readsteady.eq.0)) then ! now set to 0 --> only do it first time
+ call init_frac_and_lai(ngt)
+
+ call read_init(iread,ngt)
+
+ ! else ! Then after the first day, it will read the results of the previous year
+ ! iread = 1
+ ! ifrac_rd = 1
+ ! ilai_rd = 1
+ ! call init_frac_and_lai(ngt)
+ ! call read_init(iread,ngt)
+ ! endif
+
+
c=======================================================================
c land use change and associated carbon reservoir changes
c=======================================================================
@@ -492,7 +519,7 @@ c reads cover fractions for crops and pastures
c=======================================================================
if ((ifrac_rd.ne.1).and.(ilu.eq.1)) call read_luspecies(ngt)
-
+
c=======================================================================
c reads climatological means (for Koppen zone calculation)
c tcel_clim, tcelkop : temperature (deg.),
@@ -502,7 +529,7 @@ c=======================================================================
if (num_ncdf.ge.1) then
- cpos_lon = ilon4ngt(ngt)
+ ncpos_lon = ilon4ngt(ngt)
ncpos_lat = jlat4ngt(ngt)
istartn = (/ ncpos_lon, ncpos_lat, 1 /)
@@ -549,15 +576,14 @@ c=======================================================================
call zonepxl2(iprint)
else
if(nyear.eq.1) then
- read(22,*)xxx,yyy,izonepxl
+ read(99922,*)xxx,yyy,izonepxl
else
izonepxl=kzone(ngt)
endif
endif
kzone(ngt)=izonepxl
- endif ! endif idaily_in
-
+ endif ! endif idaily_in
end do ! end of loop on pixels
c#######################################################################
@@ -566,7 +592,7 @@ c i dont know
c=======================================================================
c#######################################################################
- if(iczon.eq.1.or.iterun.eq.1) close(22)
+ if(iczon.eq.1.or.first_iter.eq.1) close(99922)
if (idaily_in.eq.0) then
@@ -595,35 +621,35 @@ c#######################################################################
if (num_ncdf.ge.1) then
ncpos_lon = 1 ! Reset ncpos_lon and ncpos_lat
ncpos_lat = 1
- endif
+ endif
c=======================================================================
c calculates pixel neighborhood
c=======================================================================
- if (iterun.eq.1) then
+ if (first_iter.eq.1) then
do ngt = 1, n_pix
call neighborhood_old(ngt)
end do
if(imig.eq.1) then
if(isp.eq.8) then
- open(401,file='midpoints_picea.txt')
- open(402,file='probdenswind_picea.txt')
+ open(999401,file='midpoints_picea.txt')
+ open(999402,file='probdenswind_picea.txt')
do i = 1, 96
- read(401,*) igrr,midpoint_pi(i)
- read(402,*) igrr,probdenswind_pi(i)
+ read(999401,*) igrr,midpoint_pi(i)
+ read(999402,*) igrr,probdenswind_pi(i)
enddo
- close(401)
- close(402)
+ close(999401)
+ close(999402)
endif
endif
endif !(iterun)
if (((imig.eq.1).and.(isteady.eq.0)).or.
- & ((imig.eq.1).and.(isteady.eq.1).and.(iterun.ne.1)))
+ & ((imig.eq.1).and.(isteady.eq.1).and.(first_iter.ne.1)))
& call read_mig
c#######################################################################
@@ -631,8 +657,12 @@ c=======================================================================
c loop on grid cells for subroutines
c=======================================================================
c#######################################################################
-
+!!$OMP PARALLEL DO default(shared)
+!!$OMP. private(ngt, igr, iread, izonepxl, niter, ipr, nprt, ny0prv)
+!!$OMP. private(ny0max, nyr_t, nyr_t2, tbegin, tend, dn, ny0)
+ write(*,*) "second pixel loop"
do ngt = 1, n_pix
+ write(*,*) ngt
c=======================================================================
c reads climatic, vegetation and soil data
@@ -640,6 +670,7 @@ c and "sets" them into common blocks
c=======================================================================
igr = ngt
call read_in(ngt,iread)
+
if (jclonly.eq.1) go to 2455
c=======================================================================
c climate zone corresponding to the studied pixel
@@ -664,30 +695,39 @@ c=======================================================================
nprt=nstprt-1
ny0prv = 0
+c=======================================================================
+c reads or evaluates sowing dates for crops if standard values from
+c species file ("bagseas") are not used
+c=======================================================================
+
+ if (ilu.eq.1) call crop_sowing_dates
+
c=======================================================================
c Set the number of year for spin up
c If for iteration : do X loop else only do it once
c=======================================================================
- if ((itepex.eq.1).and.(isteady.eq.1).and.(readsteady.eq.0)) then
+ if ((first_iter.eq.1).and.(isteady.eq.1).and.
+ & (readsteady.eq.0)) then
nyr_t = 0
+ ! ny0max = 2
else
ny0max = 1
nyr_t = 1
- if ((readsteady.eq.0).and. (itepex.eq.1)) call read_spinup
+ if ((readsteady.eq.1).and. (first_iter.eq.1)) nyr_t = 1!call read_spinup
endif
-c=======================================================================
-c reads or evaluates sowing dates for crops if standard values from
-c species file ("bagseas") are not used
-c=======================================================================
-
- if (ilu.eq.1) call crop_sowing_dates
-
+
! if first time : do multiple year == spinup
+c#######################################################################
+c=======================================================================
+c loop on years for spinup and run
+c=======================================================================
+c#######################################################################
do ny0 = 1,ny0max
nyr_t2 = 0
- if(itepex.eq.1 .and. ny0.eq.1 .and. iread.ge.1) nyr_t2 = 1
+ if (ny0 .eq. ny0max) nyr_t = 1
+ if(first_iter.eq.1.and. ny0.eq.1 .and.iread.ge.1) nyr_t2= 1
c=====================================================================72
c Call subroutine set_frac to set vegetation fraction
@@ -783,6 +823,20 @@ c=======================================================================
c=======================================================================
c grid point results
c=======================================================================
+ ! if (ny0.eq.ny0max-1) then
+ ! call close_file(iread)
+ ! filexto = trim(CHAR(iyrrGE-1))
+ ! call open_file(iread)
+ ! call wri_res
+ ! endif
+ ! if ((ny0.eq.ny0max) .and.(jdarGE.eq.jda0GE+1)) then
+ ! call close_file(iread)
+ ! filexto = trim(CHAR(iyrrGE-1))
+ ! call open_file(iread)
+ ! call wri_res
+ ! else if ((ny0.eq.ny0max) .and.(jdarGE.ne.jda0GE+1)) then
+ ! call wri_res
+ ! endif
if (ny0.eq.ny0max) then
call wri_res
@@ -800,18 +854,22 @@ c writes net budgets for last in output nr 27
c (test of convergence file)
c=====================================================================72
- call wri_1st(y,nyear,ngt)
-
+ call wri_1st(y,nyear,ngt)
+ call fill_table(ngt)
+
2455 continue
enddo ! end of loop on ngt
+!!$OMP END PARALLEL DO
+ call write_res
+ call make_marglf
c=======================================================================
c records CO2 of previous years
c=======================================================================
- write(326,'(20(1x,f10.3))')(co2_prev(iy),iy=1,20)
-
+ !write(999326,'(20(1x,f10.3))')(co2_prev(iy),iy=1,20)
+
c=======================================================================
c dispersion module
c=======================================================================
@@ -852,11 +910,13 @@ cc idayct = idayct + nd
DEALLOCATE(rhu1year)
DEALLOCATE(win1year)
- if(iclim.eq.1) close(666)
- close(5)
- close(28)
- close(61)
- if (ifull.eq.1) close(24)
-
- stop
+ if(iclim.eq.1) close(999666)
+ close(9995)
+ close(99928)
+ close(99961)
+ if (ifull.eq.1) close(99924)
+
+ ! call system('cp ./CARAIB_couplage/results/Belgium/*
+ ! . /climato_tmp1/tdethinne/res_car/')
+
end subroutine caraib
\ No newline at end of file