Context Navigation

← Previous Revision
Next Revision →
Normal
Revision Log

etat0dyn_netcdf.F90

Last change on this file was 5186, checked in by abarral, 2 months ago
Encapsulate files in modules
File size: 22.7 KB

Rev	Line
[2293]	1	MODULE etat0dyn
[5099]	2
[5116]	3	!*******************************************************************************
	4	! Purpose: Create dynamical initial state using atmospheric fields from a
	5	! database of atmospheric to initialize the model.
	6	!-------------------------------------------------------------------------------
	7	! Comments:
[5099]	8
[5116]	9	! * This module is designed to work for Earth (and with ioipsl)
[5099]	10
[5116]	11	! * etat0dyn_netcdf routine can access to NetCDF data through the following
	12	! routine (to be called after restget):
	13	! CALL startget_dyn3d(varname, lon_in, lat_in, pls, workvar,&
	14	! champ, lon_in2, lat_in2)
[5099]	15
[5116]	16	! * Variables should have the following names in the NetCDF files:
	17	! 'U' : East ward wind (in "ECDYN.nc")
	18	! 'V' : Northward wind (in "ECDYN.nc")
	19	! 'TEMP' : Temperature (in "ECDYN.nc")
	20	! 'R' : Relative humidity (in "ECDYN.nc")
	21	! 'RELIEF' : High resolution orography (in "Relief.nc")
[5099]	22
[5116]	23	! * The land mask and corresponding weights can be:
	24	! 1) already known (in particular if etat0dyn has been called before) ;
	25	! in this case, ANY(masque(:,:)/=-99999.) = .TRUE.
	26	! 2) computed using the ocean mask from the ocean model (to ensure ocean
	27	! fractions are the same for atmosphere and ocean) for coupled runs.
	28	! File name: "o2a.nc" ; variable name: "OceMask"
	29	! 3) computed from topography file "Relief.nc" for forced runs.
[5099]	30
[5116]	31	! * There is a big mess with the longitude size. Should it be iml or iml+1 ?
	32	! I have chosen to use the iml+1 as an argument to this routine and we declare
	33	! internaly smaller fields when needed. This needs to be cleared once and for
	34	! all in LMDZ. A convention is required.
	35	!-------------------------------------------------------------------------------
	36	USE ioipsl, ONLY: flininfo, flinopen, flinget, flinclo, histclo
	37	USE lmdz_assert_eq, ONLY: assert_eq
[2597]	38	USE comconst_mod, ONLY: pi, cpp, kappa
[2600]	39	USE comvert_mod, ONLY: ap, bp, preff, pressure_exner
[2673]	40	USE temps_mod, ONLY: annee_ref, day_ref, itau_dyn, itau_phy, start_time
[5117]	41	USE lmdz_strings, ONLY: strLower
[5118]	42	USE lmdz_iniprint, ONLY: lunout, prt_level
[5134]	43	USE lmdz_comdissnew, ONLY: lstardis, nitergdiv, nitergrot, niterh, tetagdiv, &
	44	tetagrot, tetatemp, coefdis, vert_prof_dissip
[5136]	45	USE lmdz_comgeom2
[5116]	46
[5159]	47	USE lmdz_dimensions, ONLY: iim, jjm, llm, ndm
	48	USE lmdz_paramet
[2293]	49	IMPLICIT NONE
	50
	51	PRIVATE
	52	PUBLIC :: etat0dyn_netcdf
	53
[5159]	54
	55
[2293]	56	REAL, SAVE :: deg2rad
[5116]	57	INTEGER, SAVE :: iml_dyn, jml_dyn, llm_dyn, ttm_dyn, fid_dyn
	58	REAL, ALLOCATABLE, SAVE :: lon_dyn(:, :), lat_dyn(:, :), levdyn_ini(:)
	59	CHARACTER(LEN = 120), PARAMETER :: dynfname = 'ECDYN.nc'
[2293]	60
	61	CONTAINS
	62
[5116]	63	!-------------------------------------------------------------------------------
[5099]	64
[5116]	65	SUBROUTINE etat0dyn_netcdf(masque, phis)
[5099]	66
[5116]	67	!-------------------------------------------------------------------------------
	68	! Purpose: Create dynamical initial states.
	69	!-------------------------------------------------------------------------------
	70	! Notes: 1) This routine is designed to work for Earth
	71	! 2) If masque(:,:)/=-99999., masque and phis are already known.
	72	! Otherwise: compute it.
	73	!-------------------------------------------------------------------------------
	74	USE control_mod
	75	USE regr_lat_time_coefoz_m, ONLY: regr_lat_time_coefoz
	76	USE regr_pr_o3_m, ONLY: regr_pr_o3
	77	USE press_coefoz_m, ONLY: press_coefoz
	78	USE exner_hyb_m, ONLY: exner_hyb
	79	USE exner_milieu_m, ONLY: exner_milieu
[5182]	80	USE lmdz_infotrac, ONLY: nqtot, tracers
[5116]	81	USE lmdz_filtreg
	82	USE lmdz_cppkeys_wrapper, ONLY: CPPKEY_INCA
[5117]	83	USE lmdz_q_sat, ONLY: q_sat
[5186]	84	#ifndef CPP_PARA
	85	USE lmdz_dynredem, ONLY: dynredem0, dynredem1
	86	#endif
[5116]	87	IMPLICIT NONE
	88	!-------------------------------------------------------------------------------
	89	! Arguments:
	90	REAL, INTENT(INOUT) :: masque(iip1, jjp1) !--- Land-ocean mask
	91	REAL, INTENT(INOUT) :: phis (iip1, jjp1) !--- Ground geopotential
	92	!-------------------------------------------------------------------------------
	93	! Local variables:
	94	CHARACTER(LEN = 256) :: modname, fmt
	95	INTEGER :: i, j, l, ji, itau, iday, iq
	96	REAL :: xpn, xps, time, phystep
	97	REAL, DIMENSION(iip1, jjp1) :: psol
	98	REAL, DIMENSION(iip1, jjp1, llm + 1) :: p3d
	99	REAL, DIMENSION(iip1, jjp1, llm) :: uvent, t3d, tpot, qsat, qd
	100	REAL, DIMENSION(iip1, jjp1, llm) :: pk, pls, y, masse
	101	REAL, DIMENSION(iip1, jjm, llm) :: vvent
	102	REAL, DIMENSION(ip1jm, llm) :: pbarv
	103	REAL, DIMENSION(ip1jmp1, llm) :: pbaru, phi, w
	104	REAL, DIMENSION(ip1jmp1) :: pks
	105	REAL, DIMENSION(iim) :: xppn, xpps
	106	REAL, ALLOCATABLE :: q3d(:, :, :, :)
	107	!-------------------------------------------------------------------------------
	108	modname = 'etat0dyn_netcdf'
[2293]	109
[5116]	110	deg2rad = pi / 180.0
	111	y(:, :, :) = 0 !ym warning unitialized variable
[2293]	112
[5116]	113	! Compute psol AND tsol, knowing phis.
	114	!*******************************************************************************
	115	CALL start_init_dyn(rlonv, rlatu, rlonu, rlatv, phis, psol)
[2293]	116
[5116]	117	! Mid-levels pressure computation
	118	!*******************************************************************************
	119	CALL pression(ip1jmp1, ap, bp, psol, p3d) !--- Update p3d
	120	IF(pressure_exner) THEN !--- Update pk, pks
	121	CALL exner_hyb (ip1jmp1, psol, p3d, pks, pk)
	122	ELSE
	123	CALL exner_milieu(ip1jmp1, psol, p3d, pks, pk)
	124	END IF
	125	pls(:, :, :) = preff * (pk(:, :, :) / cpp)**(1. / kappa) !--- Update pls
[2293]	126
[5116]	127	! Update uvent, vvent, t3d and tpot
	128	!*******************************************************************************
	129	uvent(:, :, :) = 0.0 ; vvent(:, :, :) = 0.0 ; t3d (:, :, :) = 0.0
	130	CALL startget_dyn3d('u', rlonu, rlatu, pls, y, uvent, rlonv, rlatv)
	131	CALL startget_dyn3d('v', rlonv, rlatv, pls(:, :jjm, :), y(:, :jjm, :), vvent, &
	132	rlonu, rlatu(:jjm))
	133	CALL startget_dyn3d('t', rlonv, rlatu, pls, y, t3d, rlonu, rlatv)
	134	tpot(:, :, :) = t3d(:, :, :)
	135	CALL startget_dyn3d('tpot', rlonv, rlatu, pls, pk, tpot, rlonu, rlatv)
[2293]	136
[5116]	137	WRITE(lunout, *) 'T3D min,max:', MINVAL(t3d(:, :, :)), MAXVAL(t3d(:, :, :))
	138	WRITE(lunout, *) 'PLS min,max:', MINVAL(pls(:, :, :)), MAXVAL(pls(:, :, :))
[2293]	139
[5116]	140	! Humidity at saturation computation
	141	!*******************************************************************************
	142	WRITE(lunout, *) 'avant q_sat'
	143	CALL q_sat(llm * jjp1 * iip1, t3d, pls, qsat)
	144	WRITE(lunout, *) 'apres q_sat'
	145	WRITE(lunout, *) 'QSAT min,max:', MINVAL(qsat(:, :, :)), MAXVAL(qsat(:, :, :))
	146	! WRITE(lunout,*) 'QSAT :',qsat(10,20,:)
	147	qd (:, :, :) = 0.0
	148	CALL startget_dyn3d('q', rlonv, rlatu, pls, qsat, qd, rlonu, rlatv)
	149	ALLOCATE(q3d(iip1, jjp1, llm, nqtot)); q3d(:, :, :, :) = 0.0 ; q3d(:, :, :, 1) = qd(:, :, :)
	150	CALL flinclo(fid_dyn)
[2293]	151
[5116]	152	! Parameterization of ozone chemistry:
	153	!*******************************************************************************
	154	! Look for ozone tracer:
	155	IF (CPPKEY_INCA) THEN
	156	DO iq = 1, nqtot; IF(strLower(tracers(iq)%name)=="o3") EXIT;
	157	END DO
	158	IF(iq/=nqtot + 1) THEN
	159	CALL regr_lat_time_coefoz
	160	CALL press_coefoz
	161	CALL regr_pr_o3(p3d, q3d(:, :, :, iq))
	162	q3d(:, :, :, iq) = q3d(:, :, :, iq) * 48. / 29. !--- Mole->mass fraction
	163	END IF
	164	END IF
	165	q3d(iip1, :, :, :) = q3d(1, :, :, :)
	166
	167	! Writing
	168	!*******************************************************************************
	169	CALL inidissip(lstardis, nitergdiv, nitergrot, niterh, tetagdiv, tetagrot, &
	170	tetatemp, vert_prof_dissip)
	171	WRITE(lunout, *)'sortie inidissip'
	172	itau = 0
	173	itau_dyn = 0
	174	itau_phy = 0
	175	iday = dayref + itau / day_step
	176	time = FLOAT(itau - (iday - dayref) * day_step) / day_step
	177	IF(time>1.) THEN
	178	time = time - 1
	179	iday = iday + 1
	180	END IF
	181	day_ref = dayref
	182	annee_ref = anneeref
	183	CALL geopot(ip1jmp1, tpot, pk, pks, phis, phi)
	184	WRITE(lunout, *)'sortie geopot'
	185	CALL caldyn0(itau, uvent, vvent, tpot, psol, masse, pk, phis, &
	186	phi, w, pbaru, pbarv, time + iday - dayref)
	187	WRITE(lunout, *)'sortie caldyn0'
	188	start_time = 0.
[2336]	189	#ifdef CPP_PARA
	190	CALL dynredem0_loc( "start.nc", dayref, phis)
	191	#else
[5116]	192	CALL dynredem0("start.nc", dayref, phis)
[2336]	193	#endif
[5116]	194	WRITE(lunout, *)'sortie dynredem0'
[2336]	195	#ifdef CPP_PARA
	196	CALL dynredem1_loc( "start.nc", 0.0, vvent, uvent, tpot, q3d, masse, psol)
	197	#else
[5116]	198	CALL dynredem1("start.nc", 0.0, vvent, uvent, tpot, q3d, masse, psol)
[2336]	199	#endif
[5116]	200	WRITE(lunout, *)'sortie dynredem1'
	201	CALL histclo()
[2293]	202
[5116]	203	END SUBROUTINE etat0dyn_netcdf
[5099]	204
[5116]	205	!-------------------------------------------------------------------------------
[2293]	206
	207
[5116]	208	!-------------------------------------------------------------------------------
[5099]	209
[5116]	210	SUBROUTINE startget_dyn3d(var, lon_in, lat_in, pls, workvar, &
	211	champ, lon_in2, lat_in2)
[5128]	212
	213
[5116]	214	IMPLICIT NONE
	215	!===============================================================================
	216	! Purpose: Compute some quantities (u,v,t,q,tpot) using variables U,V,TEMP and R
	217	! (3D fields) of file dynfname.
	218	!-------------------------------------------------------------------------------
	219	! Note: An input auxilliary field "workvar" has to be specified in two cases:
	220	! * for "q": the saturated humidity.
	221	! * for "tpot": the Exner function.
	222	!===============================================================================
	223	! Arguments:
	224	CHARACTER(LEN = *), INTENT(IN) :: var
	225	REAL, INTENT(IN) :: lon_in(:) ! dim (iml)
	226	REAL, INTENT(IN) :: lat_in(:) ! dim (jml)
	227	REAL, INTENT(IN) :: pls (:, :, :) ! dim (iml, jml, lml)
	228	REAL, INTENT(IN) :: workvar(:, :, :) ! dim (iml, jml, lml)
	229	REAL, INTENT(INOUT) :: champ (:, :, :) ! dim (iml, jml, lml)
	230	REAL, INTENT(IN) :: lon_in2(:) ! dim (iml)
	231	REAL, INTENT(IN) :: lat_in2(:) ! dim (jml2)
	232	!-------------------------------------------------------------------------------
	233	! Local variables:
	234	CHARACTER(LEN = 10) :: vname
	235	CHARACTER(LEN = 256) :: msg, modname = "startget_dyn3d"
	236	INTEGER :: iml, jml, jml2, lml, il
	237	REAL :: xppn, xpps
	238	!-------------------------------------------------------------------------------
	239	iml = assert_eq([SIZE(lon_in), SIZE(pls, 1), SIZE(workvar, 1), SIZE(champ, 1), &
	240	SIZE(lon_in2)], TRIM(modname) // " iml")
	241	jml = assert_eq(SIZE(lat_in), SIZE(pls, 2), SIZE(workvar, 2), SIZE(champ, 2), &
	242	TRIM(modname) // " jml")
	243	lml = assert_eq(SIZE(pls, 3), SIZE(workvar, 3), SIZE(champ, 3), &
	244	TRIM(modname) // " lml")
	245	jml2 = SIZE(lat_in2)
[2293]	246
[5116]	247	!--- CHECK IF THE FIELD IS KNOWN
	248	SELECT CASE(var)
	249	CASE('u'); vname = 'U'
	250	CASE('v'); vname = 'V'
	251	CASE('t'); vname = 'TEMP'
	252	CASE('q'); vname = 'R'; msg = 'humidity as the saturated humidity'
	253	CASE('tpot'); msg = 'potential temperature as the Exner function'
	254	CASE DEFAULT; msg = 'No rule to extract variable ' // TRIM(var)
	255	CALL abort_gcm(modname, TRIM(msg) // ' from any data set', 1)
	256	END SELECT
[2293]	257
[5116]	258	!--- CHECK IF SOMETHING IS MISSING
	259	IF((var=='tpot'.OR.var=='q').AND.MINVAL(workvar)==MAXVAL(workvar)) THEN
	260	msg = 'Could not compute ' // TRIM(msg) // ' is missing or constant.'
	261	CALL abort_gcm(modname, TRIM(msg), 1)
	262	END IF
[2293]	263
[5116]	264	!--- INTERPOLATE 3D FIELD IF NEEDED
	265	IF(var/='tpot') CALL start_inter_3d(TRIM(vname), lon_in, lat_in, lon_in2, &
	266	lat_in2, pls, champ)
[2293]	267
[5116]	268	!--- COMPUTE THE REQUIRED FILED
	269	SELECT CASE(var)
	270	CASE('u'); DO il = 1, lml; champ(:, :, il) = champ(:, :, il) * cu(:, 1:jml);
	271	END DO
	272	champ(iml, :, :) = champ(1, :, :) !--- Eastward wind
[2293]	273
[5116]	274	CASE('v'); DO il = 1, lml; champ(:, :, il) = champ(:, :, il) * cv(:, 1:jml);
	275	END DO
	276	champ(iml, :, :) = champ(1, :, :) !--- Northward wind
[2293]	277
[5116]	278	CASE('tpot', 'q')
	279	IF(var=='tpot') THEN; champ = champ * cpp / workvar !--- Potential temperature
	280	ELSE; champ = champ * .01 * workvar !--- Relative humidity
	281	WHERE(champ<0.) champ = 1.0E-10
[2299]	282	END IF
[5116]	283	DO il = 1, lml
	284	xppn = SUM(aire(:, 1) * champ(:, 1, il)) / apoln
	285	xpps = SUM(aire(:, jml) * champ(:, jml, il)) / apols
	286	champ(:, 1, il) = xppn
	287	champ(:, jml, il) = xpps
[2299]	288	END DO
[5116]	289	END SELECT
[2293]	290
[5116]	291	END SUBROUTINE startget_dyn3d
[5099]	292
[5116]	293	!-------------------------------------------------------------------------------
[2293]	294
	295
[5116]	296	!-------------------------------------------------------------------------------
[5099]	297
[5116]	298	SUBROUTINE start_init_dyn(lon_in, lat_in, lon_in2, lat_in2, zs, psol)
[5099]	299
[5116]	300	!-------------------------------------------------------------------------------
	301	IMPLICIT NONE
	302	!===============================================================================
	303	! Purpose: Compute psol, knowing phis.
	304	!===============================================================================
	305	! Arguments:
	306	REAL, INTENT(IN) :: lon_in (:), lat_in (:) ! dim (iml) (jml)
	307	REAL, INTENT(IN) :: lon_in2(:), lat_in2(:) ! dim (iml) (jml2)
	308	REAL, INTENT(IN) :: zs (:, :) ! dim (iml,jml)
	309	REAL, INTENT(OUT) :: psol(:, :) ! dim (iml,jml)
	310	!-------------------------------------------------------------------------------
	311	! Local variables:
	312	CHARACTER(LEN = 256) :: modname = 'start_init_dyn'
	313	REAL :: date, dt
	314	INTEGER :: iml, jml, jml2, itau(1)
	315	REAL, ALLOCATABLE :: lon_rad(:), lon_ini(:), var_ana(:, :)
	316	REAL, ALLOCATABLE :: lat_rad(:), lat_ini(:)
	317	REAL, ALLOCATABLE :: z(:, :), ps(:, :), ts(:, :)
	318	!-------------------------------------------------------------------------------
	319	iml = assert_eq(SIZE(lon_in), SIZE(zs, 1), SIZE(psol, 1), SIZE(lon_in2), &
	320	TRIM(modname) // " iml")
	321	jml = assert_eq(SIZE(lat_in), SIZE(zs, 2), SIZE(psol, 2), TRIM(modname) // " jml")
	322	jml2 = SIZE(lat_in2)
[2293]	323
[5116]	324	WRITE(lunout, *) 'Opening the surface analysis'
	325	CALL flininfo(dynfname, iml_dyn, jml_dyn, llm_dyn, ttm_dyn, fid_dyn)
	326	WRITE(lunout, *) 'Values read: ', iml_dyn, jml_dyn, llm_dyn, ttm_dyn
[2293]	327
[5116]	328	ALLOCATE(lon_dyn(iml_dyn, jml_dyn), lat_dyn(iml_dyn, jml_dyn))
	329	ALLOCATE(levdyn_ini(llm_dyn))
	330	CALL flinopen(dynfname, .FALSE., iml_dyn, jml_dyn, llm_dyn, &
	331	lon_dyn, lat_dyn, levdyn_ini, ttm_dyn, itau, date, dt, fid_dyn)
[2293]	332
[5116]	333	!--- IF ANGLES ARE IN DEGREES, THEY ARE CONVERTED INTO RADIANS
	334	ALLOCATE(lon_ini(iml_dyn), lat_ini(jml_dyn))
	335	lon_ini(:) = lon_dyn(:, 1); IF(MAXVAL(lon_dyn)>pi) lon_ini = lon_ini * deg2rad
	336	lat_ini(:) = lat_dyn(1, :); IF(MAXVAL(lat_dyn)>pi) lat_ini = lat_ini * deg2rad
[2293]	337
[5116]	338	ALLOCATE(var_ana(iml_dyn, jml_dyn), lon_rad(iml_dyn), lat_rad(jml_dyn))
	339	CALL get_var_dyn('Z', z) !--- SURFACE GEOPOTENTIAL
	340	CALL get_var_dyn('SP', ps) !--- SURFACE PRESSURE
	341	CALL get_var_dyn('ST', ts) !--- SURFACE TEMPERATURE
	342	! CALL flinclo(fid_dyn)
	343	DEALLOCATE(var_ana, lon_rad, lat_rad, lon_ini, lat_ini)
[2293]	344
[5116]	345	!--- PSOL IS COMPUTED IN PASCALS
	346	psol(:iml - 1, :) = ps(:iml - 1, :) * (1.0 + (z(:iml - 1, :) - zs(:iml - 1, :)) / 287.0 &
	347	/ ts(:iml - 1, :))
	348	psol(iml, :) = psol(1, :)
	349	DEALLOCATE(z, ps, ts)
	350	psol(:, 1) = SUM(aire(1:iml - 1, 1) * psol(1:iml - 1, 1)) / apoln !--- NORTH POLE
	351	psol(:, jml) = SUM(aire(1:iml - 1, jml) * psol(1:iml - 1, jml)) / apols !--- SOUTH POLE
[2293]	352
[5116]	353	CONTAINS
[2293]	354
[5116]	355	!-------------------------------------------------------------------------------
[5099]	356
[5116]	357	SUBROUTINE get_var_dyn(title, field)
[5099]	358
[5116]	359	!-------------------------------------------------------------------------------
	360	USE conf_dat_m, ONLY: conf_dat2d
	361	IMPLICIT NONE
	362	!-------------------------------------------------------------------------------
	363	! Arguments:
	364	CHARACTER(LEN = *), INTENT(IN) :: title
	365	REAL, ALLOCATABLE, INTENT(INOUT) :: field(:, :)
	366	!-------------------------------------------------------------------------------
	367	! Local variables:
	368	CHARACTER(LEN = 256) :: msg
	369	INTEGER :: tllm
	370	!-------------------------------------------------------------------------------
	371	SELECT CASE(title)
	372	CASE('Z'); tllm = 0; msg = 'geopotential'
	373	CASE('SP'); tllm = 0; msg = 'surface pressure'
	374	CASE('ST'); tllm = llm_dyn; msg = 'temperature'
	375	END SELECT
	376	IF(.NOT.ALLOCATED(field)) THEN
	377	ALLOCATE(field(iml, jml))
	378	CALL flinget(fid_dyn, title, iml_dyn, jml_dyn, tllm, ttm_dyn, 1, 1, var_ana)
	379	CALL conf_dat2d(title, lon_ini, lat_ini, lon_rad, lat_rad, var_ana, .TRUE.)
	380	CALL interp_startvar(title, .TRUE., lon_rad, lat_rad, var_ana, &
	381	lon_in, lat_in, lon_in2, lat_in2, field)
	382	ELSE IF(SIZE(field)/=SIZE(z)) THEN
	383	msg = 'The ' // TRIM(msg) // ' field we have does not have the right size'
	384	CALL abort_gcm(TRIM(modname), msg, 1)
	385	END IF
[2293]	386
[5116]	387	END SUBROUTINE get_var_dyn
[5099]	388
[5116]	389	!-------------------------------------------------------------------------------
[2293]	390
[5116]	391	END SUBROUTINE start_init_dyn
[5099]	392
[5116]	393	!-------------------------------------------------------------------------------
[2293]	394
	395
[5116]	396	!-------------------------------------------------------------------------------
[5099]	397
[5116]	398	SUBROUTINE start_inter_3d(var, lon_in, lat_in, lon_in2, lat_in2, pls_in, var3d)
[5099]	399
[5116]	400	!-------------------------------------------------------------------------------
	401	USE conf_dat_m, ONLY: conf_dat3d
	402	USE lmdz_libmath_pch, ONLY: pchsp_95, pchfe_95
	403	USE lmdz_libmath, ONLY: minmax
[2293]	404
[5128]	405
[5116]	406	IMPLICIT NONE
	407	!-------------------------------------------------------------------------------
	408	! Arguments:
	409	CHARACTER(LEN = *), INTENT(IN) :: var
	410	REAL, INTENT(IN) :: lon_in(:), lat_in(:) ! dim (iml) (jml)
	411	REAL, INTENT(IN) :: lon_in2(:), lat_in2(:) ! dim (iml) (jml2)
	412	REAL, INTENT(IN) :: pls_in(:, :, :) ! dim (iml,jml,lml)
	413	REAL, INTENT(OUT) :: var3d (:, :, :) ! dim (iml,jml,lml)
	414	!-------------------------------------------------------------------------------
	415	! Local variables:
	416	CHARACTER(LEN = 256) :: modname = 'start_inter_3d'
	417	LOGICAL :: skip
	418	REAL :: chmin, chmax
	419	INTEGER :: iml, jml, lml, jml2, ii, ij, il, ierr
	420	INTEGER :: n_extrap ! Extrapolated points number
	421	REAL, ALLOCATABLE :: ax(:), lon_rad(:), lon_ini(:), lev_dyn(:), yder(:)
	422	REAL, ALLOCATABLE :: ay(:), lat_rad(:), lat_ini(:), var_tmp3d(:, :, :)
	423	REAL, ALLOCATABLE, SAVE :: var_ana3d(:, :, :)
	424	!-------------------------------------------------------------------------------
	425	iml = assert_eq(SIZE(lon_in), SIZE(lon_in2), SIZE(pls_in, 1), SIZE(var3d, 1), TRIM(modname) // " iml")
	426	jml = assert_eq(SIZE(lat_in), SIZE(pls_in, 2), SIZE(var3d, 2), TRIM(modname) // " jml")
	427	lml = assert_eq(SIZE(pls_in, 3), SIZE(var3d, 3), TRIM(modname) // " lml"); jml2 = SIZE(lat_in2)
[2293]	428
[5116]	429	WRITE(lunout, *)'Going into flinget to extract the 3D field.'
	430	IF(.NOT.ALLOCATED(var_ana3d)) ALLOCATE(var_ana3d(iml_dyn, jml_dyn, llm_dyn))
	431	CALL flinget(fid_dyn, var, iml_dyn, jml_dyn, llm_dyn, ttm_dyn, 1, 1, var_ana3d)
[2293]	432
[5116]	433	!--- ANGLES IN DEGREES ARE CONVERTED INTO RADIANS
	434	ALLOCATE(lon_ini(iml_dyn), lat_ini(jml_dyn))
	435	lon_ini(:) = lon_dyn(:, 1); IF(MAXVAL(lon_dyn)>pi) lon_ini = lon_ini * deg2rad
	436	lat_ini(:) = lat_dyn(1, :); IF(MAXVAL(lat_dyn)>pi) lat_ini = lat_ini * deg2rad
[2293]	437
[5116]	438	!--- FIELDS ARE PROCESSED TO BE ON STANDARD ANGULAR DOMAINS
	439	ALLOCATE(lon_rad(iml_dyn), lat_rad(jml_dyn), lev_dyn(llm_dyn))
	440	CALL conf_dat3d(var, lon_ini, lat_ini, levdyn_ini, &
	441	lon_rad, lat_rad, lev_dyn, var_ana3d, .TRUE.)
	442	DEALLOCATE(lon_ini, lat_ini)
[2293]	443
[5116]	444	!--- COMPUTE THE REQUIRED FIELDS USING ROUTINE grid_noro
	445	ALLOCATE(var_tmp3d(iml, jml, llm_dyn))
	446	DO il = 1, llm_dyn
	447	CALL interp_startvar(var, il==1, lon_rad, lat_rad, var_ana3d(:, :, il), &
	448	lon_in, lat_in, lon_in2, lat_in2, var_tmp3d(:, :, il))
[2293]	449	END DO
[5116]	450	DEALLOCATE(lon_rad, lat_rad)
[2293]	451
[5116]	452	!--- VERTICAL INTERPOLATION FROM TOP OF ATMOSPHERE TO GROUND
	453	ALLOCATE(ax(llm_dyn), ay(llm_dyn), yder(llm_dyn))
	454	ax = lev_dyn(llm_dyn:1:-1)
	455	skip = .FALSE.
	456	n_extrap = 0
	457	DO ij = 1, jml
	458	DO ii = 1, iml - 1
	459	ay = var_tmp3d(ii, ij, llm_dyn:1:-1)
	460	yder = pchsp_95(ax, ay, ibeg = 2, iend = 2, vc_beg = 0., vc_end = 0.)
	461	CALL pchfe_95(ax, ay, yder, skip, pls_in(ii, ij, lml:1:-1), &
	462	var3d(ii, ij, lml:1:-1), ierr)
	463	IF(ierr<0) CALL abort_gcm(TRIM(modname), 'error in pchfe_95', 1)
	464	n_extrap = n_extrap + ierr
	465	END DO
	466	END DO
	467	IF(n_extrap/=0) WRITE(lunout, *)TRIM(modname) // " pchfe_95: n_extrap=", n_extrap
	468	var3d(iml, :, :) = var3d(1, :, :)
[2293]	469
[5116]	470	END SUBROUTINE start_inter_3d
[2293]	471
[5116]	472	!-------------------------------------------------------------------------------
[2293]	473
[5099]	474
[5116]	475	!-------------------------------------------------------------------------------
[5099]	476
[5116]	477	SUBROUTINE interp_startvar(nam, ibeg, lon, lat, vari, lon1, lat1, lon2, lat2, varo)
[2293]	478
[5116]	479	!-------------------------------------------------------------------------------
	480	USE inter_barxy_m, ONLY: inter_barxy
	481	IMPLICIT NONE
	482	!-------------------------------------------------------------------------------
	483	! Arguments:
	484	CHARACTER(LEN = *), INTENT(IN) :: nam
	485	LOGICAL, INTENT(IN) :: ibeg
	486	REAL, INTENT(IN) :: lon(:), lat(:) ! dim (ii) (jj)
	487	REAL, INTENT(IN) :: vari(:, :) ! dim (ii,jj)
	488	REAL, INTENT(IN) :: lon1(:), lat1(:) ! dim (i1) (j1)
	489	REAL, INTENT(IN) :: lon2(:), lat2(:) ! dim (i1) (j2)
	490	REAL, INTENT(OUT) :: varo(:, :) ! dim (i1) (j1)
	491	!-------------------------------------------------------------------------------
	492	! Local variables:
	493	CHARACTER(LEN = 256) :: modname = "interp_startvar"
	494	INTEGER :: ii, jj, i1, j1, j2
	495	REAL, ALLOCATABLE :: vtmp(:, :)
	496	!-------------------------------------------------------------------------------
	497	ii = assert_eq(SIZE(lon), SIZE(vari, 1), TRIM(modname) // " ii")
	498	jj = assert_eq(SIZE(lat), SIZE(vari, 2), TRIM(modname) // " jj")
	499	i1 = assert_eq(SIZE(lon1), SIZE(lon2), SIZE(varo, 1), TRIM(modname) // " i1")
	500	j1 = assert_eq(SIZE(lat1), SIZE(varo, 2), TRIM(modname) // " j1")
	501	j2 = SIZE(lat2)
	502	ALLOCATE(vtmp(i1 - 1, j1))
	503	IF(ibeg.AND.prt_level>1) THEN
	504	WRITE(lunout, *)"---------------------------------------------------------"
	505	WRITE(lunout, *)"$$$ Interpolation barycentrique pour " // TRIM(nam) // " $$$"
	506	WRITE(lunout, *)"---------------------------------------------------------"
	507	END IF
	508	CALL inter_barxy(lon, lat(:jj - 1), vari, lon2(:i1 - 1), lat2, vtmp)
	509	CALL gr_int_dyn(vtmp, varo, i1 - 1, j1)
[5099]	510
[5116]	511	END SUBROUTINE interp_startvar
[2293]	512
[5116]	513	!-------------------------------------------------------------------------------
	514
[2293]	515	END MODULE etat0dyn
[5099]	516
[2293]	517	!*******************************************************************************

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

Download in other formats:

Original Format