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

Last change on this file since 3598 was 3598, checked in by jbclement, 11 days ago

PEM:
Correction of typos and minor bugs.
JBC
File size: 18.4 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    enddo
105    write(*,*) "Data for co2_ice downloaded."
106
107    do islope = 1,nslope
108        write(num,'(i2.2)') islope
109        call get_var3("h2o_ice_s_slope"//num,h2o_ice_s_dyn(:,:,islope,:))
110    enddo
111    write(*,*) "Data for h2o_ice_s downloaded."
112
113    do islope = 1,nslope
114        write(num,'(i2.2)') islope
115        call get_var3("tsurf_slope"//num,tsurf_dyn(:,:,islope,:))
116    enddo
117    write(*,*) "Data for tsurf downloaded."
118
119#ifndef CPP_STD
120    do islope = 1,nslope
121        write(num,'(i2.2)') islope
122        call get_var3("watercap_slope"//num,watercap(:,:,islope,:))
123    enddo
124    write(*,*) "Data for watercap downloaded."
125
126    do islope = 1,nslope
127        write(num,'(i2.2)') islope
128        call get_var3("perennial_co2ice_slope"//num,perennial_co2ice(:,:,islope,:))
129    enddo
130    write(*,*) "Data for perennial_co2ice downloaded."
131#endif
132endif
133
134! Close the NetCDF file
135write(*,*) "Closing "//filename1//"..."
136call error_msg(nf90_close(fID),"close",filename1)
137
138! Compute the minimum over the year for each point
139write(*,*) "Computing the min of h2o_ice..."
140where (h2o_ice_s_dyn < 0.) h2o_ice_s_dyn = 0.
141min_h2o_ice_dyn = minval(h2o_ice_s_dyn + watercap,4)
142write(*,*) "Computing the min of co2_ice..."
143where (co2_ice_slope_dyn < 0.) co2_ice_slope_dyn = 0.
144min_co2_ice_dyn = minval(co2_ice_slope_dyn + perennial_co2ice,4)
145
146! Compute averages over the year for each point
147write(*,*) "Computing the average of tsurf..."
148tsurf_avg_dyn = sum(tsurf_dyn,4)/timelen
149
150#ifndef CPP_1D
151    call gr_dyn_fi(nslope,iim_input + 1,jjm_input + 1,ngrid,tsurf_avg_dyn,tsurf_avg)
152    do islope = 1,nslope
153        call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,min_co2_ice_dyn(:,:,islope),min_co2_ice(:,islope,1))
154        call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,min_h2o_ice_dyn(:,:,islope),min_h2o_ice(:,islope,1))
155    enddo
156#else
157    tsurf_avg_yr1(1,:) = tsurf_avg_dyn(1,1,:)
158    min_co2_ice(1,:,1) = min_co2_ice_dyn(1,1,:)
159    min_h2o_ice(1,:,1) = min_h2o_ice_dyn(1,1,:)
160#endif
161write(*,*) "Downloading data Y1 done!"
162
163!------------------ Year 2 ------------------
164
165! Open the NetCDF file
166write(*,*) "Opening "//filename2//"..."
167call error_msg(NF90_OPEN(filename2,NF90_NOWRITE,fID),"open",filename2)
168
169! Download the data from the file
170call get_var3("ps",ps_dyn)
171write(*,*) "Data for surface pressure downloaded."
172
173call get_var3("co2_layer1",q_co2_dyn)
174write(*,*) "Data for vmr co2 downloaded."
175
176call get_var3("h2o_layer1",q_h2o_dyn)
177write(*,*) "Data for vmr h2o downloaded."
178
179if (nslope == 1) then ! There is no slope
180    call get_var3("co2ice",co2_ice_slope_dyn(:,:,1,:))
181    write(*,*) "Data for co2_ice downloaded."
182
183    call get_var3("h2o_ice_s",h2o_ice_s_dyn(:,:,1,:))
184    write(*,*) "Data for h2o_ice_s downloaded."
185
186    call get_var3("tsurf",tsurf_dyn(:,:,1,:))
187    write(*,*) "Data for tsurf downloaded."
188
189#ifndef CPP_STD
190    call get_var3("watercap",watercap(:,:,1,:))
191    write(*,*) "Data for watercap downloaded."
192
193    call get_var3("perennial_co2ice",perennial_co2ice(:,:,1,:))
194    write(*,*) "Data for perennial_co2ice downloaded."
195
196    if (soil_pem) then
197        call get_var4("soiltemp",tsoil_dyn(:,:,:,1,:))
198        write(*,*) "Data for soiltemp downloaded."
199
200        call get_var4("waterdensity_soil",watersoil_density_dyn(:,:,:,1,:))
201        write(*,*) "Data for waterdensity_soil downloaded."
202
203        call get_var3("waterdensity_surface",watersurf_density_dyn(:,:,1,:))
204        write(*,*) "Data for waterdensity_surface downloaded."
205    endif !soil_pem
206#endif
207else ! We use slopes
208    do islope = 1,nslope
209        write(num,'(i2.2)') islope
210        call get_var3("co2ice_slope"//num,co2_ice_slope_dyn(:,:,islope,:))
211    enddo
212    write(*,*) "Data for co2_ice downloaded."
213
214    do islope = 1,nslope
215        write(num,'(i2.2)') islope
216        call get_var3("h2o_ice_s_slope"//num,h2o_ice_s_dyn(:,:,islope,:))
217    enddo
218    write(*,*) "Data for h2o_ice_s downloaded."
219
220    do islope = 1,nslope
221        write(num,'(i2.2)') islope
222        call get_var3("tsurf_slope"//num,tsurf_dyn(:,:,islope,:))
223    enddo
224    write(*,*) "Data for tsurf downloaded."
225
226#ifndef CPP_STD
227    do islope = 1,nslope
228        write(num,'(i2.2)') islope
229        call get_var3("watercap_slope"//num,watercap(:,:,islope,:))
230    enddo
231    write(*,*) "Data for watercap downloaded."
232
233    do islope = 1,nslope
234        write(num,'(i2.2)') islope
235        call get_var3("perennial_co2ice_slope"//num,perennial_co2ice(:,:,islope,:))
236    enddo
237    write(*,*) "Data for perennial_co2ice downloaded."
238
239    if (soil_pem) then
240        do islope = 1,nslope
241            write(num,'(i2.2)') islope
242            call get_var4("soiltemp_slope"//num,tsoil_dyn(:,:,:,islope,:))
243        enddo
244        write(*,*) "Data for soiltemp downloaded."
245
246        do islope = 1,nslope
247            write(num,'(i2.2)') islope
248            call get_var4("waterdensity_soil_slope"//num,watersoil_density_dyn(:,:,:,islope,:))
249        enddo
250        write(*,*) "Data for waterdensity_soil downloaded."
251
252        do islope = 1,nslope
253            write(num,'(i2.2)') islope
254            call get_var3("waterdensity_surface"//num,watersurf_density_dyn(:,:,islope,:))
255        enddo
256        write(*,*) "Data for waterdensity_surface downloaded."
257    endif !soil_pem
258#endif
259endif
260
261! Close the NetCDF file
262write(*,*) "Closing "//filename2//"..."
263call error_msg(nf90_close(fID),"close",filename2)
264
265! Compute the minimum over the year for each point
266write(*,*) "Computing the min of h2o_ice..."
267where (h2o_ice_s_dyn < 0.) h2o_ice_s_dyn = 0.
268min_h2o_ice_dyn = minval(h2o_ice_s_dyn + watercap,4)
269write(*,*) "Computing the min of co2_ice..."
270where (co2_ice_slope_dyn < 0.) co2_ice_slope_dyn = 0.
271min_co2_ice_dyn = minval(co2_ice_slope_dyn + perennial_co2ice,4)
272
273! Compute averages over the year for each point
274write(*,*) "Computing the average of tsurf..."
275tsurf_avg_dyn = sum(tsurf_dyn,4)/timelen
276
277write(*,*) "Computing the average of Ps..."
278ps_avg_dyn = sum(ps_dyn,3)/timelen
279
280if (soil_pem) then
281    write(*,*) "Computing the average of tsoil..."
282    tsoil_avg_dyn = sum(tsoil_dyn,5)/timelen
283    write(*,*) "Computing the average of waterdensity_surface..."
284    watersurf_density_dyn_avg = sum(watersurf_density_dyn,4)/timelen
285endif
286
287! By definition, we get rid of the negative values
288A = (1./m_co2 - 1./m_noco2)
289B = 1./m_noco2
290do i = 1,iim_input + 1
291    do j = 1,jjm_input + 1
292        do t = 1, timelen
293            if (q_co2_dyn(i,j,t) < 0) then
294                q_co2_dyn(i,j,t) = 1.e-10
295            else if (q_co2_dyn(i,j,t) > 1) then
296                q_co2_dyn(i,j,t) = 1.
297            endif
298            if (q_h2o_dyn(i,j,t) < 0) then
299                q_h2o_dyn(i,j,t) = 1.e-10
300            else if (q_h2o_dyn(i,j,t) > 1) then
301                q_h2o_dyn(i,j,t) = 1.
302            endif
303            mmean = 1/(A*q_co2_dyn(i,j,t) + B)
304            vmr_co2_dyn(i,j,t) = q_co2_dyn(i,j,t)*mmean/m_co2
305        enddo
306    enddo
307enddo
308
309#ifndef CPP_1D
310    call gr_dyn_fi(timelen,iim_input + 1,jjm_input + 1,ngrid,vmr_co2_dyn,vmr_co2_PCM)
311    call gr_dyn_fi(timelen,iim_input + 1,jjm_input + 1,ngrid,ps_dyn,ps_timeseries)
312    call gr_dyn_fi(timelen,iim_input + 1,jjm_input + 1,ngrid,ps_avg_dyn,ps_avg)
313    call gr_dyn_fi(timelen,iim_input + 1,jjm_input + 1,ngrid,q_co2_dyn,q_co2)
314    call gr_dyn_fi(timelen,iim_input + 1,jjm_input + 1,ngrid,q_h2o_dyn,q_h2o)
315    do islope = 1,nslope
316        call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,min_co2_ice_dyn(:,:,islope),min_co2_ice(:,islope,2))
317        call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,min_h2o_ice_dyn(:,:,islope),min_h2o_ice(:,islope,2))
318        if (soil_pem) then
319            do l = 1,nsoilmx
320                call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,tsoil_avg_dyn(:,:,l,islope),tsoil_avg(:,l,islope))
321                do t = 1,timelen
322                    call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,tsoil_dyn(:,:,l,islope,t),tsoil_timeseries(:,l,islope,t))
323                    call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,watersoil_density_dyn(:,:,l,islope,t),watersoil_density_timeseries(:,l,islope,t))
324                enddo
325            enddo
326            call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,watersurf_density_dyn_avg(:,:,islope),watersurf_density_avg(:,islope))
327        endif ! soil_pem
328    enddo
329    call gr_dyn_fi(nslope,iim_input + 1,jjm_input + 1,ngrid,tsurf_avg_dyn,tsurf_avg)
330#else
331    vmr_co2_PCM(1,:) = vmr_co2_dyn(1,1,:)
332    ps_timeseries(1,:) = ps_dyn(1,1,:)
333    ps_avg(1) = ps_avg_dyn(1,1)
334    q_co2(1,:) = q_co2_dyn(1,1,:)
335    q_h2o(1,:) = q_h2o_dyn(1,1,:)
336    min_co2_ice(1,:,2) = min_co2_ice_dyn(1,1,:)
337    min_h2o_ice(1,:,2) = min_h2o_ice_dyn(1,1,:)
338    if (soil_pem) then
339        tsoil_timeseries(1,:,:,:) = tsoil_dyn(1,1,:,:,:)
340        tsoil_avg(1,:,:) = tsoil_avg_dyn(1,1,:,:)
341        watersoil_density_timeseries(1,:,:,:) = watersoil_density_dyn(1,1,:,:,:)
342        watersurf_density_avg(1,:) = watersurf_density_dyn_avg(1,1,:)
343    endif ! soil_pem
344    tsurf_avg(1,:) = tsurf_avg_dyn(1,1,:)
345#endif
346write(*,*) "Downloading data Y2 done!"
347
348END SUBROUTINE read_data_PCM
349
350!=======================================================================
351
352SUBROUTINE check_dim(n1,n2,str1,str2)
353
354implicit none
355
356integer,            intent(in) :: n1, n2
357character(len = *), intent(in) :: str1, str2
358
359character(256) :: s1, s2
360
361if (n1 /= n2) then
362    s1 = 'value of '//trim(str1)//' ='
363    s2 = ' read in starting file differs from parametrized '//trim(str2)//' ='
364    write(msg,'(10x,a,i4,2x,a,i4)')trim(s1),n1,trim(s2),n2
365    call abort_gcm(trim(modname),trim(msg),1)
366endif
367
368END SUBROUTINE check_dim
369
370!=======================================================================
371
372SUBROUTINE get_var1(var,v)
373
374implicit none
375
376character(len = *), intent(in)  :: var
377real, dimension(:), intent(out) :: v
378
379call error_msg(NF90_INQ_VARID(fID,var,vID),"inq",var)
380call error_msg(NF90_GET_VAR(fID,vID,v),"get",var)
381
382END SUBROUTINE get_var1
383
384!=======================================================================
385
386SUBROUTINE get_var3(var,v) ! on U grid
387
388implicit none
389
390character(len = *),     intent(in)  :: var
391real, dimension(:,:,:), intent(out) :: v
392
393call error_msg(NF90_INQ_VARID(fID,var,vID),"inq",var)
394call error_msg(NF90_GET_VAR(fID,vID,v),"get",var)
395
396END SUBROUTINE get_var3
397
398!=======================================================================
399
400SUBROUTINE get_var4(var,v)
401
402implicit none
403
404character(len = *),       intent(in)  :: var
405real, dimension(:,:,:,:), intent(out) :: v
406
407call error_msg(NF90_INQ_VARID(fID,var,vID),"inq",var)
408call error_msg(NF90_GET_VAR(fID,vID,v),"get",var)
409
410END SUBROUTINE get_var4
411
412!=======================================================================
413
414SUBROUTINE error_msg(ierr,typ,nam)
415
416implicit none
417
418integer,            intent(in) :: ierr !--- NetCDF ERROR CODE
419character(len = *), intent(in) :: typ  !--- TYPE OF OPERATION
420character(len = *), intent(in) :: nam  !--- FIELD/FILE NAME
421
422if (ierr == nf90_noerr) return
423select case(typ)
424    case('inq');   msg="Field <"//trim(nam)//"> is missing"
425    case('get');   msg="Reading failed for <"//trim(nam)//">"
426    case('put');   msg="Writing failed for <"//trim(nam)//">"
427    case('open');  msg="File opening failed for <"//trim(nam)//">"
428    case('close'); msg="File closing failed for <"//trim(nam)//">"
429    case default
430        write(*,*) 'There is no message for this error.'
431        error stop
432end select
433call abort_gcm(trim(modname),trim(msg),ierr)
434
435END SUBROUTINE error_msg
436
437END MODULE read_data_PCM_mod
Note: See TracBrowser for help on using the repository browser.