source: XIOS/dev/branch_openmp/src/test/toy_cmip6_omp.f90 @ 1350

Last change on this file since 1350 was 1350, checked in by yushan, 6 years ago

toy_cmip6_omp from Curie
File size: 12.5 KB
Line 
1PROGRAM toy_cmip6_omp 
2
3  USE xios
4  USE mod_wait
5  IMPLICIT NONE
6  INCLUDE "mpif.h"
7
8  INTEGER,PARAMETER :: il_unit=10
9  INTEGER :: comm, rank, size_loc, ierr
10  INTEGER :: ni,ibegin,iend,nj,jbegin,jend
11  INTEGER :: i,j,l,ts,n, nb_pt, il_run
12  DOUBLE PRECISION :: sypd, timestep_in_seconds, simulated_seconds_per_seconds, elapsed_per_timestep
13  CHARACTER(len=*),PARAMETER :: id="client"
14  CHARACTER(1000):: duration, timestep
15  INTEGER :: start_year,start_month,start_day
16  TYPE(xios_date) :: cdate, edate
17  TYPE(xios_duration)  :: dtime
18  TYPE(xios_context) :: ctx_hdl
19  REAL :: ilon,jlat
20  DOUBLE PRECISION,ALLOCATABLE :: lon_glo(:,:),lat_glo(:,:),lval(:)
21  DOUBLE PRECISION,ALLOCATABLE :: field_A_glo (:,:,:), pressure_glo (:,:,:), height_glo (:,:,:)
22  DOUBLE PRECISION,ALLOCATABLE :: bounds_lon_glo(:,:,:),bounds_lat_glo(:,:,:)
23  DOUBLE PRECISION,ALLOCATABLE :: pressure (:,:,:), height (:,:,:)
24  DOUBLE PRECISION,ALLOCATABLE :: lon(:,:),lat(:,:),lonvalue(:,:)
25  DOUBLE PRECISION,ALLOCATABLE :: bounds_lon(:,:,:),bounds_lat(:,:,:) 
26  DOUBLE PRECISION,ALLOCATABLE :: field_atm_2D(:,:),field_atm_3D(:,:,:),field_srf_2D(:),field_srf_3D(:,:)
27  DOUBLE PRECISION,ALLOCATABLE :: field_atm_2D_miss(:,:)
28  DOUBLE PRECISION,ALLOCATABLE :: field_oce_2D(:,:),field_oce_3D(:,:,:)
29  INTEGER, ALLOCATABLE :: kindex(:)
30  INTEGER :: provided
31
32  INTEGER :: ni_glo, nj_glo,llm
33
34  NAMELIST /param_toy/ ni_glo, nj_glo,llm,timestep,duration,sypd,start_year,start_month,start_day
35
36!!! MPI Initialization
37
38  CALL MPI_INIT_THREAD(3, provided, ierr)
39    if(provided .NE. 3) then
40      print*, "provided thread level = ", provided
41      call MPI_Abort()
42    endif
43
44  CALL init_wait
45
46  CALL MPI_COMM_RANK(MPI_COMM_WORLD,rank,ierr)
47  CALL MPI_COMM_SIZE(MPI_COMM_WORLD,size_loc,ierr)
48  if(rank < size_loc-2) then
49
50!!! Lecture des parametres du run
51
52  OPEN(unit=il_unit, file='param.def',status='old',iostat=ierr)
53  READ (il_unit, nml=param_toy)
54  !PRINT *, ni_glo, nj_glo,llm,duration
55
56  !$omp parallel default(firstprivate)
57
58!!! XIOS Initialization (get the local communicator)
59
60  CALL xios_initialize(id,return_comm=comm)
61
62  CALL MPI_COMM_RANK(comm,rank,ierr)
63  CALL MPI_COMM_SIZE(comm,size_loc,ierr)
64
65
66!!! Initialisation et allocation des coordonnées globales et locales pour la grille réguliÚre
67
68  ALLOCATE (lon_glo(ni_glo,nj_glo),lat_glo(ni_glo,nj_glo))
69  ALLOCATE(bounds_lon_glo(4,ni_glo,nj_glo))
70  ALLOCATE(bounds_lat_glo(4,ni_glo,nj_glo))
71  ALLOCATE (field_A_glo(ni_glo,nj_glo,llm))
72  ALLOCATE (pressure_glo(ni_glo,nj_glo,llm))
73  ALLOCATE (height_glo(ni_glo,nj_glo,llm))
74  ALLOCATE (lval(llm))
75
76  DO j=1,nj_glo
77    DO i=1,ni_glo
78
79      ilon=i-0.5
80      jlat=j-0.5
81
82      lat_glo(i,j)= 90-(jlat*180./nj_glo)
83      lon_glo(i,j)= (ilon*360./ni_glo)
84      !print*, 'i/lon=',i,'lon=',lon_glo(i,j), 'j/lat=',j,'lat=',lat_glo(i,j)
85
86      bounds_lat_glo(1,i,j)= 90-((jlat-0.5)*180./nj_glo)
87      bounds_lon_glo(1,i,j)=((ilon-0.5)*360./ni_glo)
88
89      bounds_lat_glo(2,i,j)= 90-((jlat-0.5)*180./nj_glo)
90      bounds_lon_glo(2,i,j)=((ilon+0.5)*360./ni_glo)
91
92      bounds_lat_glo(3,i,j)= 90-((jlat+0.5)*180./nj_glo)
93      bounds_lon_glo(3,i,j)=((ilon+0.5)*360./ni_glo)
94
95      bounds_lat_glo(4,i,j)= 90-((jlat+0.5)*180./nj_glo)
96      bounds_lon_glo(4,i,j)=((ilon-0.5)*360./ni_glo)
97
98      WHERE (abs(bounds_lat_glo(:,i,j) - 90) < 0.000000001) bounds_lat_glo(:,i,j) = 90
99      WHERE (abs(bounds_lat_glo(:,i,j) + 90) < 0.000000001) bounds_lat_glo(:,i,j) = -90
100
101      DO l=1,llm
102         field_A_glo(i,j,l)=(i-1)+(j-1)*ni_glo+10000*l
103         ! pressure at half levels. First index value is high altitude, low pressure
104         pressure_glo(i,j,l)=((l-0.)/llm)*100000 + (jlat -nj_glo/2.)/nj_glo * 10000
105         height_glo(i,j,l)=(llm-l+0.5)/llm * 15000 + jlat * 100
106      ENDDO
107    ENDDO
108  ENDDO
109  ni=ni_glo ; ibegin=0
110
111  jbegin=0
112  DO n=0,size_loc-1
113    nj=nj_glo/size_loc
114    IF (n<MOD(nj_glo,size_loc)) nj=nj+1
115    IF (n==rank) exit
116    jbegin=jbegin+nj
117  ENDDO
118
119  iend=ibegin+ni-1 ; jend=jbegin+nj-1
120
121  ALLOCATE(lon(ni,nj),lat(ni,nj),lonvalue(ni,nj))
122  ALLOCATE(bounds_lon(4,ni,nj))
123  ALLOCATE(bounds_lat(4,ni,nj))
124  lon(:,:)=lon_glo(ibegin+1:iend+1,jbegin+1:jend+1)
125  lat(:,:)=lat_glo(ibegin+1:iend+1,jbegin+1:jend+1)
126  bounds_lon(:,:,:)=bounds_lon_glo(:,ibegin+1:iend+1,jbegin+1:jend+1)
127  bounds_lat(:,:,:)=bounds_lat_glo(:,ibegin+1:iend+1,jbegin+1:jend+1)
128
129
130  DO i=1,llm
131    lval(i)=i
132 ENDDO
133 
134
135 
136!###########################################################################
137! Contexte ATM
138!###########################################################################
139 ALLOCATE(field_atm_2D(0:ni+1,-1:nj+2),field_atm_3D(0:ni+1,-1:nj+2,llm))
140 ALLOCATE(field_atm_2D_miss(0:ni+1,-1:nj+2))
141 ALLOCATE(pressure(0:ni+1,-1:nj+2,llm))
142 ALLOCATE(height(0:ni+1,-1:nj+2,llm))
143  field_atm_2D(1:ni,1:nj)=field_A_glo(ibegin+1:iend+1,jbegin+1:jend+1,1)
144  field_atm_2D_miss(1:ni,1:nj)=field_A_glo(ibegin+1:iend+1,jbegin+1:jend+1,1)
145  field_atm_3D(1:ni,1:nj,:)=field_A_glo(ibegin+1:iend+1,jbegin+1:jend+1,:)
146  pressure(1:ni,1:nj,:)=pressure_glo(ibegin+1:iend+1,jbegin+1:jend+1,:)
147  height(1:ni,1:nj,:)=height_glo(ibegin+1:iend+1,jbegin+1:jend+1,:)
148
149  CALL xios_context_initialize("arpsfx",comm)
150  CALL xios_define_calendar("Gregorian", &
151     start_date=xios_date(start_year,start_month,start_day,0,0,0), &
152     time_origin=xios_date(1850,1,1,0,0,0))
153
154  !write(0,*) 'atm context initialized' ; call flush(0)
155  CALL xios_get_handle("arpsfx",ctx_hdl)
156  CALL xios_set_current_context(ctx_hdl)
157
158  CALL xios_set_axis_attr("axis_atm",n_glo=llm ,value=lval) ;
159
160  CALL xios_set_domain_attr("domain_atm",ni_glo=ni_glo, nj_glo=nj_glo, ibegin=ibegin, &
161       ni=ni,jbegin=jbegin,nj=nj, type='rectilinear')
162  CALL xios_set_domain_attr("domain_atm",data_dim=2, data_ibegin=-1, &
163       data_ni=ni+2, data_jbegin=-2, data_nj=nj+4)
164  CALL xios_set_domain_attr("domain_atm",lonvalue_2D=lon,latvalue_2D=lat)
165  CALL xios_set_domain_attr("domain_atm", nvertex=4, bounds_lon_2d=bounds_lon, bounds_lat_2d=bounds_lat)
166  print *,'latmax/min=',minval(lat),maxval(lat)
167  print *,'latmax/min bounds=',minval(bounds_lat),maxval(bounds_lat)
168  print *,'lonmax/min=',minval(lon),maxval(lon)
169  print *,'lonmax/min bounds=',minval(bounds_lon),maxval(bounds_lon)
170
171
172!!! Definition du timestep
173
174  CALL xios_get_start_date(cdate)
175  edate=cdate+xios_duration_convert_from_string(duration)
176  dtime=xios_duration_convert_from_string(timestep)
177  CALL xios_set_timestep(timestep=dtime)
178
179!!! Fin de la definition du contexte
180
181  !CALL xios_close_context_definition()
182
183!!! Calcul de temps elaps par seconde pour respecter le SYPD (hyp : pas de délai d'I/O)
184 
185  timestep_in_seconds=xios_date_convert_to_seconds(cdate+dtime) - xios_date_convert_to_seconds(cdate)
186  simulated_seconds_per_seconds=sypd * 365 
187  elapsed_per_timestep=timestep_in_seconds/simulated_seconds_per_seconds ! in seconds
188
189!###########################################################################
190! Contexte SRF
191!###########################################################################
192
193!!! Initialisation des coordonnées globales et locales pour la grille indexee (1 point sur 2)
194
195    nb_pt=ni*nj/2
196    ALLOCATE(kindex(nb_pt),field_srf_2D(nb_pt),field_srf_3D(nb_pt,llm))
197    DO i=1,nb_pt
198      kindex(i)=2*i-1
199    ENDDO
200    field_srf_2D(1:nb_pt)=RESHAPE(field_A_glo(ibegin+1:iend+1:2,jbegin+1:jend+1,1),(/ nb_pt /))
201    field_srf_3D(1:nb_pt,:)=RESHAPE(field_A_glo(ibegin+1:iend+1:2,jbegin+1:jend+1,:),(/ nb_pt,llm /))
202
203  !CALL xios_context_initialize("surface",comm)
204  !CALL xios_get_handle("surface",ctx_hdl)
205  !CALL xios_set_current_context(ctx_hdl)
206
207  CALL xios_set_axis_attr("axis_srf",n_glo=llm ,value=lval)
208  CALL xios_set_domain_attr("domain_srf",ni_glo=ni_glo, nj_glo=nj_glo, ibegin=ibegin, ni=ni,jbegin=jbegin,nj=nj, type='rectilinear')
209  CALL xios_set_domain_attr("domain_srf",data_dim=1, data_ibegin=0, data_ni=nb_pt)
210  CALL xios_set_domain_attr("domain_srf",data_i_index=kindex)
211  CALL xios_set_domain_attr("domain_srf",lonvalue_2D=lon,latvalue_2D=lat)
212  CALL xios_set_domain_attr("domain_srf", nvertex=4, bounds_lon_2d=bounds_lon, bounds_lat_2d=bounds_lat)
213
214!!! Definition du timestep
215
216  !dtime%second=timestep
217  !dtime=xios_duration_convert_from_string(timestep)
218  !CALL xios_set_timestep(timestep=dtime)
219
220!!! Fin de la definition du contexte SRF
221
222  CALL xios_close_context_definition()
223  !write(0,*) 'srf context def closed' ; call flush(0)
224
225!###########################################################################
226! Contexte OCE
227!###########################################################################
228  ALLOCATE(field_oce_2D(0:ni+1,-1:nj+2),field_oce_3D(0:ni+1,-1:nj+2,llm))
229  field_oce_2D(1:ni,1:nj)=field_A_glo(ibegin+1:iend+1,jbegin+1:jend+1,1)
230  field_oce_3D(1:ni,1:nj,:)=field_A_glo(ibegin+1:iend+1,jbegin+1:jend+1,:)
231
232  CALL xios_context_initialize("nemo",comm)
233  CALL xios_define_calendar("Gregorian", &
234     start_date=xios_date(start_year,start_month,start_day,0,0,0), &
235     time_origin=xios_date(1850,1,1,0,0,0))
236  CALL xios_get_handle("nemo",ctx_hdl)
237  CALL xios_set_current_context(ctx_hdl)
238
239  CALL xios_set_axis_attr("axis_oce",n_glo=llm ,value=lval) ;
240
241  CALL xios_set_domain_attr("domain_oce",ni_glo=ni_glo, nj_glo=nj_glo, ibegin=ibegin, ni=ni,jbegin=jbegin,nj=nj, type='curvilinear')
242  CALL xios_set_domain_attr("domain_oce",data_dim=2, data_ibegin=-1, data_ni=ni+2, data_jbegin=-2, data_nj=nj+4)
243  CALL xios_set_domain_attr("domain_oce",lonvalue_2D=lon,latvalue_2D=lat)
244
245
246!!! Definition du timestep
247
248  dtime=xios_duration_convert_from_string(timestep)
249  CALL xios_set_timestep(timestep=dtime)
250
251  CALL xios_close_context_definition()
252
253!####################################################################################
254!!! Boucle temporelle
255!####################################################################################
256  ts=1
257  cdate=cdate+dtime
258  DO while ( cdate <= edate )
259
260      CALL xios_get_handle("arpsfx",ctx_hdl)
261      CALL xios_set_current_context(ctx_hdl)
262
263!!! Mise a jour du pas de temps
264
265      CALL xios_update_calendar(ts)
266
267!!! On donne la valeur du champ atm
268
269      !print *,'sending field_atm_2d at timestep',ts
270      CALL xios_send_field("field_atm_scalar",field_atm_2D(1,1)+ts)
271      CALL xios_send_field("field_atm_1D",field_atm_3D(1,1,:)+ts)
272      CALL xios_send_field("field_atm_2D",field_atm_2D+ts)
273      CALL xios_send_field("field_atm_3D",field_atm_3D+ts)
274      CALL xios_send_field("pressure" ,pressure)
275      CALL xios_send_field("height"   ,height)
276      if (mod(ts,2)==0) then
277         CALL xios_send_field("field_sub",field_atm_2D+ts)
278      endif
279      !! On crée un champ avec des missings qui bougent
280      !! dans le temps : un bande verticale de 1 à ni-3
281      field_atm_2D_miss(:,:)= field_atm_2D(:,:)+ts
282      field_atm_2D_miss(mod(ts,ni-3)+1,:)=1.e+20
283      CALL xios_send_field("field_miss",field_atm_2D_miss)
284     
285!!! On change de contexte
286
287      !CALL xios_get_handle("surface",ctx_hdl)
288      !CALL xios_set_current_context(ctx_hdl)
289
290!!! Mise a jour du pas de temps
291
292      !CALL xios_update_calendar(ts)
293
294!!! On donne la valeur du champ srf
295
296      CALL xios_send_field("field_srf_2D",field_srf_2D)
297      CALL xios_send_field("field_srf_3D",field_srf_3D)
298
299!!! On change de contexte
300
301      CALL xios_get_handle("nemo",ctx_hdl)
302      CALL xios_set_current_context(ctx_hdl)
303
304!!! Mise a jour du pas de temps
305
306      CALL xios_update_calendar(ts)
307
308!!! On donne la valeur du champ oce
309
310      CALL xios_send_field("field_oce_scalar",field_oce_2D(1,1)+ts)
311      CALL xios_send_field("field_oce_grid_2D",field_oce_2D)
312      CALL xios_send_field("field_oce_grid_3D",field_oce_3D)
313
314      CALL wait_us(int(elapsed_per_timestep*1.e6))   ! micro-secondes
315      cdate=cdate+dtime
316      ts=ts+1
317
318   ENDDO
319
320!####################################################################################
321!!! Finalisation
322!####################################################################################
323
324!!! Fin des contextes
325
326    CALL xios_context_finalize()
327    CALL xios_get_handle("arpsfx",ctx_hdl)
328    CALL xios_set_current_context(ctx_hdl)
329    CALL xios_context_finalize()
330    !CALL xios_get_handle("surface",ctx_hdl)
331    !CALL xios_set_current_context(ctx_hdl)
332    !CALL xios_context_finalize()
333
334    DEALLOCATE(lon, lat, lonvalue, field_atm_2D, field_atm_3D)
335    DEALLOCATE(pressure,height,pressure_glo,height_glo)
336    DEALLOCATE(kindex, field_srf_2D, field_srf_3D)
337    DEALLOCATE(field_oce_2D, field_oce_3D)
338    DEALLOCATE(lon_glo,lat_glo,field_A_glo,lval)
339
340!!! Fin de XIOS
341
342    CALL xios_finalize()
343    print*, "xios finalize OK", rank, size_loc
344   
345    !$omp master
346    !call MPI_Barrier(comm)
347    CALL MPI_COMM_FREE(comm, ierr)
348    !$omp end master
349
350    !$omp barrier
351   
352    !$omp end parallel
353
354
355    else
356
357    CALL xios_init_server
358    print *, "Server : xios_finalize "
359 
360    endif 
361   
362
363  CALL MPI_FINALIZE(ierr)
364
365END PROGRAM toy_cmip6_omp
366
367
368
369
370
371
Note: See TracBrowser for help on using the repository browser.