Context Navigation

startvar.F90 @ 2660

Last change on this file since 2660 was 1425, checked in by lguez, 14 years ago

Replaced Numerical Recipes procedures for spline interpolation (not in
the public domain) by procedures from the Pchip package in the Slatec
library. This only affects the program "ce0l", not the program
"gcm". Tested on Brodie SX8 with "-debug" and "-prod", "-parallel
none" and "-parallel mpi". "start.nc" and "limit.nc" are
changed. "startphy.nc" is not changed. The relative change is of order
1e-7 or less. The revision makes the program faster (tested on Brodie
with "-prod -d 144x142x39", CPU time is 38 s, instead of 54
s). Procedures from Slatec are untouched, except for
"i1mach.F". Created procedures "pchfe_95" and "pchsp_95" which are
wrappers for "pchfe" and "pchsp" from Slatec. "pchfe_95" and
"pchsp_95" have a safer and simpler interface.

Replaced "make" by "sxgmake" in "arch-SX8_BRODIE.fcm". Added files for
compilation by FCM with "g95".

In "arch-linux-32bit.fcm", replaced "pgf90" by "pgf95". There was no
difference between "dev" and "debug" so added "-O1" to "dev". Added
debugging options. Removed "-Wl,-Bstatic
-L/usr/lib/gcc-lib/i386-linux/2.95.2", which usually produces an error
at link-time.

Bash is now ubiquitous while KornShell? is not so use Bash instead of
KornShell? in FCM.

Replaced some statements "write(6,*)" by "write(lunout,*)". Replaced
"stop" by "stop 1" in the case where "abort_gcm" is called with "ierr
/= 0". Removed "stop" statements at the end of procedures
"limit_netcdf" and main program "ce0l" (why not let the program end
normally?).

Made some arrays automatic instead of allocatable in "start_inter_3d".

Zeroed "wake_pe", "fm_therm", "entr_therm" and "detr_therm" in
"dyn3dpar/etat0_netcdf.F90". The parallel and sequential results of
"ce0l" are thus identical.

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 1425 2010-09-02 13:45:23Z musat $
[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