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

Last change on this file since 3583 was 3571, checked in by jbclement, 11 months ago

PEM:

  • New way to manage the pressure: now the PEM manages only the average pressure and keeps the pressure deviation with the instantaneous pressure from the start to reconstruct the pressure at the end ('ps_avg = ps_start + ps_dev'). As a consequence, everything related to pressure in the PEM is modified accordingly.
  • Surface temperatures management is now simpler. It follows the strategy for the pressure (and soil temperature) described above.
  • Soil temperatures are now adapted to match the surface temperature changes occured during the PEM by modifying the soil temperature deviation at the end.
  • Few simplifications/optimizations: notably, the two PCM years are now read in one go in 'read_data_PCM_mod.F90' and only the needed variables are extracted.
  • Deletion of unused variables and unnecessary intermediate variables (memory saving and loop deletion in some cases).
  • Renaming of variables and subroutines to make everything clearer. In particular, the suffixes: '_avg' = average, '_start' = PCM start file, '_dev' = deviation, '_ini' or '0' = initial, '_dyn' = dynamical grid, '_timeseries' = daily average of last PCM year.
  • Cosmetic cleanings for readability.

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