source: trunk/LMDZ.MARS/util/startarchive2icosa/start_archive2icosa.f90 @ 3648

Last change on this file since 3648 was 2941, checked in by romain.vande, 3 years ago

Mars Dynamico : Adapt routines to the new version of Dynamico and the subslope parametrization.
File size: 24.7 KB
Line 
1PROGRAM start_archive2icosa
2
3  USE xios
4  USE mod_wait
5  USE netcdf
6 
7  IMPLICIT NONE
8  INCLUDE "mpif.h"
9  INTEGER :: rank
10  INTEGER :: size
11  INTEGER :: ierr
12
13  CHARACTER(len=*),PARAMETER :: id="client"
14  INTEGER :: comm
15  TYPE(xios_duration) :: dtime
16  CHARACTER(len=15) :: calendar_type
17  TYPE(xios_context) :: ctx_hdl
18
19  INTEGER :: n,l
20  INTEGER :: src_ibegin, src_iend, src_topo_ibegin, src_topo_iend
21  INTEGER :: src_jbegin, src_jend, src_topo_jbegin, src_topo_jend
22  INTEGER :: src_ni, src_ni_glo, src_topo_ni, src_topo_ni_glo
23  INTEGER :: src_nj, src_nj_glo, src_topo_nj, src_topo_nj_glo
24  INTEGER :: src_nlev ! number of vertical layers
25  INTEGER :: src_nq=11 ! number of tracers
26  INTEGER :: src_nt=1 ! number of time steps
27  INTEGER :: src_nsoil=18 ! number of soil layers
28  INTEGER :: src_nslope=1 ! number of subslope
29  DOUBLE PRECISION,ALLOCATABLE :: lev_values(:) ! vertical axis
30  DOUBLE PRECISION,ALLOCATABLE :: lev_p1_values(:) ! vertical axis
31  DOUBLE PRECISION,ALLOCATABLE :: nq_values(:) ! tracer # axis
32  DOUBLE PRECISION,ALLOCATABLE :: soil_layers_values(:) ! soil axis
33  DOUBLE PRECISION,ALLOCATABLE :: subslope_values(:) ! subslope axis
34  DOUBLE PRECISION,ALLOCATABLE :: src_lon(:) ! mesh center coordinate
35  DOUBLE PRECISION,ALLOCATABLE :: src_lat(:)
36  DOUBLE PRECISION,ALLOCATABLE :: src_ap(:)
37  DOUBLE PRECISION,ALLOCATABLE :: src_bp(:)
38  DOUBLE PRECISION,ALLOCATABLE :: src_controle(:)
39  DOUBLE PRECISION,ALLOCATABLE :: src_field_2D(:,:)
40  DOUBLE PRECISION,ALLOCATABLE :: src_field_2D_subslope(:,:,:)
41  DOUBLE PRECISION,ALLOCATABLE :: src_pk(:,:)
42  DOUBLE PRECISION,ALLOCATABLE :: src_field_3D(:,:,:)
43  DOUBLE PRECISION,ALLOCATABLE :: src_pressure(:,:,:)
44  DOUBLE PRECISION,ALLOCATABLE :: src_theta_rhodz(:,:,:)
45  DOUBLE PRECISION,ALLOCATABLE :: src_topo_lon(:) ! mesh center coordinate
46  DOUBLE PRECISION,ALLOCATABLE :: src_topo_lat(:)
47  DOUBLE PRECISION,ALLOCATABLE :: src_topo(:,:)
48  DOUBLE PRECISION,ALLOCATABLE :: src_field_3D_soil(:,:,:) ! 3D grid in soil, for tsoil et inertia
49  DOUBLE PRECISION,ALLOCATABLE :: src_field_3D_soil_subslope(:,:,:,:) ! 3D grid in soil, for tsoil et inertia
50  DOUBLE PRECISION,ALLOCATABLE :: src_soil_layers(:) ! soil_depth
51  DOUBLE PRECISION,ALLOCATABLE :: src_field_3D_p1(:,:,:) !q2 field
52 
53  CHARACTER(LEN=*),PARAMETER :: src_file="start_archive_nc4.nc"
54  CHARACTER(LEN=*),PARAMETER :: src_topo_file="surface_nc4.nc"
55  CHARACTER(LEN=*),PARAMETER :: dst_coord_file="start_icosa_ref.nc"
56  CHARACTER(LEN=*),PARAMETER :: src_controle_file="startphy_icosa_ref.nc"
57  CHARACTER(LEN=*),PARAMETER :: output_start_file="start_icosa_prefinalize.nc"
58  CHARACTER(LEN=*),PARAMETER :: output_startfi_file="startfi_prefinalize.nc"
59  DOUBLE PRECISION,ALLOCATABLE :: dst_lon(:),dst_lat(:)
60  DOUBLE PRECISION,ALLOCATABLE :: dst_boundslon(:,:) ! mesh corner coordinates
61  DOUBLE PRECISION,ALLOCATABLE :: dst_boundslat(:,:)
62  INTEGER :: dst_ibegin !, dst_iend
63  INTEGER :: dst_ni, dst_ni_glo
64  INTEGER :: dst_nvertex
65  INTEGER :: ncid
66  INTEGER :: dimids(5)
67  INTEGER :: varid
68 
69  INTEGER :: div, remain
70  INTEGER :: ts ! time step #
71  DOUBLE PRECISION,PARAMETER :: pi=acos(-1.d0)
72  DOUBLE PRECISION :: gravity,kappa,preff
73! Tracers
74  CHARACTER(LEN=11)          :: i_trac,format_string
75  INTEGER                    :: i
76
77!!! MPI Initialization
78  CALL MPI_INIT(ierr)
79  CALL init_wait
80
81!!! XIOS Initialization (get the local communicator)
82  CALL xios_initialize(id,return_comm=comm)
83! get local rank of MPI process
84  CALL MPI_COMM_RANK(comm,rank,ierr)
85! get total number of MPI processes
86  CALL MPI_COMM_SIZE(comm,size,ierr)
87
88!!! Open files and load sizes and coordinates
89  ierr=NF90_OPEN(src_topo_file, NF90_NOWRITE, ncid)
90  ierr=NF90_INQ_VARID(ncid,"zMOL",varid)
91  ierr=NF90_INQUIRE_VARIABLE(ncid, varid,dimids=dimids)
92  write(*,*) "rank=",rank,"dimids=",dimids
93  ierr=NF90_INQUIRE_DIMENSION(ncid, dimids(1), len=src_topo_ni_glo)
94  ierr=NF90_INQUIRE_DIMENSION(ncid, dimids(2), len=src_topo_nj_glo)
95  write(*,*) "rank=",rank," src_topo_ni_glo=",src_topo_ni_glo ! longitude
96  write(*,*) "rank=",rank," src_topo_nj_glo=",src_topo_nj_glo ! latitude
97
98! assume domain splitup with MPI only along latitudes
99  src_topo_ni=src_topo_ni_glo
100  src_topo_ibegin=0
101  src_topo_iend=src_topo_ibegin+src_ni-1
102  write(*,*) "rank=",rank," src_topo_ni=",src_topo_ni
103 
104  src_topo_jbegin=0
105  DO n=0,size-1
106    src_topo_nj=src_topo_nj_glo/size
107    IF (n<MOD(src_topo_nj_glo,size)) src_topo_nj=src_topo_nj+1
108    IF (n==rank) exit
109    src_topo_jbegin=src_topo_jbegin+src_topo_nj
110  ENDDO
111  src_topo_jend=src_topo_jbegin+src_topo_nj-1
112  write(*,*) "rank=",rank," src_topo_nj=",src_topo_nj, &
113             " src_topo_jbegin=",src_topo_jbegin
114
115  ALLOCATE(src_topo_lon(src_topo_ni))
116  ALLOCATE(src_topo_lat(src_topo_nj))
117  ALLOCATE(src_topo(src_topo_ni,src_topo_nj))
118
119! load src_topo_lon and src_topo_lat
120  ierr=NF90_INQ_VARID(ncid,"longitude",varid)
121  ierr=NF90_GET_VAR(ncid,varid, src_topo_lon, &
122                    start=(/src_topo_ibegin+1/),count=(/src_topo_ni/))
123  WRITE(*,*) rank,":src_topo_lon(1:2)=",src_topo_lon(1:2)
124  ierr=NF90_INQ_VARID(ncid,"latitude",varid)
125  ierr=NF90_GET_VAR(ncid,varid, src_topo_lat, &
126                    start=(/src_topo_jbegin+1/),count=(/src_topo_nj/))
127  WRITE(*,*) rank,":src_topo_lat(1:2)=",src_topo_lat(1:2)
128
129! from start_archive.nc file
130  ierr=NF90_OPEN(src_file, NF90_NOWRITE, ncid)
131  ierr=NF90_INQ_VARID(ncid,"temp",varid)
132  ierr=NF90_INQUIRE_VARIABLE(ncid, varid,dimids=dimids)
133  write(*,*) "rank=",rank,"dimids=",dimids
134  ierr=NF90_INQUIRE_DIMENSION(ncid, dimids(1), len=src_ni_glo)
135  ierr=NF90_INQUIRE_DIMENSION(ncid, dimids(2), len=src_nj_glo)
136  ierr=NF90_INQUIRE_DIMENSION(ncid, dimids(3), len=src_nlev)
137  write(*,*) "rank=",rank," src_ni_glo=",src_ni_glo ! longitude
138  write(*,*) "rank=",rank," src_nj_glo=",src_nj_glo ! latitude
139  write(*,*) "rank=",rank," src_nlev=",src_nlev ! number of vertical layers
140  write(*,*) "rank=",rank," src_nq=",src_nq ! number of tracers
141! soil_depth with tsoil variable
142  ierr=NF90_INQ_VARID(ncid,"tsoil",varid)
143  ierr=NF90_INQUIRE_VARIABLE(ncid, varid,dimids=dimids)
144  write(*,*) "rank=",rank,"dimids tsoil =",dimids
145  ierr=NF90_INQUIRE_DIMENSION(ncid, dimids(3), len=src_nsoil)
146 ierr=NF90_INQUIRE_DIMENSION(ncid, dimids(4), len=src_nslope)
147  write(*,*) "rank=",rank," src_nsoil=",src_nsoil ! number of vertical layers
148  write(*,*) "rank=",rank," src_nslope=",src_nslope ! number of subslope
149
150! assume domain splitup with MPI only along latitudes
151  src_ni=src_ni_glo
152  src_ibegin=0
153  src_iend=src_ibegin+src_ni-1
154  write(*,*) "rank=",rank," src_ni=",src_ni
155 
156  src_jbegin=0
157  DO n=0,size-1
158    src_nj=src_nj_glo/size
159    IF (n<MOD(src_nj_glo,size)) src_nj=src_nj+1
160    IF (n==rank) exit
161    src_jbegin=src_jbegin+src_nj
162  ENDDO
163  src_jend=src_jbegin+src_nj-1
164  write(*,*) "rank=",rank," src_nj=",src_nj," src_jbegin=",src_jbegin
165
166  ALLOCATE(src_lon(src_ni))
167  ALLOCATE(src_lat(src_nj))
168  ALLOCATE(src_field_2D(src_ni,src_nj))
169  ALLOCATE(src_field_2D_subslope(src_ni,src_nj,src_nslope))
170  ALLOCATE(src_pk(src_ni,src_nj))
171  ALLOCATE(src_field_3D(src_ni,src_nj,src_nlev))
172  ALLOCATE(src_field_3D_p1(src_ni,src_nj,src_nlev+1))
173  ALLOCATE(src_field_3D_soil(src_ni,src_nj,src_nsoil))
174  ALLOCATE(src_field_3D_soil_subslope(src_ni,src_nj,src_nsoil,src_nslope))
175!  ALLOCATE(src_field_soil_layers(src_nsoil))
176  ALLOCATE(src_pressure(src_ni,src_nj,src_nlev+1))
177  ALLOCATE(src_theta_rhodz(src_ni,src_nj,src_nlev))
178
179! load src_lon and src_lat
180  ierr=NF90_INQ_VARID(ncid,"rlonv",varid)
181  ierr=NF90_GET_VAR(ncid,varid, src_lon, &
182                    start=(/src_ibegin+1/),count=(/src_ni/))
183! convert rad to deg
184  src_lon(1:src_ni)=src_lon(1:src_ni)*(180.d0/pi)
185  WRITE(*,*) rank,":src_lon=",src_lon
186  ierr=NF90_INQ_VARID(ncid,"rlatu",varid)
187  ierr=NF90_GET_VAR(ncid,varid, src_lat, &
188                    start=(/src_jbegin+1/),count=(/src_nj/))
189! convert rad to deg
190  src_lat(1:src_nj)=src_lat(1:src_nj)*(180.d0/pi)
191  WRITE(*,*) rank,":src_lat=",src_lat
192
193! load ap, bp and controle
194  ALLOCATE(src_ap(src_nlev+1),src_bp(src_nlev+1))
195  ierr=NF90_INQ_VARID(ncid,"ap",varid)
196  ierr=NF90_GET_VAR(ncid,varid,src_ap)
197  WRITE(*,*) rank,":src_ap(1:5)=",src_ap(1:5)
198  ierr=NF90_INQ_VARID(ncid,"bp",varid)
199  ierr=NF90_GET_VAR(ncid,varid,src_bp)
200  WRITE(*,*) rank,":src_bp(1:5)=",src_bp(1:5)
201 
202! controle is taken in startphy_icosa_ref as start_archive_nc4 is too old
203  ierr=NF90_OPEN(src_controle_file, NF90_NOWRITE, ncid)
204  ALLOCATE(src_controle(100))
205  ierr=NF90_INQ_VARID(ncid,"controle",varid)
206  ierr=NF90_GET_VAR(ncid,varid,src_controle)
207  ! day_ini set to 0 as startphy_icosa_ref is a restart
208  src_controle(3)=0
209  WRITE(*,*) rank,":src_controle(1:5)=",src_controle(1:5)
210  gravity=src_controle(7)
211  WRITE(*,*) rank,":gravity=",gravity
212  kappa=src_controle(9)
213  WRITE(*,*) rank,":kappa=",kappa
214!  preff=src_controle(18)
215! AD: Warning preff is hardcoded because not in controle of restartfi
216  preff=610.
217  WRITE(*,*) rank,":preff=",preff
218
219! destination coordinates
220  ierr=NF90_OPEN(dst_coord_file, NF90_NOWRITE, ncid)
221  ierr=NF90_INQ_VARID(ncid,"bounds_lon_mesh",varid)
222  ierr=NF90_INQUIRE_VARIABLE(ncid, varid,dimids=dimids)
223  ierr=NF90_INQUIRE_DIMENSION(ncid, dimids(1), len=dst_nvertex)
224  ierr=NF90_INQUIRE_DIMENSION(ncid, dimids(2), len=dst_ni_glo)
225  write(*,*) "rank=",rank," dst_nvertex=",dst_nvertex ! vertex
226  write(*,*) "rank=",rank," dst_ni_glo=",dst_ni_glo ! vertex boundaries
227
228! evenly split into MPI domains
229  div    = dst_ni_glo/size
230  remain = MOD( dst_ni_glo, size )
231  IF (rank < remain) THEN
232    dst_ni=div+1 ;
233    dst_ibegin=rank*(div+1) ;
234  ELSE
235    dst_ni=div ;
236    dst_ibegin= remain * (div+1) + (rank-remain) * div ;
237  ENDIF
238  write(*,*) "rank=",rank," dst_ni=",dst_ni
239
240  ALLOCATE(dst_lon(dst_ni))
241  ALLOCATE(dst_lat(dst_ni))
242  ALLOCATE(dst_boundslon(dst_nvertex,dst_ni))
243  ALLOCATE(dst_boundslat(dst_nvertex,dst_ni))
244
245! load dst_lon, dst_lat, dst_boundslon and dst_boundslat
246  ierr=NF90_INQ_VARID(ncid,"lon_mesh",varid)
247  ierr=NF90_GET_VAR(ncid,varid, dst_lon, &
248                    start=(/dst_ibegin+1/), count=(/dst_ni/))
249  WRITE(*,*) rank,":dst_lon(1:5)=",dst_lon(1:5)
250  ierr=NF90_INQ_VARID(ncid,"lat_mesh",varid)
251  ierr=NF90_GET_VAR(ncid,varid, dst_lat, &
252                    start=(/dst_ibegin+1/), count=(/dst_ni/))
253  WRITE(*,*) rank,":dst_lat(1:5)=",dst_lat(1:5)
254  ierr=NF90_INQ_VARID(ncid,"bounds_lon_mesh",varid)
255  ierr=NF90_GET_VAR(ncid,varid,dst_boundslon, &
256                    start=(/1,dst_ibegin+1/), count=(/dst_nvertex,dst_ni/))
257  WRITE(*,*) rank,":dst_boundslon(:,1:2)=",dst_boundslon(:,1:2)
258  ierr=NF90_INQ_VARID(ncid,"bounds_lat_mesh",varid)
259  ierr=NF90_GET_VAR(ncid,varid, dst_boundslat, &
260                    start=(/1,dst_ibegin+1/), count=(/dst_nvertex,dst_ni/))
261  WRITE(*,*) rank,":dst_boundslat(:,1:2)=",dst_boundslat(:,1:2)
262
263
264! Initialize XIOS context
265  WRITE(*,*) rank,":CALL xios_context_initialize()"
266  CALL xios_context_initialize("test",comm)
267  CALL xios_get_handle("test",ctx_hdl)
268  CALL xios_set_current_context(ctx_hdl)
269
270! Set XIOS calendar and timestep
271  CALL xios_get_calendar_type(calendar_type)
272  WRITE(*,*) rank,":calendar_type = ", calendar_type
273  dtime%second = 3600
274  CALL xios_set_timestep(dtime)
275
276! Set axes
277  ! vertical atm axis
278  ALLOCATE(lev_values(src_nlev))
279  lev_values=(/ (l,l=1,src_nlev) /)
280  CALL xios_set_axis_attr("lev",n_glo=src_nlev,value=lev_values)
281  ! lev+1 for q2
282  ALLOCATE(lev_p1_values(src_nlev+1))
283  lev_p1_values=(/ (l,l=1,src_nlev+1) /)
284  CALL xios_set_axis_attr("lev_p1",n_glo=src_nlev+1,value=lev_p1_values)
285  ! tracers axis
286  ALLOCATE(nq_values(src_nq))
287  nq_values=(/(l,l=1,src_nq)/)
288  CALL xios_set_axis_attr("nq",n_glo=src_nq,value=nq_values)
289  ! soil layers axis
290  ALLOCATE(soil_layers_values(src_nsoil))
291  soil_layers_values=(/ (l,l=1,src_nsoil) /)
292  CALL xios_set_axis_attr("soil_layers",n_glo=src_nsoil,value=soil_layers_values)
293  ALLOCATE(subslope_values(src_nslope))
294  subslope_values=(/ (l,l=1,src_nslope) /)
295  CALL xios_set_axis_attr("subslope",n_glo=src_nslope,value=subslope_values)
296
297! Set domains
298  CALL xios_set_domain_attr("src_domain_regular", &
299                            ni_glo=src_ni_glo, nj_glo=src_nj_glo, &
300                            ibegin=src_ibegin, ni=src_ni, &
301                            jbegin=src_jbegin, nj=src_nj, &
302                            type='rectilinear')
303  CALL xios_set_domain_attr("src_domain_regular", &
304                             data_dim=2, &
305                             data_ibegin=0, data_ni=src_ni, &
306                             data_jbegin=0, data_nj=src_nj)
307  CALL xios_set_domain_attr("src_domain_regular", &
308                            lonvalue_1D=src_lon, &
309                            latvalue_1D=src_lat)
310
311  CALL xios_set_domain_attr("src_topo_domain_regular", &
312                            ni_glo=src_topo_ni_glo, nj_glo=src_topo_nj_glo, &
313                            ibegin=src_topo_ibegin, ni=src_topo_ni, &
314                            jbegin=src_topo_jbegin, nj=src_topo_nj, &
315                            type='rectilinear')
316  CALL xios_set_domain_attr("src_topo_domain_regular", &
317                             data_dim=2, &
318                             data_ibegin=0, data_ni=src_topo_ni, &
319                             data_jbegin=0, data_nj=src_topo_nj)
320  CALL xios_set_domain_attr("src_topo_domain_regular", &
321                            lonvalue_1D=src_topo_lon, &
322                            latvalue_1D=src_topo_lat)
323
324  CALL xios_set_domain_attr("src_domain_regular_clean", &
325                            ni_glo=src_ni_glo-1, nj_glo=src_nj_glo, &
326                            ibegin=src_ibegin, ni=src_ni-1, &
327                            jbegin=src_jbegin, nj=src_nj, &
328                            type='rectilinear')
329  CALL xios_set_domain_attr("src_domain_regular_clean", &
330                             data_dim=2, &
331                             data_ibegin=0, data_ni=src_ni-1, &
332                             data_jbegin=0, data_nj=src_nj)
333  CALL xios_set_domain_attr("src_domain_regular_clean", &
334                            lonvalue_1D=src_lon(1:src_ni-1), &
335                            latvalue_1D=src_lat)
336
337  CALL xios_set_domain_attr("dst_domain_unstructured", &
338                            ni_glo=dst_ni_glo, &
339                            ibegin=dst_ibegin, &
340                            ni=dst_ni, &
341                            type="unstructured")
342  CALL xios_set_domain_attr("dst_domain_unstructured", &
343                            lonvalue_1D=dst_lon, &
344                            latvalue_1D=dst_lat, &
345                            bounds_lon_1D=dst_boundslon, &
346                            bounds_lat_1D=dst_boundslat, &
347                            nvertex=dst_nvertex)
348
349! Finalize XIOS context definition
350  WRITE(*,*) rank,":CALL xios_close_context_definition()"
351  CALL xios_close_context_definition()
352  CALL xios_get_handle("test",ctx_hdl)
353  CALL xios_set_current_context(ctx_hdl)
354
355! Temporal loop
356  DO ts=1,src_nt
357    WRITE(*,*) rank,":ts=",ts
358    ! Update calendar
359    CALL xios_update_calendar(ts)
360
361    ! Topography
362    CALL xios_recv_field("zMOL",src_topo)
363    WRITE(*,*) rank,":topo(1:2,1:3)=",src_topo(1:2,1:3)
364    ! Send surface geopotential
365    CALL xios_send_field("topo",src_topo(:,:)*gravity*1000.)
366
367    ! Surface pressure:
368    !! get data using XIOS:
369    CALL xios_recv_field("src_ps",src_field_2D)
370    WRITE(*,*) rank,":src_ps(1:2,1:3)=",src_field_2D(1:2,1:3)
371    !! write data using XIOS
372    CALL xios_send_field("ps_clean",src_field_2D(1:src_ni-1,1:src_nj))
373
374    ! compute inter-layer pressures
375    DO l=1,src_nlev+1
376      src_pressure(:,:,l) = src_ap(l)+src_bp(l)*src_field_2D(:,:)
377    ENDDO
378   
379    ! surface temperature:
380    CALL xios_recv_field("src_tsurf",src_field_2D_subslope)
381    WRITE(*,*) rank,":src_tsurf(1:2,1:3,1)=",src_field_2D_subslope(1:2,1:3,1)
382    CALL xios_send_field("tsurf_clean",src_field_2D_subslope(1:src_ni-1,1:src_nj,:))
383
384    ! co2 ice coverage:
385!    CALL xios_recv_field("src_co2ice",src_field_2D)
386!    WRITE(*,*) rank,":src_co2ice(1:2,1:3)=",src_field_2D(1:2,1:3)
387!    CALL xios_send_field("co2ice_clean",src_field_2D(1:src_ni-1,1:src_nj))
388
389    ! emissivity:
390    CALL xios_recv_field("src_emis",src_field_2D_subslope)
391    WRITE(*,*) rank,":src_emis(1:2,1:3,1)=",src_field_2D_subslope(1:2,1:3,1)
392    CALL xios_send_field("emis_clean",src_field_2D_subslope(1:src_ni-1,1:src_nj,:))
393
394    ! q2 : q2surf for first layer of q2, the rest is q2atm
395    CALL xios_recv_field("src_q2surf",src_field_2D)
396    WRITE(*,*) rank,":src_q2(1:2,1:3)=",src_field_2D(1:2,1:3)
397    CALL xios_recv_field("src_q2atm",src_field_3D)
398    src_field_3D_p1(:,:,1)=src_field_2D(:,:)
399    src_field_3D_p1(:,:,2:src_nlev+1)=src_field_3D(:,:,:)
400
401    CALL xios_send_field("q2_clean",src_field_3D_p1(1:src_ni-1,1:src_nj,1:src_nlev+1))
402   
403    ! Temperature:
404    CALL xios_recv_field("src_temp",src_field_3D)
405    ! compute theta_rhodz
406    DO l=1,src_nlev
407      src_pk(:,:)=((.5/preff)*(src_pressure(:,:,l)+src_pressure(:,:,l+1)))**kappa
408      src_theta_rhodz(:,:,l) = src_field_3D(:,:,l) * &
409      ((src_pressure(:,:,l)-src_pressure(:,:,l+1))/gravity)/src_pk(:,:)
410    ENDDO
411    CALL xios_send_field("theta_rhodz_clean", &
412                         src_theta_rhodz(1:src_ni-1,1:src_nj,1:src_nlev))
413
414    ! zonal wind
415    CALL xios_recv_field("src_u",src_field_3D)
416    CALL xios_send_field("u_clean", &
417                          src_field_3D(1:src_ni-1,1:src_nj,1:src_nlev))
418   
419    ! meridional wind
420    CALL xios_recv_field("src_v",src_field_3D)
421    CALL xios_send_field("v_clean", &
422                          src_field_3D(1:src_ni-1,1:src_nj,1:src_nlev))
423    ! tracers
424    DO i=1,src_nq
425       if (i < 10) then
426           format_string = "(A8,I1)"
427       else
428           format_string = "(A8,I2)"
429       endif
430       write (i_trac, format_string) "src_trac", i
431       CALL xios_recv_field(i_trac,src_field_3D(:,:,:))
432       if (i==1) THEN
433          CALL xios_send_field("co2_clean", &
434                        src_field_3D(1:src_ni-1,1:src_nj,1:src_nlev))
435          CALL xios_recv_field("src_co2_surf",src_field_2D_subslope)
436          CALL xios_send_field("co2_surf_clean", &
437                        src_field_2D_subslope(1:src_ni-1,1:src_nj,:))
438       elseif (i==2) THEN
439          CALL xios_send_field("dust_number_clean", &
440                        src_field_3D(1:src_ni-1,1:src_nj,1:src_nlev))
441          CALL xios_recv_field("src_dust_number_surf",src_field_2D_subslope)
442          CALL xios_send_field("dust_number_surf_clean", &
443                        src_field_2D_subslope(1:src_ni-1,1:src_nj,:))
444       elseif (i==3) THEN
445          CALL xios_send_field("dust_mass_clean", &
446                        src_field_3D(1:src_ni-1,1:src_nj,1:src_nlev))
447          CALL xios_recv_field("src_dust_mass_surf",src_field_2D_subslope)
448          CALL xios_send_field("dust_mass_surf_clean", &
449                        src_field_2D_subslope(1:src_ni-1,1:src_nj,:))
450       elseif (i==4) THEN
451          CALL xios_send_field("ccn_number_clean", &
452                        src_field_3D(1:src_ni-1,1:src_nj,1:src_nlev))
453          CALL xios_recv_field("src_ccn_number_surf",src_field_2D_subslope)
454          CALL xios_send_field("ccn_number_surf_clean", &
455                        src_field_2D_subslope(1:src_ni-1,1:src_nj,:))
456       elseif (i==5) THEN
457          CALL xios_send_field("ccn_mass_clean", &
458                        src_field_3D(1:src_ni-1,1:src_nj,1:src_nlev))
459          CALL xios_recv_field("src_ccn_mass_surf",src_field_2D_subslope)
460          CALL xios_send_field("ccn_mass_surf_clean", &
461                        src_field_2D_subslope(1:src_ni-1,1:src_nj,:))
462       elseif (i==6) THEN
463          CALL xios_send_field("h2o_ice_clean", &
464                        src_field_3D(1:src_ni-1,1:src_nj,1:src_nlev))
465          CALL xios_recv_field("src_h2o_ice_surf",src_field_2D_subslope)
466          CALL xios_send_field("h2o_ice_surf_clean", &
467                        src_field_2D_subslope(1:src_ni-1,1:src_nj,:))
468       elseif (i==7) THEN
469          CALL xios_send_field("h2o_vap_clean", &
470                        src_field_3D(1:src_ni-1,1:src_nj,1:src_nlev))
471          CALL xios_recv_field("src_h2o_vap_surf",src_field_2D_subslope)
472          CALL xios_send_field("h2o_vap_surf_clean", &
473                        src_field_2D_subslope(1:src_ni-1,1:src_nj,:))
474        elseif (i==8) THEN
475          CALL xios_send_field("stormdust_number_clean", &
476                        src_field_3D(1:src_ni-1,1:src_nj,1:src_nlev))
477          CALL xios_recv_field("src_stormdust_number_surf",src_field_2D_subslope)
478          CALL xios_send_field("stormdust_number_surf_clean", &
479                        src_field_2D_subslope(1:src_ni-1,1:src_nj,:))
480        elseif (i==9) THEN
481          CALL xios_send_field("stormdust_mass_clean", &
482                        src_field_3D(1:src_ni-1,1:src_nj,1:src_nlev))
483          CALL xios_recv_field("src_stormdust_mass_surf",src_field_2D_subslope)
484          CALL xios_send_field("stormdust_mass_surf_clean", &
485                        src_field_2D_subslope(1:src_ni-1,1:src_nj,:))
486        elseif (i==10) THEN
487          CALL xios_send_field("topdust_number_clean", &
488                        src_field_3D(1:src_ni-1,1:src_nj,1:src_nlev))
489          CALL xios_recv_field("src_topdust_number_surf",src_field_2D_subslope)
490          CALL xios_send_field("topdust_number_surf_clean", &
491                        src_field_2D_subslope(1:src_ni-1,1:src_nj,:))
492        elseif (i==11) THEN
493          CALL xios_send_field("topdust_mass_clean", &
494                        src_field_3D(1:src_ni-1,1:src_nj,1:src_nlev))
495          CALL xios_recv_field("src_topdust_mass_surf",src_field_2D_subslope)
496          CALL xios_send_field("topdust_mass_surf_clean", &
497                        src_field_2D_subslope(1:src_ni-1,1:src_nj,:))
498       ENDIF
499    ENDDO
500
501    ! soil temperature
502    CALL xios_recv_field("src_tsoil",src_field_3D_soil_subslope)
503    CALL xios_send_field("tsoil_clean", &
504                          src_field_3D_soil_subslope(1:src_ni-1,1:src_nj,1:src_nsoil,:))
505    ! soil thermal intertia
506    CALL xios_recv_field("src_inertiedat",src_field_3D_soil)
507    CALL xios_send_field("inertiedat_clean", &
508                          src_field_3D_soil(1:src_ni-1,1:src_nj,1:src_nsoil))
509
510    ! ZMEA "zmea Orographie sous-maille"
511    CALL xios_recv_field("src_ZMEA",src_field_2D)
512    WRITE(*,*) rank,":src_ZMEA(1:2,1:3)=",src_field_2D(1:2,1:3)
513    CALL xios_send_field("ZMEA_clean",src_field_2D(1:src_ni-1,1:src_nj))
514   
515    ! ZSTD "zstd Orographie sous-maille"
516    CALL xios_recv_field("src_ZSTD",src_field_2D)
517    WRITE(*,*) rank,":src_ZSTD(1:2,1:3)=",src_field_2D(1:2,1:3)
518    CALL xios_send_field("ZSTD_clean",src_field_2D(1:src_ni-1,1:src_nj))
519   
520    ! ZSIG "zsig Orographie sous-maille"
521    CALL xios_recv_field("src_ZSIG",src_field_2D)
522    WRITE(*,*) rank,":src_ZSIG(1:2,1:3)=",src_field_2D(1:2,1:3)
523    CALL xios_send_field("ZSIG_clean",src_field_2D(1:src_ni-1,1:src_nj))
524   
525    ! ZGAM "zgam Orographie sous-maille"
526    CALL xios_recv_field("src_ZGAM",src_field_2D)
527    WRITE(*,*) rank,":src_ZGAM(1:2,1:3)=",src_field_2D(1:2,1:3)
528    CALL xios_send_field("ZGAM_clean",src_field_2D(1:src_ni-1,1:src_nj))
529   
530    ! ZTHE "zthe Orographie sous-maille"
531    CALL xios_recv_field("src_ZTHE",src_field_2D)
532    WRITE(*,*) rank,":src_ZTHE(1:2,1:3)=",src_field_2D(1:2,1:3)
533    CALL xios_send_field("ZTHE_clean",src_field_2D(1:src_ni-1,1:src_nj))
534
535    ! albedodat "albedo"
536    CALL xios_recv_field("src_albedodat",src_field_2D)
537    WRITE(*,*) rank,":src_albedodat(1:2,1:3)=",src_field_2D(1:2,1:3)
538    CALL xios_send_field("albedodat_clean",src_field_2D(1:src_ni-1,1:src_nj))
539
540    ! z0 "surface roughness"
541    CALL xios_recv_field("src_z0",src_field_2D)
542    WRITE(*,*) rank,":src_z0(1:2,1:3)=",src_field_2D(1:2,1:3)
543    CALL xios_send_field("z0_clean",src_field_2D(1:src_ni-1,1:src_nj))
544
545
546  ENDDO ! of DO ts=1,src_nt
547 
548!! Finalize
549  write(*,*) rank,":Finalize: call xios_context_finalize"
550  CALL xios_context_finalize()
551
552  write(*,*) rank,":Finalize: call MPI_COMM_FREE"
553  CALL MPI_COMM_FREE(comm, ierr)
554  write(*,*) rank,":Finalize: call xios_finalize"
555  CALL xios_finalize()
556
557  if (rank==0) then
558    ! add a couple of things in the "startphy_icosa.nc" file
559    write(*,*) rank,"Write controle() to startphy_icosa.nc"
560    ierr=NF90_OPEN(output_startfi_file,NF90_WRITE,ncid)
561    ierr=NF90_REDEF(ncid) ! switch to define mode
562    ierr=NF90_DEF_DIM(ncid,"index",100,dimids(1))
563    ierr=NF90_DEF_VAR(ncid,"controle",NF90_DOUBLE,dimids(1),varid)
564    ierr=NF90_ENDDEF(ncid) ! switch out of define mode
565    ierr=NF90_PUT_VAR(ncid,varid,src_controle(1:100))
566    if (ierr.ne.NF90_NOERR) then
567      write(*,*) "NetCDF Error:",NF90_STRERROR(ierr)
568    endif
569    ierr=NF90_CLOSE(ncid)
570    ! add a couple of things in the "start_icosa.nc" file
571    ierr=NF90_OPEN(output_start_file,NF90_WRITE,ncid)
572    ierr=NF90_REDEF(ncid)
573    ierr=NF90_DEF_DIM(ncid,"nvertex_u",2,dimids(1))
574    ierr=NF90_DEF_VAR(ncid,"iteration",NF90_FLOAT,varid)
575    ierr=NF90_ENDDEF(ncid)
576    if (ierr.ne.NF90_NOERR) then
577      write(*,*) "NetCDF Error:",NF90_STRERROR(ierr)
578    endif
579    ierr=NF90_PUT_VAR(ncid,varid,0) ! set "iteration" value to 0
580    ierr=NF90_CLOSE(ncid)
581   
582  endif ! of if (rank==0)
583
584  write(*,*) rank,":Finalize: call MPI_FINALIZE"
585  CALL MPI_FINALIZE(ierr)
586
587  write(*,*) rank,":my_remap: all is well that ends well!"
588
589END PROGRAM start_archive2icosa
Note: See TracBrowser for help on using the repository browser.