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

Last change on this file since 3149 was 3149, checked in by jbclement, 20 months ago

PEM:

  • Simplification of the algorithm managing the stopping criteria;
  • Complete rework of the ice management in the PEM (H2O & CO2);

    Subroutines to evolve the H2O and CO2 ice are now in the same module "evol_ice_mod.F90".
    Tendencies are computed from the variation of "ice + frost" between the 2 PCM runs.
    Evolving ice in the PEM is now called 'h2o_ice' or 'co2_ice' (not anymore in 'qsurf' and free of 'water_reservoir').
    Default value 'ini_h2o_bigreservoir' (= 10 m) initializes the H2O ice of the first PEM run where there is 'watercap'. For the next PEM runs, initialization is done with the value kept in "startpem.nc". CO2 ice is taken from 'perennial_co2ice' of the PCM (paleoclimate flag must be true).
    Simplification of the condition to compute the surface ice cover needed for the stopping criteria.
    Frost ('qsurf') is not evolved by the PEM and given back to the PCM.
    New default threshold value 'inf_h2oice_threshold' (= 2 m) to decide at the end of the PEM run if the H2O ice should be 'watercap' or not for the next PCM runs. If H2O ice cannot be 'watercap', then the remaining H2O ice is transferred to the frost ('qsurf').

  • Renaming of variables/subroutines for clarity;
  • Some cleanings throughout the code;
  • Small updates in files of the deftank.

JBC

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