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

Last change on this file since 3129 was 3123, checked in by llange, 2 years ago

MARS PEM
1) Adapting the Mars PEM soil grid to the one of the PCM. The first layers in the PEM follow those from the PCM (PYM grid), and then, for layers at depth, we use the classical power low grid.
2) Correction when reading the soil temperature, water density at the surface, and initialization of the ice table.
LL
File size: 15.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(filename,timelen,iim_input,jjm_input,ngrid,nslope,vmr_co2_gcm_phys,ps_timeseries,          &
18                         min_co2_ice,min_h2o_ice,tsurf_ave,tsoil_ave,tsurf_gcm,tsoil_gcm,q_co2,q_h2o,co2_ice_slope, &
19                         watersurf_density_ave,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:
37
38character(len = *), intent(in) :: filename                            ! File name
39integer,            intent(in) :: timelen                             ! number of times stored in the file
40integer                        :: iim_input, jjm_input, ngrid, nslope ! number of points in the lat x lon dynamical grid, number of subgrid slopes
41! Ouputs
42real, dimension(ngrid,nslope),         intent(out) :: min_co2_ice      ! Minimum of co2 ice  per slope of the year [kg/m^2]
43real, dimension(ngrid,nslope),         intent(out) :: min_h2o_ice      ! Minimum of h2o ice per slope of the year [kg/m^2]
44real, dimension(ngrid,timelen),        intent(out) :: vmr_co2_gcm_phys ! Physics x Times  co2 volume mixing ratio retrieve from the PCM [m^3/m^3]
45real, dimension(ngrid,timelen),        intent(out) :: ps_timeseries    ! Surface Pressure [Pa]
46real, dimension(ngrid,timelen),        intent(out) :: q_co2            ! CO2 mass mixing ratio in the first layer [kg/m^3]
47real, dimension(ngrid,timelen),        intent(out) :: q_h2o            ! H2O mass mixing ratio in the first layer [kg/m^3]
48real, dimension(ngrid,nslope,timelen), intent(out) :: co2_ice_slope    ! co2 ice amount per  slope of the year [kg/m^2]
49!SOIL
50real, dimension(ngrid,nslope),                 intent(out) :: tsurf_ave             ! Average surface temperature of the concatenated file [K]
51real, dimension(ngrid,nsoilmx,nslope),         intent(out) :: tsoil_ave             ! Average soil temperature of the concatenated file [K]
52real, dimension(ngrid,nslope,timelen),         intent(out) :: tsurf_gcm             ! Surface temperature of the concatenated file, time series [K]
53real, dimension(ngrid,nsoilmx,nslope,timelen), intent(out) :: tsoil_gcm             ! Soil temperature of the concatenated file, time series [K]
54real, dimension(ngrid,nslope),                 intent(out) :: watersurf_density_ave ! Water density at the surface [kg/m^3]
55real, dimension(ngrid,nsoilmx,nslope,timelen), intent(out) :: watersoil_density     ! Water density in the soil layer, time series [kg/m^3]
56!=======================================================================
57! Local Variables
58character(12)                   :: start_file_type = "earth" ! default start file type
59real, dimension(:), allocatable :: time      ! times stored in start
60integer                         :: indextime ! index of selected time
61integer                         :: edges(4), corner(4)
62integer                         :: i, j, l, t          ! loop variables
63real                            :: A, B, mmean         ! Molar Mass of co2 and no co2, A;B intermediate variables to compute the mean molar mass of the layer
64integer                         :: islope              ! loop for variables
65character(2)                    :: num                 ! for reading sloped variables
66real, dimension(iim_input + 1,jjm_input + 1,nslope,timelen)         :: h2o_ice_s_dyn             ! h2o ice per slope of the concatenated file [kg/m^2]
67real, dimension(iim_input + 1,jjm_input + 1,nslope,timelen)         :: watercap
68real, dimension(iim_input + 1,jjm_input + 1,timelen)                :: vmr_co2_gcm               ! CO2 volume mixing ratio in the first layer [mol/m^3]
69real, dimension(iim_input + 1,jjm_input + 1,timelen)                :: ps_GCM                    ! Surface Pressure [Pa]
70real, dimension(iim_input + 1,jjm_input + 1,nslope)                 :: min_co2_ice_dyn
71real, dimension(iim_input + 1,jjm_input + 1,nslope)                 :: min_h2o_ice_dyn
72real, dimension(iim_input + 1,jjm_input + 1,nslope)                 :: tsurf_ave_dyn             ! Average surface temperature of the concatenated file [K]
73real, dimension(iim_input + 1,jjm_input + 1,nsoilmx,nslope)         :: tsoil_ave_dyn             ! Average soil temperature of the concatenated file [K]
74real, dimension(iim_input + 1,jjm_input + 1,nslope,timelen)         :: tsurf_gcm_dyn             ! Surface temperature of the concatenated file, time series [K]
75real, dimension(iim_input + 1,jjm_input + 1,nsoilmx,nslope,timelen) :: tsoil_gcm_dyn             ! Soil temperature of the concatenated file, time series [K]
76real, dimension(iim_input + 1,jjm_input + 1,timelen)                :: q_co2_dyn                 ! CO2 mass mixing ratio in the first layer [kg/m^3]
77real, dimension(iim_input + 1,jjm_input + 1,timelen)                :: q_h2o_dyn                 ! H2O mass mixing ratio in the first layer [kg/m^3]
78real, dimension(iim_input + 1,jjm_input + 1,nslope,timelen)         :: co2_ice_slope_dyn         ! co2 ice amount per  slope of the year [kg/m^2]
79real, dimension(iim_input + 1,jjm_input + 1,nslope,timelen)         :: watersurf_density_dyn     ! Water density at the surface, time series [kg/m^3]
80real, dimension(iim_input + 1,jjm_input + 1,nslope)                 :: watersurf_density_dyn_ave ! Water density at the surface, dynamic grid, yearly averaged [kg/m^3]
81real, dimension(iim_input + 1,jjm_input + 1,nsoilmx,nslope,timelen) :: watersoil_density_dyn     ! Water density in the soil layer, time series [kg/m^3]
82real, dimension(ngrid,nslope,timelen)                               :: watersurf_density         ! Water density at the surface, physical grid, time series [kg/m^3]
83!-----------------------------------------------------------------------
84! Open the NetCDF file
85write(*,*) "Opening "//filename//"..."
86call error_msg(NF90_OPEN(filename,NF90_NOWRITE,fID),"open",filename)
87
88! Dowload the data from the file
89call get_var3("ps",ps_GCM)
90write(*,*) "Data for surface pressure downloaded!"
91
92call get_var3("co2_layer1",q_co2_dyn)
93write(*,*) "Data for vmr co2 downloaded!"
94
95call get_var3("h2o_layer1",q_h2o_dyn)
96write(*,*) "Data for vmr h2o downloaded!"
97
98if (nslope == 1) then ! There is no slope
99    call get_var3("co2ice",co2_ice_slope_dyn(:,:,1,:))
100    write(*,*) "Data for co2_ice downloaded!"
101
102    call get_var3("h2o_ice_s",h2o_ice_s_dyn(:,:,1,:))
103    write(*,*) "Data for h2o_ice_s downloaded!"
104
105#ifndef CPP_STD
106    call get_var3("watercap",watercap(:,:,1,:))
107    write(*,*) "Data for watercap downloaded!"
108#endif
109
110    call get_var3("tsurf",tsurf_gcm_dyn(:,:,1,:))
111    write(*,*) "Data for tsurf downloaded!"
112
113#ifndef CPP_STD
114    if (soil_pem) then
115        call get_var4("soiltemp",tsoil_gcm_dyn(:,:,:,1,:))
116        write(*,*) "Data for soiltemp downloaded!"
117
118        call get_var4("Waterdensity_soil",watersoil_density_dyn(:,:,:,1,:))
119        write(*,*) "Data for waterdensity_soil downloaded!"
120
121        call get_var3("Waterdensity_surface",watersurf_density_dyn(:,:,1,:))
122        write(*,*) "Data for waterdensity_surface downloaded!"
123    endif !soil_pem
124#endif
125else ! We use slopes
126    do islope = 1,nslope
127        write(num,'(i2.2)') islope
128        call get_var3("co2ice_slope"//num,co2_ice_slope_dyn(:,:,islope,:))
129    enddo
130    write(*,*) "Data for co2_ice downloaded!"
131
132    do islope = 1,nslope
133        write(num,'(i2.2)') islope
134        call get_var3("h2o_ice_s_slope"//num,h2o_ice_s_dyn(:,:,islope,:))
135    enddo
136    write(*,*) "Data for h2o_ice_s downloaded!"
137
138#ifndef CPP_STD
139    do islope = 1,nslope
140        write(num,'(i2.2)') islope
141        call get_var3("watercap_slope"//num,watercap(:,:,islope,:))
142    enddo
143    write(*,*) "Data for watercap downloaded!"
144#endif
145
146    do islope = 1,nslope
147        write(num,'(i2.2)') islope
148        call get_var3("tsurf_slope"//num,tsurf_gcm_dyn(:,:,islope,:))
149    enddo
150    write(*,*) "Data for tsurf downloaded!"
151
152#ifndef CPP_STD
153    if (soil_pem) then
154        do islope = 1,nslope
155            write(num,'(i2.2)') islope
156            call get_var4("soiltemp_slope"//num,tsoil_gcm_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,4)
182
183! Compute averages over the year for each point
184write(*,*) "Computing the average of tsurf..."
185tsurf_ave_dyn(:,:,:) = sum(tsurf_gcm_dyn(:,:,:,:),4)/timelen
186
187if (soil_pem) then
188    write(*,*) "Computing average of tsoil..."
189    tsoil_ave_dyn(:,:,:,:) = sum(tsoil_gcm_dyn(:,:,:,:,:),5)/timelen
190    write(*,*) "Computing average of waterdensity_surface..."
191    watersurf_density_dyn_ave(:,:,:) = 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-30
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_gcm(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_gcm,vmr_co2_gcm_phys)
218    call gr_dyn_fi(timelen,iim_input + 1,jjm_input + 1,ngrid,ps_GCM,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_ave_dyn(:,:,l,islope),tsoil_ave(:,l,islope))
228                do t = 1,timelen
229                    call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,tsoil_gcm_dyn(:,:,l,islope,t),tsoil_gcm(:,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_ave(:,:,islope),watersurf_density_ave(:,islope))
234        endif ! soil_pem
235        do t = 1,timelen
236            call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,tsurf_GCM_dyn(:,:,islope,t),tsurf_GCM(:,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_ave_dyn,tsurf_ave)
242#else
243    vmr_co2_gcm_phys(1,:) = vmr_co2_gcm(1,1,:)
244    ps_timeseries(1,:) = ps_GCM(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_ave(1,:,:) = tsoil_ave_dyn(1,1,:,:)
251        tsoil_gcm(1,:,:,:) = tsoil_gcm_dyn(1,1,:,:,:)
252        watersoil_density(1,:,:,:) = watersoil_density_dyn(1,1,:,:,:)
253        watersurf_density_ave(1,:) = watersurf_density_dyn_ave(1,1,:)
254    endif ! soil_pem
255    tsurf_GCM(1,:,:) = tsurf_GCM_dyn(1,1,:,:)
256    co2_ice_slope(1,:,:) = co2_ice_slope_dyn(1,1,:,:)
257    tsurf_ave(1,:) = tsurf_ave_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(filename,v)
285
286implicit none
287
288character(len = *), intent(in)  :: filename
289real, dimension(:), intent(out) :: v
290
291call error_msg(NF90_INQ_VARID(fID,filename,vID),"inq",filename)
292call error_msg(NF90_GET_VAR(fID,vID,v),"get",filename)
293
294END SUBROUTINE get_var1
295
296!=======================================================================
297
298SUBROUTINE get_var3(filename,v) ! on U grid
299
300implicit none
301
302character(len = *),     intent(in)  :: filename
303real, dimension(:,:,:), intent(out) :: v
304
305call error_msg(NF90_INQ_VARID(fID,filename,vID),"inq",filename)
306call error_msg(NF90_GET_VAR(fID,vID,v),"get",filename)
307
308END SUBROUTINE get_var3
309
310!=======================================================================
311
312SUBROUTINE get_var4(filename,v)
313
314implicit none
315
316character(len = *),       intent(in)  :: filename
317real, dimension(:,:,:,:), intent(out) :: v
318
319call error_msg(NF90_INQ_VARID(fID,filename,vID),"inq",filename)
320call error_msg(NF90_GET_VAR(fID,vID,v),"get",filename)
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.