source: trunk/LMDZ.COMMON/libf/evolution/read_data_PCM_mod.F90 @ 3620

Last change on this file since 3620 was 3620, checked in by jbclement, 28 hours ago

PEM:
Few small adjustments in the code.
JBC
File size: 14.3 KB
Line 
1MODULE read_data_PCM_mod
2
3use netcdf, only: nf90_open, NF90_NOWRITE, nf90_noerr, nf90_strerror, &
4                  nf90_get_var, nf90_inq_varid, nf90_inq_dimid,       &
5                  nf90_inquire_dimension, nf90_close
6
7implicit none
8
9character(13), parameter :: modname = 'read_data_PCM'
10character(256)           :: msg      ! for reading
11integer                  :: fID, vID ! for reading
12
13!=======================================================================
14contains
15!=======================================================================
16
17SUBROUTINE read_data_PCM(filename1,filename2,timelen,iim_input,jjm_input,ngrid,nslope,vmr_co2_PCM,ps_timeseries,ps_avg,tsurf_avg_yr1,tsurf_avg, &
18                         tsoil_avg,tsoil_timeseries,min_co2_ice,min_h2o_ice,q_co2,q_h2o,watersurf_density_avg,watersoil_density_timeseries)
19use comsoil_h,             only: nsoilmx
20use comsoil_h_PEM,         only: soil_pem
21use constants_marspem_mod, only: m_co2, m_noco2
22
23implicit none
24
25!=======================================================================
26!
27! Purpose: Read initial confitions file from the PCM
28!
29! Authors: JBC
30!=======================================================================
31
32include "dimensions.h"
33
34!=======================================================================
35! Inputs:
36character(*), intent(in) :: filename1, filename2                ! File names
37integer,      intent(in) :: timelen                             ! Number of times stored in the file
38integer,      intent(in) :: iim_input, jjm_input, ngrid, nslope ! Number of points in the lat x lon dynamical grid, number of subgrid slopes
39! Ouputs
40real, dimension(ngrid,nslope,2),               intent(out) :: min_co2_ice                  ! Minimum of co2 ice per slope of the year [kg/m^2]
41real, dimension(ngrid,nslope,2),               intent(out) :: min_h2o_ice                  ! Minimum of h2o ice per slope of the year [kg/m^2]
42real, dimension(ngrid,timelen),                intent(out) :: vmr_co2_PCM                  ! Grid points x Times co2 volume mixing ratio retrieve from the PCM [m^3/m^3]
43real, dimension(ngrid,timelen),                intent(out) :: q_co2                        ! CO2 mass mixing ratio in the first layer [kg/m^3]
44real, dimension(ngrid,timelen),                intent(out) :: q_h2o                        ! H2O mass mixing ratio in the first layer [kg/m^3]
45real, dimension(ngrid,timelen),                intent(out) :: ps_timeseries                ! Surface pressure timeseries [Pa]
46real, dimension(ngrid),                        intent(out) :: ps_avg                       ! Averaged surface pressure [K]
47real, dimension(ngrid,nslope),                 intent(out) :: tsurf_avg                    ! Averaged surface temperature [K]
48real, dimension(ngrid,nslope),                 intent(out) :: tsurf_avg_yr1                ! Averaged surface temperature for year 1 [K]
49real, dimension(ngrid,nsoilmx,nslope),         intent(out) :: tsoil_avg                    ! Averaged soil temperature for year 2 [K]
50real, dimension(ngrid,nsoilmx,nslope,timelen), intent(out) :: tsoil_timeseries             ! Soil temperature timeseries [K]
51real, dimension(ngrid,nslope),                 intent(out) :: watersurf_density_avg        ! Water density at the surface [kg/m^3]
52real, dimension(ngrid,nsoilmx,nslope,timelen), intent(out) :: watersoil_density_timeseries ! Water density timeseries in the soil layer [kg/m^3]
53! Local variables
54integer                               :: i, j, t, islope                       ! Loop variables
55real                                  :: A, B                                  ! Intermediate variables to compute the mean molar mass of the layer
56character(:), allocatable             :: num                                   ! For reading sloped variables
57real, dimension(:,:),     allocatable :: var2_read                             ! Variables for reading (2 dimensions)
58real, dimension(:,:,:),   allocatable :: var3_read_1, var3_read_2, var3_read_3 ! Variables for reading (3 dimensions)
59real, dimension(:,:,:,:), allocatable :: var4_read                             ! Variables for reading (4 dimensions)
60!-----------------------------------------------------------------------
61
62!------------------------------- Year 1 --------------------------------
63! Open the NetCDF file
64write(*,*) "Opening "//filename1//"..."
65call error_msg(NF90_OPEN(filename1,NF90_NOWRITE,fID),"open",filename1)
66
67! Download the data from the file
68allocate(var2_read(iim_input + 1,jjm_input + 1),var3_read_1(iim_input + 1,jjm_input + 1,timelen),var3_read_2(iim_input + 1,jjm_input + 1,timelen))
69
70if (nslope == 1) then ! No slope
71    allocate(character(0) :: num)
72else ! We use slopes
73    allocate(character(8) :: num)
74endif
75
76do islope = 1,nslope
77    if (nslope /= 1) then
78        num='  '
79        write(num,'(i2.2)') islope
80        num = '_slope'//num
81    endif
82
83    ! CO2 ice
84    !--------
85    call get_var3("co2ice"//num,var3_read_1)
86    where (var3_read_1 < 0.) var3_read_1 = 0.
87    write(*,*) "Data for co2_ice"//num//" downloaded."
88#ifndef CPP_STD
89    call get_var3("perennial_co2ice"//num,var3_read_2)
90    write(*,*) "Data for perennial_co2ice"//num//" downloaded."
91#endif
92
93    ! Compute the minimum over the year for each point
94    var2_read = minval(var3_read_1 + var3_read_2,3)
95#ifndef CPP_1D
96    call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,var2_read,min_co2_ice(:,islope,1))
97#else
98    min_co2_ice(1,islope,1) = var2_read(1,1)
99#endif
100    write(*,*) "Min of co2_ice"//num//" computed."
101
102    ! H2O ice
103    !--------
104    call get_var3("h2o_ice_s"//num,var3_read_1)
105    where (var3_read_1 < 0.) var3_read_1 = 0.
106    write(*,*) "Data for h2o_ice_s"//num//" downloaded."
107#ifndef CPP_STD
108    call get_var3("watercap"//num,var3_read_2)
109    write(*,*) "Data for watercap"//num//" downloaded."
110#endif
111
112    ! Compute the minimum over the year for each point
113    var2_read = minval(var3_read_1 + var3_read_2,3)
114#ifndef CPP_1D
115    call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,var2_read,min_h2o_ice(:,islope,1))
116#else
117    min_h2o_ice(1,islope,1) = var2_read(1,1)
118#endif
119    write(*,*) "Min of h2o_ice"//num//" computed."
120
121    ! Tsurf
122    !------
123    call get_var3("tsurf"//num,var3_read_1)
124    write(*,*) "Data for tsurf"//num//" downloaded."
125
126    ! Compute average over the year for each point
127    var2_read = sum(var3_read_1,3)/timelen
128#ifndef CPP_1D
129    call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,var2_read,tsurf_avg_yr1(:,islope))
130#else
131    tsurf_avg_yr1(1,islope) = var2_read(1,1)
132#endif
133    write(*,*) "Average of tsurf"//num//" computed."
134enddo
135
136! Close the NetCDF file
137call error_msg(nf90_close(fID),"close",filename1)
138write(*,*) '> "'//filename1//'" processed!'
139
140!------------------------------- Year 2 --------------------------------
141! Open the NetCDF file
142write(*,*) "Opening "//filename2//"..."
143call error_msg(NF90_OPEN(filename2,NF90_NOWRITE,fID),"open",filename2)
144
145! Download the data from the file
146allocate(var3_read_3(iim_input + 1,jjm_input + 1,nsoilmx),var4_read(iim_input + 1,jjm_input + 1,nsoilmx,timelen))
147
148! Surface pressure
149!-----------------
150call get_var3("ps",var3_read_1)
151write(*,*) "Data for surface pressure downloaded."
152
153! Compute average over the year for each point
154var2_read = sum(var3_read_1,3)/timelen
155#ifndef CPP_1D
156    call gr_dyn_fi(timelen,iim_input + 1,jjm_input + 1,ngrid,var3_read_1,ps_timeseries)
157    call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,var2_read,ps_avg)
158#else
159    ps_timeseries(1,:) = var3_read_1(1,1,:)
160    ps_avg(1) = var2_read(1,1)
161#endif
162write(*,*) "Average of surface pressure computed."
163
164! CO2 vmr
165!--------
166call get_var3("co2_layer1",var3_read_1)
167where (var3_read_1 < 0.)
168    var3_read_1 = 1.e-10
169else where (var3_read_1 > 1.)
170    var3_read_1 = 1.
171end where
172#ifndef CPP_1D
173    call gr_dyn_fi(timelen,iim_input + 1,jjm_input + 1,ngrid,var3_read_1,q_co2)
174#else
175    q_co2(1,:) = var3_read_1(1,1,:)
176#endif
177A = (1./m_co2 - 1./m_noco2)
178B = 1./m_noco2
179vmr_co2_PCM = q_co2/(A*q_co2 + B)/m_co2
180write(*,*) "Data for CO2 vmr downloaded."
181
182! H2O vmr
183!--------
184call get_var3("h2o_layer1",var3_read_1)
185where (var3_read_1 < 0.)
186    var3_read_1 = 1.e-10
187else where (var3_read_1 > 1.)
188    var3_read_1 = 1.
189end where
190#ifndef CPP_1D
191    call gr_dyn_fi(timelen,iim_input + 1,jjm_input + 1,ngrid,var3_read_1,q_h2o)
192#else
193    q_h2o(1,:) = var3_read_1(1,1,:)
194#endif
195write(*,*) "Data for H2O vmr downloaded."
196
197do islope = 1,nslope
198    if (nslope /= 1) then
199        num='  '
200        write(num,'(i2.2)') islope
201        num = '_slope'//num
202    endif
203
204    ! CO2 ice
205    !--------
206    call get_var3("co2ice"//num,var3_read_1)
207    where (var3_read_1 < 0.) var3_read_1 = 0.
208    write(*,*) "Data for co2_ice"//num//" downloaded."
209#ifndef CPP_STD
210    call get_var3("perennial_co2ice"//num,var3_read_2)
211    write(*,*) "Data for perennial_co2ice"//num//" downloaded."
212#endif
213
214    ! Compute the minimum over the year for each point
215    var2_read = minval(var3_read_1 + var3_read_2,3)
216#ifndef CPP_1D
217    call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,var2_read,min_co2_ice(:,islope,2))
218#else
219    min_co2_ice(1,islope,2) = var2_read(1,1)
220#endif
221    write(*,*) "Min of co2_ice"//num//" computed."
222
223    ! H2O ice
224    !--------
225    call get_var3("h2o_ice_s"//num,var3_read_1)
226    where (var3_read_1 < 0.) var3_read_1 = 0.
227    write(*,*) "Data for h2o_ice_s"//num//" downloaded."
228#ifndef CPP_STD
229    call get_var3("watercap"//num,var3_read_2)
230    write(*,*) "Data for watercap"//num//" downloaded."
231#endif
232
233    ! Compute the minimum over the year for each point
234    var2_read = minval(var3_read_1 + var3_read_2,3)
235#ifndef CPP_1D
236    call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,var2_read,min_h2o_ice(:,islope,2))
237#else
238    min_h2o_ice(1,islope,2) = var2_read(1,1)
239#endif
240    write(*,*) "Min of h2o_ice"//num//" computed."
241
242    ! Tsurf
243    !------
244    call get_var3("tsurf"//num,var3_read_1)
245    write(*,*) "Data for tsurf"//num//" downloaded."
246
247    ! Compute average over the year for each point
248    var2_read = sum(var3_read_1,3)/timelen
249#ifndef CPP_1D
250    call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,var2_read,tsurf_avg(:,islope))
251#else
252    tsurf_avg(1,islope) = var2_read(1,1)
253#endif
254    write(*,*) "Average of tsurf"//num//" computed."
255
256    if (soil_pem) then
257        ! Tsoil
258        !------
259        call get_var4("soiltemp"//num,var4_read)
260        write(*,*) "Data for soiltemp"//num//" downloaded."
261
262        ! Compute average over the year for each point
263        var3_read_3 = sum(var4_read,4)/timelen
264#ifndef CPP_1D
265        do t = 1,timelen
266            call gr_dyn_fi(nsoilmx,iim_input + 1,jjm_input + 1,ngrid,var4_read(:,:,:,t),tsoil_timeseries(:,:,islope,t))
267        enddo
268        call gr_dyn_fi(nsoilmx,iim_input + 1,jjm_input + 1,ngrid,var3_read_3,tsoil_avg(:,:,islope))
269#else
270        tsoil_timeseries(1,:,islope,:) = var4_read(1,1,:,:)
271        tsoil_avg(1,:,islope) = var3_read_3(1,1,:)
272#endif
273        write(*,*) "Average of tsoil"//num//" computed."
274
275        ! Soil water density
276        !-------------------
277        call get_var4("waterdensity_soil"//num,var4_read)
278#ifndef CPP_1D
279        do t = 1,timelen
280            call gr_dyn_fi(nsoilmx,iim_input + 1,jjm_input + 1,ngrid,var4_read(:,:,:,t),watersoil_density_timeseries(:,:,islope,t))
281        enddo
282#else
283        watersoil_density_timeseries(1,:,islope,:) = var4_read(1,1,:,:)
284#endif
285        write(*,*) "Data for waterdensity_soil"//num//" downloaded."
286
287        ! Surface water density
288        !----------------------
289        call get_var3("waterdensity_surface"//num,var3_read_1)
290        write(*,*) "Data for waterdensity_surface"//num//" downloaded."
291
292        ! Compute average over the year for each point
293        var2_read = sum(var3_read_1,3)/timelen
294#ifndef CPP_1D
295        call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,var2_read,watersurf_density_avg(:,islope))
296#else
297        watersurf_density_avg(1,islope) = var2_read(1,1)
298#endif
299        write(*,*) "Average of waterdensity_surface"//num//" computed."
300    endif
301enddo
302deallocate(var2_read,var3_read_1,var3_read_2,var3_read_3,var4_read,num)
303
304! Close the NetCDF file
305call error_msg(nf90_close(fID),"close",filename2)
306write(*,*) '> "'//filename2//'" processed!'
307
308END SUBROUTINE read_data_PCM
309
310!=======================================================================
311
312SUBROUTINE check_dim(n1,n2,str1,str2)
313
314implicit none
315
316integer,            intent(in) :: n1, n2
317character(len = *), intent(in) :: str1, str2
318
319character(256) :: s1, s2
320
321if (n1 /= n2) then
322    s1 = 'value of '//trim(str1)//' ='
323    s2 = ' read in starting file differs from parametrized '//trim(str2)//' ='
324    write(msg,'(10x,a,i4,2x,a,i4)')trim(s1),n1,trim(s2),n2
325    call abort_gcm(trim(modname),trim(msg),1)
326endif
327
328END SUBROUTINE check_dim
329
330!=======================================================================
331
332SUBROUTINE get_var1(var,v)
333
334implicit none
335
336character(len = *), intent(in)  :: var
337real, dimension(:), intent(out) :: v
338
339call error_msg(NF90_INQ_VARID(fID,var,vID),"inq",var)
340call error_msg(NF90_GET_VAR(fID,vID,v),"get",var)
341
342END SUBROUTINE get_var1
343
344!=======================================================================
345
346SUBROUTINE get_var3(var,v) ! on U grid
347
348implicit none
349
350character(len = *),     intent(in)  :: var
351real, dimension(:,:,:), intent(out) :: v
352
353call error_msg(NF90_INQ_VARID(fID,var,vID),"inq",var)
354call error_msg(NF90_GET_VAR(fID,vID,v),"get",var)
355
356END SUBROUTINE get_var3
357
358!=======================================================================
359
360SUBROUTINE get_var4(var,v)
361
362implicit none
363
364character(len = *),       intent(in)  :: var
365real, dimension(:,:,:,:), intent(out) :: v
366
367call error_msg(NF90_INQ_VARID(fID,var,vID),"inq",var)
368call error_msg(NF90_GET_VAR(fID,vID,v),"get",var)
369
370END SUBROUTINE get_var4
371
372!=======================================================================
373
374SUBROUTINE error_msg(ierr,typ,nam)
375
376implicit none
377
378integer,            intent(in) :: ierr !--- NetCDF ERROR CODE
379character(len = *), intent(in) :: typ  !--- TYPE OF OPERATION
380character(len = *), intent(in) :: nam  !--- FIELD/FILE NAME
381
382if (ierr == nf90_noerr) return
383select case(typ)
384    case('inq');   msg="Field <"//trim(nam)//"> is missing"
385    case('get');   msg="Reading failed for <"//trim(nam)//">"
386    case('put');   msg="Writing failed for <"//trim(nam)//">"
387    case('open');  msg="File opening failed for <"//trim(nam)//">"
388    case('close'); msg="File closing failed for <"//trim(nam)//">"
389    case default
390        write(*,*) 'There is no message for this error.'
391        error stop
392end select
393call abort_gcm(trim(modname),trim(msg),ierr)
394
395END SUBROUTINE error_msg
396
397END MODULE read_data_PCM_mod
398
Note: See TracBrowser for help on using the repository browser.