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

Last change on this file since 3137 was 3130, checked in by jbclement, 2 years ago

PEM:
The perennial co2 ice is now taken into account with co2 frost (qsurf) to compute the tendency and to make the update + Rework of how co2 frost is converted to perennial co2 ice at the end of the PEM run + Correction of the value of 'threshold_co2_frost2perennial' to correspond to 10 m + Perennial co2 ice is now handled outside 'paleoclimate' in "phyetat0_mod.F90" of the Mars PCM + Some cleanings.

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