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

Last change on this file since 3345 was 3210, checked in by jbclement, 10 months ago

PEM:

  • Addition in the deftank of a bash script "modify_startfi_var.sh" to modify the value of a variable in a "startfi.nc".
  • Small corrections due to r3206.

JBC
File size: 15.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(filename,timelen,iim_input,jjm_input,ngrid,nslope,vmr_co2_PCM_phys,ps_timeseries,          &
18                         min_co2_ice,min_h2o_ice,tsurf_avg,tsoil_avg,tsurf_PCM,tsoil_PCM,q_co2,q_h2o,co2_ice_slope, &
19                         watersurf_density_avg,watersoil_density)
20use comsoil_h,             only: nsoilmx
21use comsoil_h_PEM,         only: soil_pem
22use constants_marspem_mod, only: m_co2, m_noco2
23
24implicit none
25
26!=======================================================================
27!
28! Purpose: Read initial confitions file from the PCM
29!
30! Authors: RV & LL
31!=======================================================================
32
33include "dimensions.h"
34
35!=======================================================================
36! Arguments:
37character(*), intent(in) :: filename                            ! File name
38integer,      intent(in) :: timelen                             ! number of times stored in the file
39integer,      intent(in) :: iim_input, jjm_input, ngrid, nslope ! number of points in the lat x lon dynamical grid, number of subgrid slopes
40! Ouputs
41real, dimension(ngrid,nslope),         intent(out) :: min_co2_ice      ! Minimum of co2 ice per slope of the year [kg/m^2]
42real, dimension(ngrid,nslope),         intent(out) :: min_h2o_ice      ! Minimum of h2o ice per slope of the year [kg/m^2]
43real, dimension(ngrid,timelen),        intent(out) :: vmr_co2_PCM_phys ! Physics x Times co2 volume mixing ratio retrieve from the PCM [m^3/m^3]
44real, dimension(ngrid,timelen),        intent(out) :: ps_timeseries    ! Surface Pressure [Pa]
45real, dimension(ngrid,timelen),        intent(out) :: q_co2            ! CO2 mass mixing ratio in the first layer [kg/m^3]
46real, dimension(ngrid,timelen),        intent(out) :: q_h2o            ! H2O mass mixing ratio in the first layer [kg/m^3]
47real, dimension(ngrid,nslope,timelen), intent(out) :: co2_ice_slope    ! co2 ice amount per slope of the year [kg/m^2]
48!SOIL
49real, dimension(ngrid,nslope),                 intent(out) :: tsurf_avg             ! Average surface temperature of the concatenated file [K]
50real, dimension(ngrid,nsoilmx,nslope),         intent(out) :: tsoil_avg             ! Average soil temperature of the concatenated file [K]
51real, dimension(ngrid,nslope,timelen),         intent(out) :: tsurf_PCM             ! Surface temperature of the concatenated file, time series [K]
52real, dimension(ngrid,nsoilmx,nslope,timelen), intent(out) :: tsoil_PCM             ! Soil temperature of the concatenated file, time series [K]
53real, dimension(ngrid,nslope),                 intent(out) :: watersurf_density_avg ! Water density at the surface [kg/m^3]
54real, dimension(ngrid,nsoilmx,nslope,timelen), intent(out) :: watersoil_density     ! Water density in the soil layer, time series [kg/m^3]
55!=======================================================================
56! Local Variables
57integer                         :: i, j, l, t                                                    ! loop variables
58real                            :: A, B, mmean                                                   ! Molar Mass of co2 and no co2, A;B intermediate variables to compute the mean molar mass of the layer
59integer                         :: islope                                                        ! loop for variables
60character(2)                    :: num                                                           ! for reading sloped variables
61real, dimension(iim_input + 1,jjm_input + 1,nslope,timelen)         :: h2o_ice_s_dyn             ! h2o ice per slope of the concatenated file [kg/m^2]
62real, dimension(iim_input + 1,jjm_input + 1,nslope,timelen)         :: watercap
63real, dimension(iim_input + 1,jjm_input + 1,nslope,timelen)         :: perennial_co2ice
64real, dimension(iim_input + 1,jjm_input + 1,timelen)                :: vmr_co2_PCM               ! CO2 volume mixing ratio in the first layer [mol/m^3]
65real, dimension(iim_input + 1,jjm_input + 1,timelen)                :: ps_PCM                    ! Surface Pressure [Pa]
66real, dimension(iim_input + 1,jjm_input + 1,nslope)                 :: min_co2_ice_dyn
67real, dimension(iim_input + 1,jjm_input + 1,nslope)                 :: min_h2o_ice_dyn
68real, dimension(iim_input + 1,jjm_input + 1,nslope)                 :: tsurf_avg_dyn             ! Average surface temperature of the concatenated file [K]
69real, dimension(iim_input + 1,jjm_input + 1,nsoilmx,nslope)         :: tsoil_avg_dyn             ! Average soil temperature of the concatenated file [K]
70real, dimension(iim_input + 1,jjm_input + 1,nslope,timelen)         :: tsurf_PCM_dyn             ! Surface temperature of the concatenated file, time series [K]
71real, dimension(iim_input + 1,jjm_input + 1,nsoilmx,nslope,timelen) :: tsoil_PCM_dyn             ! Soil temperature of the concatenated file, time series [K]
72real, dimension(iim_input + 1,jjm_input + 1,timelen)                :: q_co2_dyn                 ! CO2 mass mixing ratio in the first layer [kg/m^3]
73real, dimension(iim_input + 1,jjm_input + 1,timelen)                :: q_h2o_dyn                 ! H2O mass mixing ratio in the first layer [kg/m^3]
74real, dimension(iim_input + 1,jjm_input + 1,nslope,timelen)         :: co2_ice_slope_dyn         ! co2 ice amount per  slope of the year [kg/m^2]
75real, dimension(iim_input + 1,jjm_input + 1,nslope,timelen)         :: watersurf_density_dyn     ! Water density at the surface, time series [kg/m^3]
76real, dimension(iim_input + 1,jjm_input + 1,nslope)                 :: watersurf_density_dyn_avg ! Water density at the surface, dynamic grid, yearly averaged [kg/m^3]
77real, dimension(iim_input + 1,jjm_input + 1,nsoilmx,nslope,timelen) :: watersoil_density_dyn     ! Water density in the soil layer, time series [kg/m^3]
78!-----------------------------------------------------------------------
79! Open the NetCDF file
80write(*,*) "Opening "//filename//"..."
81call error_msg(NF90_OPEN(filename,NF90_NOWRITE,fID),"open",filename)
82
83! Dowload the data from the file
84call get_var3("ps",ps_PCM)
85write(*,*) "Data for surface pressure downloaded!"
86
87call get_var3("co2_layer1",q_co2_dyn)
88write(*,*) "Data for vmr co2 downloaded!"
89
90call get_var3("h2o_layer1",q_h2o_dyn)
91write(*,*) "Data for vmr h2o downloaded!"
92
93if (nslope == 1) then ! There is no slope
94    call get_var3("co2ice",co2_ice_slope_dyn(:,:,1,:))
95    write(*,*) "Data for co2_ice downloaded!"
96
97    call get_var3("h2o_ice_s",h2o_ice_s_dyn(:,:,1,:))
98    write(*,*) "Data for h2o_ice_s downloaded!"
99
100    call get_var3("tsurf",tsurf_PCM_dyn(:,:,1,:))
101    write(*,*) "Data for tsurf downloaded!"
102
103#ifndef CPP_STD
104    call get_var3("watercap",watercap(:,:,1,:))
105    write(*,*) "Data for watercap downloaded!"
106
107    call get_var3("perennial_co2ice",perennial_co2ice(:,:,1,:))
108    write(*,*) "Data for perennial_co2ice downloaded!"
109
110    if (soil_pem) then
111        call get_var4("soiltemp",tsoil_PCM_dyn(:,:,:,1,:))
112        write(*,*) "Data for soiltemp downloaded!"
113
114        call get_var4("waterdensity_soil",watersoil_density_dyn(:,:,:,1,:))
115        write(*,*) "Data for waterdensity_soil downloaded!"
116
117        call get_var3("waterdensity_surface",watersurf_density_dyn(:,:,1,:))
118        write(*,*) "Data for waterdensity_surface downloaded!"
119    endif !soil_pem
120#endif
121else ! We use slopes
122    do islope = 1,nslope
123        write(num,'(i2.2)') islope
124        call get_var3("co2ice_slope"//num,co2_ice_slope_dyn(:,:,islope,:))
125    enddo
126    write(*,*) "Data for co2_ice downloaded!"
127
128    do islope = 1,nslope
129        write(num,'(i2.2)') islope
130        call get_var3("h2o_ice_s_slope"//num,h2o_ice_s_dyn(:,:,islope,:))
131    enddo
132    write(*,*) "Data for h2o_ice_s downloaded!"
133
134    do islope = 1,nslope
135        write(num,'(i2.2)') islope
136        call get_var3("tsurf_slope"//num,tsurf_PCM_dyn(:,:,islope,:))
137    enddo
138    write(*,*) "Data for tsurf downloaded!"
139
140#ifndef CPP_STD
141    do islope = 1,nslope
142        write(num,'(i2.2)') islope
143        call get_var3("watercap_slope"//num,watercap(:,:,islope,:))
144    enddo
145    write(*,*) "Data for watercap downloaded!"
146
147    do islope = 1,nslope
148        write(num,'(i2.2)') islope
149        call get_var3("perennial_co2ice_slope"//num,perennial_co2ice(:,:,islope,:))
150    enddo
151    write(*,*) "Data for perennial_co2ice downloaded!"
152
153    if (soil_pem) then
154        do islope = 1,nslope
155            write(num,'(i2.2)') islope
156            call get_var4("soiltemp_slope"//num,tsoil_PCM_dyn(:,:,:,islope,:))
157        enddo
158        write(*,*) "Data for soiltemp downloaded!"
159
160        do islope = 1,nslope
161            write(num,'(i2.2)') islope
162            call get_var4("waterdensity_soil_slope"//num,watersoil_density_dyn(:,:,:,islope,:))
163        enddo
164        write(*,*) "Data for waterdensity_soil downloaded!"
165
166        do islope = 1,nslope
167            write(num,'(i2.2)') islope
168            call get_var3("waterdensity_surface"//num,watersurf_density_dyn(:,:,islope,:))
169        enddo
170        write(*,*) "Data for waterdensity_surface downloaded!"
171    endif !soil_pem
172#endif
173endif
174
175! Compute the minimum over the year for each point
176write(*,*) "Computing the min of h2o_ice..."
177where (h2o_ice_s_dyn < 0.) h2o_ice_s_dyn = 0.
178min_h2o_ice_dyn = minval(h2o_ice_s_dyn + watercap,4)
179write(*,*) "Computing the min of co2_ice..."
180where (co2_ice_slope_dyn < 0.) co2_ice_slope_dyn = 0.
181min_co2_ice_dyn = minval(co2_ice_slope_dyn + perennial_co2ice,4)
182
183! Compute averages over the year for each point
184write(*,*) "Computing the average of tsurf..."
185tsurf_avg_dyn = sum(tsurf_PCM_dyn,4)/timelen
186
187if (soil_pem) then
188    write(*,*) "Computing average of tsoil..."
189    tsoil_avg_dyn = sum(tsoil_PCM_dyn,5)/timelen
190    write(*,*) "Computing average of waterdensity_surface..."
191    watersurf_density_dyn_avg = sum(watersurf_density_dyn,4)/timelen
192endif
193
194! By definition, we get rid of the negative values
195A = (1./m_co2 - 1./m_noco2)
196B = 1./m_noco2
197do i = 1,iim + 1
198    do j = 1,jjm_input + 1
199        do t = 1, timelen
200            if (q_co2_dyn(i,j,t) < 0) then
201                q_co2_dyn(i,j,t) = 1.e-10
202            else if (q_co2_dyn(i,j,t) > 1) then
203                q_co2_dyn(i,j,t) = 1.
204            endif
205            if (q_h2o_dyn(i,j,t) < 0) then
206                q_h2o_dyn(i,j,t) = 1.e-10
207            else if (q_h2o_dyn(i,j,t) > 1) then
208                q_h2o_dyn(i,j,t) = 1.
209            endif
210            mmean = 1/(A*q_co2_dyn(i,j,t) + B)
211            vmr_co2_PCM(i,j,t) = q_co2_dyn(i,j,t)*mmean/m_co2
212        enddo
213    enddo
214enddo
215
216#ifndef CPP_1D
217    call gr_dyn_fi(timelen,iim_input + 1,jjm_input + 1,ngrid,vmr_co2_PCM,vmr_co2_PCM_phys)
218    call gr_dyn_fi(timelen,iim_input + 1,jjm_input + 1,ngrid,ps_PCM,ps_timeseries)
219    call gr_dyn_fi(timelen,iim_input + 1,jjm_input + 1,ngrid,q_co2_dyn,q_co2)
220    call gr_dyn_fi(timelen,iim_input + 1,jjm_input + 1,ngrid,q_h2o_dyn,q_h2o)
221
222    do islope = 1,nslope
223        call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,min_co2_ice_dyn(:,:,islope),min_co2_ice(:,islope))
224        call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,min_h2o_ice_dyn(:,:,islope),min_h2o_ice(:,islope))
225        if (soil_pem) then
226            do l = 1,nsoilmx
227                call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,tsoil_avg_dyn(:,:,l,islope),tsoil_avg(:,l,islope))
228                do t = 1,timelen
229                    call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,tsoil_PCM_dyn(:,:,l,islope,t),tsoil_PCM(:,l,islope,t))
230                    call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,watersoil_density_dyn(:,:,l,islope,t),watersoil_density(:,l,islope,t))
231                enddo
232            enddo
233            call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,watersurf_density_dyn_avg(:,:,islope),watersurf_density_avg(:,islope))
234        endif ! soil_pem
235        do t = 1,timelen
236            call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,tsurf_PCM_dyn(:,:,islope,t),tsurf_PCM(:,islope,t))
237            call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,co2_ice_slope_dyn(:,:,islope,t),co2_ice_slope(:,islope,t))
238        enddo
239    enddo
240
241    call gr_dyn_fi(nslope,iim_input + 1,jjm_input + 1,ngrid,tsurf_avg_dyn,tsurf_avg)
242#else
243    vmr_co2_PCM_phys(1,:) = vmr_co2_PCM(1,1,:)
244    ps_timeseries(1,:) = ps_PCM(1,1,:)
245    q_co2(1,:) = q_co2_dyn(1,1,:)
246    q_h2o(1,:) = q_h2o_dyn(1,1,:)
247    min_co2_ice(1,:) = min_co2_ice_dyn(1,1,:)
248    min_h2o_ice(1,:) = min_h2o_ice_dyn(1,1,:)
249    if (soil_pem) then
250        tsoil_avg(1,:,:) = tsoil_avg_dyn(1,1,:,:)
251        tsoil_PCM(1,:,:,:) = tsoil_PCM_dyn(1,1,:,:,:)
252        watersoil_density(1,:,:,:) = watersoil_density_dyn(1,1,:,:,:)
253        watersurf_density_avg(1,:) = watersurf_density_dyn_avg(1,1,:)
254    endif ! soil_pem
255    tsurf_PCM(1,:,:) = tsurf_PCM_dyn(1,1,:,:)
256    co2_ice_slope(1,:,:) = co2_ice_slope_dyn(1,1,:,:)
257    tsurf_avg(1,:) = tsurf_avg_dyn(1,1,:)
258#endif
259
260END SUBROUTINE read_data_PCM
261
262!=======================================================================
263
264SUBROUTINE check_dim(n1,n2,str1,str2)
265
266implicit none
267
268integer,            intent(in) :: n1, n2
269character(len = *), intent(in) :: str1, str2
270
271character(256) :: s1, s2
272
273if (n1 /= n2) then
274    s1 = 'value of '//trim(str1)//' ='
275    s2 = ' read in starting file differs from parametrized '//trim(str2)//' ='
276    write(msg,'(10x,a,i4,2x,a,i4)')trim(s1),n1,trim(s2),n2
277    call abort_gcm(trim(modname),trim(msg),1)
278endif
279
280END SUBROUTINE check_dim
281
282!=======================================================================
283
284SUBROUTINE get_var1(var,v)
285
286implicit none
287
288character(len = *), intent(in)  :: var
289real, dimension(:), intent(out) :: v
290
291call error_msg(NF90_INQ_VARID(fID,var,vID),"inq",var)
292call error_msg(NF90_GET_VAR(fID,vID,v),"get",var)
293
294END SUBROUTINE get_var1
295
296!=======================================================================
297
298SUBROUTINE get_var3(var,v) ! on U grid
299
300implicit none
301
302character(len = *),     intent(in)  :: var
303real, dimension(:,:,:), intent(out) :: v
304
305call error_msg(NF90_INQ_VARID(fID,var,vID),"inq",var)
306call error_msg(NF90_GET_VAR(fID,vID,v),"get",var)
307
308END SUBROUTINE get_var3
309
310!=======================================================================
311
312SUBROUTINE get_var4(var,v)
313
314implicit none
315
316character(len = *),       intent(in)  :: var
317real, dimension(:,:,:,:), intent(out) :: v
318
319call error_msg(NF90_INQ_VARID(fID,var,vID),"inq",var)
320call error_msg(NF90_GET_VAR(fID,vID,v),"get",var)
321
322END SUBROUTINE get_var4
323
324!=======================================================================
325
326SUBROUTINE error_msg(ierr,typ,nam)
327
328implicit none
329
330integer,            intent(in) :: ierr !--- NetCDF ERROR CODE
331character(len = *), intent(in) :: typ  !--- TYPE OF OPERATION
332character(len = *), intent(in) :: nam  !--- FIELD/FILE NAME
333
334if (ierr == nf90_noerr) return
335select case(typ)
336    case('inq');   msg="Field <"//trim(nam)//"> is missing"
337    case('get');   msg="Reading failed for <"//trim(nam)//">"
338    case('open');  msg="File opening failed for <"//trim(nam)//">"
339    case('close'); msg="File closing failed for <"//trim(nam)//">"
340    case default
341        write(*,*) 'There is no message for this error.'
342        error stop
343end select
344call abort_gcm(trim(modname),trim(msg),ierr)
345
346END SUBROUTINE error_msg
347
348END MODULE read_data_PCM_mod
Note: See TracBrowser for help on using the repository browser.