Context Navigation

limit_netcdf.F90 @ 2337

Last change on this file since 2337 was 2298, checked in by Laurent Fairhead, 9 years ago
Merged trunk changes -r2237:2291 into testing branch
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
File size: 27.5 KB

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

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

Download in other formats:

Original Format