Context Navigation

startvar.F90 @ 1950

Last change on this file since 1950 was 1907, checked in by lguez, 11 years ago

Added a copyright property to every file of the distribution, except
for the fcm files (which have their own copyright). Use svn propget on
a file to see the copyright. For instance:

$ svn propget copyright libf/phylmd/physiq.F90
Name of program: LMDZ
Creation date: 1984
Version: LMDZ5
License: CeCILL version 2
Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
See the license file in the root directory

Also added the files defining the CeCILL version 2 license, in French
and English, at the top of the LMDZ tree.

Property copyright set to
Name of program: LMDZ
Creation date: 1984
Version: LMDZ5
License: CeCILL version 2
Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
See the license file in the root directory
Property svn:eol-style set to native
Property svn:keywords set to Author Date Id Revision

File size: 31.2 KB

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

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

Download in other formats:

Original Format