Context Navigation

source: branches/ORCHIDEE_2_2/ORCHIDEE/src_driver/dim2_driver.f90 @ 7264

Last change on this file since 7264 was 7264, checked in by agnes.ducharne, 3 years ago
Makes all keywords related to SPREAD_PREC include an A in spreAd, cf ticket #730, comment 5 and on. All forms with SPRED_PREC get deprecated. The file orchidee.default is adapted, and now includes info on SPREAD_PREC_SEC and SPREAD_PREC_CONT (only used by new driver). The configurations, however, have not yet been updated (the only committed configuration affected is FG3nd with SPRED_PREC_SEC = 5400).
Property svn:keywords set to `HeadURL Date Author Revision`
File size: 64.8 KB

Rev	Line
[4263]	1	! ==============================================================================================================================\n
	2	! PROGRAM : driver
	3	!
[4470]	4	! CONTACT : orchidee-help _at_ listes.ipsl.fr
[4263]	5	!
	6	! LICENCE : IPSL (2006)
	7	! This software is governed by the CeCILL licence see ORCHIDEE/ORCHIDEE_CeCILL.LIC
	8	!
	9	!>\BRIEF Main program for the offline driver dim2_driver
	10	!!
	11	!!\n DESCRIPTION : This PROGRAM is the driver of the dim2 version of ORCHIDEE.
	12	!!
	13	!! The forcing data needs to be in netCDF format and should contain the following variables :
	14	!! - Incoming SW radiation
	15	!! - Incoming LW radiation
	16	!! - Precipitation
	17	!! - Air temperature at a reference level
	18	!! - Air humidity at the same level
	19	!! - wind at the same level
	20	!! - surface pressure
	21	!!
	22	!! RECENT CHANGE(S): None
	23	!!
	24	!! REFERENCE(S) : None
	25	!!
	26	!! SVN :
	27	!! $HeadURL$
	28	!! $Date$
	29	!! $Revision$
	30	!! \n
	31	!_ ================================================================================================================================
	32
[8]	33	PROGRAM driver
[4263]	34
[8]	35	USE netcdf
[1078]	36	USE ioipsl_para
[8]	37	USE grid
[4287]	38	USE intersurf, ONLY : intersurf_main_2d, intersurf_initialize_2d, intersurf_clear
[8]	39	USE constantes
[4646]	40	USE time
[8]	41	USE readdim2
[1078]	42	USE mod_orchidee_para
[8]	43	USE timer
[4263]	44
[8]	45	IMPLICIT NONE
[4263]	46
[8]	47	INTEGER :: iim, jjm, llm
	48	INTEGER :: im, jm, lm, tm, is, force_id, itest, jtest
	49	REAL :: dt, dt_force, dt_rest, date0, date0_rest
[7256]	50	REAL :: date_cur ! Starting date of the current execution
[8]	51	REAL :: zlflu
	52	REAL :: alpha
[4263]	53
[8]	54	REAL, ALLOCATABLE, DIMENSION(:,:) :: &
[2961]	55	& swdown, coszang, precip_rain, precip_snow, tair_obs, u, v, &
[8]	56	& qair_obs, pb, lwdown, &
	57	& eair_obs, zlev_vec, zlevuv_vec, relax
	58	!- Variables which are forcings for SECHIBA
	59	REAL, ALLOCATABLE, DIMENSION(:,:) :: &
	60	& petAcoef, peqAcoef, petBcoef, peqBcoef, cdrag, &
[2961]	61	& for_u, for_v, for_swnet, for_swdown, for_coszang, for_lwdown, &
[8]	62	& for_psurf, for_qair, for_tair, for_eair, &
	63	& for_ccanopy, for_rau
[4263]	64
[8]	65	REAL, ALLOCATABLE, DIMENSION(:,:) :: &
	66	& for_contfrac, old_zlev, old_qair, old_eair, tsol_rad, vevapp, &
	67	& temp_sol_NEW, temp_sol_old, qsurf, dtdt, coastalflow, riverflow, &
	68	& fluxsens, fluxlat, emis, z0
[4263]	69
[8]	70	INTEGER, ALLOCATABLE, DIMENSION(:,:,:) :: for_neighbours
[4263]	71
[8]	72	REAL, ALLOCATABLE, DIMENSION(:,:,:) :: for_resolution
[4263]	73
[8]	74	REAL, ALLOCATABLE, DIMENSION(:,:,:) :: albedo
	75	REAL, ALLOCATABLE, DIMENSION(:,:) :: albedo_vis
	76	REAL, ALLOCATABLE, DIMENSION(:,:) :: albedo_nir
[4263]	77
[8]	78	INTEGER, ALLOCATABLE, DIMENSION(:) :: kindex
	79	REAL, ALLOCATABLE, DIMENSION(:,:) :: lon, lat
[1224]	80	REAL, ALLOCATABLE, DIMENSION(:) :: tmplev
[4263]	81
[8]	82	REAL :: old_tair
	83	REAL :: atmco2
	84	INTEGER :: nbindex
	85	REAL :: julian, ss
	86	INTEGER :: yy, mm, dd
[4263]	87
[1224]	88	LOGICAL :: relaxation
[4263]	89
[1056]	90	CHARACTER(LEN=80) :: filename, driv_restname_in, driv_restname_out
[8]	91	CHARACTER(LEN=30) :: time_str, var_name
[4263]	92
[8]	93	INTEGER :: it, istp, istp_old, rest_id, it_force
[4263]	94
[8]	95	INTEGER :: split, split_start, nb_spread, for_offset
	96	INTEGER :: itau_dep, itau_dep_rest, itau_fin, itau_skip, itau_len
	97
	98	INTEGER,DIMENSION(2) :: ml
[4263]	99
[3850]	100	LOGICAL :: lstep_init, lstep_last
[4629]	101	LOGICAL :: lwdown_cons !! Flag to conserve lwdown radiation from forcing
	102	LOGICAL :: swdown_cons !! Flag to conserve swdown radiation from forcing
[4263]	103
[8]	104	! to check variables passed to intersurf
	105	INTEGER :: ik
	106	INTEGER :: i,j
[2348]	107	INTEGER :: printlev_loc !! local write level
[8]	108
	109	REAL, ALLOCATABLE, DIMENSION(:,:) :: &
	110	& fluxsens_g,vevapp_g,old_zlev_g,old_qair_g,old_eair_g,for_rau_g, &
	111	& petAcoef_g, petBcoef_g,peqAcoef_g,peqBcoef_g,albedo_g,z0_g
	112	LOGICAL :: Flag
[1328]	113	LOGICAL :: driver_reset_time
[8]	114
	115	REAL :: fill_init
	116
	117	fill_init=REAL(nf90_fill_real,r_std)
	118	CALL ioconf_expval(val_exp)
	119	!-
	120	! Init parallelism
	121
[1693]	122	CALL Init_orchidee_para
[353]	123	CALL init_timer
[1693]	124	CALL start_timer(timer_global)
[1078]	125	CALL start_timer(timer_mpi)
[2293]	126
[2348]	127	! Set specific write level to dim2_driver using PRINTLEV_dim2_driver=[0-4]
	128	! in run.def. The global printlev is used as default value.
	129	printlev_loc=get_printlev('dim2_driver')
[2293]	130
[8]	131	!=====================================================================
	132	!- 1.0 This section defines the general set-up of the experiment :
	133	!- - Forcing data to be used with its time step
	134	!- - Restart file to be used
	135	!- - The time step that will be used
	136	!- - The starting date of the simulation
	137	!- - Length of integration
	138	!- - Basic flags for SSIPSL
[4629]	139
	140	!Config Key = LWDOWN_CONS
	141	!Config Desc = Conserve longwave downwelling radiation in the forcing
	142	!Config Def = n
	143	!Config If =
	144	!Config Help = If LWDOWN_CONS=False a non conservative interpolation is used
	145	!Config Units = [FLAG]
	146	lwdown_cons = .FALSE.
	147	CALL getin_p('LWDOWN_CONS', lwdown_cons)
	148
	149	!Config Key = SWDOWN_CONS
	150	!Config Desc = Conserve shortwave downwelling radiation in the forcing
	151	!Config Def = LWDOWN_CONS
	152	!Config If =
	153	!Config Help = If SWDOWN_CONS=False a non conservative interpolation is used
	154	!Config Units = [FLAG]
	155	swdown_cons = lwdown_cons
	156	CALL getin_p('SWDOWN_CONS', swdown_cons)
	157
[8]	158	!=====================================================================
	159	!- 1.1 Initialize the driving variables. It essentialy gets the mode
	160	!- used and the size of the driving variables.
	161	!=====================================================================
[4693]	162	IF (printlev_loc>=4) WRITE(numout,*) 'Reading name of the forcing file'
[8]	163	!-
[567]	164	!Config Key = FORCING_FILE
	165	!Config Desc = Name of file containing the forcing data
	166	!Config If = [-]
[621]	167	!Config Def = forcing_file.nc
[567]	168	!Config Help = This is the name of the file which should be opened
[8]	169	!Config for reading the forcing data of the dim0 model.
	170	!Config The format of the file has to be netCDF and COADS
	171	!Config compliant.
[844]	172	!Config Units = [FILE]
[8]	173	!-
[621]	174	filename='forcing_file.nc'
[8]	175	CALL getin_p('FORCING_FILE',filename)
	176	!-
[4693]	177	IF (printlev_loc>=4) WRITE(numout,*) 'Opening forcing file'
[8]	178	!-
	179	! We call flininfo to obtain the dimensions
	180	! of iim, jjm and others.
	181	! This information will allow us to allocate all the space needed.
	182	!-
	183	CALL forcing_info &
	184	& (filename, iim, jjm, llm, tm, date0, dt_force, force_id)
	185	!-
[4693]	186	IF (printlev>=1) WRITE(numout,*) 'Information about forcing file : date0 ', date0, &
[8]	187	'iim, jjm, llm, tm',iim,jjm,llm,tm,' dt_force ',dt_force
	188	!-
[1078]	189	CALL init_ioipsl_para
	190	!-
[4693]	191	IF (printlev_loc>=4) THEN
[8]	192	WRITE(numout,*) 'Allocate memory for the driver :', iim, jjm, llm
	193	ENDIF
	194	!-
[1224]	195	ALLOCATE (tmplev(llm))
[3019]	196	tmplev(:)=0.
	197
[8]	198	ALLOCATE &
[2961]	199	& (swdown(iim,jjm), coszang(iim,jjm), precip_rain(iim,jjm), precip_snow(iim,jjm), &
[8]	200	& tair_obs(iim,jjm), u(iim,jjm), v(iim,jjm), qair_obs(iim,jjm), &
	201	& pb(iim,jjm), lwdown(iim,jjm), &
	202	& eair_obs(iim,jjm), zlev_vec(iim,jjm), zlevuv_vec(iim,jjm), relax(iim,jjm))
	203	!-
	204	ALLOCATE &
	205	& (petAcoef(iim,jjm), peqAcoef(iim,jjm), &
	206	& petBcoef(iim,jjm), peqBcoef(iim,jjm), &
	207	& cdrag(iim,jjm), for_u(iim,jjm), for_v(iim,jjm), &
[2961]	208	& for_swnet(iim,jjm), for_swdown(iim,jjm), for_coszang(iim,jjm), for_lwdown(iim,jjm), &
[8]	209	& for_psurf(iim,jjm), for_qair(iim,jjm), for_tair(iim,jjm), &
	210	& for_eair(iim,jjm), for_ccanopy(iim,jjm), for_rau(iim,jjm))
[4263]	211
[8]	212	ALLOCATE &
	213	& (for_contfrac(iim,jjm), for_neighbours(iim,jjm,8), for_resolution(iim,jjm,2), &
	214	& old_zlev(iim,jjm), old_qair(iim,jjm), old_eair(iim,jjm), &
	215	& tsol_rad(iim,jjm), vevapp(iim,jjm), &
	216	& temp_sol_NEW(iim,jjm), temp_sol_old(iim,jjm), &
	217	& dtdt(iim,jjm), coastalflow(iim,jjm), riverflow(iim,jjm), &
	218	& fluxsens(iim,jjm), fluxlat(iim,jjm), emis(iim,jjm), &
	219	& z0(iim,jjm), qsurf(iim,jjm))
[4263]	220
[8]	221	ALLOCATE(albedo(iim,jjm,2))
	222	ALLOCATE(albedo_vis(iim,jjm),albedo_nir(iim,jjm))
	223	ALLOCATE(kindex(iim*jjm))
	224	ALLOCATE(lon(iim,jjm), lat(iim,jjm))
[4263]	225
[8]	226	swdown(:,:) = fill_init
	227	precip_rain(:,:) = 0.0
	228	precip_snow(:,:) = 0.0
	229	tair_obs(:,:) = 0.0
	230	u(:,:) = fill_init
	231	v(:,:) = fill_init
	232	qair_obs(:,:) = fill_init
	233	pb(:,:) = fill_init
	234	lwdown(:,:) = fill_init
	235	eair_obs(:,:) = fill_init
	236	zlev_vec(:,:) = 0.0
	237	zlevuv_vec(:,:) = 0.0
	238	relax(:,:) = 0.0
	239	petAcoef(:,:) = 0.0
	240	peqAcoef(:,:) = 0.0
	241	petBcoef(:,:) = 0.0
	242	peqBcoef(:,:) = 0.0
	243	cdrag(:,:) = 0.0
	244	for_u(:,:) = fill_init
	245	for_v(:,:) = fill_init
	246	for_swnet(:,:) = fill_init
	247	for_swdown(:,:) = fill_init
	248	for_lwdown(:,:) = fill_init
	249	for_psurf(:,:) = fill_init
	250	for_qair(:,:) = fill_init
	251	for_tair(:,:) = fill_init
	252	for_eair(:,:) = fill_init
	253	for_ccanopy(:,:) = 0.0
	254	for_rau(:,:) = fill_init
	255	for_contfrac(:,:) = 0.0
	256	for_neighbours(:,:,:) = 0
	257	for_resolution(:,:,:) = 0.0
	258	old_zlev(:,:) = 0.0
	259	old_qair(:,:) = 0.0
	260	old_eair(:,:) = 0.0
	261	tsol_rad(:,:) = 0.0
	262	vevapp(:,:) = 0.0
	263	temp_sol_NEW(:,:) = fill_init
	264	temp_sol_old(:,:) = fill_init
	265	dtdt(:,:) = 0.0
	266	coastalflow(:,:) = 0.0
	267	riverflow(:,:) = 0.0
[1062]	268	fluxsens(:,:) = fill_init
	269	fluxlat(:,:) = fill_init
[8]	270	emis(:,:) = 0.0
	271	z0(:,:) = fill_init
	272	qsurf(:,:) = 0.0
	273	albedo(:,:,:) = fill_init
	274	albedo_vis(:,:) = fill_init
	275	albedo_nir(:,:) = fill_init
	276	kindex(:) = 0
	277	lon(:,:) = 0.0
	278	lat(:,:) = 0.0
	279	!-
	280	! We need to know the grid.
	281	! Then we can initialize the restart files, and then we
	282	! can give read the restart files in the forcing subroutine.
	283	!-
[1224]	284	CALL forcing_grid (iim,jjm,llm,lon,lat,init_f=.FALSE.)
[8]	285	!=====================================================================
	286	!- 1.2 Time step to be used.
[4629]	287	!- This is relative to the time step of the forcing data
[8]	288	!=====================================================================
	289	IF ( .NOT. weathergen ) THEN
[2575]	290	!Config Key = DT_SECHIBA
	291	!Config Desc = Time-step of the SECHIBA component
[567]	292	!Config If = NOT(WEATHERGEN)
[2575]	293	!Config Def = 1800.
	294	!Config Help = Determines the time resolution at which
	295	!Config the calculations in the SECHIBA component
	296	!Config are done
	297	!Config Units = [seconds]
	298	dt = 1800
	299	CALL getin_p('DT_SECHIBA', dt)
	300	split = INT(dt_force/dt)
[4693]	301	IF (printlev_loc >= 1) WRITE(numout,*) 'Time step in forcing file: dt_force=',dt_force
	302	IF (printlev_loc >= 1) WRITE(numout,*) 'Time step in sechiba component: dt_sechiba=',dt
	303	IF (printlev_loc >= 1) WRITE(numout,*) 'Splitting of each forcing time step: split=',split
[2575]	304
	305	IF ( split .LT. 1. ) THEN
	306	CALL ipslerr_p ( 3, 'dim2_driver',&
	307	'Time step of the forcing file is higher than the time step in sechiba',&
	308	'Please, modify DT_SECHIBA parameter value !','')
	309	END IF
[8]	310	ELSE
[2575]	311	! Case weathergen:
	312	! The model time step in sechiba is always the same as the forcing time step
	313	dt = dt_force
[8]	314	split = 1
	315	ENDIF
[2575]	316
[8]	317	!=====================================================================
	318	!- 1.3 Initialize the restart file for the driver
	319	!=====================================================================
[567]	320	!Config Key = RESTART_FILEIN
	321	!Config Desc = Name of restart to READ for initial conditions
	322	!Config If = [-]
	323	!Config Def = NONE
	324	!Config Help = This is the name of the file which will be opened
[8]	325	!Config to extract the initial values of all prognostic
	326	!Config values of the model. This has to be a netCDF file.
	327	!Config Not truly COADS compliant. NONE will mean that
	328	!Config no restart file is to be expected.
[844]	329	!Config Units = [FILE]
[8]	330	!-
[1056]	331	driv_restname_in = 'NONE'
	332	CALL getin_p('RESTART_FILEIN',driv_restname_in)
[4693]	333	if (printlev_loc>=4) WRITE(numout,*) 'INPUT RESTART_FILE : ',TRIM(driv_restname_in)
[8]	334	!-
[567]	335	!Config Key = RESTART_FILEOUT
	336	!Config Desc = Name of restart files to be created by the driver
	337	!Config If = [-]
	338	!Config Def = driver_rest_out.nc
	339	!Config Help = This variable give the name for
	340	!Config the restart files. The restart software within
	341	!Config IOIPSL will add .nc if needed
[844]	342	!Config Units = [FILE]
[8]	343	!-
[1056]	344	driv_restname_out = 'driver_rest_out.nc'
	345	CALL getin_p('RESTART_FILEOUT', driv_restname_out)
[4693]	346	if (printlev_loc>=4) WRITE(numout,*) 'OUTPUT RESTART_FILE : ',TRIM(driv_restname_out)
[8]	347	!-
	348	! Set default values for the start and end of the simulation
	349	! in the forcing chronology.
[1328]	350
[1078]	351	CALL gather2D_mpi(lon,lon_g)
	352	CALL gather2D_mpi(lat,lat_g)
[1328]	353
	354
	355	!Config Key = DRIVER_reset_time
	356	!Config Desc = Overwrite time values from the driver restart file
	357	!Config If = [-]
[2135]	358	!Config Def = n
[1328]	359	!Config Units = [FLAG]
[8]	360
[1328]	361	driver_reset_time=.FALSE.
	362	CALL getin_p('DRIVER_reset_time', driver_reset_time)
[4693]	363	IF (printlev_loc>=4) WRITE(numout,*) 'driver_reset_time=',driver_reset_time
[1328]	364
[8]	365	IF (is_root_prc) THEN
[1328]	366	! Set default values for the time variables
	367	itau_dep_rest = 0
	368	date0_rest = date0
	369	dt_rest = dt
	370
[4693]	371	IF (printlev_loc>=4) WRITE(numout,*) &
[1328]	372	'Before driver restart file initialization : itau_dep_rest, date0_rest, dt_rest = ', &
	373	itau_dep_rest, date0_rest, dt_rest
	374
[8]	375	CALL restini &
[1328]	376	(driv_restname_in, iim_g, jjm_g, lon_g, lat_g, llm, tmplev, &
	377	driv_restname_out, itau_dep_rest, date0_rest, dt_rest, rest_id, driver_reset_time)
	378
[4693]	379	IF (printlev_loc>=4) WRITE(numout,*) &
[1328]	380	'After driver restart file initialization : itau_dep_rest, date0_rest, dt_rest = ', &
	381	itau_dep_rest, date0_rest, dt_rest
	382
	383	IF (itau_dep_rest == 0 .OR. driver_reset_time) THEN
	384	! Restart file did not exist or we decide to overwrite time values in it.
	385	! Set time values to read the begining of the forcing file.
	386	itau_dep=0
	387	itau_fin=tm
	388	date0_rest = date0
	389	ELSE
	390	! Take time values from restart file
[8]	391	itau_dep = itau_dep_rest
[254]	392	itau_fin = itau_dep+tm
[8]	393	ENDIF
[1328]	394
[4693]	395	IF (printlev_loc >= 1) WRITE(numout,*) &
[1328]	396	'Restart file initialized : itau_dep, itau_fin, date0_rest, dt_rest = ', &
	397	itau_dep, itau_fin, date0_rest, dt_rest
[8]	398	ENDIF
[1328]	399
	400	! Communicate values from root_prc to all the others
[8]	401	CALL bcast (itau_dep_rest)
	402	CALL bcast (itau_dep)
	403	CALL bcast (itau_fin)
	404	CALL bcast (date0_rest)
	405	CALL bcast (dt_rest)
[1328]	406
[8]	407	!=====================================================================
	408	!- 1.4 Here we do the first real reading of the driving. It only
	409	!- gets a few variables.
	410	!=====================================================================
[1224]	411
[8]	412	! prepares kindex table from the information obtained
	413	! from the forcing data and reads some initial values for
	414	! temperature, etc.
	415	!-
	416	kindex(1) = 1
	417	!-
	418	CALL forcing_READ &
	419	& (filename, rest_id, .TRUE., .FALSE., &
[4629]	420	& 0, itau_dep, 0, split, nb_spread, lwdown_cons, swdown_cons, &
[8]	421	& date0, dt_force, iim, jjm, lon, lat, zlev_vec, zlevuv_vec, tm, &
[2961]	422	& swdown, coszang, precip_rain, precip_snow, tair_obs, &
[8]	423	& u, v, qair_obs, pb, lwdown, for_contfrac, for_neighbours, for_resolution, &
	424	& for_swnet, eair_obs, petAcoef, peqAcoef, petBcoef, peqBcoef, cdrag, for_ccanopy, &
	425	& kindex, nbindex, force_id)
	426	!-
[3941]	427	IF (printlev_loc >= 2) WRITE (numout,*) ">> Number of land points =",nbindex
[8]	428	IF (nbindex == 0) THEN
	429	WRITE(numout,*) "Limits : (W,E / N,S)", limit_west, limit_east, &
	430	& limit_north, limit_south
[1078]	431	CALL ipslerr_p ( 3, 'dim2_driver','number of land points error.', &
[8]	432	& ' is zero !','stop driver')
	433	ENDIF
[3447]	434	!-
	435	DO ik=1,nbindex
	436	jlandindex(ik) = (((kindex(ik)-1)/iim) + 1)
	437	ilandindex(ik) = (kindex(ik) - (jlandindex(ik)-1)*iim)
	438	ENDDO
[2348]	439	IF (printlev_loc>=4) THEN
[8]	440	WRITE(numout,*) "kindex of land points : ", kindex(1:nbindex)
	441	WRITE(numout,*) "index i of land points : ", ilandindex
	442	WRITE(numout,*) "index j of land points : ", jlandindex
	443	ENDIF
	444
	445	im = iim; jm = jjm; lm = llm;
	446	IF ( (iim > 1).AND.(jjm > 1) ) THEN
	447	jtest = INT((kindex(INT(nbindex/2))-1)/iim)+1
	448	itest = MAX( 1, kindex(INT(nbindex/2))-(jtest-1)*iim )
	449	ELSE
	450	jtest = 1
	451	itest = 1
	452	ENDIF
[2348]	453	IF (printlev_loc>=3) WRITE(numout,*) "test point in dim2_driver : ",itest,jtest
[8]	454	!-
	455	IF ((im /= iim) .AND. (jm /= jjm) .AND. (lm /= llm)) THEN
	456	WRITE(numout,*) ' dimensions are not good. Verify FILE :'
	457	WRITE(numout,*) ' filename = ',filename
	458	WRITE(numout,*) ' im, jm, lm lus = ', im, jm, lm
	459	WRITE(numout,*) ' iim, jjm, llm demandes = ', iim, jjm, llm
[2373]	460	CALL ipslerr_p(3,'dim2_driver','Pb in dimensions','','')
[8]	461	ENDIF
	462	!=====================================================================
	463	!- 1.5 Configures the time-steps and other parameters
	464	!- of the run.
	465	!=====================================================================
	466	!-
	467	! If the time steping of the restart is different from the one
	468	! of the forcing we need to convert the itau_dep into the
	469	! chronology of the forcing. This ensures that the forcing will
	470	! start at the date of the restart file. Obviously the itau_fin
	471	! needs to be shifted as well !
	472	!-
	473	IF ( (dt_rest /= dt_force).AND.(itau_dep > 1) ) THEN
	474	itau_dep = NINT((itau_dep*dt_rest )/dt_force)
[254]	475	itau_fin = itau_dep+tm
[2348]	476	if (printlev_loc>=3) WRITE(numout,*) &
[8]	477	& 'The time steping of the restart is different from the one ',&
	478	& 'of the forcing we need to convert the itau_dep into the ',&
	479	& 'chronology of the forcing. This ensures that the forcing will ',&
	480	& 'start at the date of the restart file. Obviously the itau_fin ',&
	481	& 'needs to be shifted as well : itau_dep, itau_fin ', &
	482	& itau_dep, itau_fin
	483	ENDIF
	484	!-
	485	! Same things if the starting dates are not the same.
	486	! Everything should look as if we had only one forcing file !
	487	!-
	488	IF (date0 /= date0_rest) THEN
	489	WRITE(numout,*) 'date0_rest , date0 : ',date0_rest , date0
	490	for_offset = NINT((date0_rest-date0)*one_day/dt_force)
	491	ELSE
	492	for_offset = 0
	493	ENDIF
[4693]	494	IF (printlev_loc >= 3) WRITE(numout,*) 'OFFSET FOR THE data read :', for_offset
[8]	495
[999]	496	CALL ioconf_startdate(date0_rest)
[8]	497	!-
[567]	498	!Config Key = TIME_SKIP
	499	!Config Desc = Time in the forcing file at which the model is started.
	500	!Config If = [-]
	501	!Config Def = 0
	502	!Config Help = This time give the point in time at which the model
[8]	503	!Config should be started. If exists, the date of the restart file is use.
[567]	504	!Config The FORMAT of this date can be either of the following :
	505	!Config n : time step n within the forcing file
	506	!Config nS : n seconds after the first time-step in the file
	507	!Config nD : n days after the first time-step
	508	!Config nM : n month after the first time-step (year of 365 days)
	509	!Config nY : n years after the first time-step (year of 365 days)
[8]	510	!Config Or combinations :
[567]	511	!Config nYmM: n years and m month
[844]	512	!Config Units = [seconds, days, months, years]
[8]	513	!-
	514	itau_skip = 0
	515	WRITE(time_str,'(I10)') itau_skip
	516	CALL getin_p('TIME_SKIP', time_str)
	517	!-
	518	! Transform into itau
	519	!-
	520	CALL tlen2itau (time_str, dt_force, date0, itau_skip)
	521	!-
	522	itau_dep = itau_dep+itau_skip
	523	!-
	524	! We need to select the right position of the splited time steps.
	525	!-
	526	istp = itau_dep*split+1
	527	IF (MOD(istp-1,split) /= 0) THEN
	528	split_start = MOD(istp-1,split)+1
	529	ELSE
	530	split_start = 1
	531	ENDIF
	532	istp_old = itau_dep_rest
	533	itau_len = itau_fin-itau_dep
[4263]	534
[567]	535	!Config Key = TIME_LENGTH
	536	!Config Desc = Length of the integration in time.
	537	!Config If = [-]
[2861]	538	!Config Def = Full length of the forcing file
[567]	539	!Config Help = Length of integration. By default the entire length
[8]	540	!Config of the forcing is used. The FORMAT of this date can
	541	!Config be either of the following :
[567]	542	!Config n : time step n within the forcing file
	543	!Config nS : n seconds after the first time-step in the file
	544	!Config nD : n days after the first time-step
	545	!Config nM : n month after the first time-step (year of 365 days)
	546	!Config nY : n years after the first time-step (year of 365 days)
[8]	547	!Config Or combinations :
[567]	548	!Config nYmM: n years and m month
[844]	549	!Config Units = [seconds, days, months, years]
[8]	550	!-
	551	WRITE(time_str,'(I10)') itau_len
	552	CALL getin_p('TIME_LENGTH', time_str)
	553	!-
[7256]	554	! Transform time lenght into number of time step (itau_len) using date of current year
[8]	555	!-
[7256]	556	date_cur = itau2date(itau_dep, date0, dt_force)
	557	CALL tlen2itau (time_str, dt_force, date_cur, itau_len)
[8]	558	!-
	559	itau_fin = itau_dep+itau_len
	560	!-
[4693]	561	IF (printlev_loc >= 1) THEN
	562	WRITE(numout,*) '>> Time origine in the forcing file :', date0
	563	WRITE(numout,*) '>> Time origine in the restart file :', date0_rest
	564	WRITE(numout,*) '>> Simulate starts at forcing time-step : ', itau_dep
	565	WRITE(numout,*) '>> The splited time-steps start at (Sets the '
	566	WRITE(numout,*) '>> chronology for the history and restart files):',istp
	567	WRITE(numout,*) '>> The time spliting starts at :', split_start
	568	WRITE(numout,*) '>> Simulation ends at forcing time-step: ', itau_fin
	569	WRITE(numout,*) '>> Length of the simulation is thus :', itau_len
	570	WRITE(numout,*) '>> Length of the forcing data is in time-steps : ', tm
	571	WRITE(numout,*) '>> Time steps : true, forcing and restart : ', dt,dt_force,dt_rest
	572	END IF
	573
[8]	574	IF (tm < itau_len) THEN
[1078]	575	CALL ipslerr_p ( 2, 'dim2_driver','Length of the simulation is greater than.', &
[4693]	576	' Length of the forcing data is in time-steps','verify TIME_LENGTH parameter.')
[8]	577	ENDIF
	578
[4693]	579
[8]	580	!=====================================================================
	581	!- 2.0 This section is going to define the details by which
	582	!- the input data is going to be used to force the
	583	!- land-surface scheme. The tasks are the following :
	584	!- - Is the coupling going to be explicit or implicit
	585	!- - Type of interpolation to be used.
	586	!- - At which height are the atmospheric forcings going to be used ?
	587	!- - Is a relaxation method going to be used on the forcing
	588	!- - Does net radiation in the interpolated data need to be conserved
	589	!- - How do we distribute precipitation.
	590	!=====================================================================
[567]	591	!Config Key = RELAXATION
	592	!Config Desc = method of forcing
	593	!Config If = [-]
	594	!Config Def = n
	595	!Config Help = A method is proposed by which the first atmospheric
[8]	596	!Config level is not directly forced by observations but
	597	!Config relaxed with a time constant towards observations.
	598	!Config For the moment the methods tends to smooth too much
	599	!Config the diurnal cycle and introduces a time shift.
	600	!Config A more sophisticated method is needed.
[844]	601	!Config Units = [FLAG]
[8]	602	!-
	603	relaxation = .FALSE.
	604	CALL getin_p('RELAXATION', relaxation)
	605	IF ( relaxation ) THEN
	606	WRITE(numout,*) 'dim2_driver : The relaxation option is temporarily disabled as it does not'
	607	WRITE(numout,*) ' produce energy conservation in ORCHIDEE. If you intend to use it'
	608	WRITE(numout,*) ' you should validate it.'
[2373]	609	CALL ipslerr_p(3,'dim2_driver','relaxation option is not activated.','This option needs to be validated.','')
	610
[567]	611	!Config Key = RELAX_A
	612	!Config Desc = Time constant of the relaxation layer
	613	!Config If = RELAXATION
	614	!Config Def = 1.0
	615	!Config Help = The time constant associated to the atmospheric
[8]	616	!Config conditions which are going to be computed
	617	!Config in the relaxed layer. To avoid too much
	618	!Config damping the value should be larger than 1000.
[844]	619	!Config Units = [days?]
[8]	620	!-
	621	alpha = 1000.0
	622	CALL getin_p('RELAX_A', alpha)
	623	ENDIF
[847]	624
[7264]	625	!Config Key = SPREAD_PREC
[567]	626	!Config Desc = Spread the precipitation.
	627	!Config If = [-]
[2639]	628	!Config Def = Half of the forcing time step or uniform, depending on dt_force and dt_sechiba
[7264]	629	!Config Help = Spread the precipitation over SPREAD_PREC steps of the splited forcing
[8]	630	!Config time step. This ONLY applied if the forcing time step has been splited.
	631	!Config If the value indicated is greater than SPLIT_DT, SPLIT_DT is used for it.
[844]	632	!Config Units = [-]
[8]	633	!-
[2639]	634	IF ( dt_force >= 3*one_hour) THEN
	635	! Distribut the precipitations over the half of the forcing time step
	636	nb_spread = INT(0.5 * (dt_force/dt))
	637	ELSE
	638	! Distribut the precipitations uniformly over the forcing time step
	639	nb_spread = dt_force/dt
	640	END IF
	641
[7264]	642	CALL getin_p('SPREAD_PREC', nb_spread)
[8]	643	IF (nb_spread > split) THEN
	644	WRITE(numout,*) 'WARNING : nb_spread is too large it will be '
	645	WRITE(numout,*) ' set to the value of split'
	646	nb_spread = split
	647	ELSE IF (split == 1) THEN
	648	nb_spread = 1
	649	ENDIF
[2639]	650
[4629]	651
[8]	652	!=====================================================================
	653	!- 3.0 Finaly we can prepare all the variables before launching
	654	!- the simulation !
	655	!=====================================================================
	656	! Initialize LOGICAL and the length of the integration
	657	!-
[3850]	658	lstep_init = .TRUE.
	659	lstep_last = .FALSE.
	660
[8]	661	temp_sol_NEW(:,:) = tp_00
	662	!-
[567]	663	!Config Key = ATM_CO2
[4981]	664	!Config Desc = Value to precribe atmosoheric CO2
	665	!Config If = [FORCE_CO2_VEG=y or Offline mode]
[567]	666	!Config Def = 350.
[4981]	667	!Config Help = Used in offline mode or in coupled mode if FORCE_CO2_VEG=y
[844]	668	!Config Units = [ppm]
[8]	669	atmco2=350.
	670	CALL getin_p('ATM_CO2',atmco2)
	671	for_ccanopy(:,:)=atmco2
	672	!-
	673	! Preparing for the implicit scheme.
	674	! This means loading the prognostic variables from the restart file.
	675	!-
[353]	676	Flag=.FALSE.
[8]	677	IF (is_root_prc) THEN
[353]	678	ALLOCATE(fluxsens_g(iim_g,jjm_g))
[8]	679	var_name= 'fluxsens'
	680	CALL restget &
	681	& (rest_id, var_name, iim_g, jjm_g, 1, istp_old, .TRUE., fluxsens_g)
	682	IF (ALL(fluxsens_g(:,:) == val_exp)) THEN
[353]	683	Flag=.TRUE.
	684	ELSE
	685	Flag=.FALSE.
[8]	686	ENDIF
[353]	687	ELSE
	688	ALLOCATE(fluxsens_g(0,1))
[8]	689	ENDIF
[353]	690	CALL bcast(Flag)
	691	IF (.NOT. Flag) THEN
[1078]	692	CALL scatter2D_mpi(fluxsens_g,fluxsens)
[353]	693	ELSE
	694	fluxsens(:,:) = zero
	695	ENDIF
	696	DEALLOCATE(fluxsens_g)
[8]	697	!-
	698	IF (is_root_prc) THEN
[353]	699	ALLOCATE(vevapp_g(iim_g,jjm_g))
[8]	700	var_name= 'vevapp'
	701	CALL restget &
	702	& (rest_id, var_name, iim_g, jjm_g, 1, istp_old, .TRUE., vevapp_g)
	703	IF (ALL(vevapp_g(:,:) == val_exp)) THEN
[353]	704	Flag=.TRUE.
	705	ELSE
	706	Flag=.FALSE.
[8]	707	ENDIF
[353]	708	ELSE
	709	ALLOCATE(vevapp_g(0,1))
[8]	710	ENDIF
[353]	711	CALL bcast(Flag)
	712	IF (.NOT. Flag) THEN
[1078]	713	CALL scatter2D_mpi(vevapp_g,vevapp)
[353]	714	ELSE
	715	vevapp(:,:) = zero
	716	ENDIF
	717	DEALLOCATE(vevapp_g)
[8]	718	!-
	719	IF (is_root_prc) THEN
[353]	720	ALLOCATE(old_zlev_g(iim_g,jjm_g))
[8]	721	var_name= 'zlev_old'
	722	CALL restget &
	723	& (rest_id, var_name, iim_g, jjm_g, 1, istp_old, .TRUE., old_zlev_g)
	724	IF (ALL(old_zlev_g(:,:) == val_exp)) THEN
	725	Flag=.TRUE.
	726	ELSE
	727	Flag=.FALSE.
	728	ENDIF
[353]	729	ELSE
	730	ALLOCATE(old_zlev_g(0,1))
[8]	731	ENDIF
	732	CALL bcast(Flag)
[353]	733	IF ( .NOT. Flag ) THEN
[1078]	734	CALL scatter2D_mpi(old_zlev_g,old_zlev)
[353]	735	ELSE
	736	old_zlev(:,:)=zlev_vec(:,:)
	737	ENDIF
	738	DEALLOCATE(old_zlev_g)
[8]	739	!-
	740	IF (is_root_prc) THEN
[353]	741	ALLOCATE(old_qair_g(iim_g,jjm_g))
[8]	742	var_name= 'qair_old'
	743	CALL restget &
	744	& (rest_id, var_name, iim_g, jjm_g, 1, istp_old, .TRUE., old_qair_g)
	745	IF (ALL(old_qair_g(:,:) == val_exp)) THEN
	746	Flag=.TRUE.
	747	ELSE
	748	Flag=.FALSE.
	749	ENDIF
[353]	750	ELSE
	751	ALLOCATE(old_qair_g(0,1))
[8]	752	ENDIF
	753	CALL bcast(Flag)
[353]	754	IF ( .NOT. Flag ) THEN
[1078]	755	CALL scatter2D_mpi(old_qair_g,old_qair)
[353]	756	ELSE
	757	old_qair(:,:) = qair_obs(:,:)
	758	ENDIF
	759	DEALLOCATE(old_qair_g)
[8]	760	!-
	761	IF (is_root_prc) THEN
[353]	762	ALLOCATE(old_eair_g(iim_g,jjm_g))
[8]	763	var_name= 'eair_old'
	764	CALL restget &
	765	& (rest_id, var_name, iim_g, jjm_g, 1, istp_old, .TRUE., old_eair_g)
	766	IF (ALL(old_eair_g(:,:) == val_exp)) THEN
	767	Flag=.TRUE.
	768	ELSE
	769	Flag=.FALSE.
	770	ENDIF
[353]	771	ELSE
	772	ALLOCATE(old_eair_g(0,1))
[8]	773	ENDIF
	774	CALL bcast(Flag)
[353]	775	IF ( .NOT. Flag ) THEN
[1078]	776	CALL scatter2D_mpi(old_eair_g,old_eair)
[353]	777	ELSE
[8]	778	DO ik=1,nbindex
	779	i=ilandindex(ik)
	780	j=jlandindex(ik)
	781	old_eair(i,j) = cp_air * tair_obs(i,j) + cte_grav*zlev_vec(i,j)
	782	ENDDO
	783	ENDIF
[353]	784	DEALLOCATE(old_eair_g)
[8]	785	!-
	786	! old density is also needed because we do not yet have the right pb
	787	!-
	788	!=> obsolète ??!! (tjrs calculé après forcing_read)
	789	IF (is_root_prc) THEN
[353]	790	ALLOCATE(for_rau_g(iim_g,jjm_g))
[8]	791	var_name= 'rau_old'
	792	CALL restget &
	793	& (rest_id, var_name, iim_g, jjm_g, 1, istp_old, .TRUE., for_rau_g)
	794	IF (ALL(for_rau_g(:,:) == val_exp)) THEN
	795	Flag=.TRUE.
	796	ELSE
	797	Flag=.FALSE.
	798	ENDIF
[353]	799	ELSE
	800	ALLOCATE(for_rau_g(0,1))
[8]	801	ENDIF
	802	CALL bcast(Flag)
[353]	803	IF ( .NOT. Flag ) THEN
[1078]	804	CALL scatter2D_mpi(for_rau_g,for_rau)
[353]	805	ELSE
[8]	806	DO ik=1,nbindex
	807	i=ilandindex(ik)
	808	j=jlandindex(ik)
	809	for_rau(i,j) = pb(i,j)/(cte_molr*(tair_obs(i,j)))
	810	ENDDO
	811	ENDIF
[353]	812	DEALLOCATE(for_rau_g)
[8]	813	!-
	814	! For this variable the restart is extracted by SECHIBA
	815	!-
	816	temp_sol_NEW(:,:) = tair_obs(:,:)
	817	!-
	818	if (.NOT. is_watchout) THEN
	819	!-
	820	! This does not yield a correct restart in the case of relaxation
	821	!-
	822	IF (is_root_prc) THEN
[353]	823	ALLOCATE(petAcoef_g(iim_g,jjm_g))
[8]	824	var_name= 'petAcoef'
	825	CALL restget &
	826	& (rest_id, var_name, iim_g, jjm_g, 1, istp_old, .TRUE., petAcoef_g)
	827	IF (ALL(petAcoef_g(:,:) == val_exp)) THEN
	828	Flag=.TRUE.
	829	ELSE
	830	Flag=.FALSE.
	831	ENDIF
[353]	832	ELSE
	833	ALLOCATE(petAcoef_g(0,1))
[8]	834	ENDIF
	835	CALL bcast(Flag)
[353]	836	IF ( .NOT. Flag ) THEN
[1078]	837	CALL scatter2D_mpi(petAcoef_g,petAcoef)
[353]	838	ELSE
	839	petAcoef(:,:) = zero
	840	ENDIF
	841	DEALLOCATE(petAcoef_g)
[8]	842	!--
	843	IF (is_root_prc) THEN
[353]	844	ALLOCATE(petBcoef_g(iim_g,jjm_g))
[8]	845	var_name= 'petBcoef'
	846	CALL restget &
	847	& (rest_id, var_name, iim_g, jjm_g, 1, istp_old, .TRUE., petBcoef_g)
	848	IF (ALL(petBcoef_g(:,:) == val_exp)) THEN
	849	Flag=.TRUE.
	850	ELSE
	851	Flag=.FALSE.
	852	ENDIF
[353]	853	ELSE
	854	ALLOCATE(petBcoef_g(0,1))
[8]	855	ENDIF
	856	CALL bcast(Flag)
[353]	857	IF ( .NOT. Flag ) THEN
[1078]	858	CALL scatter2D_mpi(petBcoef_g,petBcoef)
[353]	859	ELSE
	860	petBcoef(:,:) = old_eair(:,:)
	861	ENDIF
	862	DEALLOCATE(petBcoef_g)
[8]	863	!--
	864	IF (is_root_prc) THEN
[353]	865	ALLOCATE(peqAcoef_g(iim_g,jjm_g))
[8]	866	var_name= 'peqAcoef'
	867	CALL restget &
	868	& (rest_id, var_name, iim_g, jjm_g, 1, istp_old, .TRUE., peqAcoef_g)
	869	IF (ALL(peqAcoef_g(:,:) == val_exp)) THEN
	870	Flag=.TRUE.
	871	ELSE
	872	Flag=.FALSE.
	873	ENDIF
[353]	874	ELSE
	875	ALLOCATE(peqAcoef_g(0,1))
[8]	876	ENDIF
	877	CALL bcast(Flag)
[353]	878	IF ( .NOT. Flag ) THEN
[1078]	879	CALL scatter2D_mpi(peqAcoef_g,peqAcoef)
[353]	880	ELSE
	881	peqAcoef(:,:) = zero
	882	ENDIF
	883	DEALLOCATE(peqAcoef_g)
[8]	884	!--
	885	IF (is_root_prc) THEN
[353]	886	ALLOCATE(peqBcoef_g(iim_g,jjm_g))
[8]	887	var_name= 'peqBcoef'
	888	CALL restget &
	889	& (rest_id, var_name, iim_g, jjm_g, 1, istp_old, .TRUE., peqBcoef_g)
	890	IF (ALL(peqBcoef_g(:,:) == val_exp)) THEN
	891	Flag=.TRUE.
	892	ELSE
	893	Flag=.FALSE.
	894	ENDIF
[353]	895	ELSE
	896	ALLOCATE(peqBcoef_g(0,1))
[8]	897	ENDIF
	898	CALL bcast(Flag)
[353]	899	IF ( .NOT. Flag ) THEN
[1078]	900	CALL scatter2D_mpi(peqBcoef_g,peqBcoef)
[353]	901	ELSE
	902	peqBcoef(:,:) = old_qair(:,:)
	903	ENDIF
	904	DEALLOCATE(peqBcoef_g)
[8]	905	ENDIF
	906	!-
	907	! And other variables which need initial variables. These variables
	908	! will get properly initialized by ORCHIDEE when it is called for
	909	! the first time.
	910	!-
	911	albedo(:,:,:) = 0.13
	912	emis(:,:) = 1.0
	913	z0(:,:) = 0.1
	914	!--
	915	!=====================================================================
	916	!- 4.0 START THE TIME LOOP
	917	!=====================================================================
	918
	919	it = itau_dep+1
	920	DO WHILE ( it <= itau_fin )
	921	!----
	922	it_force = it+for_offset
	923	IF (it_force < 0) THEN
	924	WRITE(numout,*) 'The day is not in the forcing file :', &
	925	& it_force, it, for_offset
[2373]	926	CALL ipslerr_p(3,'dim2_driver','Pb in forcing file','','')
[8]	927	ENDIF
	928	is=split_start
	929	DO WHILE ( is <= split )
	930	!-----
[999]	931	julian = itau2date(istp, date0_rest, dt)
[8]	932	CALL ju2ymds(julian, yy, mm, dd, ss)
[2348]	933	IF (printlev_loc>=3) THEN
[8]	934	WRITE(numout,*) "=============================================================="
	935	WRITE(numout,"('New iteration at date : ',I4,'-',I2.2,'-',I2.2,':',F8.4)") &
	936	& yy,mm,dd,ss/3600.
	937	#ifdef CPP_PARA
	938	IF (is_root_prc) THEN
	939	WRITE(,) "=============================================================="
	940	WRITE(*,"('New iteration at date : ',I4,'-',I2.2,'-',I2.2,':',F8.4)") &
	941	& yy,mm,dd,ss/3600.
	942	ENDIF
	943	#endif
	944	ENDIF
	945	!-----
	946	IF ( (it == itau_fin).AND.(is == split) ) THEN
[3850]	947	lstep_last = .TRUE.
[8]	948	ENDIF
	949	!-----
[2348]	950	IF (printlev_loc>=3) WRITE(numout,*) 'Into forcing_read'
[8]	951	!-----
	952	CALL forcing_READ &
[3850]	953	& (filename, rest_id, .FALSE., lstep_last, &
[4629]	954	& it_force, istp, is, split, nb_spread, lwdown_cons, swdown_cons, &
[999]	955	& date0_rest, dt_force, iim, jjm, lon, lat, zlev_vec, zlevuv_vec, tm, &
[2961]	956	& swdown, coszang, precip_rain, precip_snow, tair_obs, &
[8]	957	& u, v, qair_obs, pb, for_lwdown, for_contfrac, for_neighbours, for_resolution, &
	958	& for_swnet, eair_obs, petAcoef, peqAcoef, petBcoef, peqBcoef, cdrag, for_ccanopy, &
	959	& kindex, nbindex, force_id)
	960
	961	!-----
	962	!---- SECHIBA expects surface pressure in hPa
	963	!-----
	964	for_psurf(:,:) = pb(:,:)/100.
	965
[2348]	966	IF (printlev_loc>=4) THEN
[8]	967	WRITE(numout,*) "dim2_driver 0 ",it_force
[353]	968	WRITE(numout,*) ">> Index of land points =",kindex(1:nbindex)
[8]	969	WRITE(numout,*) "Lowest level wind speed North = ", &
	970	& (/ ( u(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
	971	WRITE(numout,*) "Lowest level wind speed East = ", &
	972	& (/ ( v(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
	973	WRITE(numout,*) "z0 ; Surface roughness = ", &
	974	& (/ ( z0(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
	975	WRITE(numout,*) "Height of first layer = ", &
	976	& (/ ( zlev_vec(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
	977	WRITE(numout,*) "Lowest level specific humidity = ", &
	978	& (/ ( qair_obs(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
	979	WRITE(numout,*) "Rain precipitation = ", &
	980	& (/ ( precip_rain(ilandindex(ik), jlandindex(ik))*dt,ik=1,nbindex ) /)
	981	WRITE(numout,*) "Snow precipitation = ", &
	982	& (/ ( precip_snow(ilandindex(ik), jlandindex(ik))*dt,ik=1,nbindex ) /)
	983	WRITE(numout,*) "Down-welling long-wave flux = ", &
	984	& (/ ( for_lwdown(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
	985	WRITE(numout,*) "Net surface short-wave flux = ", &
	986	& (/ ( for_swnet(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
	987	WRITE(numout,*) "Downwelling surface short-wave flux = ", &
	988	& (/ ( swdown(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
	989	WRITE(numout,*) "Air temperature in Kelvin = ", &
	990	& (/ ( tair_obs(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
	991	WRITE(numout,*) "Air potential energy = ", &
	992	& (/ ( eair_obs(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
	993	WRITE(numout,*) "CO2 concentration in the canopy = ", &
	994	& (/ ( for_ccanopy(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
	995	WRITE(numout,*) "Coeficients A from the PBL resolution = ", &
	996	& (/ ( petAcoef(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
	997	WRITE(numout,*) "One for T and another for q = ", &
	998	& (/ ( peqAcoef(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
	999	WRITE(numout,*) "Coeficients B from the PBL resolution = ", &
	1000	& (/ ( petBcoef(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
	1001	WRITE(numout,*) "One for T and another for q = ", &
	1002	& (/ ( peqBcoef(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
	1003	WRITE(numout,*) "Cdrag = ", &
	1004	& (/ ( cdrag(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
	1005	WRITE(numout,*) "CO2 concentration in the canopy = ", &
	1006	& (/ ( for_ccanopy(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
	1007	WRITE(numout,*) "Lowest level pressure = ", &
	1008	& (/ ( for_psurf(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
	1009	WRITE(numout,*) "Geographical coordinates lon = ", &
	1010	& (/ ( lon(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
	1011	WRITE(numout,*) "Geographical coordinates lat = ", &
	1012	& (/ ( lat(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
	1013	WRITE(numout,*) "Fraction of continent in the grid = ", &
	1014	& (/ ( for_contfrac(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
	1015	ENDIF
	1016	!-----
	1017	!---- Prepare : tmp_qair, tmp_eair, tmp_tair, tmp_pb
	1018	!---- and : for_u, for_v, for_lwdown, for_swnet, for_swdown
	1019	!---- All the work is done in forcing_read
[2348]	1020	IF (printlev_loc>=3) WRITE(numout,*) 'Prepare the atmospheric forcing'
[8]	1021	!-----
	1022	IF (.NOT. is_watchout) THEN
	1023	DO ik=1,nbindex
	1024	i=ilandindex(ik)
	1025	j=jlandindex(ik)
	1026	eair_obs(i,j) = cp_airtair_obs(i,j)+cte_gravzlev_vec(i,j)
	1027	for_swnet(i,j) = (1.-(albedo(i,j,1)+albedo(i,j,2))/2.)*swdown(i,j)
	1028	ENDDO
	1029	ENDIF
[48]	1030	DO ik=1,nbindex
	1031	i=ilandindex(ik)
	1032	j=jlandindex(ik)
	1033	for_swdown(i,j) = swdown(i,j)
[2961]	1034	for_coszang(i,j) = coszang(i,j)
[48]	1035	ENDDO
[8]	1036	!-----
	1037	!---- Computing the buffer zone !
	1038	!-----
	1039	IF (relaxation) THEN
	1040	DO ik=1,nbindex
	1041	i=ilandindex(ik)
	1042	j=jlandindex(ik)
	1043	for_qair(i,j) = peqAcoef(i,j)(-1.) vevapp(i,j)*dt+peqBcoef(i,j)
	1044	!-------
	1045	for_eair(i,j) = petAcoef(i,j)(-1.) fluxsens(i,j)+petBcoef(i,j)
	1046	!-------
	1047	ENDDO
	1048	DO ik=1,nbindex
	1049	i=ilandindex(ik)
	1050	j=jlandindex(ik)
	1051	for_tair(i,j) = (for_eair(i,j) - cte_grav*zlev_vec(i,j))/cp_air
	1052	!-------
	1053	!!$ if (.NOT. is_watchout) &
	1054	!!$ epot_sol(:,:) = cp_air*temp_sol_NEW(:,:)
	1055	!-------
	1056	ENDDO
	1057	DO ik=1,nbindex
	1058	i=ilandindex(ik)
	1059	j=jlandindex(ik)
	1060	for_rau(i,j) = pb(i,j) / (cte_molr*for_tair(i,j))
	1061	!-------
	1062	relax(i,j) = for_rau(i,j)*alpha
	1063	ENDDO
[1224]	1064
[8]	1065	DO ik=1,nbindex
	1066	i=ilandindex(ik)
	1067	j=jlandindex(ik)
[1224]	1068	zlflu = zlev_vec(i,j)/2.0*dt
[8]	1069	peqAcoef(i,j) = 1.0/(zlflu+relax(i,j))
	1070	peqBcoef(i,j) = (relax(i,j) * qair_obs(i,j)/(zlflu+relax(i,j))) + &
	1071	& for_qair(i,j)/(1.0+relax(i,j)/zlflu)
	1072	ENDDO
	1073	!-------
	1074	! relax(:,:) = for_rau(:,:)*alpha
	1075	DO ik=1,nbindex
	1076	i=ilandindex(ik)
	1077	j=jlandindex(ik)
	1078	petAcoef(i,j) = 1.0/(zlflu+relax(i,j))
	1079	petBcoef(i,j) = ( relax(i,j) * eair_obs(i,j) / (zlflu+relax(i,j)) ) &
	1080	& + for_eair(i,j)/(1.0+relax(i,j)/zlflu)
	1081	ENDDO
	1082	ELSE
	1083	for_qair(:,:) = fill_init
	1084	for_eair(:,:) = fill_init
	1085	for_tair(:,:) = fill_init
	1086	DO ik=1,nbindex
	1087	i=ilandindex(ik)
	1088	j=jlandindex(ik)
	1089	for_qair(i,j) = qair_obs(i,j)
	1090	for_eair(i,j) = eair_obs(i,j)
	1091	for_tair(i,j) = tair_obs(i,j)
	1092	ENDDO
	1093	!-------
	1094	!!$ if (.NOT. is_watchout) &
	1095	!!$ epot_sol(:,:) = cp_air*temp_sol_NEW(:,:)
	1096	!-------
	1097	DO ik=1,nbindex
	1098	i=ilandindex(ik)
	1099	j=jlandindex(ik)
	1100	for_rau(i,j) = pb(i,j) / (cte_molr*for_tair(i,j))
	1101	ENDDO
	1102	!-------
	1103	IF (.NOT. is_watchout) THEN
	1104	petAcoef(:,:) = 0.0
	1105	peqAcoef(:,:) = 0.0
	1106	DO ik=1,nbindex
	1107	i=ilandindex(ik)
	1108	j=jlandindex(ik)
	1109	petBcoef(i,j) = eair_obs(i,j)
	1110	peqBcoef(i,j) = qair_obs(i,j)
	1111	ENDDO
	1112	ENDIF
	1113	ENDIF
	1114	!-----
	1115	IF (.NOT. is_watchout) &
	1116	cdrag(:,:) = 0.0
	1117	for_ccanopy(:,:)=atmco2
	1118	!-----
	1119	!---- SECHIBA expects wind, temperature and humidity at the same height.
	1120	!---- If this is not the case then we need to correct for that.
	1121	!-----
[1224]	1122	DO ik=1,nbindex
	1123	i=ilandindex(ik)
	1124	j=jlandindex(ik)
	1125	for_u(i,j) = u(i,j)*LOG(zlev_vec(i,j)/z0(i,j)) / &
	1126	LOG(zlevuv_vec(i,j)/z0(i,j))
	1127	for_v(i,j) = v(i,j)*LOG(zlev_vec(i,j)/z0(i,j)) / &
	1128	LOG(zlevuv_vec(i,j)/z0(i,j))
	1129	END DO
	1130
[8]	1131	!-----
	1132	!---- Prepare the other variables WITH the special CASE
	1133	!---- of splited time steps
	1134	!----
[2348]	1135	!---- PRINT input value for printlev_loc>=3
[8]	1136	!-----
[2348]	1137	IF (printlev_loc>=3) THEN
[8]	1138	WRITE(numout,*) ' >>>>>> time step it_force = ',it_force
	1139	WRITE(numout,*) &
	1140	& ' tair, qair, eair = ', &
	1141	& for_tair(itest,jtest),for_qair(itest,jtest), &
	1142	& for_eair(itest,jtest)
	1143	WRITE(numout,*) &
	1144	& ' OBS : tair, qair, eair = ', &
	1145	& tair_obs(itest,jtest),qair_obs(itest,jtest), &
	1146	& eair_obs(itest,jtest)
	1147	WRITE(numout,*) ' u et v = ',for_u(itest,jtest),for_v(itest,jtest)
	1148	WRITE(numout,*) ' precip rain et snow = ', &
	1149	& precip_rain(itest,jtest),precip_snow(itest,jtest)
	1150	WRITE(numout,*) ' lwdown et swnet = ', &
	1151	& for_lwdown(itest,jtest),for_swnet(itest,jtest)
	1152	WRITE(numout,*) ' petAcoef et peqAcoef = ', &
	1153	& petAcoef(itest,jtest), peqAcoef(itest,jtest)
	1154	WRITE(numout,*) ' petBcoef et peqAcoef = ', &
	1155	& petBcoef(itest,jtest),peqBcoef(itest,jtest)
	1156	WRITE(numout,*) ' zlev = ',zlev_vec(itest,jtest)
	1157	ENDIF
	1158	!-----
[4263]	1159
[3850]	1160	IF (lstep_init) THEN
[8]	1161
	1162	DO ik=1,nbindex
	1163	i=ilandindex(ik)
	1164	j=jlandindex(ik)
	1165	for_swdown(i,j) = swdown(i,j)
[2961]	1166	for_coszang(i,j) = coszang(i,j)
[8]	1167	ENDDO
[2348]	1168	IF (printlev_loc>=4) THEN
[3850]	1169	WRITE(numout,*) "dim2_driver lstep_init ",it_force
[353]	1170	WRITE(numout,*) ">> Index of land points =",kindex(1:nbindex)
[8]	1171	WRITE(numout,*) "Lowest level wind speed North = ", &
	1172	& (/ ( for_u(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
	1173	WRITE(numout,*) "Lowest level wind speed East = ", &
	1174	& (/ ( for_v(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
	1175	WRITE(numout,*) "z0 ; Surface roughness = ", &
	1176	& (/ ( z0(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
	1177	WRITE(numout,*) "Height of first layer = ", &
	1178	& (/ ( zlev_vec(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
	1179	WRITE(numout,*) "Lowest level specific humidity = ", &
	1180	& (/ ( for_qair(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
	1181	WRITE(numout,*) "Rain precipitation = ", &
	1182	& (/ ( precip_rain(ilandindex(ik), jlandindex(ik))*dt,ik=1,nbindex ) /)
	1183	WRITE(numout,*) "Snow precipitation = ", &
	1184	& (/ ( precip_snow(ilandindex(ik), jlandindex(ik))*dt,ik=1,nbindex ) /)
	1185	WRITE(numout,*) "Down-welling long-wave flux = ", &
	1186	& (/ ( for_lwdown(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
	1187	WRITE(numout,*) "Net surface short-wave flux = ", &
	1188	& (/ ( for_swnet(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
	1189	WRITE(numout,*) "Downwelling surface short-wave flux = ", &
	1190	& (/ ( for_swdown(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
	1191	WRITE(numout,*) "Air temperature in Kelvin = ", &
	1192	& (/ ( for_tair(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
	1193	WRITE(numout,*) "Air potential energy = ", &
	1194	& (/ ( for_eair(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
	1195	WRITE(numout,*) "CO2 concentration in the canopy = ", &
	1196	& (/ ( for_ccanopy(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
	1197	WRITE(numout,*) "Coeficients A from the PBL resolution = ", &
	1198	& (/ ( petAcoef(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
	1199	WRITE(numout,*) "One for T and another for q = ", &
	1200	& (/ ( peqAcoef(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
	1201	WRITE(numout,*) "Coeficients B from the PBL resolution = ", &
	1202	& (/ ( petBcoef(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
	1203	WRITE(numout,*) "One for T and another for q = ", &
	1204	& (/ ( peqBcoef(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
	1205	WRITE(numout,*) "Cdrag = ", &
	1206	& (/ ( cdrag(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
	1207	WRITE(numout,*) "Lowest level pressure = ", &
	1208	& (/ ( for_psurf(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
	1209	WRITE(numout,*) "Geographical coordinates lon = ", &
	1210	& (/ ( lon(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
	1211	WRITE(numout,*) "Geographical coordinates lat = ", &
	1212	& (/ ( lat(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
	1213	WRITE(numout,*) "Fraction of continent in the grid = ", &
	1214	& (/ ( for_contfrac(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
	1215	ENDIF
	1216	!-------
	1217	!------ CALL sechiba to initialize fields
	1218	!------ and have some initial results: emis, albedo, z0
	1219	!-------
[2605]	1220	CALL intersurf_initialize_2d &
[8]	1221	& (istp_old, iim, jjm, nbindex, kindex, dt, &
[3850]	1222	& lstep_init, .FALSE., lon, lat, for_contfrac, for_resolution, date0_rest, &
[8]	1223	! first level conditions
	1224	& zlev_vec, for_u, for_v, &
	1225	& for_qair, for_tair, for_eair, for_ccanopy, &
	1226	! Variables for the implicit coupling
	1227	& cdrag, petAcoef, peqAcoef, petBcoef, peqBcoef, &
	1228	! Rain, snow, radiation and surface pressure
	1229	& precip_rain, precip_snow, &
[898]	1230	& for_lwdown, for_swnet, for_swdown, for_psurf, &
[8]	1231	! Output : Fluxes
	1232	& vevapp, fluxsens, fluxlat, coastalflow, riverflow, &
	1233	! Surface temperatures and surface properties
[2605]	1234	& tsol_rad, temp_sol_NEW, qsurf, albedo, emis, z0 )
[8]	1235
[1078]	1236	CALL Stop_timer(timer_global)
	1237	CALL Stop_timer(timer_mpi)
	1238	CALL Start_timer(timer_global)
	1239	CALL Start_timer(timer_mpi)
[8]	1240	!
[3850]	1241	lstep_init = .FALSE.
[8]	1242	!
	1243	! Get Restart values for albedo and z0,
	1244	! as they modify forcing variables swnet and wind.
	1245	!-------
	1246	! albedo
[353]	1247	IF (is_root_prc) THEN
	1248	ALLOCATE(albedo_g(iim_g,jjm_g))
	1249	ELSE
	1250	ALLOCATE(albedo_g(0,1))
	1251	ENDIF
[8]	1252	!
	1253	IF (is_root_prc) THEN
	1254	var_name= 'albedo_vis'
	1255	CALL restget &
	1256	& (rest_id, var_name, iim_g, jjm_g, 1, istp_old, .TRUE., albedo_g)
	1257	IF (ALL(albedo_g(:,:) == val_exp)) THEN
	1258	Flag=.TRUE.
	1259	ELSE
	1260	Flag=.FALSE.
	1261	ENDIF
	1262	ENDIF
	1263	CALL bcast(Flag)
[353]	1264	IF ( .NOT. Flag ) THEN
[1078]	1265	CALL scatter2D_mpi(albedo_g,albedo_vis)
[353]	1266	albedo(:,:,1)=albedo_vis(:,:)
	1267	ELSE
	1268	albedo_vis(:,:)=albedo(:,:,1)
	1269	ENDIF
[8]	1270	!
	1271	IF (is_root_prc) THEN
	1272	var_name= 'albedo_nir'
	1273	CALL restget &
	1274	& (rest_id, var_name, iim_g, jjm_g, 1, istp_old, .TRUE., albedo_g)
	1275	IF (ALL(albedo_g(:,:) == val_exp)) THEN
	1276	Flag=.TRUE.
	1277	ELSE
	1278	Flag=.FALSE.
	1279	ENDIF
	1280	ENDIF
	1281	CALL bcast(Flag)
[353]	1282	IF ( .NOT. Flag ) THEN
[1078]	1283	CALL scatter2D_mpi(albedo_g,albedo_nir)
[353]	1284	albedo(:,:,2)=albedo_nir(:,:)
	1285	ELSE
	1286	albedo_nir(:,:)=albedo(:,:,2)
	1287	ENDIF
[8]	1288	!
[353]	1289	DEALLOCATE(albedo_g)
[8]	1290	!--
	1291	! z0
	1292	IF (is_root_prc) THEN
[353]	1293	ALLOCATE(z0_g(iim_g,jjm_g))
[8]	1294	var_name= 'z0'
	1295	CALL restget &
	1296	& (rest_id, var_name, iim_g, jjm_g, 1, istp_old, .TRUE., z0_g)
	1297	IF (ALL(z0_g(:,:) == val_exp)) THEN
	1298	Flag=.TRUE.
	1299	ELSE
	1300	Flag=.FALSE.
	1301	ENDIF
[353]	1302	ELSE
	1303	ALLOCATE(z0_g(0,1))
[8]	1304	ENDIF
	1305	CALL bcast(Flag)
[353]	1306	IF (.NOT. Flag) &
[1078]	1307	CALL scatter2D_mpi(z0_g,z0)
[353]	1308	DEALLOCATE(z0_g)
[8]	1309	!-------
	1310	DO ik=1,nbindex
	1311	i=ilandindex(ik)
	1312	j=jlandindex(ik)
	1313	temp_sol_old(i,j) = temp_sol_NEW(i,j)
	1314	for_swnet(i,j) = (1.- (albedo(i,j,1)+albedo(i,j,2))/2.)*swdown(i,j)
	1315	for_swdown(i,j) = swdown(i,j)
[2961]	1316	for_coszang(i,j) = coszang(i,j)
[8]	1317	ENDDO
	1318	!
	1319	! MM : z0 have been modified then we must lower the wind again
	1320	!-----
	1321	!---- SECHIBA expects wind, temperature and humidity at the same height.
	1322	!---- If this is not the case then we need to correct for that.
	1323	!-----
[1224]	1324	DO ik=1,nbindex
	1325	i=ilandindex(ik)
	1326	j=jlandindex(ik)
	1327	for_u(i,j) = u(i,j) * LOG(zlev_vec(i,j)/z0(i,j)) / &
	1328	LOG(zlevuv_vec(i,j)/z0(i,j))
	1329	for_v(i,j) = v(i,j) * LOG(zlev_vec(i,j)/z0(i,j)) / &
	1330	LOG(zlevuv_vec(i,j)/z0(i,j))
	1331	END DO
	1332
[8]	1333	!-----
[3850]	1334	!---- PRINT input value after lstep_init for printlev_loc>=3
[8]	1335	!-----
[2348]	1336	IF (printlev_loc>=3) THEN
[3850]	1337	WRITE(numout,*) ' >>>>>> after lstep_init = ',lstep_init
[8]	1338	WRITE(numout,*) ' u et v = ',for_u(itest,jtest),for_v(itest,jtest)
	1339	WRITE(numout,*) ' swnet = ', for_swnet(itest,jtest)
	1340	ENDIF
	1341	!-------
[2348]	1342	IF (printlev_loc>=4) THEN
[3850]	1343	WRITE(numout,*) "dim2_driver lstep_init outputs"
[8]	1344	! Output : Fluxes
	1345	WRITE(numout,*) "vevapp ; Total of evaporation = ", &
	1346	& (/ ( vevapp(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
	1347	WRITE(numout,*) "Sensible heat flux = ", &
	1348	& (/ ( fluxsens(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
	1349	WRITE(numout,*) "Latent heat flux = ", &
	1350	& (/ ( fluxlat(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
	1351	WRITE(numout,*) "coastalflow ; Diffuse flow of water into the ocean (m^3/dt) = ", &
	1352	& (/ ( coastalflow(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
	1353	WRITE(numout,*) "riverflow ; Largest rivers flowing into the ocean (m^3/dt) = ", &
	1354	& (/ ( riverflow(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
	1355	! Surface temperatures and surface properties
	1356	WRITE(numout,*) "tsol_rad ; Radiative surface temperature = ", &
	1357	& (/ ( tsol_rad(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
	1358	WRITE(numout,*) "temp_sol_new ; New soil temperature = ", &
	1359	& (/ ( temp_sol_NEW(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
	1360	WRITE(numout,*) "qsurf ; Surface specific humidity = ", &
	1361	& (/ ( qsurf(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
	1362	WRITE(numout,*) "albedoVIS = ", &
	1363	& (/ ( albedo(ilandindex(ik), jlandindex(ik), 1),ik=1,nbindex ) /)
	1364	WRITE(numout,*) "albedoNIR = ", &
	1365	& (/ ( albedo(ilandindex(ik), jlandindex(ik), 2),ik=1,nbindex ) /)
	1366	WRITE(numout,*) "emis ; Emissivity = ", &
	1367	& (/ ( emis(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
	1368	WRITE(numout,*) "z0 ; Surface roughness = ", &
	1369	& (/ ( z0(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
	1370	ENDIF
	1371	!-------
[2348]	1372	IF (printlev_loc>=3) THEN
[8]	1373	WRITE(numout,*) &
	1374	& ' OUT rest : z0, albedoVIS, albedoNIR, emis = ', &
	1375	& z0(itest,jtest),albedo(itest,jtest,1), &
	1376	& albedo(itest,jtest,2),emis(itest,jtest)
	1377	WRITE(numout,*) ' OUT rest : coastal and river flow = ', &
	1378	& coastalflow(itest,jtest), riverflow(itest,jtest)
	1379	WRITE(numout,*) ' OUT rest : tsol_rad, vevapp = ', &
	1380	& tsol_rad(itest,jtest), vevapp(itest,jtest)
	1381	WRITE(numout,*) ' OUT rest : temp_sol_new =', &
	1382	& temp_sol_NEW(itest,jtest)
	1383	ENDIF
	1384
[4263]	1385	ENDIF ! lstep_init
[8]	1386	!-----
	1387	!---- Calling SECHIBA and doing the number crunching.
	1388	!---- Note that for the first time step SECHIBA is called twice.
	1389	!----
	1390	!---- All H_2O fluxes are now in Kg/m^2s
	1391	!-----
[2348]	1392	IF (printlev_loc>=4) THEN
[8]	1393	WRITE(numout,*) "dim2_driver ",it_force
[353]	1394	WRITE(numout,*) ">> Index of land points =",kindex(1:nbindex)
[8]	1395	WRITE(numout,*) "Lowest level wind speed North = ", &
	1396	& (/ ( for_u(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
	1397	WRITE(numout,*) "Lowest level wind speed East = ", &
	1398	& (/ ( for_v(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
	1399	WRITE(numout,*) "z0 ; Surface roughness = ", &
	1400	& (/ ( z0(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
	1401	WRITE(numout,*) "Height of first layer = ", &
	1402	& (/ ( zlev_vec(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
	1403	WRITE(numout,*) "Lowest level specific humidity = ", &
	1404	& (/ ( for_qair(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
	1405	WRITE(numout,*) "Rain precipitation = ", &
	1406	& (/ ( precip_rain(ilandindex(ik), jlandindex(ik))*dt,ik=1,nbindex ) /)
	1407	WRITE(numout,*) "Snow precipitation = ", &
	1408	& (/ ( precip_snow(ilandindex(ik), jlandindex(ik))*dt,ik=1,nbindex ) /)
	1409	WRITE(numout,*) "Down-welling long-wave flux = ", &
	1410	& (/ ( for_lwdown(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
	1411	WRITE(numout,*) "Net surface short-wave flux = ", &
	1412	& (/ ( for_swnet(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
	1413	WRITE(numout,*) "Downwelling surface short-wave flux = ", &
	1414	& (/ ( for_swdown(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
	1415	WRITE(numout,*) "Air temperature in Kelvin = ", &
	1416	& (/ ( for_tair(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
	1417	WRITE(numout,*) "Air potential energy = ", &
	1418	& (/ ( for_eair(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
	1419	WRITE(numout,*) "CO2 concentration in the canopy = ", &
	1420	& (/ ( for_ccanopy(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
	1421	WRITE(numout,*) "Coeficients A from the PBL resolution = ", &
	1422	& (/ ( petAcoef(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
	1423	WRITE(numout,*) "One for T and another for q = ", &
	1424	& (/ ( peqAcoef(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
	1425	WRITE(numout,*) "Coeficients B from the PBL resolution = ", &
	1426	& (/ ( petBcoef(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
	1427	WRITE(numout,*) "One for T and another for q = ", &
	1428	& (/ ( peqBcoef(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
	1429	WRITE(numout,*) "Cdrag = ", &
	1430	& (/ ( cdrag(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
	1431	WRITE(numout,*) "Lowest level pressure = ", &
	1432	& (/ ( for_psurf(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
	1433	WRITE(numout,*) "Geographical coordinates lon = ", &
	1434	& (/ ( lon(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
	1435	WRITE(numout,*) "Geographical coordinates lat = ", &
	1436	& (/ ( lat(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
	1437	WRITE(numout,*) "Fraction of continent in the grid = ", &
	1438	& (/ ( for_contfrac(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
	1439	ENDIF
	1440
[2596]	1441	CALL intersurf_main_2d &
[8]	1442	& (istp, iim, jjm, nbindex, kindex, dt, &
[3850]	1443	& lstep_init, lstep_last, lon, lat, for_contfrac, for_resolution, date0_rest, &
[8]	1444	! first level conditions
	1445	& zlev_vec, for_u, for_v, &
	1446	& for_qair, for_tair, for_eair, for_ccanopy, &
	1447	! Variables for the implicit coupling
	1448	& cdrag, petAcoef, peqAcoef, petBcoef, peqBcoef, &
	1449	! Rain, snow, radiation and surface pressure
	1450	& precip_rain, precip_snow, &
[898]	1451	& for_lwdown, for_swnet, for_swdown, for_psurf, &
[8]	1452	! Output : Fluxes
	1453	& vevapp, fluxsens, fluxlat, coastalflow, riverflow, &
	1454	! Surface temperatures and surface properties
[898]	1455	& tsol_rad, temp_sol_NEW, qsurf, albedo, emis, z0, &
	1456	! VOC : radiation
[2961]	1457	& for_coszang)
[898]	1458
[8]	1459	!-------
[2348]	1460	IF (printlev_loc>=4) THEN
[8]	1461	WRITE(numout,*) "dim2_driver outputs"
	1462	! Output : Fluxes
	1463	WRITE(numout,*) "vevapp ; Total of evaporation = ", &
	1464	& (/ ( vevapp(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
	1465	WRITE(numout,*) "Sensible heat flux = ", &
	1466	& (/ ( fluxsens(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
	1467	WRITE(numout,*) "Latent heat flux = ", &
	1468	& (/ ( fluxlat(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
	1469	WRITE(numout,*) "coastalflow ; Diffuse flow of water into the ocean (m^3/dt) = ", &
	1470	& (/ ( coastalflow(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
	1471	WRITE(numout,*) "riverflow ; Largest rivers flowing into the ocean (m^3/dt) = ", &
	1472	& (/ ( riverflow(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
	1473	! Surface temperatures and surface properties
	1474	WRITE(numout,*) "tsol_rad ; Radiative surface temperature = ", &
	1475	& (/ ( tsol_rad(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
	1476	WRITE(numout,*) "temp_sol_new ; New soil temperature = ", &
	1477	& (/ ( temp_sol_NEW(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
	1478	WRITE(numout,*) "qsurf ; Surface specific humidity = ", &
	1479	& (/ ( qsurf(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
	1480	WRITE(numout,*) "albedoVIS = ", &
	1481	& (/ ( albedo(ilandindex(ik), jlandindex(ik), 1),ik=1,nbindex ) /)
	1482	WRITE(numout,*) "albedoNIR = ", &
	1483	& (/ ( albedo(ilandindex(ik), jlandindex(ik), 2),ik=1,nbindex ) /)
	1484	WRITE(numout,*) "emis ; Emissivity = ", &
	1485	& (/ ( emis(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
	1486	WRITE(numout,*) "z0 ; Surface roughness = ", &
	1487	& (/ ( z0(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
	1488	ENDIF
	1489	!-----
	1490	dtdt(:,:) = zero
	1491	DO ik=1,nbindex
	1492	i=ilandindex(ik)
	1493	j=jlandindex(ik)
	1494	dtdt(i,j) = ABS(temp_sol_NEW(i,j)-temp_sol_old(i,j))/dt
	1495	ENDDO
	1496	!-----
	1497	!---- Test if the point with the largest change has more than 5K per dt
	1498	!-----
[4693]	1499	IF (printlev_loc >=3) THEN
[3941]	1500	IF (MAXVAL(dtdt(:,:)) > 5./dt) THEN
	1501	ml = MAXLOC(dtdt)
	1502	CALL ju2ymds(julian, yy, mm, dd, ss)
	1503	WRITE(numout,"('ATT :',I5,' big temperature jumps on ', &
	1504	I4,'-',I2.2,'-',I2.2,':',F8.4)") &
	1505	COUNT(dtdt(:,:) > 5./dt),yy,mm,dd,ss/3600.
	1506	WRITE(numout,*) &
	1507	'Maximum change of surface temperature located at :', &
	1508	lon(ml(1),ml(2)),lat(ml(1),ml(2))
	1509	WRITE(numout,*) 'Coordinates in grid space: ',ml(1),ml(2)
	1510	WRITE(numout,*) 'Change from ',temp_sol_old(ml(1),ml(2)), &
	1511	' to ',temp_sol_new(ml(1),ml(2)),&
	1512	'with sw_in = ',for_swnet(ml(1),ml(2))
	1513	old_tair = &
	1514	(old_eair(ml(1),ml(2))-cte_grav*old_zlev(ml(1),ml(2)))/cp_air
	1515	WRITE(numout,*) 'Air temperature change from ',old_tair, &
	1516	' to ',for_tair(ml(1),ml(2))
	1517	WRITE(numout,*) 'Max of dtdt : ',dtdt(ml(1),ml(2)),' with dt = ',dt
	1518	ENDIF
	1519	END IF
	1520
[8]	1521	temp_sol_old(:,:) = temp_sol_NEW(:,:)
	1522	!-----
[2348]	1523	!---- PRINT output value for printlev_loc>=3
[8]	1524	!-----
[2348]	1525	IF (printlev_loc>=3) THEN
[8]	1526	WRITE(numout,*) ' OUT : z0, albedoVIS, albedoNIR, emis = ', &
	1527	& z0(itest,jtest),albedo(itest,jtest,1), &
	1528	& albedo(itest,jtest,2),emis(itest,jtest)
	1529	WRITE(numout,*) ' OUT : coastal and river flow = ',&
	1530	& coastalflow(itest,jtest), riverflow(itest,jtest)
	1531	WRITE(numout,*) ' OUT : tsol_rad, vevapp = ', &
	1532	& tsol_rad(itest,jtest), vevapp(itest,jtest)
	1533	WRITE(numout,*) ' OUT : temp_sol_new =', temp_sol_NEW(itest,jtest)
	1534	ENDIF
	1535	!-----
	1536	!---- Give some variables to the output package
	1537	!---- for saving on the history tape
	1538	!-----
[2348]	1539	IF (printlev_loc>=3) WRITE(numout,*) 'history written for ', istp
[8]	1540	!-----
	1541	istp_old = istp
	1542	istp = istp+1
	1543	!-----
	1544	old_zlev(:,:) = zlev_vec(:,:)
	1545	old_qair(:,:) = for_qair(:,:)
	1546	old_eair(:,:) = for_eair(:,:)
	1547	!-----
	1548	is = is + 1
[4263]	1549	ENDDO ! DO WHILE (is <= split)
[8]	1550	split_start = 1
	1551	IF (it==itau_fin-1) THEN
[1078]	1552
[1101]	1553	CALL Write_Load_Balance(REAL(Get_cpu_time(timer_mpi),r_std))
[1078]	1554
[8]	1555	ENDIF
	1556	it = it + 1
[4263]	1557	ENDDO ! DO WHILE (it <= itau_fin)
[8]	1558	!-
	1559	! Archive in restart file the prognostic variables
	1560	!-
[2348]	1561	IF (printlev_loc>=3) WRITE(numout,*) 'Write the restart for the driver', istp_old
[8]	1562	!-
	1563	var_name = 'fluxsens'
[353]	1564	IF (is_root_prc) THEN
	1565	ALLOCATE(fluxsens_g(iim_g,jjm_g))
	1566	ELSE
	1567	ALLOCATE(fluxsens_g(0,1))
	1568	ENDIF
[1078]	1569	CALL gather2D_mpi(fluxsens , fluxsens_g)
[8]	1570	IF(is_root_prc) CALL restput (rest_id, var_name, iim_g, jjm_g, 1, istp_old, fluxsens_g)
[353]	1571	DEALLOCATE(fluxsens_g)
[8]	1572
	1573	var_name = 'vevapp'
[353]	1574	IF (is_root_prc) THEN
	1575	ALLOCATE(vevapp_g(iim_g,jjm_g))
	1576	ELSE
	1577	ALLOCATE(vevapp_g(0,1))
	1578	ENDIF
[1078]	1579	CALL gather2D_mpi( vevapp, vevapp_g)
[8]	1580	IF(is_root_prc) CALL restput (rest_id, var_name, iim_g, jjm_g, 1, istp_old, vevapp_g)
[353]	1581	DEALLOCATE(vevapp_g)
[8]	1582
	1583	var_name = 'zlev_old'
[353]	1584	IF (is_root_prc) THEN
	1585	ALLOCATE(old_zlev_g(iim_g,jjm_g))
	1586	ELSE
	1587	ALLOCATE(old_zlev_g(0,1))
	1588	ENDIF
[1078]	1589	CALL gather2D_mpi( old_zlev, old_zlev_g)
[8]	1590	IF(is_root_prc) CALL restput (rest_id, var_name, iim_g, jjm_g, 1, istp_old, old_zlev_g)
[353]	1591	DEALLOCATE(old_zlev_g)
[8]	1592
	1593	var_name = 'qair_old'
[353]	1594	IF (is_root_prc) THEN
	1595	ALLOCATE(old_qair_g(iim_g,jjm_g))
	1596	ELSE
	1597	ALLOCATE(old_qair_g(0,1))
	1598	ENDIF
[1078]	1599	CALL gather2D_mpi( old_qair, old_qair_g)
[8]	1600	IF(is_root_prc) CALL restput (rest_id, var_name, iim_g, jjm_g, 1, istp_old, old_qair_g)
[353]	1601	DEALLOCATE(old_qair_g)
[8]	1602
	1603	var_name = 'eair_old'
[353]	1604	IF (is_root_prc) THEN
	1605	ALLOCATE(old_eair_g(iim_g,jjm_g))
	1606	ELSE
	1607	ALLOCATE(old_eair_g(0,1))
	1608	ENDIF
[1078]	1609	CALL gather2D_mpi( old_eair, old_eair_g)
[8]	1610	IF(is_root_prc) CALL restput (rest_id, var_name, iim_g, jjm_g, 1, istp_old, old_eair_g)
[353]	1611	DEALLOCATE(old_eair_g)
[8]	1612
	1613	var_name = 'rau_old'
[353]	1614	IF (is_root_prc) THEN
	1615	ALLOCATE(for_rau_g(iim_g,jjm_g))
	1616	ELSE
	1617	ALLOCATE(for_rau_g(0,1))
	1618	ENDIF
[1078]	1619	CALL gather2D_mpi( for_rau, for_rau_g)
[8]	1620	IF(is_root_prc) CALL restput (rest_id, var_name, iim_g, jjm_g, 1, istp_old, for_rau_g)
[353]	1621	DEALLOCATE(for_rau_g)
[8]	1622
[353]	1623	IF (is_root_prc) THEN
	1624	ALLOCATE(albedo_g(iim_g,jjm_g))
	1625	ELSE
	1626	ALLOCATE(albedo_g(0,1))
	1627	ENDIF
[8]	1628	var_name= 'albedo_vis'
	1629	albedo_vis(:,:)=albedo(:,:,1)
[1078]	1630	CALL gather2D_mpi(albedo_vis,albedo_g)
[8]	1631	IF(is_root_prc) CALL restput (rest_id, var_name, iim_g, jjm_g, 1, istp_old, albedo_g)
	1632	!
	1633	var_name= 'albedo_nir'
	1634	albedo_nir(:,:)=albedo(:,:,2)
[1078]	1635	CALL gather2D_mpi(albedo_nir,albedo_g)
[8]	1636	IF(is_root_prc) CALL restput (rest_id, var_name, iim_g, jjm_g, 1, istp_old, albedo_g)
[353]	1637	DEALLOCATE(albedo_g)
[8]	1638
[353]	1639	IF (is_root_prc) THEN
	1640	ALLOCATE(z0_g(iim_g,jjm_g))
	1641	ELSE
	1642	ALLOCATE(z0_g(0,1))
	1643	ENDIF
[8]	1644	var_name= 'z0'
[1078]	1645	CALL gather2D_mpi(z0,z0_g)
[8]	1646	IF(is_root_prc) CALL restput (rest_id, var_name, iim_g, jjm_g, 1, istp_old, z0_g)
[353]	1647	DEALLOCATE(z0_g)
[8]	1648
	1649	if (.NOT. is_watchout) THEN
[353]	1650	var_name = 'petAcoef'
	1651	IF (is_root_prc) THEN
	1652	ALLOCATE(petAcoef_g(iim_g,jjm_g))
	1653	ELSE
	1654	ALLOCATE(petAcoef_g(0,1))
	1655	ENDIF
[1078]	1656	CALL gather2D_mpi( petAcoef, petAcoef_g)
[8]	1657	IF(is_root_prc) CALL restput (rest_id, var_name, iim_g, jjm_g, 1, istp_old, petAcoef_g)
[353]	1658	DEALLOCATE(petAcoef_g)
[8]	1659
[353]	1660	var_name = 'petBcoef'
	1661	IF (is_root_prc) THEN
	1662	ALLOCATE(petBcoef_g(iim_g,jjm_g))
	1663	ELSE
	1664	ALLOCATE(petBcoef_g(0,1))
	1665	ENDIF
[1078]	1666	CALL gather2D_mpi( petBcoef, petBcoef_g)
[8]	1667	IF(is_root_prc) CALL restput (rest_id, var_name, iim_g, jjm_g, 1, istp_old, petBcoef_g)
[353]	1668	DEALLOCATE(petBcoef_g)
[8]	1669
[353]	1670	var_name = 'peqAcoef'
	1671	IF (is_root_prc) THEN
	1672	ALLOCATE(peqAcoef_g(iim_g,jjm_g))
	1673	ELSE
	1674	ALLOCATE(peqAcoef_g(0,1))
	1675	ENDIF
[1078]	1676	CALL gather2D_mpi( peqAcoef, peqAcoef_g)
[8]	1677	IF(is_root_prc) CALL restput (rest_id, var_name, iim_g, jjm_g, 1, istp_old, peqAcoef_g)
[353]	1678	DEALLOCATE(peqAcoef_g)
[8]	1679
[353]	1680	var_name = 'peqBcoef'
	1681	IF (is_root_prc) THEN
	1682	ALLOCATE(peqBcoef_g(iim_g,jjm_g))
	1683	ELSE
	1684	ALLOCATE(peqBcoef_g(0,1))
	1685	ENDIF
[1078]	1686	CALL gather2D_mpi( peqBcoef, peqBcoef_g)
[8]	1687	IF(is_root_prc) CALL restput (rest_id, var_name, iim_g, jjm_g, 1, istp_old, peqBcoef_g)
[353]	1688	DEALLOCATE(peqBcoef_g)
[8]	1689	ENDIF
	1690	!-
[2348]	1691	IF (printlev_loc>=3) WRITE(numout,*) 'Restart for the driver written'
[8]	1692	!=====================================================================
	1693	!- 5.0 Closing all files
	1694	!=====================================================================
	1695	CALL flinclo(force_id)
[2348]	1696	IF ( printlev_loc>=3 ) WRITE(numout,*) 'FLIN CLOSED'
[8]	1697	CALL histclo
[2348]	1698	IF ( printlev_loc>=3 ) WRITE(numout,*) 'HIST CLOSED'
[8]	1699
	1700	IF(is_root_prc) THEN
	1701	CALL restclo
[2348]	1702	IF ( printlev_loc>=3 ) WRITE(numout,*) 'REST CLOSED'
[8]	1703	CALL getin_dump
[2348]	1704	IF ( printlev_loc>=3 ) WRITE(numout,*) 'GETIN CLOSED'
[8]	1705	ENDIF
[1078]	1706
	1707
	1708
	1709	WRITE(numout,*) '-------------------------------------------'
	1710	WRITE(numout,*) '------> CPU Time Global ',Get_cpu_Time(timer_global)
	1711	WRITE(numout,*) '------> CPU Time without mpi ',Get_cpu_Time(timer_mpi)
	1712	WRITE(numout,*) '------> Real Time Global ',Get_real_Time(timer_global)
	1713	WRITE(numout,*) '------> real Time without mpi ',Get_real_Time(timer_mpi)
	1714	WRITE(numout,*) '-------------------------------------------'
[4287]	1715
	1716	! Call driver_clear for deallocation and reset of initialization variables
	1717	CALL driver_clear
	1718
[1078]	1719	CALL Finalize_mpi
	1720
[4263]	1721
[1099]	1722	WRITE(numout,*) 'END of dim2_driver'
[4263]	1723
[4287]	1724
	1725	CONTAINS
	1726
	1727	!! ================================================================================================================================
	1728	!! SUBROUTINE : driver_clear
	1729	!!
	1730	!>\BRIEF Clear driver main program and call clear funcions for underlaying module intersurf
	1731	!!
	1732	!! DESCRIPTION : Deallocate memory and reset initialization variables to there original values
	1733	!! This subroutine call intersurf_clear which will call sechiba_clear.
	1734	!!
	1735	!_ ================================================================================================================================
	1736	SUBROUTINE driver_clear
	1737
	1738	CALL intersurf_clear
	1739
	1740	END SUBROUTINE driver_clear
	1741
[8]	1742	END PROGRAM driver

Note: See TracBrowser for help on using the repository browser.

Download in other formats: