Context Navigation

limit_netcdf.F90 @ 1523

Last change on this file since 1523 was 1441, checked in by Ehouarn Millour, 14 years ago

minor fix: enable routine to be compilable when compiling without Earth physics.

Property svn:eol-style set to native
Property svn:keywords set to Author Date Id Revision

File size: 26.0 KB

Rev	Line
[1000]	1	!
[1279]	2	! $Id: limit_netcdf.F90 1441 2010-10-13 13:06:56Z idelkadi $
[1319]	3	!-------------------------------------------------------------------------------
[1000]	4	!
[1319]	5	SUBROUTINE limit_netcdf(interbar, extrap, oldice, masque)
	6	!
	7	!-------------------------------------------------------------------------------
	8	! Author : L. Fairhead, 27/01/94
	9	!-------------------------------------------------------------------------------
	10	! Purpose: Boundary conditions files building for new model using climatologies.
	11	! Both grids have to be regular.
	12	!-------------------------------------------------------------------------------
	13	! Note: This routine is designed to work for Earth
	14	!-------------------------------------------------------------------------------
	15	! Modification history:
	16	! * 23/03/1994: Z. X. Li
	17	! * 09/1999: L. Fairhead (netcdf reading in LMDZ.3.3)
	18	! * 07/2001: P. Le Van
	19	! * 11/2009: L. Guez (ozone day & night climatos, see etat0_netcdf.F90)
	20	! * 12/2009: D. Cugnet (f77->f90, calendars, files from coupled runs)
	21	!-------------------------------------------------------------------------------
[1425]	22	USE control_mod
[1279]	23	#ifdef CPP_EARTH
[1319]	24	USE dimphy
	25	USE ioipsl, ONLY : ioget_year_len
	26	USE phys_state_var_mod, ONLY : pctsrf
	27	USE netcdf, ONLY : NF90_OPEN, NF90_CREATE, NF90_CLOSE, &
	28	NF90_DEF_DIM, NF90_DEF_VAR, NF90_PUT_VAR, NF90_PUT_ATT, &
	29	NF90_NOERR, NF90_NOWRITE, NF90_DOUBLE, NF90_GLOBAL, &
[1425]	30	NF90_CLOBBER, NF90_ENDDEF, NF90_UNLIMITED, NF90_FLOAT
[1323]	31	USE inter_barxy_m, only: inter_barxy
[1319]	32	#endif
	33	IMPLICIT NONE
	34	!-------------------------------------------------------------------------------
	35	! Arguments:
[1000]	36	#include "dimensions.h"
	37	#include "paramet.h"
[1319]	38	#include "iniprint.h"
	39	LOGICAL, INTENT(IN) :: interbar ! barycentric interpolation
	40	LOGICAL, INTENT(IN) :: extrap ! SST extrapolation flag
	41	LOGICAL, INTENT(IN) :: oldice ! old way ice computation
	42	REAL, DIMENSION(iip1,jjp1), INTENT(IN) :: masque ! land mask
	43	#ifndef CPP_EARTH
	44	WRITE(lunout,*)'limit_netcdf: Earth-specific routine, needs Earth physics'
	45	#else
	46	!-------------------------------------------------------------------------------
	47	! Local variables:
[1000]	48	#include "logic.h"
	49	#include "comvert.h"
	50	#include "comgeom2.h"
	51	#include "comconst.h"
	52	#include "indicesol.h"
	53
[1319]	54	!--- For fractionary sub-cell use (old coding used soil index: 0,1,2,3) --------
	55	LOGICAL, PARAMETER :: fracterre=.TRUE.
[1000]	56
[1319]	57	!--- INPUT NETCDF FILES NAMES --------------------------------------------------
	58	CHARACTER(LEN=25) :: icefile, sstfile, dumstr
	59	CHARACTER(LEN=25), PARAMETER :: famipsst='amipbc_sst_1x1.nc ', &
	60	famipsic='amipbc_sic_1x1.nc ', &
	61	fclimsst='amipbc_sst_1x1_clim.nc ', &
	62	fclimsic='amipbc_sic_1x1_clim.nc ', &
	63	fcpldsst='cpl_atm_sst.nc ', &
	64	fcpldsic='cpl_atm_sic.nc ', &
	65	frugo ='Rugos.nc ', &
	66	falbe ='Albedo.nc '
[1000]	67
[1319]	68	!--- OUTPUT VARIABLES FOR NETCDF FILE ------------------------------------------
	69	REAL, DIMENSION(klon) :: fi_ice, verif
	70	REAL, DIMENSION(:,:), POINTER :: phy_rug=>NULL(), phy_ice=>NULL()
	71	REAL, DIMENSION(:,:), POINTER :: phy_sst=>NULL(), phy_alb=>NULL()
	72	REAL, DIMENSION(:,:), ALLOCATABLE :: phy_bil
	73	REAL, DIMENSION(:,:,:), ALLOCATABLE :: pctsrf_t
	74	INTEGER :: nbad
[1000]	75
[1319]	76	!--- VARIABLES FOR OUTPUT FILE WRITING -----------------------------------------
	77	INTEGER :: ierr, nid, ndim, ntim, k
	78	INTEGER, DIMENSION(2) :: dims
	79	INTEGER :: id_tim, id_SST, id_BILS, id_RUG, id_ALB
	80	INTEGER :: id_FOCE, id_FSIC, id_FTER, id_FLIC
	81	INTEGER :: NF90_FORMAT
	82	LOGICAL :: lCPL !--- T: IPCC-IPSL cpl model output files
[1328]	83	INTEGER :: ndays !--- Depending on the output calendar
[1000]	84
[1319]	85	!--- INITIALIZATIONS -----------------------------------------------------------
	86	#ifdef NC_DOUBLE
	87	NF90_FORMAT=NF90_DOUBLE
	88	#else
	89	NF90_FORMAT=NF90_FLOAT
	90	#endif
[1000]	91
[1319]	92	pi = 4.*ATAN(1.)
	93	rad = 6371229.
	94	daysec= 86400.
	95	omeg = 2.*pi/daysec
	96	g = 9.8
	97	kappa = 0.2857143
	98	cpp = 1004.70885
[1441]	99	dtvr = daysec/REAL(day_step)
[1319]	100	CALL inigeom
[1000]	101
[1319]	102	!--- Beware: anneeref (from gcm.def) is used to determine output time sampling
	103	ndays=ioget_year_len(anneeref)
[1000]	104
[1319]	105	!--- RUGOSITY TREATMENT --------------------------------------------------------
	106	WRITE(lunout,*) 'Traitement de la rugosite'
	107	CALL get_2Dfield(frugo,'RUG',interbar,ndays,phy_rug,mask=masque(1:iim,:))
[1000]	108
[1319]	109	!--- OCEAN TREATMENT -----------------------------------------------------------
	110	PRINT*, 'Traitement de la glace oceanique' ; icefile=''; lCPL=.FALSE.
[1000]	111
[1319]	112	! Input SIC file selection
	113	icefile='fake'
	114	IF(NF90_OPEN(famipsic,NF90_NOWRITE,nid)==NF90_NOERR) icefile=TRIM(famipsic)
	115	IF(NF90_OPEN(fclimsic,NF90_NOWRITE,nid)==NF90_NOERR) icefile=TRIM(fclimsic)
	116	IF(NF90_OPEN(fcpldsic,NF90_NOWRITE,nid)==NF90_NOERR) icefile=TRIM(fcpldsic)
	117	SELECT CASE(icefile)
	118	CASE(famipsic); dumstr='Amip.'
	119	CASE(fclimsic); dumstr='Amip climatologique.'
	120	CASE(fcpldsic); dumstr='de sortie du modele couplé IPSL/IPCC.';lCPL=.TRUE.
	121	CASE('fake'); CALL abort_gcm('limit_netcdf','Fichier SIC non reconnu.',1)
	122	END SELECT
	123	ierr=NF90_CLOSE(nid)
	124	WRITE(lunout,*)'Pour la glace de mer a ete choisi un fichier '//TRIM(dumstr)
[1000]	125
[1328]	126	CALL get_2Dfield(icefile,'SIC',interbar,ndays,phy_ice,flag=oldice,lCPL=lCPL)
[1000]	127
[1319]	128	ALLOCATE(pctsrf_t(klon,nbsrf,ndays))
	129	DO k=1,ndays
	130	fi_ice=phy_ice(:,k)
	131	WHERE(fi_ice>=1.0 ) fi_ice=1.0
	132	WHERE(fi_ice<EPSFRA) fi_ice=0.0
	133	IF(fracterre) THEN
	134	pctsrf_t(:,is_ter,k)=pctsrf(:,is_ter) ! land soil
	135	pctsrf_t(:,is_lic,k)=pctsrf(:,is_lic) ! land ice
[1328]	136	IF(lCPL) THEN ! SIC=pICE*(1-LIC-TER)
[1319]	137	pctsrf_t(:,is_sic,k)=fi_ice*(1-pctsrf(:,is_lic)-pctsrf(:,is_ter))
	138	ELSE ! SIC=pICE-LIC
	139	pctsrf_t(:,is_sic,k)=fi_ice-pctsrf_t(:,is_lic,k)
	140	END IF
	141	WHERE(pctsrf_t(:,is_sic,k)<=0) pctsrf_t(:,is_sic,k)=0.
	142	WHERE(1.0-zmasq<EPSFRA)
	143	pctsrf_t(:,is_sic,k)=0.0
	144	pctsrf_t(:,is_oce,k)=0.0
	145	ELSEWHERE
	146	WHERE(pctsrf_t(:,is_sic,k)>=1.0-zmasq)
	147	pctsrf_t(:,is_sic,k)=1.0-zmasq
	148	pctsrf_t(:,is_oce,k)=0.0
	149	ELSEWHERE
	150	pctsrf_t(:,is_oce,k)=1.0-zmasq-pctsrf_t(:,is_sic,k)
	151	WHERE(pctsrf_t(:,is_oce,k)<EPSFRA)
	152	pctsrf_t(:,is_oce,k)=0.0
	153	pctsrf_t(:,is_sic,k)=1.0-zmasq
	154	END WHERE
	155	END WHERE
	156	END WHERE
	157	nbad=COUNT(pctsrf_t(:,is_oce,k)<0.0)
	158	IF(nbad>0) WRITE(lunout,*) 'pb sous maille pour nb point = ',nbad
	159	nbad=COUNT(abs(sum(pctsrf_t(:,:,k),dim=2)-1.0)>EPSFRA)
	160	IF(nbad>0) WRITE(lunout,*) 'pb sous surface pour nb points = ',nbad
	161	ELSE
	162	pctsrf_t(:,is_ter,k)=pctsrf(:,is_ter)
	163	WHERE(NINT(pctsrf(:,is_ter))==1)
	164	pctsrf_t(:,is_sic,k)=0.
	165	pctsrf_t(:,is_oce,k)=0.
	166	WHERE(fi_ice>=1.e-5)
	167	pctsrf_t(:,is_lic,k)=fi_ice
	168	pctsrf_t(:,is_ter,k)=pctsrf_t(:,is_ter,k)-pctsrf_t(:,is_lic,k)
	169	ELSEWHERE
	170	pctsrf_t(:,is_lic,k)=0.0
	171	END WHERE
	172	ELSEWHERE
	173	pctsrf_t(:,is_lic,k) = 0.0
	174	WHERE(fi_ice>=1.e-5)
	175	pctsrf_t(:,is_ter,k)=0.0
	176	pctsrf_t(:,is_sic,k)=fi_ice
	177	pctsrf_t(:,is_oce,k)=1.0-pctsrf_t(:,is_sic,k)
	178	ELSEWHERE
	179	pctsrf_t(:,is_sic,k)=0.0
	180	pctsrf_t(:,is_oce,k)=1.0
	181	END WHERE
	182	END WHERE
	183	verif=sum(pctsrf_t(:,:,k),dim=2)
	184	nbad=COUNT(verif<1.0-1.e-5.OR.verif>1.0+1.e-5)
	185	IF(nbad>0) WRITE(lunout,*) 'pb sous maille pour nb point = ',nbad
	186	END IF
	187	END DO
	188	DEALLOCATE(phy_ice)
[1000]	189
[1319]	190	!--- SST TREATMENT -------------------------------------------------------------
	191	WRITE(lunout,*) 'Traitement de la sst' ; sstfile=''; lCPL=.FALSE.
[1000]	192
[1319]	193	! Input SST file selection
	194	sstfile='fake'
	195	IF(NF90_OPEN(famipsst,NF90_NOWRITE,nid)==NF90_NOERR) sstfile=TRIM(famipsst)
	196	IF(NF90_OPEN(fclimsst,NF90_NOWRITE,nid)==NF90_NOERR) sstfile=TRIM(fclimsst)
	197	IF(NF90_OPEN(fcpldsst,NF90_NOWRITE,nid)==NF90_NOERR) sstfile=TRIM(fcpldsst)
	198	SELECT CASE(icefile)
	199	CASE(famipsic); dumstr='Amip.'
	200	CASE(fclimsic); dumstr='Amip climatologique.'
	201	CASE(fcpldsic); dumstr='de sortie du modele couplé IPSL/IPCC.';lCPL=.TRUE.
	202	CASE('fake'); CALL abort_gcm('limit_netcdf','Fichier SST non reconnu',1)
	203	END SELECT
	204	ierr=NF90_CLOSE(nid)
	205	WRITE(lunout,*)'Pour la temperature de mer a ete choisi un fichier '//TRIM(dumstr)
[1000]	206
[1319]	207	CALL get_2Dfield(trim(sstfile),'SST',interbar,ndays,phy_sst,flag=extrap)
[1000]	208
[1319]	209	!--- ALBEDO TREATMENT ----------------------------------------------------------
	210	WRITE(lunout,*) 'Traitement de l albedo'
	211	CALL get_2Dfield(falbe,'ALB',interbar,ndays,phy_alb)
[1000]	212
[1319]	213	!--- REFERENCE GROUND HEAT FLUX TREATMENT --------------------------------------
	214	ALLOCATE(phy_bil(klon,ndays)); phy_bil=0.0
[1000]	215
[1319]	216	!--- OUTPUT FILE WRITING -------------------------------------------------------
	217	WRITE(lunout,*) 'Ecriture du fichier limit'
[1000]	218
[1319]	219	!--- File creation
	220	ierr=NF90_CREATE("limit.nc",NF90_CLOBBER,nid)
	221	ierr=NF90_PUT_ATT(nid,NF90_GLOBAL,"title","Fichier conditions aux limites")
[1000]	222
[1319]	223	!--- Dimensions creation
	224	ierr=NF90_DEF_DIM(nid,"points_physiques",klon,ndim)
	225	ierr=NF90_DEF_DIM(nid,"time",NF90_UNLIMITED,ntim)
[1000]	226
[1319]	227	dims=(/ndim,ntim/)
[1000]	228
[1319]	229	!--- Variables creation
[1328]	230	ierr=NF90_DEF_VAR(nid,"TEMPS",NF90_FORMAT,(/ntim/),id_tim)
[1319]	231	ierr=NF90_DEF_VAR(nid,"FOCE", NF90_FORMAT,dims,id_FOCE)
	232	ierr=NF90_DEF_VAR(nid,"FSIC", NF90_FORMAT,dims,id_FSIC)
	233	ierr=NF90_DEF_VAR(nid,"FTER", NF90_FORMAT,dims,id_FTER)
	234	ierr=NF90_DEF_VAR(nid,"FLIC", NF90_FORMAT,dims,id_FLIC)
	235	ierr=NF90_DEF_VAR(nid,"SST", NF90_FORMAT,dims,id_SST)
	236	ierr=NF90_DEF_VAR(nid,"BILS", NF90_FORMAT,dims,id_BILS)
	237	ierr=NF90_DEF_VAR(nid,"ALB", NF90_FORMAT,dims,id_ALB)
	238	ierr=NF90_DEF_VAR(nid,"RUG", NF90_FORMAT,dims,id_RUG)
[1000]	239
[1319]	240	!--- Attributes creation
	241	ierr=NF90_PUT_ATT(nid,id_tim, "title","Jour dans l annee")
	242	ierr=NF90_PUT_ATT(nid,id_FOCE,"title","Fraction ocean")
	243	ierr=NF90_PUT_ATT(nid,id_FSIC,"title","Fraction glace de mer")
	244	ierr=NF90_PUT_ATT(nid,id_FTER,"title","Fraction terre")
	245	ierr=NF90_PUT_ATT(nid,id_FLIC,"title","Fraction land ice")
	246	ierr=NF90_PUT_ATT(nid,id_SST ,"title","Temperature superficielle de la mer")
	247	ierr=NF90_PUT_ATT(nid,id_BILS,"title","Reference flux de chaleur au sol")
	248	ierr=NF90_PUT_ATT(nid,id_ALB, "title","Albedo a la surface")
	249	ierr=NF90_PUT_ATT(nid,id_RUG, "title","Rugosite")
[1000]	250
[1319]	251	ierr=NF90_ENDDEF(nid)
	252
	253	!--- Variables saving
[1328]	254	ierr=NF90_PUT_VAR(nid,id_tim,(/(REAL(k),k=1,ndays)/))
[1319]	255	ierr=NF90_PUT_VAR(nid,id_FOCE,pctsrf_t(:,is_oce,:),(/1,1/),(/klon,ndays/))
	256	ierr=NF90_PUT_VAR(nid,id_FSIC,pctsrf_t(:,is_sic,:),(/1,1/),(/klon,ndays/))
	257	ierr=NF90_PUT_VAR(nid,id_FTER,pctsrf_t(:,is_ter,:),(/1,1/),(/klon,ndays/))
	258	ierr=NF90_PUT_VAR(nid,id_FLIC,pctsrf_t(:,is_lic,:),(/1,1/),(/klon,ndays/))
	259	ierr=NF90_PUT_VAR(nid,id_SST ,phy_sst(:,:),(/1,1/),(/klon,ndays/))
	260	ierr=NF90_PUT_VAR(nid,id_BILS,phy_bil(:,:),(/1,1/),(/klon,ndays/))
	261	ierr=NF90_PUT_VAR(nid,id_ALB ,phy_alb(:,:),(/1,1/),(/klon,ndays/))
	262	ierr=NF90_PUT_VAR(nid,id_RUG ,phy_rug(:,:),(/1,1/),(/klon,ndays/))
	263
	264	ierr=NF90_CLOSE(nid)
	265
	266	DEALLOCATE(pctsrf_t,phy_sst,phy_bil,phy_alb,phy_rug)
[1000]	267
	268
[1319]	269	!===============================================================================
	270	!
	271	CONTAINS
	272	!
	273	!===============================================================================
[1000]	274
	275
[1319]	276	!-------------------------------------------------------------------------------
	277	!
[1328]	278	SUBROUTINE get_2Dfield(fnam, mode, ibar, ndays, champo, flag, mask, lCPL)
[1319]	279	!
[1425]	280	!-----------------------------------------------------------------------------
[1319]	281	! Comments:
	282	! There are two assumptions concerning the NetCDF files, that are satisfied
	283	! with files that are conforming NC convention:
	284	! 1) The last dimension of the variables used is the time record.
	285	! 2) Dimensional variables have the same names as corresponding dimensions.
[1425]	286	!-----------------------------------------------------------------------------
	287	USE netcdf, ONLY: NF90_OPEN, NF90_INQ_VARID, NF90_INQUIRE_VARIABLE, &
	288	NF90_CLOSE, NF90_INQ_DIMID, NF90_INQUIRE_DIMENSION, NF90_GET_VAR, &
	289	NF90_GET_ATT
[1319]	290	USE dimphy, ONLY : klon
	291	USE phys_state_var_mod, ONLY : pctsrf
[1403]	292	USE control_mod
[1425]	293	use pchsp_95_m, only: pchsp_95
	294	use pchfe_95_m, only: pchfe_95
	295	use arth_m, only: arth
	296
[1319]	297	IMPLICIT NONE
	298	#include "dimensions.h"
	299	#include "paramet.h"
	300	#include "comgeom2.h"
	301	#include "indicesol.h"
	302	#include "iniprint.h"
[1425]	303	!-----------------------------------------------------------------------------
[1319]	304	! Arguments:
	305	CHARACTER(LEN=*), INTENT(IN) :: fnam ! NetCDF file name
	306	CHARACTER(LEN=3), INTENT(IN) :: mode ! RUG, SIC, SST or ALB
	307	LOGICAL, INTENT(IN) :: ibar ! interp on pressure levels
	308	INTEGER, INTENT(IN) :: ndays ! current year number of days
[1425]	309	REAL, POINTER, DIMENSION(:, :) :: champo ! output field = f(t)
	310	LOGICAL, OPTIONAL, INTENT(IN) :: flag ! extrapol. (SST) old ice (SIC)
	311	REAL, OPTIONAL, DIMENSION(iim, jjp1), INTENT(IN) :: mask
[1328]	312	LOGICAL, OPTIONAL, INTENT(IN) :: lCPL ! Coupled model flag (for ICE)
[1425]	313	!------------------------------------------------------------------------------
[1319]	314	! Local variables:
	315	!--- NetCDF
[1425]	316	INTEGER :: ncid, varid ! NetCDF identifiers
[1319]	317	CHARACTER(LEN=30) :: dnam ! dimension name
	318	CHARACTER(LEN=80) :: varname ! NetCDF variable name
	319	!--- dimensions
	320	INTEGER, DIMENSION(4) :: dids ! NetCDF dimensions identifiers
	321	REAL, ALLOCATABLE, DIMENSION(:) :: dlon_ini ! initial longitudes vector
	322	REAL, ALLOCATABLE, DIMENSION(:) :: dlat_ini ! initial latitudes vector
	323	REAL, POINTER, DIMENSION(:) :: dlon, dlat ! reordered lon/lat vectors
	324	!--- fields
	325	INTEGER :: imdep, jmdep, lmdep ! dimensions of 'champ'
[1425]	326	REAL, ALLOCATABLE, DIMENSION(:, :) :: champ ! wanted field on initial grid
[1319]	327	REAL, ALLOCATABLE, DIMENSION(:) :: yder, timeyear
[1425]	328	REAL, DIMENSION(iim, jjp1) :: champint ! interpolated field
	329	REAL, ALLOCATABLE, DIMENSION(:, :, :) :: champtime
	330	REAL, ALLOCATABLE, DIMENSION(:, :, :) :: champan
[1319]	331	!--- input files
	332	CHARACTER(LEN=20) :: cal_in ! calendar
	333	INTEGER :: ndays_in ! number of days
	334	!--- misc
	335	INTEGER :: i, j, k, l ! loop counters
[1425]	336	REAL, ALLOCATABLE, DIMENSION(:, :) :: work ! used for extrapolation
[1319]	337	CHARACTER(LEN=25) :: title ! for messages
	338	LOGICAL :: extrp ! flag for extrapolation
	339	REAL :: chmin, chmax
[1425]	340	INTEGER ierr
	341	integer n_extrap ! number of extrapolated points
	342	logical skip
	343	!------------------------------------------------------------------------------
	344	!---Variables depending on keyword 'mode' -------------------------------------
[1319]	345	NULLIFY(champo)
	346	SELECT CASE(mode)
	347	CASE('RUG'); varname='RUGOS'; title='Rugosite'
	348	CASE('SIC'); varname='sicbcs'; title='Sea-ice'
	349	CASE('SST'); varname='tosbcs'; title='SST'
	350	CASE('ALB'); varname='ALBEDO'; title='Albedo'
	351	END SELECT
[1425]	352	extrp=.FALSE.
	353	IF ( PRESENT(flag) ) THEN
	354	IF ( flag .AND. mode=='SST' ) extrp=.TRUE.
	355	END IF
[1000]	356
[1425]	357	!--- GETTING SOME DIMENSIONAL VARIABLES FROM FILE -----------------------------
	358	ierr=NF90_OPEN(fnam, NF90_NOWRITE, ncid); CALL ncerr(ierr, fnam)
	359	ierr=NF90_INQ_VARID(ncid, varname, varid); CALL ncerr(ierr, fnam)
	360	ierr=NF90_INQUIRE_VARIABLE(ncid, varid, dimids=dids); CALL ncerr(ierr, fnam)
[1404]	361
[1319]	362	!--- Longitude
[1425]	363	ierr=NF90_INQUIRE_DIMENSION(ncid, dids(1), name=dnam, len=imdep)
	364	CALL ncerr(ierr, fnam); ALLOCATE(dlon_ini(imdep), dlon(imdep))
	365	ierr=NF90_INQ_VARID(ncid, dnam, varid); CALL ncerr(ierr, fnam)
	366	ierr=NF90_GET_VAR(ncid, varid, dlon_ini); CALL ncerr(ierr, fnam)
	367	WRITE(lunout, *) 'variable ', dnam, 'dimension ', imdep
[1000]	368
[1319]	369	!--- Latitude
[1425]	370	ierr=NF90_INQUIRE_DIMENSION(ncid, dids(2), name=dnam, len=jmdep)
	371	CALL ncerr(ierr, fnam); ALLOCATE(dlat_ini(jmdep), dlat(jmdep))
	372	ierr=NF90_INQ_VARID(ncid, dnam, varid); CALL ncerr(ierr, fnam)
	373	ierr=NF90_GET_VAR(ncid, varid, dlat_ini); CALL ncerr(ierr, fnam)
	374	WRITE(lunout, *) 'variable ', dnam, 'dimension ', jmdep
[1000]	375
[1319]	376	!--- Time (variable is not needed - it is rebuilt - but calendar is)
[1425]	377	ierr=NF90_INQUIRE_DIMENSION(ncid, dids(3), name=dnam, len=lmdep)
	378	CALL ncerr(ierr, fnam); ALLOCATE(timeyear(lmdep))
	379	ierr=NF90_INQ_VARID(ncid, dnam, varid); CALL ncerr(ierr, fnam)
[1319]	380	cal_in=' '
[1425]	381	ierr=NF90_GET_ATT(ncid, varid, 'calendar', cal_in)
[1319]	382	IF(ierr/=NF90_NOERR) THEN
	383	SELECT CASE(mode)
[1425]	384	CASE('RUG', 'ALB'); cal_in='360d'
	385	CASE('SIC', 'SST'); cal_in='gregorian'
[1319]	386	END SELECT
[1425]	387	WRITE(lunout, *)'ATTENTION: variable "time" sans attribut "calendrier" ' &
	388	// 'dans '//TRIM(fnam)//'. On choisit la valeur par defaut.'
[1319]	389	END IF
[1425]	390	WRITE(lunout, *) 'variable ', dnam, 'dimension ', lmdep, 'calendrier ', &
	391	cal_in
[1000]	392
[1425]	393	!--- CONSTRUCTING THE INPUT TIME VECTOR FOR INTERPOLATION --------------------
[1319]	394	!--- Determining input file number of days, depending on calendar
[1425]	395	ndays_in=year_len(anneeref, cal_in)
[1000]	396
[1319]	397	!--- Time vector reconstruction (time vector from file is not trusted)
	398	!--- If input records are not monthly, time sampling has to be constant !
[1425]	399	timeyear=mid_months(anneeref, cal_in, lmdep)
	400	IF (lmdep /= 12) WRITE(lunout, '(a, i3, a)') 'Note : les fichiers de ' &
	401	// TRIM(mode) // ' ne comportent pas 12, mais ', lmdep, &
	402	' enregistrements.'
[1000]	403
[1425]	404	!--- GETTING THE FIELD AND INTERPOLATING IT ----------------------------------
	405	ALLOCATE(champ(imdep, jmdep), champtime(iim, jjp1, lmdep))
	406	IF(extrp) ALLOCATE(work(imdep, jmdep))
[1000]	407
[1425]	408	WRITE(lunout, *)
	409	WRITE(lunout, '(a, i3, a)')'LECTURE ET INTERPOLATION HORIZ. DE ', lmdep, &
	410	' CHAMPS.'
	411	ierr=NF90_INQ_VARID(ncid, varname, varid); CALL ncerr(ierr, fnam)
	412	DO l=1, lmdep
	413	ierr=NF90_GET_VAR(ncid, varid, champ, (/1, 1, l/), (/imdep, jmdep, 1/))
	414	CALL ncerr(ierr, fnam)
	415	CALL conf_dat2d(title, imdep, jmdep, dlon_ini, dlat_ini, dlon, dlat, &
	416	champ, ibar)
[1000]	417
[1425]	418	IF (extrp) CALL extrapol(champ, imdep, jmdep, 999999., .TRUE., .TRUE., 2, &
	419	work)
[1000]	420
[1319]	421	IF(ibar.AND..NOT.(mode=='SIC'.AND.flag)) THEN
	422	IF(l==1) THEN
[1425]	423	WRITE(lunout, *) &
[1319]	424	'-------------------------------------------------------------------------'
[1425]	425	WRITE(lunout, *) &
	426	'Utilisation de l''interpolation barycentrique pour '//TRIM(title)//' $$$'
	427	WRITE(lunout, *) &
[1319]	428	'-------------------------------------------------------------------------'
	429	END IF
	430	IF(mode=='RUG') champ=LOG(champ)
[1323]	431	CALL inter_barxy(dlon, dlat(:jmdep-1), champ, rlonu(:iim), &
	432	rlatv, champint)
[1319]	433	IF(mode=='RUG') THEN
	434	champint=EXP(champint)
	435	WHERE(NINT(mask)/=1) champint=0.001
	436	END IF
	437	ELSE
	438	SELECT CASE(mode)
[1323]	439	CASE('RUG'); CALL rugosite(imdep, jmdep, dlon, dlat, champ, &
[1319]	440	iim, jjp1, rlonv, rlatu, champint, mask)
[1323]	441	CASE('SIC'); CALL sea_ice (imdep, jmdep, dlon, dlat, champ, &
[1319]	442	iim, jjp1, rlonv, rlatu, champint)
[1425]	443	CASE('SST', 'ALB'); CALL grille_m(imdep, jmdep, dlon, dlat, champ, &
[1319]	444	iim, jjp1, rlonv, rlatu, champint)
	445	END SELECT
	446	END IF
[1425]	447	champtime(:, :, l)=champint
[1319]	448	END DO
[1425]	449	ierr=NF90_CLOSE(ncid); CALL ncerr(ierr, fnam)
[1000]	450
[1425]	451	DEALLOCATE(dlon_ini, dlat_ini, dlon, dlat, champ)
[1319]	452	IF(extrp) DEALLOCATE(work)
[1000]	453
[1425]	454	!--- TIME INTERPOLATION ------------------------------------------------------
	455	WRITE(lunout, *)
	456	WRITE(lunout, *)'INTERPOLATION TEMPORELLE.'
	457	WRITE(lunout, "(2x, ' Vecteur temps en entree: ', 10f6.1)") timeyear
	458	WRITE(lunout, "(2x, ' Vecteur temps en sortie de 0 a ', i3)") ndays
	459	ALLOCATE(yder(lmdep), champan(iip1, jjp1, ndays))
	460	skip = .false.
	461	n_extrap = 0
	462	DO j=1, jjp1
	463	DO i=1, iim
	464	yder = pchsp_95(timeyear, champtime(i, j, :), ibeg=2, iend=2, &
	465	vc_beg=0., vc_end=0.)
	466	CALL pchfe_95(timeyear, champtime(i, j, :), yder, skip, &
	467	arth(0., real(ndays_in) / ndays, ndays), champan(i, j, :), ierr)
	468	if (ierr < 0) stop 1
	469	n_extrap = n_extrap + ierr
[1319]	470	END DO
	471	END DO
[1425]	472	if (n_extrap /= 0) then
	473	print *, "get_2Dfield pchfe_95: n_extrap = ", n_extrap
	474	end if
	475	champan(iip1, :, :)=champan(1, :, :)
	476	DEALLOCATE(yder, champtime, timeyear)
[1000]	477
[1319]	478	!--- Checking the result
[1425]	479	DO j=1, jjp1
	480	CALL minmax(iip1, champan(1, j, 10), chmin, chmax)
	481	WRITE(lunout, *)' '//TRIM(title)//' au temps 10 ', chmin, chmax, j
[1319]	482	END DO
[1000]	483
[1425]	484	!--- SPECIAL FILTER FOR SST: SST>271.38 --------------------------------------
[1319]	485	IF(mode=='SST') THEN
[1425]	486	WRITE(lunout, *) 'Filtrage de la SST: SST >= 271.38'
[1319]	487	WHERE(champan<271.38) champan=271.38
	488	END IF
[1000]	489
[1425]	490	!--- SPECIAL FILTER FOR SIC: 0.0<SIC<1.0 -------------------------------------
[1319]	491	IF(mode=='SIC') THEN
[1425]	492	WRITE(lunout, *) 'Filtrage de la SIC: 0.0 < Sea-ice < 1.0'
	493	IF(.NOT.lCPL) champan(:, :, :)=champan(:, :, :)/100.
	494	champan(iip1, :, :)=champan(1, :, :)
[1319]	495	WHERE(champan>1.0) champan=1.0
	496	WHERE(champan<0.0) champan=0.0
	497	END IF
[1000]	498
[1425]	499	!--- DYNAMICAL TO PHYSICAL GRID ----------------------------------------------
	500	ALLOCATE(champo(klon, ndays))
	501	DO k=1, ndays
	502	CALL gr_dyn_fi(1, iip1, jjp1, klon, champan(1, 1, k), champo(1, k))
[1319]	503	END DO
	504	DEALLOCATE(champan)
[1328]	505
[1319]	506	END SUBROUTINE get_2Dfield
	507	!
	508	!-------------------------------------------------------------------------------
[1000]	509
	510
	511
[1319]	512	!-------------------------------------------------------------------------------
	513	!
	514	FUNCTION year_len(y,cal_in)
	515	!
	516	!-------------------------------------------------------------------------------
	517	USE ioipsl, ONLY : ioget_calendar,ioconf_calendar,lock_calendar,ioget_year_len
	518	IMPLICIT NONE
	519	!-------------------------------------------------------------------------------
	520	! Arguments:
	521	INTEGER :: year_len
	522	INTEGER, INTENT(IN) :: y
	523	CHARACTER(LEN=*), INTENT(IN) :: cal_in
	524	!-------------------------------------------------------------------------------
	525	! Local variables:
	526	CHARACTER(LEN=20) :: cal_out ! calendar (for outputs)
	527	!-------------------------------------------------------------------------------
	528	!--- Getting the input calendar to reset at the end of the function
	529	CALL ioget_calendar(cal_out)
[1000]	530
[1319]	531	!--- Unlocking calendar and setting it to wanted one
	532	CALL lock_calendar(.FALSE.); CALL ioconf_calendar(TRIM(cal_in))
[1000]	533
[1319]	534	!--- Getting the number of days in this year
	535	year_len=ioget_year_len(y)
[1000]	536
[1319]	537	!--- Back to original calendar
	538	CALL lock_calendar(.FALSE.); CALL ioconf_calendar(TRIM(cal_out))
[1000]	539
[1319]	540	END FUNCTION year_len
	541	!
	542	!-------------------------------------------------------------------------------
[1000]	543
	544
[1319]	545	!-------------------------------------------------------------------------------
	546	!
	547	FUNCTION mid_months(y,cal_in,nm)
	548	!
	549	!-------------------------------------------------------------------------------
	550	USE ioipsl, ONLY : ioget_calendar,ioconf_calendar,lock_calendar,ioget_mon_len
	551	IMPLICIT NONE
	552	!-------------------------------------------------------------------------------
	553	! Arguments:
	554	INTEGER, INTENT(IN) :: y ! year
	555	CHARACTER(LEN=*), INTENT(IN) :: cal_in ! calendar
	556	INTEGER, INTENT(IN) :: nm ! months/year number
	557	REAL, DIMENSION(nm) :: mid_months ! mid-month times
	558	!-------------------------------------------------------------------------------
	559	! Local variables:
	560	CHARACTER(LEN=99) :: mess ! error message
	561	CHARACTER(LEN=20) :: cal_out ! calendar (for outputs)
	562	INTEGER, DIMENSION(nm) :: mnth ! months lengths (days)
	563	INTEGER :: m ! months counter
	564	INTEGER :: nd ! number of days
	565	!-------------------------------------------------------------------------------
	566	nd=year_len(y,cal_in)
[1000]	567
[1319]	568	IF(nm==12) THEN
[1000]	569
[1319]	570	!--- Getting the input calendar to reset at the end of the function
	571	CALL ioget_calendar(cal_out)
[1000]	572
[1319]	573	!--- Unlocking calendar and setting it to wanted one
	574	CALL lock_calendar(.FALSE.); CALL ioconf_calendar(TRIM(cal_in))
[1000]	575
[1319]	576	!--- Getting the length of each month
	577	DO m=1,nm; mnth(m)=ioget_mon_len(y,m); END DO
[1000]	578
[1319]	579	!--- Back to original calendar
	580	CALL lock_calendar(.FALSE.); CALL ioconf_calendar(TRIM(cal_out))
[1000]	581
[1319]	582	ELSE IF(MODULO(nd,nm)/=0) THEN
	583	WRITE(mess,'(a,i3,a,i3,a)')'Unconsistent calendar: ',nd,' days/year, but ',&
	584	nm,' months/year. Months number should divide days number.'
	585	CALL abort_gcm('mid_months',TRIM(mess),1)
[1000]	586
[1319]	587	ELSE
	588	mnth=(/(m,m=1,nm,nd/nm)/)
	589	END IF
[1000]	590
[1319]	591	!--- Mid-months times
[1441]	592	mid_months(1)=0.5*REAL(mnth(1))
[1319]	593	DO k=2,nm
[1441]	594	mid_months(k)=mid_months(k-1)+0.5*REAL(mnth(k-1)+mnth(k))
[1319]	595	END DO
[1000]	596
[1319]	597	END FUNCTION mid_months
	598	!
	599	!-------------------------------------------------------------------------------
[1000]	600
	601
[1319]	602	!-------------------------------------------------------------------------------
	603	!
	604	SUBROUTINE ncerr(ncres,fnam)
	605	!
	606	!-------------------------------------------------------------------------------
	607	! Purpose: NetCDF errors handling.
	608	!-------------------------------------------------------------------------------
	609	USE netcdf, ONLY : NF90_NOERR, NF90_STRERROR
	610	IMPLICIT NONE
	611	!-------------------------------------------------------------------------------
	612	! Arguments:
	613	INTEGER, INTENT(IN) :: ncres
	614	CHARACTER(LEN=*), INTENT(IN) :: fnam
	615	!-------------------------------------------------------------------------------
	616	#include "iniprint.h"
	617	IF(ncres/=NF90_NOERR) THEN
	618	WRITE(lunout,*)'Problem with file '//TRIM(fnam)//' in routine limit_netcdf.'
	619	CALL abort_gcm('limit_netcdf',NF90_STRERROR(ncres),1)
	620	END IF
[1000]	621
[1319]	622	END SUBROUTINE ncerr
	623	!
	624	!-------------------------------------------------------------------------------
	625
[1441]	626	#endif
	627	! of #ifdef CPP_EARTH
[1319]	628
	629	END SUBROUTINE limit_netcdf

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

Download in other formats:

Original Format