source: trunk/LMDZ.COMMON/libf/evolution/read_XIOS_data.F90 @ 3985

Last change on this file since 3985 was 3983, checked in by jbclement, 10 days ago

PEM:

  • Removing completely the ice metamorphism computed by a threshold at the end of the PCM (which was commented).
  • Addition of a module "metamorphism" to compute the PCM frost at the PEM beginning and give it back to the PCM at the PEM end. The frost is considered as the ice given by the PCM "startfi.nc" which is above the yearly minimum. Thereby, metamorphism is performed through this operation.
  • Ice reservoirs representation in the PEM is modernized.

JBC
File size: 10.1 KB
Line 
1MODULE read_XIOS_data
2
3use netcdf, only: nf90_open, nf90_close, nf90_inquire_dimension, nf90_inq_dimid, nf90_noerr, nf90_nowrite, nf90_get_var, nf90_inq_varid
4
5implicit none
6
7character(19), parameter :: file1_daily  = "Xoutdaily4pem_Y1.nc"
8character(19), parameter :: file2_daily  = "Xoutdaily4pem_Y2.nc"
9character(20), parameter :: file1_yearly = "Xoutyearly4pem_Y1.nc"
10character(20), parameter :: file2_yearly = "Xoutyearly4pem_Y2.nc"
11character(256)           :: msg      ! For reading
12integer                  :: fID, vID ! For reading
13
14interface get_var
15    module procedure get_var_1d, get_var_2d, get_var_3d, get_var_4d
16end interface get_var
17
18!=======================================================================
19contains
20!=======================================================================
21SUBROUTINE read_PCM_data(ngrid,nslope,nsoil_PCM,nsol,h2ofrost_PCM,co2frost_PCM,ps_avg,tsurf_avg,tsurf_avg_y1,tsoil_avg,tsoil_ts,watersurf_density_avg,d_h2oice,d_co2ice, &
22                         ps_ts,q_h2o_ts,q_co2_ts,watersoil_density_ts,min_h2oice,min_co2ice)
23
24use grid_conversion,  only: lonlat2vect
25use comsoil_h_PEM,    only: soil_pem
26use compute_tend_mod, only: compute_tend
27use metamorphism,     only: compute_frost
28
29implicit none
30
31! Arguments
32!----------
33integer,                       intent(in) :: ngrid, nslope, nsoil_PCM, nsol
34real, dimension(ngrid,nslope), intent(in) :: h2ofrost_PCM, co2frost_PCM
35real, dimension(ngrid),                       intent(out) :: ps_avg
36real, dimension(ngrid,nslope),                intent(out) :: tsurf_avg, tsurf_avg_y1, watersurf_density_avg, d_h2oice, d_co2ice, min_h2oice, min_co2ice
37real, dimension(ngrid,nsoil_PCM,nslope),      intent(out) :: tsoil_avg
38real, dimension(ngrid,nsol),                  intent(out) :: ps_ts, q_h2o_ts, q_co2_ts
39real, dimension(ngrid,nsoil_PCM,nslope,nsol), intent(out) :: tsoil_ts, watersoil_density_ts
40
41! Local variables
42!----------------
43integer                               :: islope, isoil, isol, nlon, nlat
44real, dimension(:,:),     allocatable :: var_read_2d
45real, dimension(:,:,:),   allocatable :: var_read_3d
46real, dimension(:,:,:,:), allocatable :: var_read_4d
47character(:),             allocatable :: num ! For reading slope variables
48real, dimension(ngrid,nslope,2)       :: min_h2operice, min_co2perice, min_h2ofrost, min_co2frost
49
50! Code
51!-----
52! Initialization
53min_h2operice = 0.
54min_co2perice = 0.
55min_h2ofrost = 0.
56min_co2frost = 0.
57if (nslope == 1) then ! No slope
58    allocate(character(0) :: num)
59else ! Using slopes
60    allocate(character(8) :: num)
61endif
62
63!~~~~~~~~~~~~~~~~~~~~~~~~ Year 1 - Yearly data ~~~~~~~~~~~~~~~~~~~~~~~~~
64! Open the NetCDF file of XIOS outputs
65write(*,*) "Opening "//file1_yearly//"..."
66call error_msg(nf90_open(file1_yearly,nf90_nowrite,fID),'open',file1_yearly)
67
68! Get the dimensions
69call error_msg(nf90_inq_dimid(fID,'lon',vID),'inq',file1_yearly)
70call error_msg(nf90_inquire_dimension(fID,vID,len = nlon),'inq',file1_yearly)
71call error_msg(nf90_inq_dimid(fID,'lat',vID),'inq',file1_yearly)
72call error_msg(nf90_inquire_dimension(fID,vID,len = nlat),'inq',file1_yearly)
73
74! Allocate and read the variables
75allocate(var_read_2d(nlon,nlat),var_read_3d(nlon,nlat,nsoil_PCM))
76do islope = 1,nslope
77    if (nslope /= 1) then
78        num='  '
79        write(num,'(i2.2)') islope
80        num = '_slope'//num
81    endif
82    call get_var('co2ice'//num,var_read_2d)          ; call lonlat2vect(nlon,nlat,ngrid,var_read_2d,min_co2frost(:,islope,1))
83    call get_var('h2o_ice_s'//num,var_read_2d)       ; call lonlat2vect(nlon,nlat,ngrid,var_read_2d,min_h2ofrost(:,islope,1))
84    call get_var('watercap'//num,var_read_2d)        ; call lonlat2vect(nlon,nlat,ngrid,var_read_2d,min_h2operice(:,islope,1))
85    call get_var('perennial_co2ice'//num,var_read_2d); call lonlat2vect(nlon,nlat,ngrid,var_read_2d,min_co2perice(:,islope,1))
86    call get_var('tsurf'//num,var_read_2d)           ; call lonlat2vect(nlon,nlat,ngrid,var_read_2d,tsurf_avg_y1(:,islope))
87enddo
88
89! Close the NetCDF file of XIOS outputs
90call error_msg(nf90_close(fID),"close",file1_yearly)
91write(*,*) '> Data from "'//file1_yearly//'" downloaded!'
92
93!~~~~~~~~~~~~~~~~~~~~~~~~ Year 2 - Yearly data ~~~~~~~~~~~~~~~~~~~~~~~~~
94! Open the NetCDF file of XIOS outputs
95write(*,*) "Opening "//file2_yearly//"..."
96call error_msg(nf90_open(file2_yearly,nf90_nowrite,fID),'open',file2_yearly)
97
98! Allocate and read the variables
99call get_var('ps',var_read_2d); call lonlat2vect(nlon,nlat,ngrid,var_read_2d,ps_avg)
100do islope = 1,nslope
101    if (nslope /= 1) then
102        num='  '
103        write(num,'(i2.2)') islope
104        num = '_slope'//num
105    endif
106    call get_var('tsurf'//num,var_read_2d)           ; call lonlat2vect(nlon,nlat,ngrid,var_read_2d,tsurf_avg(:,islope))
107    call get_var('co2ice'//num,var_read_2d)          ; call lonlat2vect(nlon,nlat,ngrid,var_read_2d,min_co2frost(:,islope,2))
108    call get_var('h2o_ice_s'//num,var_read_2d)       ; call lonlat2vect(nlon,nlat,ngrid,var_read_2d,min_h2ofrost(:,islope,2))
109    call get_var('watercap'//num,var_read_2d)        ; call lonlat2vect(nlon,nlat,ngrid,var_read_2d,min_h2operice(:,islope,2))
110    call get_var('perennial_co2ice'//num,var_read_2d); call lonlat2vect(nlon,nlat,ngrid,var_read_2d,min_co2perice(:,islope,2))
111    if (soil_pem) then
112        call get_var('soiltemp'//num,var_read_3d)
113        do isoil = 1,nsoil_PCM
114            call lonlat2vect(nlon,nlat,ngrid,var_read_3d(:,:,isoil),tsoil_avg(:,isoil,islope))
115        enddo
116        call get_var('waterdensity_surface'//num,var_read_2d); call lonlat2vect(nlon,nlat,ngrid,var_read_2d,watersurf_density_avg(:,islope))
117    endif
118enddo
119! Close the NetCDF file of XIOS outputs
120call error_msg(nf90_close(fID),"close",file1_yearly)
121write(*,*) '> Data from "'//file1_yearly//'" downloaded!'
122
123!~~~~~~~~~~~~~~~~~~~~~~~~~ Year 2 - Daily data ~~~~~~~~~~~~~~~~~~~~~~~~~
124! Open the NetCDF file of XIOS outputs
125write(*,*) "Opening "//file2_daily//"..."
126call error_msg(nf90_open(file2_daily,nf90_nowrite,fID),'open',file2_daily)
127
128! Allocate and read the variables
129deallocate(var_read_2d,var_read_3d)
130allocate(var_read_3d(nlon,nlat,nsol),var_read_4d(nlon,nlat,nsoil_PCM,nsol))
131call get_var('ps',var_read_3d)
132do isol = 1,nsoil_PCM
133    call lonlat2vect(nlon,nlat,ngrid,var_read_3d(:,:,isol),ps_ts(:,isol))
134enddo
135call get_var('h2o_layer1',var_read_3d)
136do isol = 1,nsoil_PCM
137    call lonlat2vect(nlon,nlat,ngrid,var_read_3d(:,:,isol),q_h2o_ts(:,isol))
138enddo
139call get_var('co2_layer1',var_read_3d)
140do isol = 1,nsoil_PCM
141    call lonlat2vect(nlon,nlat,ngrid,var_read_3d(:,:,isol),q_co2_ts(:,isol))
142enddo
143if (soil_pem) then
144    do islope = 1,nslope
145        if (nslope /= 1) then
146            num='  '
147            write(num,'(i2.2)') islope
148            num = '_slope'//num
149        endif
150        call get_var('soiltemp'//num,var_read_4d)
151        do isol = 1,nsol
152            do isoil = 1,nsoil_PCM
153                call lonlat2vect(nlon,nlat,ngrid,var_read_4d(:,:,isoil,isol),tsoil_ts(:,isoil,islope,isol))
154            enddo
155        enddo
156        call get_var('waterdensity_soil'//num,var_read_4d)
157        do isol = 1,nsol
158            do isoil = 1,nsoil_PCM
159                call lonlat2vect(nlon,nlat,ngrid,var_read_4d(:,:,isoil,isol),watersoil_density_ts(:,isoil,islope,isol))
160            enddo
161        enddo
162    enddo
163endif
164deallocate(var_read_3d,var_read_4d,num)
165
166! Close the NetCDF file of XIOS outputs
167call error_msg(nf90_close(fID),"close",file1_daily)
168write(*,*) '> Data from "'//file2_daily//'" downloaded!'
169
170!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
171! Compute frost from yearly minima
172call compute_frost(ngrid,nslope,h2ofrost_PCM,min_h2ofrost,co2frost_PCM,min_co2frost)
173
174!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
175! Compute ice tendencies from yearly minima
176write(*,*) '> Computing surface ice tendencies'
177call compute_tend(ngrid,nslope,min_h2operice + min_h2ofrost,d_h2oice)
178write(*,*) 'H2O ice tendencies (min/max):', minval(d_h2oice), maxval(d_h2oice)
179call compute_tend(ngrid,nslope,min_co2perice + min_co2frost,d_co2ice)
180write(*,*) 'CO2 ice tendencies (min/max):', minval(d_co2ice), maxval(d_co2ice)
181
182END SUBROUTINE read_PCM_data
183
184!=======================================================================
185SUBROUTINE error_msg(ierr,typ,nam)
186
187implicit none
188
189integer,      intent(in) :: ierr !--- NetCDF error code
190character(*), intent(in) :: typ  !--- Type of operation
191character(*), intent(in) :: nam  !--- Field/File name
192
193if (ierr == nf90_noerr) return
194select case(typ)
195    case('inq')  ; msg = "Dim/Field <"//trim(nam)//"> is missing"
196    case('get')  ; msg = "Reading failed for <"//trim(nam)//">"
197    case('put')  ; msg = "Writing failed for <"//trim(nam)//">"
198    case('open') ; msg = "File opening failed for <"//trim(nam)//">"
199    case('close'); msg = "File closing failed for <"//trim(nam)//">"
200    case default
201        write(*,*) 'There is no message for this error.'
202        error stop
203end select
204call abort_gcm(__FILE__,trim(msg),ierr)
205
206END SUBROUTINE error_msg
207
208!=======================================================================
209SUBROUTINE get_var_1d(var,v)
210
211implicit none
212
213character(*), intent(in) :: var
214real, dimension(:), intent(out) :: v
215
216call error_msg(nf90_inq_varid(fID,var,vID),"inq",var)
217call error_msg(nf90_get_var(fID,vID,v),"get",var)
218
219END SUBROUTINE get_var_1d
220
221!=======================================================================
222SUBROUTINE get_var_2d(var,v)
223
224implicit none
225
226character(*), intent(in) :: var
227real, dimension(:,:), intent(out) :: v
228
229call error_msg(nf90_inq_varid(fID,var,vID),"inq",var)
230call error_msg(nf90_get_var(fID,vID,v),"get",var)
231
232END SUBROUTINE get_var_2d
233
234!=======================================================================
235SUBROUTINE get_var_3d(var,v)
236
237implicit none
238
239character(*), intent(in) :: var
240real, dimension(:,:,:), intent(out) :: v
241
242call error_msg(nf90_inq_varid(fID,var,vID),"inq",var)
243call error_msg(nf90_get_var(fID,vID,v),"get",var)
244
245END SUBROUTINE get_var_3d
246
247!=======================================================================
248SUBROUTINE get_var_4d(var,v)
249
250implicit none
251
252character(*), intent(in) :: var
253real, dimension(:,:,:,:), intent(out) :: v
254
255call error_msg(nf90_inq_varid(fID,var,vID),"inq",var)
256call error_msg(nf90_get_var(fID,vID,v),"get",var)
257
258END SUBROUTINE get_var_4d
259
260END MODULE read_XIOS_data
Note: See TracBrowser for help on using the repository browser.