Context Navigation

create_limit_unstruct.F90 @ 5018

Last change on this file since 5018 was 4361, checked in by lguez, 2 years ago

Change calendar attribute "360d" to "360_day"

"360_day" is the correct attribute according to CF convention. "360d"
leads to an error when opening a history file with xarray.

File size: 12.0 KB

Line
1	MODULE create_limit_unstruct_mod
2	PRIVATE
3	INTEGER, PARAMETER :: lmdep=12
4
5	PUBLIC create_limit_unstruct
6
7	CONTAINS
8
9
10	SUBROUTINE create_limit_unstruct
11	USE dimphy
12	#ifdef CPP_XIOS
13	USE xios
14	USE ioipsl, ONLY : ioget_year_len
15	USE time_phylmdz_mod, ONLY : annee_ref
16	USE indice_sol_mod
17	USE phys_state_var_mod
18	USE mod_phys_lmdz_para
19	IMPLICIT NONE
20	INCLUDE "iniprint.h"
21	REAL, DIMENSION(:,:),ALLOCATABLE :: sic
22	REAL, DIMENSION(:,:),ALLOCATABLE :: sst
23	REAL, DIMENSION(klon,lmdep) :: rugos
24	REAL, DIMENSION(klon,lmdep) :: albedo
25	REAL, DIMENSION(:,:),ALLOCATABLE :: sic_mpi
26	REAL, DIMENSION(:,:),ALLOCATABLE :: sst_mpi
27	REAL, DIMENSION(klon_mpi,lmdep) :: rugos_mpi
28	REAL, DIMENSION(klon_mpi,lmdep) :: albedo_mpi
29	INTEGER :: ndays
30	REAL :: fi_ice(klon)
31	REAL, ALLOCATABLE :: sic_year(:,:)
32	REAL, ALLOCATABLE :: sst_year(:,:)
33	REAL, ALLOCATABLE :: rugos_year(:,:)
34	REAL, ALLOCATABLE :: albedo_year(:,:)
35	REAL, ALLOCATABLE :: pctsrf_t(:,:,:)
36	REAL, ALLOCATABLE :: phy_bil(:,:)
37	REAL, ALLOCATABLE :: sst_year_mpi(:,:)
38	REAL, ALLOCATABLE :: rugos_year_mpi(:,:)
39	REAL, ALLOCATABLE :: albedo_year_mpi(:,:)
40	REAL, ALLOCATABLE :: pctsrf_t_mpi(:,:,:)
41	REAL, ALLOCATABLE :: phy_bil_mpi(:,:)
42	INTEGER :: l,k
43	INTEGER :: nbad
44	INTEGER :: sic_time_axis_size
45	INTEGER :: sst_time_axis_size
46	CHARACTER(LEN=99) :: mess ! error message
47
48
49	ndays=ioget_year_len(annee_ref)
50
51	IF (is_omp_master) CALL xios_get_axis_attr("time_sic",n_glo=sic_time_axis_size)
52	CALL bcast_omp(sic_time_axis_size)
53	ALLOCATE(sic_mpi(klon_mpi,sic_time_axis_size))
54	ALLOCATE(sic(klon,sic_time_axis_size))
55
56
57	IF (is_omp_master) CALL xios_get_axis_attr("time_sst",n_glo=sst_time_axis_size)
58	CALL bcast_omp(sst_time_axis_size)
59	ALLOCATE(sst_mpi(klon_mpi,sst_time_axis_size))
60	ALLOCATE(sst(klon,sst_time_axis_size))
61
62	IF (is_omp_master) THEN
63	CALL xios_recv_field("sic_limit",sic_mpi)
64	CALL xios_recv_field("sst_limit",sst_mpi)
65	CALL xios_recv_field("rugos_limit",rugos_mpi)
66	CALL xios_recv_field("albedo_limit",albedo_mpi)
67	ENDIF
68	CALL scatter_omp(sic_mpi,sic)
69	CALL scatter_omp(sst_mpi,sst)
70	CALL scatter_omp(rugos_mpi,rugos)
71	CALL scatter_omp(albedo_mpi,albedo)
72
73	ALLOCATE(sic_year(klon,ndays))
74	ALLOCATE(sst_year(klon,ndays))
75	ALLOCATE(rugos_year(klon,ndays))
76	ALLOCATE(albedo_year(klon,ndays))
77	ALLOCATE(pctsrf_t(klon,nbsrf,ndays))
78	ALLOCATE(phy_bil(klon,ndays)); phy_bil=0.0
79
80
81	! sic
82	IF (sic_time_axis_size==lmdep) THEN
83	CALL time_interpolation(ndays,sic,'gregorian',sic_year)
84	ELSE IF (sic_time_axis_size==ndays) THEN
85	sic_year=sic
86	ELSE
87	WRITE(mess,*) 'sic time axis is nor montly, nor daily. sic time interpolation ',&
88	'is requiered but is not currently managed'
89	CALL abort_physic('create_limit_unstruct',TRIM(mess),1)
90	ENDIF
91
92	sic_year(:,:)=sic_year(:,:)/100. ! convert percent to fraction
93	WHERE(sic_year(:,:)>1.0) sic_year(:,:)=1.0 ! Some fractions have some time large negative values
94	WHERE(sic_year(:,:)<0.0) sic_year(:,:)=0.0 ! probably better to apply alse this filter before horizontal interpolation
95
96	! sst
97	IF (sst_time_axis_size==lmdep) THEN
98	CALL time_interpolation(ndays,sst,'gregorian',sst_year)
99	ELSE IF (sst_time_axis_size==ndays) THEN
100	sst_year=sst
101	ELSE
102	WRITE(mess,*)'sic time axis is nor montly, nor daily. sic time interpolation ',&
103	'is requiered but is not currently managed'
104	CALL abort_physic('create_limit_unstruct',TRIM(mess),1)
105	ENDIF
106	WHERE(sst_year(:,:)<271.38) sst_year(:,:)=271.38
107
108
109	! rugos
110	DO l=1, lmdep
111	WHERE(NINT(zmasq(:))/=1) rugos(:,l)=0.001
112	ENDDO
113	CALL time_interpolation(ndays,rugos,'360_day',rugos_year)
114
115	! albedo
116	CALL time_interpolation(ndays,albedo,'360_day',albedo_year)
117
118
119	DO k=1,ndays
120	fi_ice=sic_year(:,k)
121	WHERE(fi_ice>=1.0 ) fi_ice=1.0
122	WHERE(fi_ice<EPSFRA) fi_ice=0.0
123	pctsrf_t(:,is_ter,k)=pctsrf(:,is_ter) ! land soil
124	pctsrf_t(:,is_lic,k)=pctsrf(:,is_lic) ! land ice
125
126	!! IF (icefile==trim(fcpldsic)) THEN ! SIC=pICE*(1-LIC-TER)
127	!! pctsrf_t(:,is_sic,k)=fi_ice(:)*(1.-pctsrf(:,is_lic)-pctsrf(:,is_ter))
128	!! ELSE IF (icefile==trim(fhistsic)) THEN ! SIC=pICE
129	!! pctsrf_t(:,is_sic,k)=fi_ice(:)
130	!! ELSE ! icefile==famipsic ! SIC=pICE-LIC
131	pctsrf_t(:,is_sic,k)=fi_ice-pctsrf_t(:,is_lic,k)
132	! END IF
133	WHERE(pctsrf_t(:,is_sic,k)<=0) pctsrf_t(:,is_sic,k)=0.
134	WHERE(1.0-zmasq<EPSFRA)
135	pctsrf_t(:,is_sic,k)=0.0
136	pctsrf_t(:,is_oce,k)=0.0
137	ELSEWHERE
138	WHERE(pctsrf_t(:,is_sic,k)>=1.0-zmasq)
139	pctsrf_t(:,is_sic,k)=1.0-zmasq
140	pctsrf_t(:,is_oce,k)=0.0
141	ELSEWHERE
142	pctsrf_t(:,is_oce,k)=1.0-zmasq-pctsrf_t(:,is_sic,k)
143	WHERE(pctsrf_t(:,is_oce,k)<EPSFRA)
144	pctsrf_t(:,is_oce,k)=0.0
145	pctsrf_t(:,is_sic,k)=1.0-zmasq
146	END WHERE
147	END WHERE
148	END WHERE
149	nbad=COUNT(pctsrf_t(:,is_oce,k)<0.0)
150	IF(nbad>0) WRITE(lunout,*) 'pb sous maille pour nb point = ',nbad
151	nbad=COUNT(abs(sum(pctsrf_t(:,:,k),dim=2)-1.0)>EPSFRA)
152	IF(nbad>0) WRITE(lunout,*) 'pb sous surface pour nb points = ',nbad
153	END DO
154
155	ALLOCATE(sst_year_mpi(klon_mpi,ndays))
156	ALLOCATE(rugos_year_mpi(klon_mpi,ndays))
157	ALLOCATE(albedo_year_mpi(klon_mpi,ndays))
158	ALLOCATE(pctsrf_t_mpi(klon_mpi,nbsrf,ndays))
159	ALLOCATE(phy_bil_mpi(klon_mpi,ndays))
160
161	CALL gather_omp(pctsrf_t , pctsrf_t_mpi)
162	CALL gather_omp(sst_year , sst_year_mpi)
163	CALL gather_omp(phy_bil , phy_bil_mpi)
164	CALL gather_omp(albedo_year, albedo_year_mpi)
165	CALL gather_omp(rugos_year , rugos_year_mpi)
166
167	IF (is_omp_master) THEN
168	CALL xios_send_field("foce_limout",pctsrf_t_mpi(:,is_oce,:))
169	CALL xios_send_field("fsic_limout",pctsrf_t_mpi(:,is_sic,:))
170	CALL xios_send_field("fter_limout",pctsrf_t_mpi(:,is_ter,:))
171	CALL xios_send_field("flic_limout",pctsrf_t_mpi(:,is_lic,:))
172	CALL xios_send_field("sst_limout", sst_year_mpi)
173	CALL xios_send_field("bils_limout",phy_bil_mpi)
174	CALL xios_send_field("alb_limout", albedo_year_mpi)
175	CALL xios_send_field("rug_limout", rugos_year_mpi)
176	ENDIF
177	#endif
178	END SUBROUTINE create_limit_unstruct
179
180
181	SUBROUTINE time_interpolation(ndays,field_in,calendar,field_out)
182	USE pchsp_95_m, only: pchsp_95
183	USE pchfe_95_m, only: pchfe_95
184	USE arth_m, only: arth
185	USE dimphy, ONLY : klon
186	USE ioipsl, ONLY : ioget_year_len
187	USE time_phylmdz_mod, ONLY : annee_ref
188	USE mod_phys_lmdz_para
189	IMPLICIT NONE
190	INCLUDE "iniprint.h"
191
192	INTEGER, INTENT(IN) :: ndays
193	REAL, INTENT(IN) :: field_in(klon,lmdep)
194	CHARACTER(LEN=*),INTENT(IN) :: calendar
195	REAL, INTENT(OUT) :: field_out(klon,ndays)
196
197	INTEGER :: ndays_in
198	REAL :: timeyear(lmdep)
199	REAL :: yder(lmdep)
200	INTEGER :: ij,ierr, n_extrap
201	LOGICAL :: skip
202
203	CHARACTER (len = 50) :: modname = 'create_limit_unstruct.time_interpolation'
204	CHARACTER (len = 80) :: abort_message
205
206
207	IF (is_omp_master) ndays_in=year_len(annee_ref, calendar)
208	CALL bcast_omp(ndays_in)
209	IF (is_omp_master) timeyear=mid_months(annee_ref, calendar, lmdep)
210	CALL bcast_omp(timeyear)
211
212	n_extrap = 0
213	skip=.FALSE.
214	DO ij=1,klon
215	yder = pchsp_95(timeyear, field_in(ij, :), ibeg=2, iend=2, vc_beg=0., vc_end=0.)
216	CALL pchfe_95(timeyear, field_in(ij, :), yder, skip, arth(0., real(ndays_in) / ndays, ndays), field_out(ij, :), ierr)
217	if (ierr < 0) then
218	abort_message='error in pchfe_95'
219	CALL abort_physic(modname,abort_message,1)
220	endif
221	n_extrap = n_extrap + ierr
222	END DO
223
224	IF (n_extrap /= 0) then
225	WRITE(lunout,*) "get_2Dfield pchfe_95: n_extrap = ", n_extrap
226	ENDIF
227
228
229	END SUBROUTINE time_interpolation
230	!-------------------------------------------------------------------------------
231	!
232	FUNCTION year_len(y,cal_in)
233	!
234	!-------------------------------------------------------------------------------
235	USE ioipsl, ONLY : ioget_calendar,ioconf_calendar,lock_calendar,ioget_year_len
236	IMPLICIT NONE
237	!-------------------------------------------------------------------------------
238	! Arguments:
239	INTEGER :: year_len
240	INTEGER, INTENT(IN) :: y
241	CHARACTER(LEN=*), INTENT(IN) :: cal_in
242	!-------------------------------------------------------------------------------
243	! Local variables:
244	CHARACTER(LEN=20) :: cal_out ! calendar (for outputs)
245	!-------------------------------------------------------------------------------
246	!--- Getting the input calendar to reset at the end of the function
247	CALL ioget_calendar(cal_out)
248
249	!--- Unlocking calendar and setting it to wanted one
250	CALL lock_calendar(.FALSE.); CALL ioconf_calendar(TRIM(cal_in))
251
252	!--- Getting the number of days in this year
253	year_len=ioget_year_len(y)
254
255	!--- Back to original calendar
256	CALL lock_calendar(.FALSE.); CALL ioconf_calendar(TRIM(cal_out))
257
258	END FUNCTION year_len
259	!
260	!-------------------------------------------------------------------------------
261
262
263	!-------------------------------------------------------------------------------
264	!
265	FUNCTION mid_months(y,cal_in,nm)
266	!
267	!-------------------------------------------------------------------------------
268	USE ioipsl, ONLY : ioget_calendar,ioconf_calendar,lock_calendar,ioget_mon_len
269	IMPLICIT NONE
270	!-------------------------------------------------------------------------------
271	! Arguments:
272	INTEGER, INTENT(IN) :: y ! year
273	CHARACTER(LEN=*), INTENT(IN) :: cal_in ! calendar
274	INTEGER, INTENT(IN) :: nm ! months/year number
275	REAL, DIMENSION(nm) :: mid_months ! mid-month times
276	!-------------------------------------------------------------------------------
277	! Local variables:
278	CHARACTER(LEN=99) :: mess ! error message
279	CHARACTER(LEN=20) :: cal_out ! calendar (for outputs)
280	INTEGER, DIMENSION(nm) :: mnth ! months lengths (days)
281	INTEGER :: m ! months counter
282	INTEGER :: nd ! number of days
283	INTEGER :: k
284	!-------------------------------------------------------------------------------
285	nd=year_len(y,cal_in)
286
287	IF(nm==12) THEN
288
289	!--- Getting the input calendar to reset at the end of the function
290	CALL ioget_calendar(cal_out)
291
292	!--- Unlocking calendar and setting it to wanted one
293	CALL lock_calendar(.FALSE.); CALL ioconf_calendar(TRIM(cal_in))
294
295	!--- Getting the length of each month
296	DO m=1,nm; mnth(m)=ioget_mon_len(y,m); END DO
297
298	!--- Back to original calendar
299	CALL lock_calendar(.FALSE.); CALL ioconf_calendar(TRIM(cal_out))
300
301	ELSE IF(MODULO(nd,nm)/=0) THEN
302	WRITE(mess,'(a,i3,a,i3,a)')'Unconsistent calendar: ',nd,' days/year, but ',&
303	nm,' months/year. Months number should divide days number.'
304	CALL abort_physic('mid_months',TRIM(mess),1)
305
306	ELSE
307	mnth=(/(m,m=1,nm,nd/nm)/)
308	END IF
309
310	!--- Mid-months times
311	mid_months(1)=0.5*REAL(mnth(1))
312	DO k=2,nm
313	mid_months(k)=mid_months(k-1)+0.5*REAL(mnth(k-1)+mnth(k))
314	END DO
315
316	END FUNCTION mid_months
317
318
319	END MODULE create_limit_unstruct_mod

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

Download in other formats:

Original Format