Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
Revision Log

startvar.F90 @ 937

Last change on this file since 937 was 808, checked in by slebonnois, 12 years ago
SL: Many changes for VENUS (related to newstart) and TITAN (related to clouds). Please read DOC/chantiers/commit_importants.log (cf v808).
File size: 30.6 KB

Rev	Line
[1]	1	!
	2	! $Id: startvar.F90 1425 2010-09-02 13:45:23Z lguez $
	3	!
	4	!*******************************************************************************
	5	!
	6	MODULE startvar
	7	!
	8	!*******************************************************************************
	9	!
	10	!-------------------------------------------------------------------------------
	11	! Purpose: Access data from the database of atmospheric to initialize the model.
	12	!-------------------------------------------------------------------------------
	13	! Comments:
	14	!
	15	! * This module is designed to work for Earth (and with ioipsl)
	16	!
	17	! * There are three ways to acces data, depending on the type of field
	18	! which needs to be extracted. In any case the call should come after a restget
	19	! and should be of the type : CALL startget(...)
	20	!
	21	! - A 1D variable on the physical grid :
	22	! CALL startget_phys1d((varname, iml, jml, lon_in, lat_in, nbindex, &
	23	! champ, val_exp, jml2, lon_in2, lat_in2, ibar )
	24	!
	25	! - A 2D variable on the dynamical grid :
	26	! CALL startget_phys2d(varname, iml, jml, lon_in, lat_in, &
	27	! champ, val_exp, jml2, lon_in2, lat_in2, ibar )
	28	!
	29	! - A 3D variable on the dynamical grid :
	30	! CALL startget_dyn((varname, iml, jml, lon_in, lat_in, lml, pls, workvar, &
	31	! champ, val_exp, jml2, lon_in2, lat_in2, ibar )
	32	!
	33	! * Data needs to be in NetCDF format
	34	!
	35	! * Variables should have the following names in the files:
	36	! 'RELIEF' : High resolution orography
	37	! 'ST' : Surface temperature
	38	! 'CDSW' : Soil moisture
	39	! 'Z' : Surface geopotential
	40	! 'SP' : Surface pressure
	41	! 'U' : East ward wind
	42	! 'V' : Northward wind
	43	! 'TEMP' : Temperature
	44	! 'R' : Relative humidity
	45	!
	46	! * There is a big mess with the longitude size. Should it be iml or iml+1 ?
	47	! I have chosen to use the iml+1 as an argument to this routine and we declare
	48	! internaly smaller fields when needed. This needs to be cleared once and for
	49	! all in LMDZ. A convention is required.
	50	!-------------------------------------------------------------------------------
	51	USE ioipsl
	52	IMPLICIT NONE
	53
	54	PRIVATE
	55	PUBLIC startget_phys2d, startget_phys1d, startget_dyn
	56	! INTERFACE startget
	57	! MODULE PROCEDURE startget_phys1d, startget_phys2d, startget_dyn
	58	! END INTERFACE
	59
	60	REAL, SAVE :: deg2rad, pi
	61	INTEGER, SAVE :: iml_rel, jml_rel
	62	INTEGER, SAVE :: fid_phys, iml_phys, jml_phys
	63	INTEGER, SAVE :: fid_dyn, iml_dyn, jml_dyn, llm_dyn, ttm_dyn
	64	REAL, DIMENSION(:,:), ALLOCATABLE, TARGET, SAVE :: lon_phys, lon_dyn
	65	REAL, DIMENSION(:,:), ALLOCATABLE, TARGET, SAVE :: lat_phys, lat_dyn
	66	REAL, DIMENSION(:,:), ALLOCATABLE, TARGET, SAVE :: lon_rug, lon_alb, lon_rel
	67	REAL, DIMENSION(:,:), ALLOCATABLE, TARGET, SAVE :: lat_rug, lat_alb, lat_rel
	68	REAL, DIMENSION(:), ALLOCATABLE, TARGET, SAVE :: levdyn_ini
	69	REAL, DIMENSION(:,:), ALLOCATABLE, TARGET, SAVE :: relief, zstd, zsig, zgam
[808]	70	REAL, DIMENSION(:,:), ALLOCATABLE, TARGET, SAVE :: zthe, zpic, zval
[1]	71	REAL, DIMENSION(:,:), ALLOCATABLE, TARGET, SAVE :: rugo, phis, tsol, qsol
	72	REAL, DIMENSION(:,:), ALLOCATABLE, TARGET, SAVE :: psol_dyn
	73	REAL, DIMENSION(:,:,:), ALLOCATABLE, TARGET, SAVE :: var_ana3d
	74
	75	CONTAINS
	76
	77	!-------------------------------------------------------------------------------
	78	!
	79	SUBROUTINE startget_phys1d(varname, iml, jml, lon_in, lat_in, nbindex, champ, &
	80	val_exp ,jml2, lon_in2, lat_in2, ibar)
	81	!
	82	!-------------------------------------------------------------------------------
	83	! Comment:
	84	! This routine only works if the variable does not exist or is constant.
	85	!-------------------------------------------------------------------------------
	86	! Arguments:
	87	CHARACTER(LEN=*), INTENT(IN) :: varname
	88	INTEGER, INTENT(IN) :: iml, jml
	89	REAL, DIMENSION(iml), INTENT(IN) :: lon_in
	90	REAL, DIMENSION(jml), INTENT(IN) :: lat_in
	91	INTEGER, INTENT(IN) :: nbindex
	92	REAL, DIMENSION(nbindex), INTENT(INOUT) :: champ
	93	REAL, INTENT(IN) :: val_exp
	94	INTEGER, INTENT(IN) :: jml2
	95	REAL, DIMENSION(iml), INTENT(IN) :: lon_in2
	96	REAL, DIMENSION(jml2), INTENT(IN) :: lat_in2
	97	LOGICAL, INTENT(IN) :: ibar
	98	!-------------------------------------------------------------------------------
	99	! Local variables:
	100	#include "iniprint.h"
	101	REAL, DIMENSION(:,:), POINTER :: v2d
	102	!-------------------------------------------------------------------------------
	103	v2d=>NULL()
	104	IF(MINVAL(champ)==MAXVAL(champ).AND.MINVAL(champ)==val_exp) THEN
	105
	106	!--- CHECKING IF THE FIELD IS KNOWN ; READING UNALLOCATED FILES
	107	SELECT CASE(varname)
	108	CASE('tsol')
	109	IF(.NOT.ALLOCATED(tsol)) &
	110	CALL start_init_phys(iml,jml,lon_in,lat_in,jml2,lon_in2,lat_in2,ibar)
	111	CASE('qsol')
	112	IF(.NOT.ALLOCATED(qsol)) &
	113	CALL start_init_phys(iml,jml,lon_in,lat_in,jml2,lon_in2,lat_in2,ibar)
	114	CASE('psol')
	115	IF(.NOT.ALLOCATED(psol_dyn)) &
	116	CALL start_init_dyn (iml,jml,lon_in,lat_in,jml2,lon_in2,lat_in2,ibar)
	117	CASE('zmea','zstd','zsig','zgam','zthe','zpic','zval')
	118	IF(.NOT.ALLOCATED(relief)) &
	119	CALL start_init_orog(iml,jml,lon_in,lat_in)
	120	CASE('rads','snow','tslab','seaice','rugmer','agsno')
	121	CASE DEFAULT
	122	WRITE(lunout,*)'startget_phys1d'
	123	WRITE(lunout,*)'No rule is present to extract variable '//TRIM(varname)&
	124	//' from any data set'; STOP
	125	END SELECT
	126
	127	!--- SELECTING v2d FOR WANTED VARIABLE AND CHEKING ITS SIZE
	128	SELECT CASE(varname)
	129	CASE('rads','snow','tslab','seaice'); champ=0.0
	130	CASE('rugmer'); champ(:)=0.001
	131	CASE('agsno'); champ(:)=50.0
	132	CASE DEFAULT
	133	SELECT CASE(varname)
	134	CASE('tsol'); v2d=>tsol
	135	CASE('qsol'); v2d=>qsol
	136	CASE('psol'); v2d=>psol_dyn
	137	CASE('zmea'); v2d=>relief
	138	CASE('zstd'); v2d=>zstd
	139	CASE('zsig'); v2d=>zsig
	140	CASE('zgam'); v2d=>zgam
	141	CASE('zthe'); v2d=>zthe
	142	CASE('zpic'); v2d=>zpic
	143	CASE('zval'); v2d=>zval
	144	END SELECT
	145	IF(SIZE(v2d)/=SIZE(lon_in)*SIZE(lat_in)) THEN
	146	WRITE(lunout,*)'STARTVAR module has been initialized to the wrong size'
	147	STOP
	148	END IF
	149	CALL gr_dyn_fi(1,iml,jml,nbindex,v2d,champ)
	150	END SELECT
	151
	152	ELSE
	153
	154	!--- SOME FIELDS ARE CAUGHT: MAY BE NEEDED FOR A 3D INTEPROLATION
	155	SELECT CASE(varname)
	156	CASE('tsol')
	157	IF(.NOT.ALLOCATED(tsol)) ALLOCATE(tsol(iml,jml))
	158	CALL gr_fi_dyn(1,iml,jml,nbindex,champ,tsol)
	159	END SELECT
	160
	161	END IF
	162
	163	END SUBROUTINE startget_phys1d
	164	!
	165	!-------------------------------------------------------------------------------
	166
	167
	168	!-------------------------------------------------------------------------------
	169	!
	170	SUBROUTINE startget_phys2d(varname, iml, jml, lon_in, lat_in, champ, val_exp, &
[808]	171	jml2, lon_in2, lat_in2 , ibar)
[1]	172	!
	173	!-------------------------------------------------------------------------------
	174	! Comment:
	175	! This routine only works if the variable does not exist or is constant.
	176	!-------------------------------------------------------------------------------
	177	! Arguments:
	178	CHARACTER(LEN=*), INTENT(IN) :: varname
	179	INTEGER, INTENT(IN) :: iml, jml
	180	REAL, DIMENSION(iml), INTENT(IN) :: lon_in
	181	REAL, DIMENSION(jml), INTENT(IN) :: lat_in
	182	REAL, DIMENSION(iml,jml), INTENT(INOUT) :: champ
	183	REAL, INTENT(IN) :: val_exp
	184	INTEGER, INTENT(IN) :: jml2
	185	REAL, DIMENSION(iml), INTENT(IN) :: lon_in2
	186	REAL, DIMENSION(jml2), INTENT(IN) :: lat_in2
	187	LOGICAL, INTENT(IN) :: ibar
	188	!-------------------------------------------------------------------------------
	189	! Local variables:
	190	#include "iniprint.h"
	191	REAL, DIMENSION(:,:), POINTER :: v2d=>NULL()
	192	!-------------------------------------------------------------------------------
	193	v2d=>NULL()
	194	IF(MINVAL(champ)==MAXVAL(champ).AND.MINVAL(champ)==val_exp) THEN
	195
	196	!--- CHECKING IF THE FIELD IS KNOWN ; READING UNALLOCATED FILES
	197	SELECT CASE(varname)
	198	CASE('psol')
	199	IF(.NOT.ALLOCATED(psol_dyn)) &
	200	CALL start_init_dyn (iml,jml,lon_in,lat_in,jml2,lon_in2,lat_in2,ibar)
	201	CASE('relief')
[808]	202	IF(.NOT.ALLOCATED(relief)) CALL start_init_orog(iml,jml,lon_in,lat_in)
[1]	203	CASE('surfgeo')
	204	IF(.NOT.ALLOCATED(phis)) CALL start_init_orog(iml,jml,lon_in,lat_in)
[808]	205	CASE('rugosite')
[1]	206	IF(.NOT.ALLOCATED(rugo)) CALL start_init_orog(iml,jml,lon_in,lat_in)
	207	CASE DEFAULT
	208	WRITE(lunout,*)'startget_phys2d'
	209	WRITE(lunout,*)'No rule is present to extract variable '//TRIM(varname)&
	210	//' from any data set'; STOP
	211	END SELECT
	212
	213	!--- SELECTING v2d FOR WANTED VARIABLE AND CHEKING ITS SIZE
	214	SELECT CASE(varname)
	215	CASE('psol'); v2d=>psol_dyn
	216	CASE('relief'); v2d=>relief
	217	CASE('rugosite'); v2d=>rugo
	218	CASE('surfgeo'); v2d=>phis
	219	END SELECT
	220	IF(SIZE(champ)/=SIZE(v2d)) THEN
	221	WRITE(lunout,*) 'STARTVAR module has been initialized to the wrong size'
	222	STOP
	223	END IF
	224
	225	champ(:,:)=v2d(:,:)
	226
	227	ELSE
	228
	229	!--- SOME FIELDS ARE CAUGHT: MAY BE NEEDED FOR A 3D INTEPROLATION
	230	SELECT CASE(varname)
	231	CASE ('surfgeo')
	232	IF(.NOT.ALLOCATED(phis)) ALLOCATE(phis(iml,jml))
	233	IF(SIZE(phis)/=SIZE(champ)) THEN
	234	WRITE(lunout,*)'STARTVAR module has been initialized to the wrong size'
	235	STOP
	236	END IF
	237	phis(:,:)=champ(:,:)
	238	END SELECT
	239
	240	END IF
	241
	242	END SUBROUTINE startget_phys2d
	243	!
	244	!-------------------------------------------------------------------------------
	245
	246
	247	!-------------------------------------------------------------------------------
	248	!
	249	SUBROUTINE startget_dyn(varname, lon_in, lat_in, pls,workvar,&
	250	champ, val_exp, lon_in2, lat_in2, ibar)
	251
	252	use assert_eq_m, only: assert_eq
	253
	254
	255	!-------------------------------------------------------------------------------
	256	! Comment:
	257	! This routine only works if the variable does not exist or is constant.
	258	!-------------------------------------------------------------------------------
	259	! Arguments:
	260	CHARACTER(LEN=*), INTENT(IN) :: varname
	261	REAL, INTENT(IN) :: lon_in(:) ! dim(iml)
	262	REAL, INTENT(IN) :: lat_in(:) ! dim(jml)
	263	REAL, INTENT(IN) :: pls(:, :, :) ! dim(iml, jml, lml)
	264	REAL, INTENT(IN) :: workvar(:, :, :) ! dim(iml, jml, lml)
	265	REAL, INTENT(INOUT) :: champ(:, :, :) ! dim(iml, jml, lml)
	266	REAL, INTENT(IN) :: val_exp
	267	REAL, INTENT(IN) :: lon_in2(:) ! dim(iml)
	268	REAL, INTENT(IN) :: lat_in2(:) ! dim(jml2)
	269	LOGICAL, INTENT(IN) :: ibar
	270	!-------------------------------------------------------------------------------
	271	! Local variables:
	272	#include "iniprint.h"
	273	#include "dimensions.h"
	274	#include "comconst.h"
	275	#include "paramet.h"
	276	#include "comgeom2.h"
	277	INTEGER :: iml, jml
	278	INTEGER :: lml
	279	INTEGER :: jml2
	280	REAL, DIMENSION(:,:,:), POINTER :: v3d=>NULL()
	281	CHARACTER(LEN=10) :: vname
	282	INTEGER :: il
	283	REAL :: xppn, xpps
	284	!-------------------------------------------------------------------------------
	285	NULLIFY(v3d)
	286	IF(MINVAL(champ)==MAXVAL(champ).AND.MINVAL(champ)==val_exp) THEN
	287
	288	iml = assert_eq((/size(lon_in), size(pls, 1), size(workvar, 1), &
	289	& size(champ, 1), size(lon_in2)/), "startget_dyn iml")
	290	jml = assert_eq(size(lat_in), size(pls, 2), size(workvar, 2), &
	291	& size(champ, 2), "startget_dyn jml")
	292	lml = assert_eq(size(pls, 3), size(workvar, 3), size(champ, 3), &
	293	& "startget_dyn lml")
	294	jml2 = size(lat_in2)
	295
	296	!--- READING UNALLOCATED FILES
	297	IF(.NOT.ALLOCATED(psol_dyn)) &
	298	CALL start_init_dyn(iml,jml,lon_in,lat_in,jml2,lon_in2,lat_in2,ibar)
	299
	300	!--- CHECKING IF THE FIELD IS KNOWN AND INTERPOLATING 3D FIELDS
	301	SELECT CASE(varname)
	302	CASE('u'); vname='U'
	303	CASE('v'); vname='V'
	304	CASE('t','tpot'); vname='TEMP'
	305	CASE('q'); vname='R'
	306	CASE DEFAULT
	307	WRITE(lunout,*)'startget_dyn'
	308	WRITE(lunout,*)'No rule is present to extract variable '//TRIM(varname)&
	309	//' from any data set'; STOP
	310	END SELECT
	311	CALL start_inter_3d(TRIM(vname), iml, jml, lml, lon_in, lat_in, jml2, &
	312	lon_in2, lat_in2, pls, champ,ibar )
	313
	314	!--- COMPUTING THE REQUIRED FILED
	315	SELECT CASE(varname)
	316	CASE('u') !--- Eastward wind
	317	DO il=1,lml; champ(:,:,il)=champ(:,:,il)*cu(:,1:jml); END DO
	318	champ(iml,:,:)=champ(1,:,:)
	319
	320	CASE('v') !--- Northward wind
	321	DO il=1,lml; champ(:,:,il)=champ(:,:,il)*cv(:,1:jml); END DO
	322	champ(iml,:,:)=champ(1,:,:)
	323
	324	CASE('tpot') !--- Temperature
	325	IF(MINVAL(workvar)/=MAXVAL(workvar)) THEN
	326	champ=champ*cpp/workvar
	327	DO il=1,lml
	328	xppn = SUM(aire(:,1 )*champ(:,1 ,il))/apoln
	329	xpps = SUM(aire(:,jml)*champ(:,jml,il))/apols
	330	champ(:,1 ,il) = xppn
	331	champ(:,jml,il) = xpps
	332	END DO
	333	ELSE
	334	WRITE(lunout,*)'Could not compute potential temperature as the'
	335	WRITE(lunout,*)'Exner function is missing or constant.'; STOP
	336	END IF
	337
	338	CASE('q') !--- Relat. humidity
	339	IF(MINVAL(workvar)/=MAXVAL(workvar)) THEN
	340	champ=0.01champworkvar
	341	WHERE(champ<0.) champ=1.0E-10
	342	DO il=1,lml
	343	xppn = SUM(aire(:,1 )*champ(:,1 ,il))/apoln
	344	xpps = SUM(aire(:,jml)*champ(:,jml,il))/apols
	345	champ(:,1 ,il) = xppn
	346	champ(:,jml,il) = xpps
	347	END DO
	348	ELSE
	349	WRITE(lunout,*)'Could not compute specific humidity as the'
	350	WRITE(lunout,*)'saturated humidity is missing or constant.'; STOP
	351	END IF
	352
	353	END SELECT
	354
	355	END IF
	356
	357	END SUBROUTINE startget_dyn
	358	!
	359	!-------------------------------------------------------------------------------
	360
	361
	362	!-------------------------------------------------------------------------------
	363	!
[808]	364	SUBROUTINE start_init_orog(iml,jml,lon_in,lat_in)
[1]	365	!
	366	!-------------------------------------------------------------------------------
	367	! Arguments:
	368	INTEGER, INTENT(IN) :: iml, jml
	369	REAL, DIMENSION(iml), INTENT(IN) :: lon_in
	370	REAL, DIMENSION(jml), INTENT(IN) :: lat_in
	371	!-------------------------------------------------------------------------------
	372	! Local variables:
	373	#include "iniprint.h"
[808]	374	#include "comconst.h"
[1]	375	CHARACTER(LEN=25) :: title
	376	CHARACTER(LEN=120) :: orofname
	377	LOGICAL :: check=.TRUE.
	378	REAL, DIMENSION(1) :: lev
	379	REAL :: date, dt
	380	INTEGER, DIMENSION(1) :: itau
	381	INTEGER :: fid, llm_tmp, ttm_tmp
	382	REAL, DIMENSION(:,:), ALLOCATABLE :: relief_hi, tmp_var
	383	REAL, DIMENSION(:), ALLOCATABLE :: lon_rad, lat_rad, lon_ini, lat_ini
	384	!-------------------------------------------------------------------------------
	385	pi=2.0*ASIN(1.0); deg2rad=pi/180.0
	386
	387	orofname = 'Relief.nc'; title='RELIEF'
	388	IF(check) WRITE(lunout,*)'Reading the high resolution orography'
	389	CALL flininfo(orofname, iml_rel, jml_rel, llm_tmp, ttm_tmp, fid)
	390
	391	ALLOCATE(lat_rel(iml_rel,jml_rel),lon_rel(iml_rel,jml_rel))
	392	CALL flinopen(orofname, .FALSE., iml_rel, jml_rel, llm_tmp, lon_rel, lat_rel,&
	393	lev, ttm_tmp, itau, date, dt, fid)
	394	ALLOCATE(relief_hi(iml_rel,jml_rel))
	395	CALL flinget(fid, title, iml_rel, jml_rel, llm_tmp, ttm_tmp, 1, 1, relief_hi)
	396	CALL flinclo(fid)
	397
	398	!--- IF ANGLES ARE IN DEGREES, THEY ARE CONVERTED INTO RADIANS
	399	ALLOCATE(lon_ini(iml_rel),lat_ini(jml_rel))
	400	lon_ini(:)=lon_rel(:,1); IF(MAXVAL(lon_rel)>pi) lon_ini=lon_ini*deg2rad
	401	lat_ini(:)=lat_rel(1,:); IF(MAXVAL(lat_rel)>pi) lat_ini=lat_ini*deg2rad
	402
	403	!--- FIELDS ARE PROCESSED TO BE ON STANDARD ANGULAR DOMAINS
	404	ALLOCATE(lon_rad(iml_rel),lat_rad(jml_rel))
	405	CALL conf_dat2d(title, iml_rel, jml_rel, lon_ini, lat_ini, lon_rad, lat_rad, &
	406	relief_hi, .FALSE.)
	407	DEALLOCATE(lon_ini,lat_ini)
	408
	409	!--- COMPUTING THE REQUIRED FIELDS USING ROUTINE grid_noro
	410	IF(check) WRITE(lunout,*)'Computes all parameters needed for gravity wave dra&
	411	&g code'
	412
	413	ALLOCATE(phis(iml,jml)) ! Geopotentiel au sol
	414	ALLOCATE(zstd(iml,jml)) ! Deviation standard de l'orographie sous-maille
	415	ALLOCATE(zsig(iml,jml)) ! Pente de l'orographie sous-maille
	416	ALLOCATE(zgam(iml,jml)) ! Anisotropie de l'orographie sous maille
	417	ALLOCATE(zthe(iml,jml)) ! Orientation axe +grande pente d'oro sous maille
	418	ALLOCATE(zpic(iml,jml)) ! Hauteur pics de la SSO
	419	ALLOCATE(zval(iml,jml)) ! Hauteur vallees de la SSO
	420	ALLOCATE(relief(iml,jml)) ! Orographie moyenne
	421
	422	CALL grid_noro(iml_rel, jml_rel, lon_rad, lat_rad, relief_hi, iml-1, jml, &
[808]	423	lon_in, lat_in, phis, relief, zstd, zsig, zgam, zthe, zpic, zval)
	424	phis = phis * g
[1]	425
	426	!--- SURFACE ROUGHNESS COMPUTATION (UNUSED FOR THE MOMENT !!! )
	427	IF(check) WRITE(lunout,*)'Compute surface roughness induced by the orography'
	428	ALLOCATE(rugo (iml ,jml))
	429	ALLOCATE(tmp_var(iml-1,jml))
	430	CALL rugsoro(iml_rel, jml_rel, lon_rad, lat_rad, relief_hi, iml-1, jml, &
	431	lon_in, lat_in, tmp_var)
	432	rugo(1:iml-1,:)=tmp_var; rugo(iml,:)=tmp_var(1,:)
	433	DEALLOCATE(relief_hi,tmp_var,lon_rad,lat_rad)
	434	RETURN
	435
	436	END SUBROUTINE start_init_orog
	437	!
	438	!-------------------------------------------------------------------------------
	439
	440
	441	!-------------------------------------------------------------------------------
	442	!
	443	SUBROUTINE start_init_phys(iml,jml,lon_in,lat_in,jml2,lon_in2,lat_in2,ibar)
	444	!
	445	!-------------------------------------------------------------------------------
	446	! Arguments:
	447	INTEGER, INTENT(IN) :: iml, jml
	448	REAL, DIMENSION(iml), INTENT(IN) :: lon_in
	449	REAL, DIMENSION(jml), INTENT(IN) :: lat_in
	450	INTEGER, INTENT(IN) :: jml2
	451	REAL, DIMENSION(iml), INTENT(IN) :: lon_in2
	452	REAL, DIMENSION(jml2), INTENT(IN) :: lat_in2
	453	LOGICAL, INTENT(IN) :: ibar
	454	!-------------------------------------------------------------------------------
	455	! Local variables:
	456	#include "iniprint.h"
	457	CHARACTER(LEN=25) :: title
	458	CHARACTER(LEN=120) :: physfname
	459	LOGICAL :: check=.TRUE.
	460	REAL :: date, dt
	461	INTEGER, DIMENSION(1) :: itau
	462	INTEGER :: llm_tmp, ttm_tmp
	463	REAL, DIMENSION(:,:), ALLOCATABLE :: var_ana
	464	REAL, DIMENSION(:), ALLOCATABLE :: lon_rad, lat_rad, lon_ini, lat_ini
	465	REAL, DIMENSION(:), ALLOCATABLE :: levphys_ini
	466	!-------------------------------------------------------------------------------
	467	physfname = 'ECPHY.nc'; pi=2.0*ASIN(1.0); deg2rad=pi/180.0
	468	IF(check) WRITE(lunout,*)'Opening the surface analysis'
	469	CALL flininfo(physfname, iml_phys, jml_phys, llm_tmp, ttm_tmp, fid_phys)
	470
	471	ALLOCATE(lat_phys(iml_phys,jml_phys))
	472	ALLOCATE(lon_phys(iml_phys,jml_phys))
	473	ALLOCATE(levphys_ini(llm_tmp))
	474	CALL flinopen(physfname, .FALSE., iml_phys, jml_phys, llm_tmp, lon_phys, &
	475	lat_phys, levphys_ini, ttm_tmp, itau, date, dt, fid_phys)
	476	DEALLOCATE(levphys_ini)
	477
	478	!--- IF ANGLES ARE IN DEGREES, THEY ARE CONVERTED INTO RADIANS
	479	ALLOCATE(lon_ini(iml_phys),lat_ini(jml_phys))
	480	lon_ini(:)=lon_phys(:,1); IF(MAXVAL(lon_phys)>pi) lon_ini=lon_ini*deg2rad
	481	lat_ini(:)=lat_phys(1,:); IF(MAXVAL(lat_phys)>pi) lat_ini=lat_ini*deg2rad
	482
	483	ALLOCATE(var_ana(iml_phys,jml_phys),lon_rad(iml_phys),lat_rad(jml_phys))
	484
	485	!--- SURFACE TEMPERATURE
	486	title='ST'
	487	ALLOCATE(tsol(iml,jml))
	488	CALL flinget(fid_phys,title,iml_phys,jml_phys,llm_tmp,ttm_tmp,1,1,var_ana)
	489	CALL conf_dat2d(title,iml_phys, jml_phys, lon_ini, lat_ini, lon_rad, lat_rad,&
	490	var_ana , ibar )
	491	CALL interp_startvar(title, ibar, .TRUE., &
	492	iml_phys, jml_phys, lon_rad, lat_rad, var_ana, iml, jml, jml-1, &
	493	lon_in, lat_in, lon_in2, lat_in2, tsol)
	494
	495	!--- SOIL MOISTURE
	496	title='CDSW'
	497	ALLOCATE(qsol(iml,jml))
	498	CALL flinget(fid_phys,title,iml_phys,jml_phys,llm_tmp,ttm_tmp,1,1,var_ana)
	499	CALL conf_dat2d(title,iml_phys, jml_phys, lon_ini, lat_ini, lon_rad, lat_rad,&
	500	var_ana, ibar )
	501	CALL interp_startvar(title, ibar, .TRUE., &
	502	iml_phys, jml_phys, lon_rad, lat_rad, var_ana, iml, jml, jml-1, &
	503	lon_in, lat_in, lon_in2, lat_in2, qsol)
	504
	505	CALL flinclo(fid_phys)
	506
	507	DEALLOCATE(var_ana,lon_rad,lat_rad,lon_ini,lat_ini)
	508
	509	END SUBROUTINE start_init_phys
	510	!
	511	!-------------------------------------------------------------------------------
	512
	513
	514	!-------------------------------------------------------------------------------
	515	!
	516	SUBROUTINE start_init_dyn(iml,jml,lon_in,lat_in,jml2,lon_in2,lat_in2,ibar)
	517	!
	518	!-------------------------------------------------------------------------------
	519	! Arguments:
	520	INTEGER, INTENT(IN) :: iml, jml
	521	REAL, DIMENSION(iml), INTENT(IN) :: lon_in
	522	REAL, DIMENSION(jml), INTENT(IN) :: lat_in
	523	INTEGER, INTENT(IN) :: jml2
	524	REAL, DIMENSION(iml), INTENT(IN) :: lon_in2
	525	REAL, DIMENSION(jml2), INTENT(IN) :: lat_in2
	526	LOGICAL, INTENT(IN) :: ibar
	527	!-------------------------------------------------------------------------------
	528	! Local variables:
	529	#include "iniprint.h"
	530	#include "dimensions.h"
	531	#include "paramet.h"
	532	#include "comgeom2.h"
	533	CHARACTER(LEN=25) :: title
	534	CHARACTER(LEN=120) :: physfname
	535	LOGICAL :: check=.TRUE.
	536	REAL :: date, dt
	537	INTEGER, DIMENSION(1) :: itau
	538	INTEGER :: i, j
	539	REAL, DIMENSION(:,:), ALLOCATABLE :: var_ana, z
	540	REAL, DIMENSION(:), ALLOCATABLE :: lon_rad, lat_rad, lon_ini, lat_ini
	541	REAL, DIMENSION(:), ALLOCATABLE :: xppn, xpps
	542	!-------------------------------------------------------------------------------
	543
	544	!--- KINETIC ENERGY
	545	physfname = 'ECDYN.nc'; pi=2.0*ASIN(1.0); deg2rad=pi/180.0
	546	IF(check) WRITE(lunout,*) 'Opening the surface analysis'
	547	CALL flininfo(physfname, iml_dyn, jml_dyn, llm_dyn, ttm_dyn, fid_dyn)
	548	IF(check) WRITE(lunout,*) 'Values read: ', iml_dyn, jml_dyn, llm_dyn, ttm_dyn
	549
	550	ALLOCATE(lat_dyn(iml_dyn,jml_dyn))
	551	ALLOCATE(lon_dyn(iml_dyn,jml_dyn))
	552	ALLOCATE(levdyn_ini(llm_dyn))
	553	CALL flinopen(physfname, .FALSE., iml_dyn, jml_dyn, llm_dyn, lon_dyn,lat_dyn,&
	554	levdyn_ini, ttm_dyn, itau, date, dt, fid_dyn)
	555
	556	!--- IF ANGLES ARE IN DEGREES, THEY ARE CONVERTED INTO RADIANS
	557	ALLOCATE(lon_ini(iml_dyn),lat_ini(jml_dyn))
	558	lon_ini(:)=lon_dyn(:,1); IF(MAXVAL(lon_dyn)>pi) lon_ini=lon_ini*deg2rad
	559	lat_ini(:)=lat_dyn(1,:); IF(MAXVAL(lat_dyn)>pi) lat_ini=lat_ini*deg2rad
	560
	561	ALLOCATE(var_ana(iml_dyn,jml_dyn),lon_rad(iml_dyn),lat_rad(jml_dyn))
	562
	563	!--- SURFACE GEOPOTENTIAL
	564	title='Z'
	565	ALLOCATE(z(iml,jml))
	566	CALL flinget(fid_dyn, title, iml_dyn, jml_dyn, 0, ttm_dyn, 1, 1, var_ana)
	567	CALL conf_dat2d(title, iml_dyn, jml_dyn, lon_ini, lat_ini, lon_rad, lat_rad, &
	568	var_ana, ibar)
	569	CALL interp_startvar(title, ibar, .TRUE., &
	570	iml_dyn, jml_dyn, lon_rad, lat_rad, var_ana, iml, jml, jml-1, &
	571	lon_in, lat_in, lon_in2, lat_in2, z)
	572
	573	!--- SURFACE PRESSURE
	574	title='SP'
	575	ALLOCATE(psol_dyn(iml,jml))
	576	CALL flinget(fid_dyn, title, iml_dyn, jml_dyn, 0, ttm_dyn, 1, 1, var_ana)
	577	CALL conf_dat2d(title, iml_dyn, jml_dyn, lon_ini, lat_ini, lon_rad, lat_rad, &
	578	var_ana, ibar)
	579	CALL interp_startvar(title, ibar, .TRUE., &
	580	iml_dyn, jml_dyn, lon_rad, lat_rad, var_ana, iml, jml, jml-1, &
	581	lon_in, lat_in, lon_in2, lat_in2, psol_dyn)
	582
	583	DEALLOCATE(var_ana,lon_rad,lat_rad,lon_ini,lat_ini)
	584
	585	!--- ALLOCATION OF VARIABLES CREATED IN OR COMING FROM RESTART FILE
	586	IF(.NOT.ALLOCATED(tsol)) THEN
	587	CALL start_init_phys(iml,jml,lon_in,lat_in,jml2,lon_in2,lat_in2,ibar)
	588	ELSE
	589	IF(SIZE(tsol)/=SIZE(psol_dyn)) THEN
	590	WRITE(lunout,*)'start_init_dyn :'
	591	WRITE(lunout,*)'The temperature field we have does not have the right size'
	592	STOP
	593	END IF
	594	END IF
	595
	596	IF(.NOT.ALLOCATED(phis)) THEN
	597	CALL start_init_orog(iml,jml,lon_in,lat_in)
	598	ELSE
	599	IF(SIZE(phis)/=SIZE(psol_dyn)) THEN
	600	WRITE(lunout,*)'start_init_dyn :'
	601	WRITE(lunout,*)'The orography field we have does not have the right size'
	602	STOP
	603	END IF
	604	END IF
	605
	606	!--- PSOL IS COMPUTED IN PASCALS
	607	DO j = 1, jml
	608	DO i = 1, iml-1
	609	psol_dyn(i,j) = psol_dyn(i,j)*(1.0+(z(i,j)-phis(i,j))/287.0/tsol(i,j))
	610	END DO
	611	psol_dyn(iml,j) = psol_dyn(1,j)
	612	END DO
	613	DEALLOCATE(z)
	614
	615	ALLOCATE(xppn(iml-1),xpps(iml-1))
	616	DO i = 1, iml-1
	617	xppn(i) = aire( i,1) * psol_dyn( i,1)
	618	xpps(i) = aire( i,jml) * psol_dyn( i,jml)
	619	END DO
	620	DO i = 1, iml
	621	psol_dyn(i,1 ) = SUM(xppn)/apoln
	622	psol_dyn(i,jml) = SUM(xpps)/apols
	623	END DO
	624	DEALLOCATE(xppn,xpps)
	625
	626	RETURN
	627
	628	END SUBROUTINE start_init_dyn
	629	!
	630	!-------------------------------------------------------------------------------
	631
	632
	633	!-------------------------------------------------------------------------------
	634	!
	635	SUBROUTINE start_inter_3d(varname, iml, jml, lml, lon_in, lat_in, jml2, &
	636	lon_in2, lat_in2, pls_in, var3d, ibar)
	637
	638	use pchsp_95_m, only: pchsp_95
	639	use pchfe_95_m, only: pchfe_95
	640
	641	! Arguments:
	642	CHARACTER(LEN=*), INTENT(IN) :: varname
	643	INTEGER, INTENT(IN) :: iml, jml, lml
	644	REAL, DIMENSION(iml), INTENT(IN) :: lon_in
	645	REAL, DIMENSION(jml), INTENT(IN) :: lat_in
	646	INTEGER, INTENT(IN) :: jml2
	647	REAL, DIMENSION(iml), INTENT(IN) :: lon_in2
	648	REAL, DIMENSION(jml2), INTENT(IN) :: lat_in2
	649	REAL, DIMENSION(iml, jml, lml), INTENT(IN) :: pls_in
	650	REAL, DIMENSION(iml, jml, lml), INTENT(OUT) :: var3d
	651	LOGICAL, INTENT(IN) :: ibar
	652	!----------------------------------------------------------------------------
	653	! Local variables:
	654	#include "iniprint.h"
	655	LOGICAL:: check=.TRUE., skip
	656	REAL chmin, chmax
	657	INTEGER ii, ij, il, ierr
	658	integer n_extrap ! number of extrapolated points
	659	REAL, DIMENSION(iml, jml, llm_dyn):: var_tmp3d
	660	REAL, DIMENSION(:), ALLOCATABLE :: lon_rad, lat_rad, lon_ini, lat_ini
	661	REAL, DIMENSION(llm_dyn):: lev_dyn, ax, ay, yder
	662
	663	!---------------------------------------------------------------------------
	664	IF(check) WRITE(lunout, *)'Going into flinget to extract the 3D field.'
	665	IF(check) WRITE(lunout, *) fid_dyn, varname, iml_dyn, jml_dyn, llm_dyn, &
	666	ttm_dyn
	667	IF(check) WRITE(lunout, *) 'Allocating space for interpolation', iml, jml, &
	668	llm_dyn
	669
	670	IF(.NOT.ALLOCATED(var_ana3d)) ALLOCATE(var_ana3d(iml_dyn, jml_dyn, llm_dyn))
	671	CALL flinget(fid_dyn, varname, iml_dyn, jml_dyn, llm_dyn, ttm_dyn, 1, 1, &
	672	var_ana3d)
	673
	674	!--- IF ANGLES ARE IN DEGREES, THEY ARE CONVERTED INTO RADIANS
	675	ALLOCATE(lon_ini(iml_dyn), lat_ini(jml_dyn))
	676	lon_ini(:)=lon_dyn(:, 1); IF(MAXVAL(lon_dyn)>pi) lon_ini=lon_ini*deg2rad
	677	lat_ini(:)=lat_dyn(1, :); IF(MAXVAL(lat_dyn)>pi) lat_ini=lat_ini*deg2rad
	678
	679	!--- FIELDS ARE PROCESSED TO BE ON STANDARD ANGULAR DOMAINS
	680	ALLOCATE(lon_rad(iml_dyn), lat_rad(jml_dyn))
	681	CALL conf_dat3d (varname, iml_dyn, jml_dyn, llm_dyn, lon_ini, lat_ini, &
	682	levdyn_ini, lon_rad, lat_rad, lev_dyn, var_ana3d, ibar)
	683	DEALLOCATE(lon_ini, lat_ini)
	684
	685	!--- COMPUTING THE REQUIRED FIELDS USING ROUTINE grid_noro
	686	DO il=1, llm_dyn
	687	CALL interp_startvar(varname, ibar, il==1, iml_dyn, jml_dyn, lon_rad, &
	688	lat_rad, var_ana3d(:, :, il), iml, jml, jml2, lon_in, lat_in, &
	689	lon_in2, lat_in2, var_tmp3d(:, :, il))
	690	END DO
	691	DEALLOCATE(lon_rad, lat_rad)
	692
	693	!--- VERTICAL INTERPOLATION IS PERFORMED FROM TOP OF ATMOSPHERE TO GROUND
	694	ax = lev_dyn(llm_dyn:1:-1)
	695	skip = .false.
	696	n_extrap = 0
	697	DO ij=1, jml
	698	DO ii=1, iml-1
	699	ay = var_tmp3d(ii, ij, llm_dyn:1:-1)
	700	yder = pchsp_95(ax, ay, ibeg=2, iend=2, vc_beg=0., vc_end=0.)
	701	CALL pchfe_95(ax, ay, yder, skip, pls_in(ii, ij, lml:1:-1), &
	702	var3d(ii, ij, lml:1:-1), ierr)
	703	if (ierr < 0) stop 1
	704	n_extrap = n_extrap + ierr
	705	END DO
	706	END DO
	707	if (n_extrap /= 0) then
	708	print *, "start_inter_3d pchfe_95: n_extrap = ", n_extrap
	709	end if
	710	var3d(iml, :, :) = var3d(1, :, :)
	711
	712	DO il=1, lml
	713	CALL minmax(iml*jml, var3d(1, 1, il), chmin, chmax)
	714	WRITE(lunout, *)' '//TRIM(varname)//' min max l ', il, chmin, chmax
	715	END DO
	716
	717	END SUBROUTINE start_inter_3d
	718	!
	719	!-------------------------------------------------------------------------------
	720
	721
	722	!-------------------------------------------------------------------------------
	723	!
	724	SUBROUTINE interp_startvar(vname, ibar, ibeg, ii, jj, lon, lat, vari, &
	725	i1, j1, j2, lon1, lat1, lon2, lat2, varo)
	726	!
	727	!-------------------------------------------------------------------------------
	728
	729	USE inter_barxy_m, only: inter_barxy
	730
	731	! Arguments:
	732	CHARACTER(LEN=*), INTENT(IN) :: vname
	733	LOGICAL, INTENT(IN) :: ibar, ibeg
	734	INTEGER, INTENT(IN) :: ii, jj
	735	REAL, DIMENSION(ii), INTENT(IN) :: lon
	736	REAL, DIMENSION(jj), INTENT(IN) :: lat
	737	REAL, DIMENSION(ii,jj), INTENT(IN) :: vari
	738	INTEGER, INTENT(IN) :: i1, j1, j2
	739	REAL, DIMENSION(i1), INTENT(IN) :: lon1
	740	REAL, DIMENSION(j1), INTENT(IN) :: lat1
	741	REAL, DIMENSION(i1), INTENT(IN) :: lon2
	742	REAL, DIMENSION(j2), INTENT(IN) :: lat2
	743	REAL, DIMENSION(i1,j1), INTENT(OUT) :: varo
	744	!-------------------------------------------------------------------------------
	745	! Local variables:
	746	#include "iniprint.h"
	747	REAL, DIMENSION(i1-1,j1) :: vtmp
	748	!-------------------------------------------------------------------------------
	749	IF(ibar) THEN
	750	IF(ibeg) THEN
	751	WRITE(lunout,*) &
	752	'---------------------------------------------------------------'
	753	WRITE(lunout,*) &
	754	'$$$ Utilisation de l interpolation barycentrique pour '//TRIM(vname)//' $$$'
	755	WRITE(lunout,*) &
	756	'---------------------------------------------------------------'
	757	END IF
	758	CALL inter_barxy(lon, lat(:jj-1), vari, lon2(:i1-1), lat2(:j2), vtmp)
	759	ELSE
	760	CALL grille_m (ii, jj, lon, lat, vari, i1-1, j1, lon1, lat1, vtmp)
	761	END IF
	762	CALL gr_int_dyn(vtmp, varo, i1-1, j1)
	763
	764	END SUBROUTINE interp_startvar
	765	!
	766	!-------------------------------------------------------------------------------
	767
	768	END MODULE startvar
	769	!
	770	!*******************************************************************************

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

Download in other formats:

Original Format