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

Last change on this file since 3599 was 3599, checked in by jbclement, 15 hours ago

PEM:
Few small optimizations.
JBC

File size: 17.8 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: RV & LL
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, l, t, islope        ! Loop variables
55real                                                                :: A, B, mmean               ! Molar Mass of co2 and no co2, A;B intermediate variables to compute the mean molar mass of the layer
56character(2)                                                        :: num                       ! Hor reading sloped variables
57real, dimension(iim_input + 1,jjm_input + 1,nslope,timelen)         :: h2o_ice_s_dyn             ! H2O ice per slope of the concatenated file [kg/m^2]
58real, dimension(iim_input + 1,jjm_input + 1,nslope,timelen)         :: watercap
59real, dimension(iim_input + 1,jjm_input + 1,nslope,timelen)         :: perennial_co2ice
60real, dimension(iim_input + 1,jjm_input + 1,timelen)                :: vmr_co2_dyn               ! CO2 volume mixing ratio in the first layer [mol/m^3]
61real, dimension(iim_input + 1,jjm_input + 1,timelen)                :: ps_dyn                   ! Surface pressure [Pa]
62real, dimension(iim_input + 1,jjm_input + 1)                        :: ps_avg_dyn                ! Averaged surface pressure of the concatenated file [Pa]
63real, dimension(iim_input + 1,jjm_input + 1,nslope)                 :: min_co2_ice_dyn
64real, dimension(iim_input + 1,jjm_input + 1,nslope)                 :: min_h2o_ice_dyn
65real, dimension(iim_input + 1,jjm_input + 1,nslope)                 :: tsurf_avg_dyn             ! Averaged surface temperature of the concatenated file [K]
66real, dimension(iim_input + 1,jjm_input + 1,nsoilmx,nslope)         :: tsoil_avg_dyn             ! Averaged soil temperature of the concatenated file [K]
67real, dimension(iim_input + 1,jjm_input + 1,nslope,timelen)         :: tsurf_dyn                ! Surface temperature of the concatenated file, time series [K]
68real, dimension(iim_input + 1,jjm_input + 1,nsoilmx,nslope,timelen) :: tsoil_dyn                ! Soil temperature of the concatenated file, time series [K]
69real, dimension(iim_input + 1,jjm_input + 1,timelen)                :: q_co2_dyn                 ! CO2 mass mixing ratio in the first layer [kg/m^3]
70real, dimension(iim_input + 1,jjm_input + 1,timelen)                :: q_h2o_dyn                 ! H2O mass mixing ratio in the first layer [kg/m^3]
71real, dimension(iim_input + 1,jjm_input + 1,nslope,timelen)         :: co2_ice_slope_dyn         ! co2 ice amount per slope of the year [kg/m^2]
72real, dimension(iim_input + 1,jjm_input + 1,nslope,timelen)         :: watersurf_density_dyn     ! Water density at the surface, time series [kg/m^3]
73real, dimension(iim_input + 1,jjm_input + 1,nslope)                 :: watersurf_density_dyn_avg ! Water density at the surface, dynamic grid, yearly averaged [kg/m^3]
74real, dimension(iim_input + 1,jjm_input + 1,nsoilmx,nslope,timelen) :: watersoil_density_dyn     ! Water density in the soil layer, time series [kg/m^3]
75!-----------------------------------------------------------------------
76
77!------------------ Year 1 ------------------
78! Open the NetCDF file
79write(*,*) "Opening "//filename1//"..."
80call error_msg(NF90_OPEN(filename1,NF90_NOWRITE,fID),"open",filename1)
81
82! Download the data from the file
83if (nslope == 1) then ! There is no slope
84    call get_var3("co2ice",co2_ice_slope_dyn(:,:,1,:))
85    write(*,*) "Data for co2_ice downloaded."
86
87    call get_var3("h2o_ice_s",h2o_ice_s_dyn(:,:,1,:))
88    write(*,*) "Data for h2o_ice_s downloaded."
89
90    call get_var3("tsurf",tsurf_dyn(:,:,1,:))
91    write(*,*) "Data for tsurf downloaded."
92
93#ifndef CPP_STD
94    call get_var3("watercap",watercap(:,:,1,:))
95    write(*,*) "Data for watercap downloaded."
96
97    call get_var3("perennial_co2ice",perennial_co2ice(:,:,1,:))
98    write(*,*) "Data for perennial_co2ice downloaded."
99#endif
100else ! We use slopes
101    do islope = 1,nslope
102        write(num,'(i2.2)') islope
103        call get_var3("co2ice_slope"//num,co2_ice_slope_dyn(:,:,islope,:))
104        write(*,*) "Data for co2_ice_slope"//num//" downloaded."
105        call get_var3("h2o_ice_s_slope"//num,h2o_ice_s_dyn(:,:,islope,:))
106        write(*,*) "Data for h2o_ice_s_slope"//num//" downloaded."
107        call get_var3("tsurf_slope"//num,tsurf_dyn(:,:,islope,:))
108        write(*,*) "Data for tsurf_slope"//num//" downloaded."
109#ifndef CPP_STD
110        call get_var3("watercap_slope"//num,watercap(:,:,islope,:))
111        write(*,*) "Data for watercap_slope"//num//" downloaded."
112        call get_var3("perennial_co2ice_slope"//num,perennial_co2ice(:,:,islope,:))
113        write(*,*) "Data for perennial_co2ice_slope"//num//" downloaded."
114#endif
115    enddo
116endif
117
118! Close the NetCDF file
119write(*,*) "Closing "//filename1//"..."
120call error_msg(nf90_close(fID),"close",filename1)
121
122! Compute the minimum over the year for each point
123write(*,*) "Computing the min of h2o_ice..."
124where (h2o_ice_s_dyn < 0.) h2o_ice_s_dyn = 0.
125min_h2o_ice_dyn = minval(h2o_ice_s_dyn + watercap,4)
126write(*,*) "Computing the min of co2_ice..."
127where (co2_ice_slope_dyn < 0.) co2_ice_slope_dyn = 0.
128min_co2_ice_dyn = minval(co2_ice_slope_dyn + perennial_co2ice,4)
129
130! Compute averages over the year for each point
131write(*,*) "Computing the average of tsurf..."
132tsurf_avg_dyn = sum(tsurf_dyn,4)/timelen
133
134#ifndef CPP_1D
135    call gr_dyn_fi(nslope,iim_input + 1,jjm_input + 1,ngrid,tsurf_avg_dyn,tsurf_avg)
136    do islope = 1,nslope
137        call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,min_co2_ice_dyn(:,:,islope),min_co2_ice(:,islope,1))
138        call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,min_h2o_ice_dyn(:,:,islope),min_h2o_ice(:,islope,1))
139    enddo
140#else
141    tsurf_avg_yr1(1,:) = tsurf_avg_dyn(1,1,:)
142    min_co2_ice(1,:,1) = min_co2_ice_dyn(1,1,:)
143    min_h2o_ice(1,:,1) = min_h2o_ice_dyn(1,1,:)
144#endif
145write(*,*) "Downloading data Y1 done!"
146
147!------------------ Year 2 ------------------
148
149! Open the NetCDF file
150write(*,*) "Opening "//filename2//"..."
151call error_msg(NF90_OPEN(filename2,NF90_NOWRITE,fID),"open",filename2)
152
153! Download the data from the file
154call get_var3("ps",ps_dyn)
155write(*,*) "Data for surface pressure downloaded."
156
157call get_var3("co2_layer1",q_co2_dyn)
158write(*,*) "Data for vmr co2 downloaded."
159
160call get_var3("h2o_layer1",q_h2o_dyn)
161write(*,*) "Data for vmr h2o downloaded."
162
163if (nslope == 1) then ! There is no slope
164    call get_var3("co2ice",co2_ice_slope_dyn(:,:,1,:))
165    write(*,*) "Data for co2_ice downloaded."
166
167    call get_var3("h2o_ice_s",h2o_ice_s_dyn(:,:,1,:))
168    write(*,*) "Data for h2o_ice_s downloaded."
169
170    call get_var3("tsurf",tsurf_dyn(:,:,1,:))
171    write(*,*) "Data for tsurf downloaded."
172
173#ifndef CPP_STD
174    call get_var3("watercap",watercap(:,:,1,:))
175    write(*,*) "Data for watercap downloaded."
176
177    call get_var3("perennial_co2ice",perennial_co2ice(:,:,1,:))
178    write(*,*) "Data for perennial_co2ice downloaded."
179
180    if (soil_pem) then
181        call get_var4("soiltemp",tsoil_dyn(:,:,:,1,:))
182        write(*,*) "Data for soiltemp downloaded."
183
184        call get_var4("waterdensity_soil",watersoil_density_dyn(:,:,:,1,:))
185        write(*,*) "Data for waterdensity_soil downloaded."
186
187        call get_var3("waterdensity_surface",watersurf_density_dyn(:,:,1,:))
188        write(*,*) "Data for waterdensity_surface downloaded."
189    endif !soil_pem
190#endif
191else ! We use slopes
192    do islope = 1,nslope
193        write(num,'(i2.2)') islope
194        call get_var3("co2ice_slope"//num,co2_ice_slope_dyn(:,:,islope,:))
195        write(*,*) "Data for co2_ice_slope"//num//" downloaded."
196        call get_var3("h2o_ice_s_slope"//num,h2o_ice_s_dyn(:,:,islope,:))
197        write(*,*) "Data for h2o_ice_s_slope"//num//" downloaded."
198        call get_var3("tsurf_slope"//num,tsurf_dyn(:,:,islope,:))
199        write(*,*) "Data for tsurf_slope"//num//" downloaded."
200#ifndef CPP_STD
201        call get_var3("watercap_slope"//num,watercap(:,:,islope,:))
202        write(*,*) "Data for watercap_slope"//num//" downloaded."
203        call get_var3("perennial_co2ice_slope"//num,perennial_co2ice(:,:,islope,:))
204        write(*,*) "Data for perennial_co2ice_slope"//num//" downloaded."
205        if (soil_pem) then
206            call get_var4("soiltemp_slope"//num,tsoil_dyn(:,:,:,islope,:))
207            write(*,*) "Data for soiltemp_slope"//num//" downloaded."
208            call get_var4("waterdensity_soil_slope"//num,watersoil_density_dyn(:,:,:,islope,:))
209            write(*,*) "Data for waterdensity_soil_slope"//num//" downloaded."
210            call get_var3("waterdensity_surface"//num,watersurf_density_dyn(:,:,islope,:))
211            write(*,*) "Data for waterdensity_surface_slope"//num//" downloaded."
212        endif
213#endif
214    enddo
215endif
216
217! Close the NetCDF file
218write(*,*) "Closing "//filename2//"..."
219call error_msg(nf90_close(fID),"close",filename2)
220
221! Compute the minimum over the year for each point
222write(*,*) "Computing the min of h2o_ice..."
223where (h2o_ice_s_dyn < 0.) h2o_ice_s_dyn = 0.
224min_h2o_ice_dyn = minval(h2o_ice_s_dyn + watercap,4)
225write(*,*) "Computing the min of co2_ice..."
226where (co2_ice_slope_dyn < 0.) co2_ice_slope_dyn = 0.
227min_co2_ice_dyn = minval(co2_ice_slope_dyn + perennial_co2ice,4)
228
229! Compute averages over the year for each point
230write(*,*) "Computing the average of tsurf..."
231tsurf_avg_dyn = sum(tsurf_dyn,4)/timelen
232
233write(*,*) "Computing the average of Ps..."
234ps_avg_dyn = sum(ps_dyn,3)/timelen
235
236if (soil_pem) then
237    write(*,*) "Computing the average of tsoil..."
238    tsoil_avg_dyn = sum(tsoil_dyn,5)/timelen
239    write(*,*) "Computing the average of waterdensity_surface..."
240    watersurf_density_dyn_avg = sum(watersurf_density_dyn,4)/timelen
241endif
242
243! By definition, we get rid of the negative values
244A = (1./m_co2 - 1./m_noco2)
245B = 1./m_noco2
246do i = 1,iim_input + 1
247    do j = 1,jjm_input + 1
248        do t = 1, timelen
249            if (q_co2_dyn(i,j,t) < 0) then
250                q_co2_dyn(i,j,t) = 1.e-10
251            else if (q_co2_dyn(i,j,t) > 1) then
252                q_co2_dyn(i,j,t) = 1.
253            endif
254            if (q_h2o_dyn(i,j,t) < 0) then
255                q_h2o_dyn(i,j,t) = 1.e-10
256            else if (q_h2o_dyn(i,j,t) > 1) then
257                q_h2o_dyn(i,j,t) = 1.
258            endif
259            mmean = 1/(A*q_co2_dyn(i,j,t) + B)
260            vmr_co2_dyn(i,j,t) = q_co2_dyn(i,j,t)*mmean/m_co2
261        enddo
262    enddo
263enddo
264
265#ifndef CPP_1D
266    call gr_dyn_fi(timelen,iim_input + 1,jjm_input + 1,ngrid,vmr_co2_dyn,vmr_co2_PCM)
267    call gr_dyn_fi(timelen,iim_input + 1,jjm_input + 1,ngrid,ps_dyn,ps_timeseries)
268    call gr_dyn_fi(timelen,iim_input + 1,jjm_input + 1,ngrid,ps_avg_dyn,ps_avg)
269    call gr_dyn_fi(timelen,iim_input + 1,jjm_input + 1,ngrid,q_co2_dyn,q_co2)
270    call gr_dyn_fi(timelen,iim_input + 1,jjm_input + 1,ngrid,q_h2o_dyn,q_h2o)
271    do islope = 1,nslope
272        call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,min_co2_ice_dyn(:,:,islope),min_co2_ice(:,islope,2))
273        call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,min_h2o_ice_dyn(:,:,islope),min_h2o_ice(:,islope,2))
274        if (soil_pem) then
275            do l = 1,nsoilmx
276                call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,tsoil_avg_dyn(:,:,l,islope),tsoil_avg(:,l,islope))
277                do t = 1,timelen
278                    call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,tsoil_dyn(:,:,l,islope,t),tsoil_timeseries(:,l,islope,t))
279                    call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,watersoil_density_dyn(:,:,l,islope,t),watersoil_density_timeseries(:,l,islope,t))
280                enddo
281            enddo
282            call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,watersurf_density_dyn_avg(:,:,islope),watersurf_density_avg(:,islope))
283        endif ! soil_pem
284    enddo
285    call gr_dyn_fi(nslope,iim_input + 1,jjm_input + 1,ngrid,tsurf_avg_dyn,tsurf_avg)
286#else
287    vmr_co2_PCM(1,:) = vmr_co2_dyn(1,1,:)
288    ps_timeseries(1,:) = ps_dyn(1,1,:)
289    ps_avg(1) = ps_avg_dyn(1,1)
290    q_co2(1,:) = q_co2_dyn(1,1,:)
291    q_h2o(1,:) = q_h2o_dyn(1,1,:)
292    min_co2_ice(1,:,2) = min_co2_ice_dyn(1,1,:)
293    min_h2o_ice(1,:,2) = min_h2o_ice_dyn(1,1,:)
294    if (soil_pem) then
295        tsoil_timeseries(1,:,:,:) = tsoil_dyn(1,1,:,:,:)
296        tsoil_avg(1,:,:) = tsoil_avg_dyn(1,1,:,:)
297        watersoil_density_timeseries(1,:,:,:) = watersoil_density_dyn(1,1,:,:,:)
298        watersurf_density_avg(1,:) = watersurf_density_dyn_avg(1,1,:)
299    endif ! soil_pem
300    tsurf_avg(1,:) = tsurf_avg_dyn(1,1,:)
301#endif
302write(*,*) "Downloading data Y2 done!"
303
304END SUBROUTINE read_data_PCM
305
306!=======================================================================
307
308SUBROUTINE check_dim(n1,n2,str1,str2)
309
310implicit none
311
312integer,            intent(in) :: n1, n2
313character(len = *), intent(in) :: str1, str2
314
315character(256) :: s1, s2
316
317if (n1 /= n2) then
318    s1 = 'value of '//trim(str1)//' ='
319    s2 = ' read in starting file differs from parametrized '//trim(str2)//' ='
320    write(msg,'(10x,a,i4,2x,a,i4)')trim(s1),n1,trim(s2),n2
321    call abort_gcm(trim(modname),trim(msg),1)
322endif
323
324END SUBROUTINE check_dim
325
326!=======================================================================
327
328SUBROUTINE get_var1(var,v)
329
330implicit none
331
332character(len = *), intent(in)  :: var
333real, dimension(:), intent(out) :: v
334
335call error_msg(NF90_INQ_VARID(fID,var,vID),"inq",var)
336call error_msg(NF90_GET_VAR(fID,vID,v),"get",var)
337
338END SUBROUTINE get_var1
339
340!=======================================================================
341
342SUBROUTINE get_var3(var,v) ! on U grid
343
344implicit none
345
346character(len = *),     intent(in)  :: var
347real, dimension(:,:,:), intent(out) :: v
348
349call error_msg(NF90_INQ_VARID(fID,var,vID),"inq",var)
350call error_msg(NF90_GET_VAR(fID,vID,v),"get",var)
351
352END SUBROUTINE get_var3
353
354!=======================================================================
355
356SUBROUTINE get_var4(var,v)
357
358implicit none
359
360character(len = *),       intent(in)  :: var
361real, dimension(:,:,:,:), intent(out) :: v
362
363call error_msg(NF90_INQ_VARID(fID,var,vID),"inq",var)
364call error_msg(NF90_GET_VAR(fID,vID,v),"get",var)
365
366END SUBROUTINE get_var4
367
368!=======================================================================
369
370SUBROUTINE error_msg(ierr,typ,nam)
371
372implicit none
373
374integer,            intent(in) :: ierr !--- NetCDF ERROR CODE
375character(len = *), intent(in) :: typ  !--- TYPE OF OPERATION
376character(len = *), intent(in) :: nam  !--- FIELD/FILE NAME
377
378if (ierr == nf90_noerr) return
379select case(typ)
380    case('inq');   msg="Field <"//trim(nam)//"> is missing"
381    case('get');   msg="Reading failed for <"//trim(nam)//">"
382    case('put');   msg="Writing failed for <"//trim(nam)//">"
383    case('open');  msg="File opening failed for <"//trim(nam)//">"
384    case('close'); msg="File closing failed for <"//trim(nam)//">"
385    case default
386        write(*,*) 'There is no message for this error.'
387        error stop
388end select
389call abort_gcm(trim(modname),trim(msg),ierr)
390
391END SUBROUTINE error_msg
392
393END MODULE read_data_PCM_mod
Note: See TracBrowser for help on using the repository browser.