Context Navigation

← Previous Revision
Next Revision →
Normal
Revision Log

regr_horiz_time_climoz_m.f90

Last change on this file was 5353, checked in by yann meurdesoif, 3 weeks ago
Nvidia compiler has some difficulties to compile correctly some complex array constructor. This commit decompose it the several phases in order to achieve the compilation. Please Lionel and David, have a look to this in order to validate. Probably, in future, when compiler heuristic will be improved, this commit can be reversed. YM
File size: 32.4 KB

Rev	Line
[2788]	1	MODULE regr_horiz_time_climoz_m
	2
	3	USE interpolation, ONLY: locate
[3435]	4	USE mod_grid_phy_lmdz, ONLY: nbp_lon, nbp_lat, grid_type, unstructured
[2788]	5	USE nrtype, ONLY: pi
[5084]	6	USE netcdf, ONLY: NF90_CLOBBER, NF90_FLOAT, NF90_OPEN, &
[2788]	7	NF90_NOWRITE, NF90_NOERR, NF90_GET_ATT, NF90_GLOBAL
	8	USE netcdf95, ONLY: NF95_DEF_DIM, NF95_INQ_DIMID, NF95_INQUIRE_DIMENSION, &
[4489]	9	NF95_DEF_VAR, NF95_INQ_VARID, NF95_INQUIRE_VARIABLE, &
	10	NF95_OPEN, NF95_CREATE, NF95_GET_ATT, NF95_GW_VAR, nf95_get_var, &
	11	NF95_CLOSE, NF95_ENDDEF, NF95_PUT_ATT, NF95_PUT_VAR, NF95_COPY_ATT
[2788]	12	USE print_control_mod, ONLY: lunout
[3435]	13	USE dimphy
[2788]	14	IMPLICIT NONE
	15	PRIVATE
	16	PUBLIC :: regr_horiz_time_climoz
	17	REAL, PARAMETER :: deg2rad=pi/180.
	18	CHARACTER(LEN=13), PARAMETER :: vars_in(2)=['tro3 ','tro3_daylight']
	19
[3435]	20	INTEGER :: nlat_ou, nlon_ou
	21	REAL, ALLOCATABLE :: latitude_glo(:)
	22	!$OMP THREADPRIVATE(latitude_glo)
	23	INTEGER, ALLOCATABLE :: ind_cell_glo_glo(:)
	24	!$OMP THREADPRIVATE(ind_cell_glo_glo)
	25
[2788]	26	CONTAINS
	27
	28	!-------------------------------------------------------------------------------
	29	!
	30	SUBROUTINE regr_horiz_time_climoz(read_climoz,interpt)
	31	!
	32	!-------------------------------------------------------------------------------
	33	! Purpose: Regrid horizontally and in time zonal or 3D ozone climatologies.
	34	! * Read ozone climatology from netcdf file
	35	! * Regrid it horizontaly to LMDZ grid (quasi-conservative method)
	36	! * If interpt=T, interpolate linearly in time (one record each day)
	37	! If interpt=F, keep original time sampling (14 months).
	38	! * Save it to a new netcdf file.
	39	!-------------------------------------------------------------------------------
	40	! Remarks:
	41	! * Up to 2 variables treated: "tro3" and "tro3_daylight" (if read_climoz=2)
	42	! * Input fields coordinates: (longitudes, latitudes, pressure_levels, time)
	43	! * Output grid cells centers coordinates given by [rlonv,] rlatu.
	44	! * Output grid cells edges coordinates given by [rlonu,] rlatv.
	45	! * Input file [longitudes and] latitudes given in degrees.
	46	! * Input file pressure levels are given in Pa or hPa.
	47	! * All coordinates variables are stricly monotonic.
	48	! * Monthly fields are interpolated linearly in time to get daily values.
	49	! * Fields are known at the middle of the months, so interpolation requires an
	50	! additional record both for 1st half of january and 2nd half of december:
	51	! - For a 14-records "climoz.nc": records 1 and 14.
	52	! - For 12-records files:
	53	! record 12 of "climoz_m.nc" if available, or record 1 of "climoz.nc".
	54	! record 1 of "climoz_p.nc" if available, or record 12 of "climoz.nc".
	55	! * Calendar is taken into account to get one record each day (not 360 always).
	56	! * Missing values are filled in from sky to ground by copying lowest valid one.
	57	! Attribute "missing_value" or "_FillValue" must be present in input file.
	58	!-------------------------------------------------------------------------------
	59	USE assert_m, ONLY: assert
	60	USE cal_tools_m, ONLY: year_len, mid_month
[3435]	61	!! USE control_mod, ONLY: anneeref
	62	USE time_phylmdz_mod, ONLY: annee_ref
[2788]	63	USE ioipsl, ONLY: ioget_year_len, ioget_calendar
	64	USE regr_conserv_m, ONLY: regr_conserv
	65	USE regr_lint_m, ONLY: regr_lint
	66	USE regular_lonlat_mod, ONLY: boundslon_reg, boundslat_reg, south, west, east
	67	USE slopes_m, ONLY: slopes
[4619]	68	USE lmdz_xios
[3435]	69	USE mod_phys_lmdz_para, ONLY: is_mpi_root, is_master, is_omp_master, gather, gather_mpi, bcast_mpi, klon_mpi
	70	USE geometry_mod, ONLY : latitude_deg, ind_cell_glo
	71	USE mod_grid_phy_lmdz, ONLY: klon_glo
	72
[2788]	73	!-------------------------------------------------------------------------------
	74	! Arguments:
	75	INTEGER, INTENT(IN) :: read_climoz ! read ozone climatology, 1 or 2
	76	! 1: read a single ozone climatology used day and night
	77	! 2: same + read also a daylight climatology
	78	LOGICAL, INTENT(IN) :: interpt ! TRUE => daily interpolation
	79	! FALSE => no interpolation (14 months)
	80	!-------------------------------------------------------------------------------
	81	! Local variables:
	82
	83	!--- Input files variables
	84	INTEGER :: nlon_in ! Number of longitudes
	85	INTEGER :: nlat_in ! Number of latitudes
	86	INTEGER :: nlev_in ! Number of pressure levels
	87	INTEGER :: nmth_in ! Number of months
[4489]	88	REAL, ALLOCATABLE:: lon_in(:) ! Longitudes (ascending order, rad)
	89	REAL, ALLOCATABLE:: lat_in(:) ! Latitudes (ascending order, rad)
	90	REAL, ALLOCATABLE:: lev_in(:) ! Pressure levels (ascen. order, hPa)
[2788]	91	REAL, ALLOCATABLE :: lon_in_edge(:) ! Longitude intervals edges
	92	! (ascending order, / )
	93	REAL, ALLOCATABLE :: sinlat_in_edge(:) ! Sinus of latitude intervals edges
	94	! (ascending order, / )
	95	LOGICAL :: ldec_lon, ldec_lat, ldec_lev ! Decreasing order in input file
	96	CHARACTER(LEN=20) :: cal_in ! Calendar
	97	REAL, ALLOCATABLE :: o3_in3(:,:,:,:,:) ! Ozone climatologies
[3435]	98	REAL, ALLOCATABLE :: o3_in3bis(:,:,:,:,:) ! Ozone climatologies
[2788]	99	REAL, ALLOCATABLE :: o3_in2 (:,:,:,:) ! Ozone climatologies
[3435]	100	REAL, ALLOCATABLE :: o3_in2bis(:,:,:,:,:) ! Ozone climatologies
[2788]	101	! last index: 1 for the day-night average, 2 for the daylight field.
	102	REAL :: NaN
	103
	104	!--- Partially or totally regridded variables (:,:,nlev_in,:,read_climoz)
	105	REAL, ALLOCATABLE :: o3_regr_lon (:,:,:,:,:) ! (nlon_ou,nlat_in,:,0:13 ,:)
	106	REAL, ALLOCATABLE :: o3_regr_lonlat(:,:,:,:,:) ! (nlon_ou,nlat_ou,:,0:13 ,:)
	107	REAL, ALLOCATABLE :: o3_out3 (:,:,:,:,:) ! (nlon_ou,nlat_ou,:,ntim_ou,:)
[3435]	108	REAL, ALLOCATABLE :: o3_out3_glo (:,:,:,:) ! (nbp_lat,:,ntim_ou,:)
[2788]	109	REAL, ALLOCATABLE :: o3_regr_lat (:,:,:,:) ! (nlat_in,:,0:13 ,:)
	110	REAL, ALLOCATABLE :: o3_out2 (:,:,:,:) ! (nlat_ou,:,ntim_ou,:)
[3435]	111	REAL, ALLOCATABLE :: o3_out2_glo (:,:,:,:) ! (nbp_lat,:,ntim_ou,:)
	112	REAL, ALLOCATABLE :: o3_out (:,:,:,:) ! (nbp_lat,:,ntim_ou,:)
[2788]	113	! Dimension number \| Interval \| Contains \| For variables:
	114	! 1 (longitude) \| [rlonu(i-1), rlonu(i)] \| rlonv(i) \| all
	115	! 2 (latitude) \| [rlatv(j), rlatv(j-1)] \| rlatu(j) \| all but o3_regr_lon
	116	! 3 (press level) \| \| lev(k) \| all
	117	! Note that rlatv(0)=pi/2 and rlatv(nlat_ou)=-pi/2.
	118	! Dimension 4 is: month number (all vars but o3_out)
	119	! days elapsed since Jan. 1st 0h at mid-day (o3_out only)
	120	REAL, ALLOCATABLE :: v1(:)
	121
	122	!--- For NetCDF:
	123	INTEGER :: fID_in_m, fID_in, levID_ou, dimid, vID_in(read_climoz), ntim_ou
	124	INTEGER :: fID_in_p, fID_ou, timID_ou, varid, vID_ou(read_climoz), ndims, ncerr
[4489]	125	INTEGER, ALLOCATABLE :: dIDs(:)
[2788]	126	CHARACTER(LEN=20) :: cal_ou !--- Calendar; no time inter => same as input
	127	CHARACTER(LEN=80) :: press_unit !--- Pressure unit
	128	REAL :: tmidmonth(0:13) !--- Elapsed days since Jan-1 0h at mid-months
	129	! Additional records 0, 13 for interpolation
	130	REAL, ALLOCATABLE :: tmidday(:) !--- Output times (mid-days since Jan 1st 0h)
	131	LOGICAL :: lprev, lnext !--- Flags: previous/next files are present
	132	LOGICAL :: l3D, l2D !--- Flag: input fields are 3D or zonal
	133	INTEGER :: ii, i, j, k, l, m, dln, ib, ie, iv, dx1, dx2
	134	INTEGER, ALLOCATABLE :: sta(:), cnt(:)
	135	CHARACTER(LEN=80) :: sub, dim_nam, msg
[3435]	136	REAL :: null_array(0)
	137	LOGICAL,SAVE :: first=.TRUE.
	138	!$OMP THREADPRIVATE(first)
	139	REAL, ALLOCATABLE :: test_o3_in(:,:)
	140	REAL, ALLOCATABLE :: test_o3_out(:)
[5353]	141	REAL,ALLOCATABLE :: tmp(:)
[2788]	142
	143
[3435]	144	IF (grid_type==unstructured) THEN
	145	IF (first) THEN
	146	IF (is_master) THEN
	147	ALLOCATE(latitude_glo(klon_glo))
	148	ALLOCATE(ind_cell_glo_glo(klon_glo))
	149	ELSE
	150	ALLOCATE(latitude_glo(0))
	151	ALLOCATE(ind_cell_glo_glo(0))
	152	ENDIF
	153	CALL gather(latitude_deg, latitude_glo)
	154	CALL gather(ind_cell_glo, ind_cell_glo_glo)
	155	ENDIF
	156	ENDIF
	157
	158	IF (is_omp_master) THEN
	159	nlat_ou=nbp_lat
	160	nlon_ou=nbp_lon
	161
	162	!-------------------------------------------------------------------------------
	163	IF (is_mpi_root) THEN
	164	sub="regr_horiz_time_climoz"
	165	WRITE(lunout,*)"Call sequence information: "//TRIM(sub)
	166	CALL assert(read_climoz == 1 .OR. read_climoz == 2, "regr_lat_time_climoz")
[2788]	167
[3435]	168	CALL NF95_OPEN("climoz.nc" , NF90_NOWRITE, fID_in)
	169	lprev=NF90_OPEN("climoz_m.nc", NF90_NOWRITE, fID_in_m)==NF90_NOERR
	170	lnext=NF90_OPEN("climoz_p.nc", NF90_NOWRITE, fID_in_p)==NF90_NOERR
[2788]	171
[3435]	172	!--- Get coordinates from the input file. Converts lon/lat in radians.
	173	! Few inversions because "regr_conserv" and gcm need ascending vectors.
	174	CALL NF95_INQ_VARID(fID_in, vars_in(1), varid)
	175	CALL NF95_INQUIRE_VARIABLE(fID_in, varid, dimids=dIDs, ndims=ndims)
	176	l3D=ndims==4; l2D=ndims==3
	177	IF(l3D) WRITE(lunout,*)"Input files contain full 3D ozone fields."
	178	IF(l2D) WRITE(lunout,*)"Input files contain zonal 2D ozone fields."
	179	DO i=1,ndims
	180	CALL NF95_INQUIRE_DIMENSION(fID_in, dIDs(i), name=dim_nam, nclen=dln)
	181	CALL NF95_INQ_VARID(fID_in, dim_nam, varid)
	182	ii=i; IF(l2D) ii=i+1 !--- ndims==3:NO LONGITUDE
	183	SELECT CASE(ii)
	184	CASE(1) !--- LONGITUDE
	185	CALL NF95_GW_VAR(fID_in, varid, lon_in)
	186	ldec_lon=lon_in(1)>lon_in(dln); IF(ldec_lon) lon_in=lon_in(dln:1:-1)
	187	nlon_in=dln; lon_in=lon_in*deg2rad
	188	CASE(2) !--- LATITUDE
	189	CALL NF95_GW_VAR(fID_in, varid, lat_in)
	190	ldec_lat=lat_in(1)>lat_in(dln); IF(ldec_lat) lat_in=lat_in(dln:1:-1)
	191	nlat_in=dln; lat_in=lat_in*deg2rad
	192	CASE(3) !--- PRESSURE LEVELS
	193	CALL NF95_GW_VAR(fID_in, varid, lev_in)
	194	ldec_lev=lev_in(1)>lev_in(dln); IF(ldec_lev) lev_in=lev_in(dln:1:-1)
	195	nlev_in=dln
	196	CALL NF95_GET_ATT(fID_in, varid, "units", press_unit)
	197	k=LEN_TRIM(press_unit)
	198	DO WHILE(ICHAR(press_unit(k:k))==0)
	199	press_unit(k:k)=' '; k=LEN_TRIM(press_unit) !--- REMOVE NULL END CHAR
	200	END DO
	201	IF(press_unit == "Pa") THEN
	202	lev_in = lev_in/100. !--- CONVERT TO hPa
	203	ELSE IF(press_unit /= "hPa") THEN
	204	CALL abort_physic(sub, "the only recognized units are Pa and hPa.",1)
	205	END IF
	206	CASE(4) !--- TIME
	207	CALL NF95_INQUIRE_DIMENSION(fID_in, dIDs(i), nclen=nmth_in)
	208	cal_in='gregorian'
	209	IF(NF90_GET_ATT(fID_in, varid, 'calendar', cal_in)/=NF90_NOERR) &
	210	WRITE(lunout,*)'WARNING: missing "calendar" attribute for "'// &
	211	TRIM(dim_nam)//'" in "climoz.nc". Choosing default: "gregorian".'
	212	k=LEN_TRIM(cal_in)
	213	DO WHILE(ICHAR(cal_in(k:k))==0)
	214	cal_in(k:k)=' '; k=LEN_TRIM(cal_in) !--- REMOVE NULL END CHAR
	215	END DO
	216	END SELECT
	217	END DO
[2788]	218
[3465]	219	!--- Prepare quantities for time interpolation
	220	tmidmonth=mid_month(annee_ref, cal_in)
	221	IF(interpt) THEN
	222	ntim_ou=ioget_year_len(annee_ref)
	223	ALLOCATE(tmidday(ntim_ou))
	224	tmidday=[(REAL(k)-0.5,k=1,ntim_ou)]
	225	CALL ioget_calendar(cal_ou)
	226	ELSE
	227	ntim_ou=14
	228	cal_ou=cal_in
	229	END IF
	230	ENDIF
[2788]	231
[3465]	232	IF (grid_type==unstructured) THEN
	233	CALL bcast_mpi(nlon_in)
	234	CALL bcast_mpi(nlat_in)
	235	CALL bcast_mpi(nlev_in)
	236	CALL bcast_mpi(l3d)
	237	CALL bcast_mpi(tmidmonth)
[4847]	238	IF(interpt) CALL bcast_mpi(tmidday)
[3465]	239	CALL bcast_mpi(ntim_ou)
	240
	241	IF (is_mpi_root) THEN
	242	CALL xios_set_domain_attr("domain_climoz",nj_glo=nlat_in, nj=nlat_in, jbegin=0, latvalue_1d=lat_in/deg2rad)
	243	IF (l3D) THEN
	244	CALL xios_set_domain_attr("domain_climoz",ni_glo=nlon_in, ni=nlon_in, ibegin=0, lonvalue_1d=lon_in/deg2rad)
	245	ELSE
	246	CALL xios_set_domain_attr("domain_climoz",ni_glo=8, ni=8, ibegin=0, lonvalue_1d = (/ 0.,45.,90.,135.,180.,225.,270., 315. /))
	247	ENDIF
	248	ELSE
	249	CALL xios_set_domain_attr("domain_climoz",nj_glo=nlat_in, nj=0, jbegin=0, latvalue_1d=null_array )
	250	IF (l3D) THEN
	251	CALL xios_set_domain_attr("domain_climoz",ni_glo=nlon_in, ni=0, ibegin=0, lonvalue_1d=null_array)
	252	ELSE
	253	CALL xios_set_domain_attr("domain_climoz",ni_glo=8, ni=0, ibegin=0, lonvalue_1d=null_array)
	254	ENDIF
	255	ENDIF
	256	CALL xios_set_axis_attr("axis_climoz", n_glo=nlev_in)
	257	CALL xios_set_axis_attr("time_axis_climoz", n_glo=ntim_ou)
	258	CALL xios_set_axis_attr("time_axis_climoz", n_glo=ntim_ou)
	259	CALL xios_set_axis_attr("tr_climoz", n_glo=read_climoz)
	260	CALL xios_set_field_attr("tro3_out", enabled=.TRUE.)
	261	CALL xios_set_field_attr("tro3_out", enabled=.TRUE.)
	262
[3435]	263	IF (first) THEN
	264	first=.FALSE.
	265	RETURN
	266	ENDIF
	267	ENDIF
	268
	269
	270	IF (is_mpi_root) THEN
	271	!--- Longitudes management:
	272	! * Need to shift data if the origin of input file longitudes /= -pi
	273	! * Need to add some margin in longitude to ensure input interval contains
	274	! all the output intervals => at least one longitudes slice has to be
	275	! duplicated, possibly more for undersampling.
	276	IF(l3D) THEN
	277	IF (grid_type==unstructured) THEN
	278	dx2=0
	279	ELSE
	280	!--- Compute input edges longitudes vector (no end point yet)
	281	ALLOCATE(v1(nlon_in+1))
	282	v1(1)=(lon_in(nlon_in)+lon_in(1))/2.-pi
	283	FORALL(i=2:nlon_in) v1(i)=(lon_in(i-1)+lon_in(i))/2.
	284	v1(nlon_in+1)=v1(1)+2.*pi
	285	DEALLOCATE(lon_in)
[2788]	286
[3435]	287	!--- Shift input longitudes vector until it contains first output point boundslon_reg(1,west)
	288	v1=v1+2piREAL(FLOOR((boundslon_reg(1,west)-v1(1))/(2.*pi)))
[2788]	289
[3435]	290	!--- Ensure first input longitudes interval contains first output point boundslon_reg(1,west)
	291	dx1=locate(v1,boundslon_reg(1,west))-1
	292	v1=CSHIFT(v1,SHIFT=dx1,DIM=1)
[3465]	293	v1(nlon_in-dx1+2:)=v1(nlon_in-dx1+2:)+2.*pi
[3435]	294
	295	!--- Extend input longitudes vector until last interval contains boundslon_reg(nlat_ou,east)
[3465]	296	dx2=0; DO WHILE(v1(1+dx2)+2.*pi<=boundslon_reg(nlon_ou,east)); dx2=dx2+1; END DO
[2788]	297
[3435]	298	!--- Final edges longitudes vector (with margin and end point)
	299	ALLOCATE(lon_in_edge(nlon_in+dx2+1)); lon_in_edge=[v1,v1(2:1+dx2)+2.*pi]
	300	DEALLOCATE(v1)
	301	ENDIF
	302	END IF
[2788]	303
[3435]	304	!--- Compute sinus of intervals edges latitudes:
	305	ALLOCATE(sinlat_in_edge(nlat_in+1))
	306	sinlat_in_edge(1) = -1. ; sinlat_in_edge(nlat_in+1) = 1.
	307	FORALL(j=2:nlat_in) sinlat_in_edge(j)=SIN((lat_in(j-1)+lat_in(j))/2.)
	308	DEALLOCATE(lat_in)
[2788]	309
	310
	311
[3435]	312	!--- Check for contiguous years:
	313	ib=0; ie=13
	314	IF(nmth_in == 14) THEN; lprev=.FALSE.; lnext=.FALSE.
	315	WRITE(lunout,*)'Using 14 months ozone climatology "climoz.nc"...'
	316	ELSE
	317	IF( lprev) WRITE(lunout,*)'Using "climoz_m.nc" last record (previous year).'
	318	IF(.NOT.lprev) WRITE(lunout,*)"No previous year file ; assuming periodicity."
	319	IF( lnext) WRITE(lunout,*)'Using "climoz_p.nc" first record (next year).'
	320	IF(.NOT.lnext) WRITE(lunout,*)"No next year file ; assuming periodicity."
	321	IF(.NOT.lprev) ib=1
	322	IF(.NOT.lnext) ie=12
[2788]	323	END IF
[3435]	324	ALLOCATE(sta(ndims),cnt(ndims)); sta(:)=1
	325	IF(l3D) cnt=[nlon_in,nlat_in,nlev_in,1]
	326	IF(l2D) cnt=[ nlat_in,nlev_in,1]
	327	IF(l3D) ALLOCATE(o3_in3(nlon_in+dx2,nlat_in,nlev_in,ib:ie,read_climoz))
	328	IF(l2D) ALLOCATE(o3_in2( nlat_in,nlev_in,ib:ie,read_climoz))
[2788]	329
[3435]	330	!--- Read full current file and one record each available contiguous file
	331	DO iv=1,read_climoz
	332	CALL NF95_INQ_VARID(fID_in, vars_in(1), vID_in(iv))
[4489]	333	IF(l3D) call NF95_GET_VAR(fID_in, vID_in(iv), o3_in3(1:nlon_in,:,:,1:12,iv))
	334	IF(l2D) call NF95_GET_VAR(fID_in, vID_in(iv), o3_in2( :,:,1:12,iv))
[3435]	335	IF(lprev) THEN; sta(ndims)=12
	336	CALL NF95_INQ_VARID(fID_in_m, vars_in(1), vID_in(iv))
[4489]	337	IF(l3D) call NF95_GET_VAR(fID_in_m,vID_in(iv),o3_in3(1:nlon_in,:,:, 0,iv),sta,cnt)
	338	IF(l2d) call NF95_GET_VAR(fID_in_m,vID_in(iv),o3_in2( :,:, 0,iv),sta,cnt)
[3435]	339	END IF
	340	IF(lnext) THEN; sta(ndims)=1
	341	CALL NF95_INQ_VARID(fID_in_p, vars_in(1), vID_in(iv))
[4489]	342	IF(l3D) call NF95_GET_VAR(fID_in_p,vID_in(iv),o3_in3(1:nlon_in,:,:,13,iv),sta,cnt)
	343	IF(l2D) call NF95_GET_VAR(fID_in_p,vID_in(iv),o3_in2( :,:,13,iv),sta,cnt)
[3435]	344	END IF
[2788]	345	END DO
[3435]	346	IF(lprev.OR.lnext) DEALLOCATE(sta,cnt)
	347	IF(lprev) CALL NF95_CLOSE(fID_in_m)
	348	IF(lnext) CALL NF95_CLOSE(fID_in_p)
	349
	350	!--- Revert decreasing coordinates vector
	351	IF(l3D) THEN
	352	IF(ldec_lon) o3_in3(1:nlon_in,:,:,:,:) = o3_in3(nlon_in:1:-1,:,:,:,:)
	353	IF(ldec_lat) o3_in3 = o3_in3(:,nlat_in:1:-1,:,:,:)
	354	IF(ldec_lev) o3_in3 = o3_in3(:,:,nlev_in:1:-1,:,:)
	355
	356	IF (grid_type /= unstructured) THEN
	357	!--- Shift values for longitude and duplicate some longitudes slices
	358	o3_in3(1:nlon_in,:,:,:,:)=CSHIFT(o3_in3(1:nlon_in,:,:,:,:),SHIFT=dx1,DIM=1)
	359	o3_in3(nlon_in+1:nlon_in+dx2,:,:,:,:)=o3_in3(1:dx2,:,:,:,:)
	360	ENDIF
	361	ELSE
	362	IF(ldec_lat) o3_in2 = o3_in2( nlat_in:1:-1,:,:,:)
	363	IF(ldec_lev) o3_in2 = o3_in2( :,nlev_in:1:-1,:,:)
[2788]	364	END IF
	365
[3435]	366	!--- Deal with missing values
	367	DO m=1, read_climoz
	368	WRITE(msg,'(a,i0)')"regr_lat_time_climoz: field Nr.",m
	369	IF(NF90_GET_ATT(fID_in,vID_in(m),"missing_value",NaN)/= NF90_NOERR) THEN
	370	IF(NF90_GET_ATT(fID_in, vID_in(m),"_FillValue",NaN)/= NF90_NOERR) THEN
	371	WRITE(lunout,*)TRIM(msg)//": no missing value attribute found."; CYCLE
	372	END IF
	373	END IF
	374	WRITE(lunout,*)TRIM(msg)//": missing value attribute found."
	375	WRITE(lunout,*)"Trying to fill in NaNs ; a full field would be better."
	376
	377	!--- Check top layer contains no NaNs & search NaNs from top to ground
	378	msg=TRIM(sub)//": NaNs in top layer !"
	379	IF(l3D) THEN
	380	IF(ANY(o3_in3(:,:,1,:,m)==NaN)) CALL abort_physic(sub,msg,1)
	381	DO k = 2,nlev_in
	382	WHERE(o3_in3(:,:,k,:,m)==NaN) o3_in3(:,:,k,:,m)=o3_in3(:,:,k-1,:,m)
	383	END DO
	384	ELSE
	385	IF(ANY(o3_in2( :,1,:,m)==NaN)) THEN
	386	WRITE(lunout,*)msg
	387	!--- Fill in latitudes where all values are missing
	388	DO l=1,nmth_in
	389	!--- Next to south pole
	390	j=1; DO WHILE(o3_in2(j,1,l,m)==NaN); j=j+1; END DO
	391	IF(j>1) &
	392	o3_in2(:j-1,:,l,m)=SPREAD(o3_in2(j,:,l,m),DIM=1,ncopies=j-1)
	393	!--- Next to north pole
	394	j=nlat_in; DO WHILE(o3_in2(j,1,l,m)==NaN); j=j+1; END DO
	395	IF(j<nlat_in) &
	396	o3_in2(j+1:,:,l,m)=SPREAD(o3_in2(j,:,l,m),DIM=1,ncopies=nlat_in-j)
	397	END DO
	398	END IF
	399
	400	!--- Fill in high latitudes missing values
	401	!--- Highest level been filled-in, so has always valid values.
	402	DO k = 2,nlev_in
	403	WHERE(o3_in2(:,k,:,m)==NaN) o3_in2(:,k,:,m)=o3_in2(:,k-1,:,m)
	404	END DO
	405	END IF
[2788]	406	END DO
	407
[3435]	408	ENDIF
	409
	410	!=============================================================================
	411	IF(l3D) THEN !=== 3D FIELDS
	412	!=============================================================================
	413	IF (grid_type==unstructured) THEN
	414	nlat_ou=klon_mpi
	415
	416	IF (is_mpi_root) THEN
	417	ALLOCATE(o3_in3bis(nlon_in,nlat_in,nlev_in,0:13,read_climoz))
	418	o3_in3bis(:,:,:,ib:ie,:)=o3_in3(1:nlon_in,:,:,ib:ie,:)
	419	ELSE
	420	ALLOCATE(o3_in3bis(0,0,0,0,read_climoz))
	421	ENDIF
	422	ALLOCATE(o3_regr_lonlat(1, nlat_ou, nlev_in, 0:13, read_climoz))
	423
	424	CALL xios_send_field("tro3_in",o3_in3bis(:,:,:,:,:))
	425	CALL xios_recv_field("tro3_out",o3_regr_lonlat(1,:,:,:,:))
	426	ELSE
	427
	428	!--- Regrid in longitude
	429	ALLOCATE(o3_regr_lon(nlon_ou, nlat_in, nlev_in, ie-ib+1, read_climoz))
[5353]	430	tmp = [boundslon_reg(1,west),boundslon_reg(:,east)]
[3435]	431	CALL regr_conserv(1, o3_in3, xs = lon_in_edge, &
[5353]	432	xt = tmp , &
	433	vt = o3_regr_lon, slope = slopes(1,o3_in3, lon_in_edge))
[3435]	434	DEALLOCATE(o3_in3)
[2788]	435
[3435]	436	!--- Regrid in latitude: averaging with respect to SIN(lat) is
	437	! equivalent to weighting by COS(lat)
	438	!--- (inverted indices in "o3_regr_lonlat" because "rlatu" is decreasing)
	439	ALLOCATE(o3_regr_lonlat(nlon_ou, nlat_ou, nlev_in, 0:13, read_climoz))
[5353]	440	tmp = [- 1., SIN(boundslat_reg(nlat_ou-1:1:-1,south)), 1.]
	441	CALL regr_conserv(2, o3_regr_lon, xs = sinlat_in_edge, &
	442	xt = tmp, &
	443	vt = o3_regr_lonlat(:,nlat_ou:1:- 1,:,ib:ie,:), &
	444	slope = slopes(2,o3_regr_lon, sinlat_in_edge))
[3435]	445	DEALLOCATE(o3_regr_lon)
[2788]	446
[3435]	447	ENDIF
[2788]	448
[3435]	449	!--- Duplicate previous/next record(s) if they are not available
	450	IF(.NOT.lprev) o3_regr_lonlat(:,:,:, 0,:) = o3_regr_lonlat(:,:,:,12,:)
	451	IF(.NOT.lnext) o3_regr_lonlat(:,:,:,13,:) = o3_regr_lonlat(:,:,:, 1,:)
	452
	453	!--- Regrid in time by linear interpolation:
	454	ALLOCATE(o3_out3(nlon_ou, nlat_ou, nlev_in, ntim_ou, read_climoz))
	455	IF( interpt) CALL regr_lint(4,o3_regr_lonlat,tmidmonth,tmidday,o3_out3)
	456	IF(.NOT.interpt) o3_out3=o3_regr_lonlat
	457	DEALLOCATE(o3_regr_lonlat)
[2788]	458
[3435]	459	nlat_ou=nbp_lat
	460	IF (grid_type==unstructured) THEN
	461	CALL xios_send_field('o3_out',o3_out3)
	462	ndims=3
	463	ALLOCATE(o3_out3_glo(nlat_ou, nlev_in, ntim_ou, read_climoz))
	464	CALL gather_mpi(o3_out3(1,:,:,:,:), o3_out3_glo)
	465	ENDIF
[2788]	466
[3435]	467	!--- Create the output file and get the variable IDs:
	468	CALL prepare_out(fID_in,nlev_in,ntim_ou, fID_ou,levID_ou,timID_ou,vID_ou, &
	469	ndims, cal_ou)
[2788]	470
[3435]	471	IF (is_mpi_root) THEN
	472	!--- Write remaining coordinate variables:
	473	CALL NF95_PUT_VAR(fID_ou, levID_ou, lev_in); DEALLOCATE(lev_in)
	474	IF( interpt) CALL NF95_PUT_VAR(fID_ou, timID_ou, tmidday)
	475	IF(.NOT.interpt) CALL NF95_PUT_VAR(fID_ou, timID_ou, tmidmonth)
[2788]	476
[3435]	477	!--- Write to file (the order of "rlatu" is inverted in the output file):
	478	IF (grid_type==unstructured) THEN
[2788]	479
[3435]	480	ALLOCATE(o3_out(nlat_ou, nlev_in, ntim_ou, read_climoz))
	481	DO i=1,klon_glo
	482	o3_out(ind_cell_glo_glo(i),:,:,:)=o3_out3_glo(i,:,:,:)
	483	ENDDO
[2788]	484
[3435]	485	DO m = 1, read_climoz
	486	CALL NF95_PUT_VAR(fID_ou, vID_ou(m), o3_out(nlat_ou:1:-1,:,:,m))
	487	END DO
	488
	489	ELSE
	490	DO m = 1, read_climoz
	491	CALL NF95_PUT_VAR(fID_ou, vID_ou(m), o3_out3(:,nlat_ou:1:-1,:,:,m))
	492	END DO
	493	ENDIF
	494	CALL NF95_CLOSE(fID_ou)
[2788]	495
[3435]	496
	497	ENDIF
	498
	499
	500	!=============================================================================
	501	ELSE !=== ZONAL FIELDS
	502	!=============================================================================
	503
	504	IF (grid_type==unstructured) THEN
	505	nlat_ou=klon_mpi
	506
	507	IF (is_mpi_root) THEN
	508	ALLOCATE(o3_in2bis(8,nlat_in,nlev_in,0:13,read_climoz))
	509	o3_in2bis(:,:,:,ib:ie,:)=SPREAD(o3_in2,1,8)
	510	ELSE
	511	ALLOCATE(o3_in2bis(0,0,0,0,read_climoz))
	512	ENDIF
	513	ALLOCATE(o3_regr_lat(nlat_ou, nlev_in, 0:13, read_climoz))
	514	CALL xios_send_field("tro3_in",o3_in2bis(:,:,:,:,:))
	515	CALL xios_recv_field("tro3_out",o3_regr_lat(:,:,:,:))
	516	IF(.NOT.lprev) o3_regr_lat(:,:, 0, :) = o3_regr_lat(:,:,12,:)
	517	IF(.NOT.lnext) o3_regr_lat(:,:,13, :) = o3_regr_lat(:,:, 1,:)
	518
	519	ELSE
	520	!--- Regrid in latitude: averaging with respect to SIN(lat) is
	521	! equivalent to weighting by COS(lat)
	522	!--- (inverted indices in "o3_regr_lat" because "rlatu" is decreasing)
	523	ALLOCATE(o3_regr_lat(nlat_ou, nlev_in, 0:13, read_climoz))
[5353]	524	tmp = [- 1., SIN(boundslat_reg(nlat_ou-1:1:-1,south)), 1.]
[3435]	525	CALL regr_conserv(1, o3_in2, xs = sinlat_in_edge, &
[5353]	526	xt = tmp, &
[3435]	527	vt = o3_regr_lat(nlat_ou:1:- 1,:,ib:ie,:), &
	528	slope = slopes(1,o3_in2, sinlat_in_edge))
	529	DEALLOCATE(o3_in2)
	530
	531	!--- Duplicate previous/next record(s) if they are not available
	532	IF(.NOT.lprev) o3_regr_lat(:,:, 0,:) = o3_regr_lat(:,:,12,:)
	533	IF(.NOT.lnext) o3_regr_lat(:,:,13,:) = o3_regr_lat(:,:, 1,:)
	534
	535	ENDIF
	536
	537	!--- Regrid in time by linear interpolation:
	538	ALLOCATE(o3_out2(nlat_ou, nlev_in, ntim_ou, read_climoz))
	539	IF( interpt) CALL regr_lint(3,o3_regr_lat, tmidmonth, tmidday, o3_out2)
	540	IF(.NOT.interpt) o3_out2=o3_regr_lat
	541	DEALLOCATE(o3_regr_lat)
	542
	543	nlat_ou=nbp_lat
	544
	545	IF (grid_type==unstructured) THEN
	546	ndims=3
	547	ALLOCATE(o3_out2_glo(nlat_ou, nlev_in, ntim_ou, read_climoz))
	548	CALL gather_mpi(o3_out2, o3_out2_glo)
	549	ENDIF
	550
	551	!--- Create the output file and get the variable IDs:
	552	CALL prepare_out(fID_in,nlev_in,ntim_ou, fID_ou,levID_ou,timID_ou,vID_ou, &
	553	ndims, cal_ou)
	554
	555	IF (is_mpi_root) THEN
	556
	557	!--- Write remaining coordinate variables:
	558	CALL NF95_PUT_VAR(fID_ou, levID_ou, lev_in); DEALLOCATE(lev_in)
	559	IF( interpt) CALL NF95_PUT_VAR(fID_ou, timID_ou, tmidday)
	560	IF(.NOT.interpt) CALL NF95_PUT_VAR(fID_ou, timID_ou, tmidmonth)
	561
	562	IF (grid_type==unstructured) THEN
	563
	564	ALLOCATE(o3_out3_glo(nlat_ou, nlev_in, ntim_ou, read_climoz))
	565	DO i=1,klon_glo
	566	o3_out(ind_cell_glo_glo(i),:,:,:)=o3_out2_glo(i,:,:,:)
	567	ENDDO
	568
	569
	570	DO m = 1, read_climoz
	571	CALL NF95_PUT_VAR(fID_ou, vID_ou(m), o3_out(nlat_ou:1:-1,:,:,m))
	572	END DO
	573	ELSE
	574	!--- Write to file (the order of "rlatu" is inverted in the output file):
	575	DO m = 1, read_climoz
	576	CALL NF95_PUT_VAR(fID_ou, vID_ou(m), o3_out2(nlat_ou:1:-1,:,:,m))
	577	END DO
	578	ENDIF
	579
	580	CALL NF95_CLOSE(fID_ou)
	581
	582	ENDIF
	583
	584	!=============================================================================
	585	END IF
	586	!=============================================================================
	587
	588	IF (is_mpi_root) CALL NF95_CLOSE(fID_in)
	589
	590	ENDIF ! is_omp_master
	591
	592	first=.FALSE.
[2788]	593	END SUBROUTINE regr_horiz_time_climoz
	594	!
	595	!-------------------------------------------------------------------------------
	596
	597
	598	!-------------------------------------------------------------------------------
	599	!
	600	SUBROUTINE prepare_out(fID_in, nlev_in, ntim_ou, fID_ou, vlevID, vtimID, &
	601	vID_ou, ndims, cal_ou)
	602	!-------------------------------------------------------------------------------
	603	! Purpose: This subroutine creates the NetCDF output file, defines
	604	! dimensions and variables, and writes some of the coordinate variables.
	605	!-------------------------------------------------------------------------------
	606	USE regular_lonlat_mod, ONLY: lon_reg, lat_reg
[3435]	607	USE regular_lonlat_mod, ONLY: lon_reg, lat_reg
	608	USE mod_phys_lmdz_para, ONLY: is_mpi_root
	609	USE mod_grid_phy_lmdz, ONLY: klon_glo
	610	!
[2788]	611	!-------------------------------------------------------------------------------
	612	! Arguments:
	613	INTEGER, INTENT(IN) :: fID_in, nlev_in, ntim_ou
	614	INTEGER, INTENT(OUT) :: fID_ou, vlevID, vtimID
	615	INTEGER, INTENT(OUT) :: vID_ou(:) ! dim(1/2) 1: O3day&night 2: O3daylight
	616	INTEGER, INTENT(IN) :: ndims ! fields rank (3 or 4)
	617	CHARACTER(LEN=*), INTENT(IN) :: cal_ou ! calendar
	618	!-------------------------------------------------------------------------------
	619	! Local variables:
	620	INTEGER :: dlonID, dlatID, dlevID, dtimID, dIDs(4)
	621	INTEGER :: vlonID, vlatID, ncerr, is
[3435]	622	REAL,ALLOCATABLE :: latitude_glo_(:)
[2788]	623	CHARACTER(LEN=80) :: sub
[3435]	624	INTEGER :: i
	625
	626
[2788]	627	!-------------------------------------------------------------------------------
[3435]	628
	629	IF (is_mpi_root) THEN
	630	sub="prepare_out"
	631	WRITE(lunout,*)"CALL sequence information: "//TRIM(sub)
	632	CALL NF95_CREATE("climoz_LMDZ.nc", NF90_clobber, fID_ou)
[2788]	633
	634	!--- Dimensions:
[3435]	635	IF(ndims==4) &
	636	CALL NF95_DEF_DIM(fID_ou, "rlonv", nlon_ou, dlonID)
	637	CALL NF95_DEF_DIM(fID_ou, "rlatu", nlat_ou, dlatID)
	638	CALL NF95_DEF_DIM(fID_ou, "plev", nlev_in, dlevID)
	639	CALL NF95_DEF_DIM(fID_ou, "time", ntim_ou, dtimID)
[2788]	640
[3435]	641	!--- Define coordinate variables:
	642	IF(ndims==4) &
	643	CALL NF95_DEF_VAR(fID_ou, "rlonv", NF90_FLOAT, dlonID, vlonID)
	644	CALL NF95_DEF_VAR(fID_ou, "rlatu", NF90_FLOAT, dlatID, vlatID)
	645	CALL NF95_DEF_VAR(fID_ou, "plev", NF90_FLOAT, dlevID, vlevID)
	646	CALL NF95_DEF_VAR(fID_ou, "time", NF90_FLOAT, dtimID, vtimID)
	647	IF(ndims==4) &
	648	CALL NF95_PUT_ATT(fID_ou, vlonID, "units", "degrees_east")
	649	CALL NF95_PUT_ATT(fID_ou, vlatID, "units", "degrees_north")
	650	CALL NF95_PUT_ATT(fID_ou, vlevID, "units", "millibar")
	651	CALL NF95_PUT_ATT(fID_ou, vtimID, "units", "days since 2000-1-1")
	652	IF(ndims==4) &
	653	CALL NF95_PUT_ATT(fID_ou, vlonID, "standard_name", "longitude")
	654	CALL NF95_PUT_ATT(fID_ou, vlatID, "standard_name", "latitude")
	655	CALL NF95_PUT_ATT(fID_ou, vlevID, "standard_name", "air_pressure")
	656	CALL NF95_PUT_ATT(fID_ou, vtimID, "standard_name", "time")
	657	CALL NF95_PUT_ATT(fID_ou, vlevID, "long_name", "air pressure")
	658	CALL NF95_PUT_ATT(fID_ou, vtimID, "calendar", cal_ou)
[2788]	659
	660	!--- Define the main variables:
[3435]	661	IF(ndims==3) dIDs(1:3) = [ dlatID, dlevID, dtimID]
	662	IF(ndims==4) dIDs=[dlonID, dlatID, dlevID, dtimID]
	663	CALL NF95_DEF_VAR(fID_ou, vars_in(1), NF90_FLOAT, dIDs(1:ndims), vID_ou(1))
	664	CALL NF95_PUT_ATT(fID_ou, vID_ou(1), "long_name", "ozone mole fraction")
	665	CALL NF95_PUT_ATT(fID_ou, vID_ou(1), "standard_name", "mole_fraction_of_ozone&
[2788]	666	&_in_air")
[3435]	667	IF(SIZE(vID_ou) == 2) THEN
	668	CALL NF95_DEF_VAR(fID_ou, vars_in(2), NF90_FLOAT, dIDs(1:ndims), vID_ou(2))
	669	CALL NF95_PUT_ATT(fID_ou, vID_ou(2), "long_name","ozone mole fraction in da&
	670	&ylight")
	671	END IF
[2788]	672
	673	!--- Global attributes:
	674	! The following commands, copying attributes, may fail. That is OK.
	675	! It should just mean that the attribute is not defined in the input file.
[3435]	676	CALL NF95_COPY_ATT(fID_in,NF90_GLOBAL,"Conventions",fID_ou,NF90_GLOBAL, ncerr)
	677	CALL handle_err_copy_att("Conventions")
	678	CALL NF95_COPY_ATT(fID_in,NF90_GLOBAL,"title", fID_ou,NF90_GLOBAL, ncerr)
	679	CALL handle_err_copy_att("title")
	680	CALL NF95_COPY_ATT(fID_in,NF90_GLOBAL,"institution",fID_ou,NF90_GLOBAL, ncerr)
	681	CALL handle_err_copy_att("institution")
	682	CALL NF95_COPY_ATT(fID_in,NF90_GLOBAL,"source", fID_ou,NF90_GLOBAL, ncerr)
	683	CALL handle_err_copy_att("source")
	684	CALL NF95_PUT_ATT (fID_ou,NF90_GLOBAL,"comment", "Regridded for LMDZ")
	685	CALL NF95_ENDDEF(fID_ou)
[2788]	686
[3435]	687	IF (grid_type==unstructured) THEN
	688	ALLOCATE(latitude_glo_(klon_glo))
	689	DO i=1,klon_glo
	690	latitude_glo_(ind_cell_glo_glo(i))=latitude_glo(i)
	691	ENDDO
	692	CALL NF95_PUT_VAR(fID_ou, vlatID, latitude_glo_)
	693	ELSE
	694	!--- Write one of the coordinate variables:
	695	IF(ndims==4) CALL NF95_PUT_VAR(fID_ou, vlonID, lon_reg/deg2rad)
	696	CALL NF95_PUT_VAR(fID_ou, vlatID, lat_reg(nlat_ou:1:-1)/deg2rad)
	697	! (convert from rad to degrees and sort in ascending order)
	698	ENDIF
	699	ENDIF
	700
[2788]	701	CONTAINS
	702
	703	!-------------------------------------------------------------------------------
	704	!
	705	SUBROUTINE handle_err_copy_att(att_name)
	706	!
	707	!-------------------------------------------------------------------------------
[5084]	708	USE netcdf, ONLY: NF90_NOERR, NF90_strerror
[2788]	709	!-------------------------------------------------------------------------------
	710	! Arguments:
	711	CHARACTER(LEN=*), INTENT(IN) :: att_name
	712	!-------------------------------------------------------------------------------
	713	IF(ncerr /= NF90_NOERR) &
	714	WRITE(lunout,*)TRIM(sub)//" prepare_out NF95_COPY_ATT "//TRIM(att_name)// &
	715	" -- "//TRIM(NF90_strerror(ncerr))
	716
	717	END SUBROUTINE handle_err_copy_att
	718	!
	719	!-------------------------------------------------------------------------------
	720
	721	END SUBROUTINE prepare_out
	722	!
	723	!-------------------------------------------------------------------------------
	724
	725	END MODULE regr_horiz_time_climoz_m
	726	!
	727	!-------------------------------------------------------------------------------

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

Download in other formats:

Original Format