Context Navigation

create_limit_unstruct.F90 @ 3469

Last change on this file since 3469 was 3469, checked in by yann meurdesoif, 5 years ago

Fix problem in limit.nc creation for unstructured case in OpenMP.

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

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

Download in other formats:

Original Format