source: trunk/LMDZ.COMMON/libf/evolution/read_data_GCM_mod.F90 @ 3083

Last change on this file since 3083 was 3076, checked in by jbclement, 2 years ago

PEM:
Big cleaning/improvements of the PEM:

  • Conversion of "abort_pem.F" and "soil_settings_PEM.F" into Fortran 90;
  • Transformation of every PEM subroutines into module;
  • Rewriting of many subroutines with modern Fortran syntax;
  • Correction of a bug in "pem.F90" when calling 'recomp_tend_co2_slope'. The arguments were given in disorder and emissivity was missing;
  • Update of "launch_pem.sh" in deftank.

JBC
File size: 15.3 KB
Line 
1MODULE read_data_GCM_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(256) :: msg, var, modname ! for reading
10integer        :: fID, vID          ! for reading
11
12!=======================================================================
13contains
14!=======================================================================
15
16SUBROUTINE read_data_GCM(fichnom,timelen,iim_input,jjm_input,ngrid,nslope,vmr_co2_gcm_phys,ps_timeseries,           &
17                         min_co2_ice,min_h2o_ice,tsurf_ave,tsoil_ave,tsurf_gcm,tsoil_gcm,q_co2,q_h2o,co2_ice_slope, &
18                         watersurf_density_ave,watersoil_density)
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 GCM
28!
29! Authors: RV & LL
30!=======================================================================
31
32include "dimensions.h"
33
34!=======================================================================
35! Arguments:
36
37character(len = *), intent(in) :: fichnom                             !--- FILE NAME
38integer,            intent(in) :: timelen                             ! number of times stored in the file
39integer                        :: 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_gcm_phys ! Physics x Times  co2 volume mixing ratio retrieve from the gcm [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_ave             ! Average surface temperature of the concatenated file [K]
50real, dimension(ngrid,nsoilmx,nslope),         intent(out) :: tsoil_ave             ! Average soil temperature of the concatenated file [K]
51real, dimension(ngrid,nslope,timelen),         intent(out) :: tsurf_gcm             ! Surface temperature of the concatenated file, time series [K]
52real, dimension(ngrid,nsoilmx,nslope,timelen), intent(out) :: tsoil_gcm             ! Soil temperature of the concatenated file, time series [K]
53real, dimension(ngrid,nslope),                 intent(out) :: watersurf_density_ave ! 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
57character(12)                   :: start_file_type = "earth" ! default start file type
58real, dimension(:), allocatable :: time      ! times stored in start
59integer                         :: indextime ! index of selected time
60integer                         :: edges(4), corner(4)
61integer                         :: i, j, l, t          ! loop variables
62real                            ::  A , B, mmean       ! Molar Mass of co2 and no co2, A;B intermediate variables to compute the mean molar mass of the layer
63integer                         :: islope              ! loop for variables
64character(2)                    :: num                 ! for reading sloped variables
65real, dimension(iim_input + 1,jjm_input + 1,nslope,timelen)         :: h2o_ice_s_dyn         ! h2o ice per slope of the concatenated file [kg/m^2]
66real, dimension(iim_input + 1,jjm_input + 1,nslope,timelen)         :: watercap_slope
67real, dimension(iim_input + 1,jjm_input + 1,timelen)                :: vmr_co2_gcm           ! CO2 volume mixing ratio in the first layer  [mol/m^3]
68real, dimension(iim_input + 1,jjm_input + 1,timelen)                :: ps_GCM                ! Surface Pressure [Pa]
69real, dimension(iim_input + 1,jjm_input + 1,nslope)                 :: min_co2_ice_dyn
70real, dimension(iim_input + 1,jjm_input + 1,nslope)                 :: min_h2o_ice_dyn
71real, dimension(iim_input + 1,jjm_input + 1,nslope)                 :: tsurf_ave_dyn         ! Average surface temperature of the concatenated file [K]
72real, dimension(iim_input + 1,jjm_input + 1,nsoilmx,nslope)         :: tsoil_ave_dyn         ! Average soil temperature of the concatenated file [K]
73real, dimension(iim_input + 1,jjm_input + 1,nslope,timelen)         :: tsurf_gcm_dyn         ! Surface temperature of the concatenated file, time series [K]
74real, dimension(iim_input + 1,jjm_input + 1,nsoilmx,nslope,timelen) :: tsoil_gcm_dyn         ! Soil temperature of the concatenated file, time series [K]
75real, dimension(iim_input + 1,jjm_input + 1,timelen)                :: q_co2_dyn             ! CO2 mass mixing ratio in the first layer [kg/m^3]
76real, dimension(iim_input + 1,jjm_input + 1,timelen)                :: q_h2o_dyn             ! H2O mass mixing ratio in the first layer [kg/m^3]
77real, dimension(iim_input + 1,jjm_input + 1,nslope,timelen)         :: co2_ice_slope_dyn     ! co2 ice amount per  slope of the year [kg/m^2]
78real, dimension(iim_input + 1,jjm_input + 1,nslope,timelen)         :: watersurf_density_dyn ! Water density at the surface, time series [kg/m^3]
79real, dimension(iim_input + 1,jjm_input + 1,nsoilmx,nslope,timelen) :: watersoil_density_dyn ! Water density in the soil layer, time series [kg/m^3]
80real, dimension(ngrid,nslope,timelen)                               :: watersurf_density     ! Water density at the surface, time series [kg/m^3]
81
82!-----------------------------------------------------------------------
83modname="read_data_gcm"
84
85A = (1/m_co2 - 1/m_noco2)
86B = 1/m_noco2
87
88write(*,*) "Opening ", fichnom, "..."
89
90!  Open initial state NetCDF file
91var = fichnom
92call error_msg(NF90_OPEN(var,NF90_NOWRITE,fID),"open",var)
93
94write(*,*) "Downloading data for vmr co2..."
95
96call get_var3("co2_layer1"   ,q_co2_dyn)
97
98write(*,*) "Downloading data for vmr co2 done"
99write(*,*) "Downloading data for vmr h20..."
100
101call get_var3("h2o_layer1"   ,q_h2o_dyn)
102
103write(*,*) "Downloading data for vmr h2o done"
104write(*,*) "Downloading data for surface pressure ..."
105
106call get_var3("ps"   ,ps_GCM)
107
108write(*,*) "Downloading data for surface pressure done"
109write(*,*) "nslope=", nslope
110
111if (nslope > 1) then
112    write(*,*) "Downloading data for co2ice_slope ..."
113
114    do islope = 1,nslope
115        write(num,'(i2.2)') islope
116        call get_var3("co2ice_slope"//num,co2_ice_slope_dyn(:,:,islope,:))
117    enddo
118   
119    write(*,*) "Downloading data for co2ice_slope done"
120    write(*,*) "Downloading data for h2o_ice_s_slope ..."
121   
122    do islope = 1,nslope
123        write(num,'(i2.2)') islope
124        call get_var3("h2o_ice_s_slope"//num,h2o_ice_s_dyn(:,:,islope,:))
125    enddo
126
127    write(*,*) "Downloading data for h2o_ice_s_slope done"
128
129#ifndef CPP_STD
130    write(*,*) "Downloading data for watercap_slope ..."
131    do islope = 1,nslope
132        write(num,'(i2.2)') islope
133        call get_var3("watercap_slope"//num,watercap_slope(:,:,islope,:))
134!        watercap_slope(:,:,:,:) = 0.
135    enddo
136    write(*,*) "Downloading data for watercap_slope done"
137#endif
138
139    write(*,*) "Downloading data for tsurf_slope ..."
140
141    do islope=1,nslope
142        write(num,'(i2.2)') islope
143        call get_var3("tsurf_slope"//num,tsurf_gcm_dyn(:,:,islope,:))
144    enddo
145
146    write(*,*) "Downloading data for tsurf_slope done"
147
148#ifndef CPP_STD
149    if (soil_pem) then
150        write(*,*) "Downloading data for soiltemp_slope ..."
151
152        do islope = 1,nslope
153            write(num,'(i2.2)') islope
154            call get_var4("soiltemp_slope"//num,tsoil_gcm_dyn(:,:,:,islope,:))
155        enddo
156
157        write(*,*) "Downloading data for soiltemp_slope done"
158        write(*,*) "Downloading data for watersoil_density ..."
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
165        write(*,*) "Downloading data for  watersoil_density  done"
166        write(*,*) "Downloading data for  watersurf_density  ..."
167
168        do islope = 1,nslope
169            write(num,'(i2.2)') islope
170            call get_var3("Waterdensity_surface"//num,watersurf_density_dyn(:,:,islope,:))
171        enddo
172
173        write(*,*) "Downloading data for  watersurf_density  done"
174
175    endif !soil_pem
176#endif
177else !nslope = 1 no slope, we copy all the values
178    call get_var3("h2o_ice_s", h2o_ice_s_dyn(:,:,1,:))
179    call get_var3("co2ice", co2_ice_slope_dyn(:,:,1,:))
180    call get_var3("tsurf", tsurf_gcm_dyn(:,:,1,:))
181   
182#ifndef CPP_STD
183!    call get_var3("watercap", watercap_slope(:,:,1,:))
184    watercap_slope(:,:,1,:) = 0.
185    watersurf_density_dyn(:,:,:,:) = 0.
186    watersoil_density_dyn(:,:,:,:,:) = 0.
187#endif
188
189    if (soil_pem) call get_var4("soiltemp",tsoil_gcm_dyn(:,:,:,1,:))
190endif !nslope=1
191
192! Compute the minimum over the year for each point
193write(*,*) "Computing the min of h2o_ice_slope"
194min_h2o_ice_dyn(:,:,:) = minval(h2o_ice_s_dyn + watercap_slope,4)
195write(*,*) "Computing the min of co2_ice_slope"
196min_co2_ice_dyn(:,:,:) = minval(co2_ice_slope_dyn,4)
197
198!Compute averages
199write(*,*) "Computing average of tsurf"
200tsurf_ave_dyn(:,:,:) = sum(tsurf_gcm_dyn(:,:,:,:),4)/timelen
201
202#ifndef CPP_1D
203    do islope = 1,nslope
204        do t = 1,timelen
205            call gr_dyn_fi(1,iim_input+1,jjm_input+1,ngrid,watersurf_density_dyn(:,:,islope,t),watersurf_density(:,islope,t))
206        enddo
207    enddo
208#endif
209
210if (soil_pem) then
211    write(*,*) "Computing average of tsoil"
212    tsoil_ave_dyn(:,:,:,:) = sum(tsoil_gcm_dyn(:,:,:,:,:),5)/timelen
213    write(*,*) "Computing average of watersurf_density"
214    watersurf_density_ave(:,:) = sum(watersurf_density(:,:,:),3)/timelen
215endif
216
217! By definition, a density is positive, we get rid of the negative values
218where (min_co2_ice_dyn < 0.) min_co2_ice_dyn = 0.
219where (min_h2o_ice_dyn < 0.) min_h2o_ice_dyn = 0.
220
221do i = 1,iim + 1
222    do j = 1,jjm_input + 1
223        do t = 1, timelen
224            if (q_co2_dyn(i,j,t) < 0) then
225                q_co2_dyn(i,j,t) = 1.e-10
226            else if (q_co2_dyn(i,j,t) > 1) then
227                q_co2_dyn(i,j,t) = 1.
228            endif
229            if (q_h2o_dyn(i,j,t) < 0) then
230                q_h2o_dyn(i,j,t) = 1.e-30
231            else if (q_h2o_dyn(i,j,t) > 1) then
232                q_h2o_dyn(i,j,t) = 1.
233            endif
234            mmean = 1/(A*q_co2_dyn(i,j,t) + B)
235            vmr_co2_gcm(i,j,t) = q_co2_dyn(i,j,t)*mmean/m_co2
236        enddo
237    enddo
238enddo
239
240#ifndef CPP_1D
241    call gr_dyn_fi(timelen,iim_input + 1,jjm_input + 1,ngrid,vmr_co2_gcm,vmr_co2_gcm_phys)
242    call gr_dyn_fi(timelen,iim_input + 1,jjm_input + 1,ngrid,ps_GCM,ps_timeseries)
243    call gr_dyn_fi(timelen,iim_input + 1,jjm_input + 1,ngrid,q_co2_dyn,q_co2)
244    call gr_dyn_fi(timelen,iim_input + 1,jjm_input + 1,ngrid,q_h2o_dyn,q_h2o)
245
246    do islope = 1,nslope
247        call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,min_co2_ice_dyn(:,:,islope),min_co2_ice(:,islope))
248        call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,min_h2o_ice_dyn(:,:,islope),min_h2o_ice(:,islope))
249        if (soil_pem) then
250            do l = 1,nsoilmx
251                call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,tsoil_ave_dyn(:,:,l,islope),tsoil_ave(:,l,islope))
252                do t=1,timelen
253                    call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,tsoil_gcm_dyn(:,:,l,islope,t),tsoil_gcm(:,l,islope,t))
254                    call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,watersoil_density_dyn(:,:,l,islope,t),watersoil_density(:,l,islope,t))
255                enddo
256            enddo
257        endif !soil_pem
258        do t = 1,timelen
259            call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,tsurf_GCM_dyn(:,:,islope,t),tsurf_GCM(:,islope,t))
260            call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,co2_ice_slope_dyn(:,:,islope,t),co2_ice_slope(:,islope,t))
261        enddo
262    enddo
263
264    call gr_dyn_fi(nslope,iim_input + 1,jjm_input + 1,ngrid,tsurf_ave_dyn,tsurf_ave)
265#else
266    vmr_co2_gcm_phys(1,:) = vmr_co2_gcm(1,1,:)
267    ps_timeseries(1,:) = ps_GCM(1,1,:)
268    q_co2(1,:) = q_co2_dyn(1,1,:)
269    q_h2o(1,:) = q_h2o_dyn(1,1,:)
270    min_co2_ice(1,:) = min_co2_ice_dyn(1,1,:)
271    min_h2o_ice(1,:) = min_h2o_ice_dyn(1,1,:)
272    if (soil_pem) then
273        tsoil_ave(1,:,:) = tsoil_ave_dyn(1,1,:,:)
274        tsoil_gcm(1,:,:,:) = tsoil_gcm_dyn(1,1,:,:,:)
275        watersoil_density(1,:,:,:) = watersoil_density_dyn(1,1,:,:,:)
276    endif !soil_pem
277    tsurf_GCM(1,:,:) = tsurf_GCM_dyn(1,1,:,:)
278    co2_ice_slope(1,:,:) = co2_ice_slope_dyn(1,1,:,:)
279    tsurf_ave(1,:) = tsurf_ave_dyn(1,1,:)
280#endif
281
282END SUBROUTINE read_data_gcm
283
284SUBROUTINE check_dim(n1,n2,str1,str2)
285
286implicit none
287
288integer,            intent(in) :: n1, n2
289character(len = *), intent(in) :: str1, str2
290
291character(256) :: s1, s2
292
293if (n1 /= n2) then
294    s1 = 'value of '//trim(str1)//' ='
295    s2 = ' read in starting file differs from parametrized '//trim(str2)//' ='
296    write(msg,'(10x,a,i4,2x,a,i4)')trim(s1),n1,trim(s2),n2
297    call abort_gcm(trim(modname),trim(msg),1)
298endif
299
300END SUBROUTINE check_dim
301
302!=======================================================================
303
304SUBROUTINE get_var1(var,v)
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_var1
315
316!=======================================================================
317
318SUBROUTINE get_var3(var,v) ! on U grid
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_var3
329
330!=======================================================================
331
332SUBROUTINE get_var4(var,v)
333
334implicit none
335
336character(len = *),       intent(in)  :: var
337real, dimension(:,:,:,:), intent(out) :: v
338
339call error_msg(NF90_INQ_VARID(fID,var,vID),"inq",var)
340call error_msg(NF90_GET_VAR(fID,vID,v),"get",var)
341
342END SUBROUTINE get_var4
343
344!=======================================================================
345
346SUBROUTINE error_msg(ierr,typ,nam)
347
348implicit none
349
350integer,            intent(in) :: ierr !--- NetCDF ERROR CODE
351character(len = *), intent(in) :: typ  !--- TYPE OF OPERATION
352character(len = *), intent(in) :: nam  !--- FIELD/FILE NAME
353
354if (ierr == nf90_noerr) return
355select case(typ)
356    case('inq');   msg="Field <"//trim(nam)//"> is missing"
357    case('get');   msg="Reading failed for <"//trim(nam)//">"
358    case('open');  msg="File opening failed for <"//trim(nam)//">"
359    case('close'); msg="File closing failed for <"//trim(nam)//">"
360    case default
361        write(*,*) 'There is no message for this error.'
362        error stop
363end select
364call abort_gcm(trim(modname),trim(msg),ierr)
365
366END SUBROUTINE error_msg
367
368END MODULE read_data_GCM_mod
Note: See TracBrowser for help on using the repository browser.