Context Navigation

etat0dyn_netcdf.F90 @ 5285

Last change on this file since 5285 was 5285, checked in by abarral, 7 weeks ago
As discussed internally, remove generic ONLY: ... for new _mod_h modules
File size: 21.9 KB

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

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

Download in other formats:

Original Format