source: trunk/LMDZ.COMMON/libf/evolution/read_data_GCM.F90 @ 2909

Last change on this file since 2909 was 2897, checked in by romain.vande, 2 years ago

Mars PEM:
Deep cleaning of variables name and allocate.
All the "dyn to phys" grid change is done in subroutines and not in the main program.
File size: 14.2 KB
Line 
1!
2! $Id $
3!
4SUBROUTINE read_data_GCM(fichnom,timelen, iim_input,jjm_input,ngrid,nslope,vmr_co2_gcm_phys,ps_timeseries, &
5             min_co2_ice,min_h2o_ice,tsurf_ave,tsoil_ave,tsurf_gcm,tsoil_gcm,TI_ave,q_co2,q_h2o,co2_ice_slope, &
6             watersurf_density_ave,watersoil_density)
7
8      use netcdf, only: nf90_open,NF90_NOWRITE,nf90_noerr,nf90_strerror, &
9                        nf90_get_var, nf90_inq_varid, nf90_inq_dimid, &
10                        nf90_inquire_dimension,nf90_close
11      use comsoil_h, only: nsoilmx
12      USE comsoil_h_PEM, ONLY: soil_pem
13
14      IMPLICIT NONE
15
16!=======================================================================
17!
18! Purpose: Read initial confitions file from the GCM
19!
20! Authors: RV & LL
21!=======================================================================
22
23  include "dimensions.h"
24
25!===============================================================================
26! Arguments:
27
28  CHARACTER(LEN=*), INTENT(IN) :: fichnom          !--- FILE NAME
29  INTEGER, INTENT(IN) :: timelen                   ! number of times stored in the file
30  INTEGER :: iim_input,jjm_input,ngrid,nslope            ! number of points in the lat x lon dynamical grid, number of subgrid slopes
31! Ouputs
32  REAL, INTENT(OUT) ::  min_co2_ice(ngrid,nslope) ! Minimum of co2 ice  per slope of the year [kg/m^2]
33  REAL, INTENT(OUT) ::  min_h2o_ice(ngrid,nslope) ! Minimum of h2o ice per slope of the year [kg/m^2]
34  REAL, INTENT(OUT)  :: vmr_co2_gcm_phys(ngrid,timelen) ! Physics x Times  co2 volume mixing ratio retrieve from the gcm [m^3/m^3]
35  REAL, INTENT(OUT) ::  ps_timeseries(ngrid,timelen)! Surface Pressure [Pa]
36  REAL, INTENT(OUT) ::  q_co2(ngrid,timelen)        ! CO2 mass mixing ratio in the first layer [kg/m^3]
37  REAL, INTENT(OUT) ::  q_h2o(ngrid,timelen)        ! H2O mass mixing ratio in the first layer [kg/m^3]
38  REAL, INTENT(OUT) ::  co2_ice_slope(ngrid,nslope,timelen) ! co2 ice amount per  slope of the year [kg/m^2]
39!SOIL
40  REAL, INTENT(OUT) ::  tsurf_ave(ngrid,nslope)         ! Average surface temperature of the concatenated file [K]
41  REAL, INTENT(OUT) ::  tsoil_ave(ngrid,nsoilmx,nslope) ! Average soil temperature of the concatenated file [K]
42  REAL ,INTENT(OUT) ::  tsurf_gcm(ngrid,nslope,timelen)                  ! Surface temperature of the concatenated file, time series [K]
43  REAL , INTENT(OUT) ::  tsoil_gcm(ngrid,nsoilmx,nslope,timelen)         ! Soil temperature of the concatenated file, time series [K]
44  REAL , INTENT(OUT) ::  watersurf_density_ave(ngrid,nslope)             ! Water density at the surface [kg/m^3]
45  REAL , INTENT(OUT) ::  watersoil_density(ngrid,nsoilmx,nslope,timelen) ! Water density in the soil layer, time series [kg/m^3]
46  REAL, INTENT(OUT) ::  TI_ave(ngrid,nsoilmx,nslope)                     ! Average Thermal Inertia  of the concatenated file [SI]
47!===============================================================================
48!   Local Variables
49  CHARACTER(LEN=256) :: msg, var, modname               ! for reading
50  INTEGER :: iq, fID, vID, idecal                       ! for reading
51  INTEGER :: ierr                                       ! for reading
52  CHARACTER(len=12) :: start_file_type="earth" ! default start file type
53
54  REAL,ALLOCATABLE :: time(:) ! times stored in start
55  INTEGER :: indextime ! index of selected time
56
57  INTEGER :: edges(4),corner(4)
58  INTEGER :: i,j,l,t                                                     ! loop variables
59  real,save :: m_co2, m_noco2, A , B, mmean                            ! Molar Mass of co2 and no co2, A;B intermediate variables to compute the mean molar mass of the layer
60
61  INTEGER :: islope                                                    ! loop for variables
62  CHARACTER*2 :: num                                                   ! for reading sloped variables
63  REAL ::  h2o_ice_s_dyn(iim_input+1,jjm_input+1,nslope,timelen)       ! h2o ice per slope of the concatenated file [kg/m^2]
64  REAL ::  watercap_slope(iim_input+1,jjm_input+1,nslope,timelen)
65  REAL ::  vmr_co2_gcm(iim_input+1,jjm_input+1,timelen)                ! CO2 volume mixing ratio in the first layer  [mol/m^3]
66  REAL ::  ps_GCM(iim_input+1,jjm_input+1,timelen)                     ! Surface Pressure [Pa]
67  REAL ::  min_co2_ice_dyn(iim_input+1,jjm_input+1,nslope)
68  REAL ::  min_h2o_ice_dyn(iim_input+1,jjm_input+1,nslope)
69  REAL ::  tsurf_ave_dyn(iim_input+1,jjm_input+1,nslope)               ! Average surface temperature of the concatenated file [K]
70  REAL ::  tsoil_ave_dyn(iim_input+1,jjm_input+1,nsoilmx,nslope)       ! Average soil temperature of the concatenated file [K]
71  REAL ::  tsurf_gcm_dyn(iim_input+1,jjm_input+1,nslope,timelen)       ! Surface temperature of the concatenated file, time series [K]
72  REAL ::  tsoil_gcm_dyn(iim_input+1,jjm_input+1,nsoilmx,nslope,timelen)! Soil temperature of the concatenated file, time series [K]
73  REAL ::  TI_gcm(iim_input+1,jjm_input+1,nsoilmx,nslope,timelen)      ! Thermal Inertia  of the concatenated file, times series [SI]
74  REAL ::  TI_ave_dyn(iim_input+1,jjm_input+1,nsoilmx,nslope)          ! Average Thermal Inertia  of the concatenated file [SI]
75  REAL ::  q_co2_dyn(iim_input+1,jjm_input+1,timelen)                  ! CO2 mass mixing ratio in the first layer [kg/m^3]
76  REAL ::  q_h2o_dyn(iim_input+1,jjm_input+1,timelen)                  ! H2O mass mixing ratio in the first layer [kg/m^3]
77  REAL ::  co2_ice_slope_dyn(iim_input+1,jjm_input+1,nslope,timelen)  ! co2 ice amount per  slope of the year [kg/m^2]
78  REAL ::  watersurf_density_dyn(iim_input+1,jjm_input+1,nslope,timelen)! Water density at the surface, time series [kg/m^3]
79  REAL ::  watersurf_density(ngrid,nslope,timelen)                     ! Water density at the surface, time series [kg/m^3]
80  REAL ::  watersoil_density_dyn(iim_input+1,jjm_input+1,nsoilmx,nslope,timelen) ! Water density in the soil layer, time series [kg/m^3]
81
82!-----------------------------------------------------------------------
83  modname="read_data_gcm"
84
85      m_co2 = 44.01E-3  ! CO2 molecular mass (kg/mol)   
86      m_noco2 = 33.37E-3  ! Non condensible mol mass (kg/mol)   
87      A =(1/m_co2 - 1/m_noco2)
88      B=1/m_noco2
89
90  print *, "Opening ", fichnom, "..."
91
92!  Open initial state NetCDF file
93  var=fichnom
94  CALL err(NF90_OPEN(var,NF90_NOWRITE,fID),"open",var)
95
96     print *, "Downloading data for vmr co2..."
97
98  CALL get_var3("co2_cropped"   ,q_co2_dyn)
99
100     print *, "Downloading data for vmr co2 done"
101     print *, "Downloading data for vmr h20..."
102
103  CALL get_var3("h2o_cropped"   ,q_h2o_dyn)
104
105     print *, "Downloading data for vmr h2o done"
106     print *, "Downloading data for surface pressure ..."
107
108  CALL get_var3("ps"   ,ps_GCM)
109
110     print *, "Downloading data for surface pressure done"
111     print *, "nslope=", nslope
112     print *, "Downloading data for co2ice_slope ..."
113
114if(nslope.gt.1) then
115
116DO islope=1,nslope
117  write(num,fmt='(i2.2)') islope
118  call get_var3("co2ice_slope"//num,co2_ice_slope_dyn(:,:,islope,:))
119ENDDO
120
121     print *, "Downloading data for co2ice_slope done"
122     print *, "Downloading data for h2o_ice_s_slope ..."
123
124DO islope=1,nslope
125  write(num,fmt='(i2.2)') islope
126  call get_var3("h2o_ice_s_slope"//num,h2o_ice_s_dyn(:,:,islope,:))
127ENDDO
128
129     print *, "Downloading data for h2o_ice_s_slope done"
130
131     print *, "Downloading data for watercap_slope ..."
132DO islope=1,nslope
133       write(num,fmt='(i2.2)') islope
134       call get_var3("watercap_slope"//num,watercap_slope(:,:,islope,:))
135!        watercap_slope(:,:,:,:)= 0.
136ENDDO           
137     print *, "Downloading data for watercap_slope done"
138   
139 print *, "Downloading data for tsurf_slope ..."
140
141DO islope=1,nslope
142  write(num,fmt='(i2.2)') islope
143  call get_var3("tsurf_slope"//num,tsurf_gcm_dyn(:,:,islope,:))
144ENDDO
145
146     print *, "Downloading data for tsurf_slope done"
147
148     if(soil_pem) then
149
150     print *, "Downloading data for tsoil_slope ..."
151
152DO islope=1,nslope
153  write(num,fmt='(i2.2)') islope
154  call get_var4("tsoil_slope"//num,tsoil_gcm_dyn(:,:,:,islope,:))
155ENDDO
156
157     print *, "Downloading data for tsoil_slope done"
158     print *, "Downloading data for inertiesoil_slope ..."
159
160DO islope=1,nslope
161  write(num,fmt='(i2.2)') islope
162  call get_var4("inertiesoil_slope"//num,TI_gcm(:,:,:,islope,:))
163ENDDO
164
165     print *, "Downloading data for inertiesoil_slope done"
166
167     print *, "Downloading data for watersoil_density ..."
168
169DO islope=1,nslope
170  write(num,fmt='(i2.2)') islope
171  call get_var4("Waterdensity_soil_slope"//num,watersoil_density_dyn(:,:,:,islope,:))
172ENDDO
173
174     print *, "Downloading data for  watersoil_density  done"
175
176     print *, "Downloading data for  watersurf_density  ..."
177
178DO islope=1,nslope
179  write(num,fmt='(i2.2)') islope
180  call get_var3("Waterdensity_surface"//num,watersurf_density_dyn(:,:,islope,:))
181ENDDO
182
183     print *, "Downloading data for  watersurf_density  done"
184
185  endif !soil_pem
186
187  else !nslope=1 no slope, we copy all the values
188
189    CALL get_var3("h2o_ice_s", h2o_ice_s_dyn(:,:,1,:))
190    CALL get_var3("co2ice", co2_ice_slope_dyn(:,:,1,:))
191    call get_var3("tsurf", tsurf_gcm_dyn(:,:,1,:))
192#ifndef CPP_STD
193    call get_var3("watercap", watercap_slope(:,:,1,:))
194#endif
195
196    if(soil_pem) then
197      call get_var4("tsoil",tsoil_gcm_dyn(:,:,:,1,:))
198      call get_var4("inertiesoil",TI_gcm(:,:,:,1,:))
199    endif !soil_pem
200  endif !nslope=1
201
202! Compute the minimum over the year for each point
203  print *, "Computing the min of h2o_ice_slope"
204  min_h2o_ice_dyn(:,:,:)=minval(h2o_ice_s_dyn+watercap_slope,4)
205!  min_h2o_ice_dyn(:,:,:)=minval(h2o_ice_s_dyn,4)
206  print *, "Computing the min of co2_ice_slope"
207  min_co2_ice_dyn(:,:,:)=minval(co2_ice_slope_dyn,4)
208
209!Compute averages
210
211    print *, "Computing average of tsurf"
212    tsurf_ave_dyn(:,:,:)=SUM(tsurf_gcm_dyn(:,:,:,:),4)/timelen
213
214  DO islope = 1,nslope
215    DO t=1,timelen
216      CALL gr_dyn_fi(1,iim_input+1,jjm_input+1,ngrid,watersurf_density_dyn(:,:,islope,t),watersurf_density(:,islope,t))
217    ENDDO
218  ENDDO
219
220  if(soil_pem) then
221    print *, "Computing average of tsoil"
222    tsoil_ave_dyn(:,:,:,:)=SUM(tsoil_gcm_dyn(:,:,:,:,:),5)/timelen
223    print *, "Computing average of TI"
224    TI_ave_dyn(:,:,:,:)=SUM(TI_gcm(:,:,:,:,:),5)/timelen
225    print *, "Computing average of watersurf_density"
226    watersurf_density_ave(:,:) = SUM(watersurf_density(:,:,:),3)/timelen
227  endif
228
229! By definition, a density is positive, we get rid of the negative values
230  DO i=1,iim+1
231    DO j = 1, jjm+1
232       DO islope=1,nslope
233          if (min_co2_ice_dyn(i,j,islope).LT.0) then
234            min_co2_ice_dyn(i,j,islope)  = 0.
235          endif
236          if (min_h2o_ice_dyn(i,j,islope).LT.0) then
237            min_h2o_ice_dyn(i,j,islope)  = 0.
238          endif
239       ENDDO
240    ENDDO
241  ENDDO
242
243  DO i=1,iim+1
244    DO j = 1, jjm+1
245      DO t = 1, timelen
246         if (q_co2_dyn(i,j,t).LT.0) then
247              q_co2_dyn(i,j,t)=1E-10
248         elseif (q_co2_dyn(i,j,t).GT.1) then
249              q_co2_dyn(i,j,t)=1.
250         endif
251         if (q_h2o_dyn(i,j,t).LT.0) then
252              q_h2o_dyn(i,j,t)=1E-30
253         elseif (q_h2o_dyn(i,j,t).GT.1) then
254              q_h2o_dyn(i,j,t)=1.
255         endif
256         mmean=1/(A*q_co2_dyn(i,j,t) +B)
257         vmr_co2_gcm(i,j,t) = q_co2_dyn(i,j,t)*mmean/m_co2
258      ENDDO
259    ENDDO
260  ENDDO
261
262     CALL gr_dyn_fi(timelen,iim_input+1,jjm_input+1,ngrid,vmr_co2_gcm,vmr_co2_gcm_phys)
263     call gr_dyn_fi(timelen,iim_input+1,jjm_input+1,ngrid,ps_GCM,ps_timeseries)
264     CALL gr_dyn_fi(timelen,iim_input+1,jjm_input+1,ngrid,q_co2_dyn,q_co2)
265     CALL gr_dyn_fi(timelen,iim_input+1,jjm_input+1,ngrid,q_h2o_dyn,q_h2o)
266
267     DO islope = 1,nslope
268       CALL gr_dyn_fi(1,iim_input+1,jjm_input+1,ngrid,min_co2_ice_dyn(:,:,islope),min_co2_ice(:,islope))
269       CALL gr_dyn_fi(1,iim_input+1,jjm_input+1,ngrid,min_h2o_ice_dyn(:,:,islope),min_h2o_ice(:,islope))
270       if(soil_pem) then
271         CALL gr_dyn_fi(nsoilmx,iim_input+1,jjm_input+1,ngrid,TI_ave_dyn(:,:,:,islope),TI_ave(:,:,islope))
272       DO l=1,nsoilmx
273         CALL gr_dyn_fi(1,iim_input+1,jjm_input+1,ngrid,tsoil_ave_dyn(:,:,l,islope),tsoil_ave(:,l,islope))
274         DO t=1,timelen
275           CALL gr_dyn_fi(1,iim_input+1,jjm_input+1,ngrid,tsoil_gcm_dyn(:,:,l,islope,t),tsoil_gcm(:,l,islope,t))
276           CALL gr_dyn_fi(1,iim_input+1,jjm_input+1,ngrid,watersoil_density_dyn(:,:,l,islope,t),watersoil_density(:,l,islope,t))
277         ENDDO
278       ENDDO
279       endif !soil_pem
280       DO t=1,timelen
281         CALL gr_dyn_fi(1,iim_input+1,jjm_input+1,ngrid,tsurf_GCM_dyn(:,:,islope,t),tsurf_GCM(:,islope,t))
282         CALL gr_dyn_fi(1,iim_input+1,jjm_input+1,ngrid,co2_ice_slope_dyn(:,:,islope,t),co2_ice_slope(:,islope,t))
283       ENDDO
284     ENDDO
285
286     CALL gr_dyn_fi(nslope,iim_input+1,jjm_input+1,ngrid,tsurf_ave_dyn,tsurf_ave)
287
288  CONTAINS
289
290SUBROUTINE check_dim(n1,n2,str1,str2)
291  INTEGER,          INTENT(IN) :: n1, n2
292  CHARACTER(LEN=*), INTENT(IN) :: str1, str2
293  CHARACTER(LEN=256) :: s1, s2
294  IF(n1/=n2) THEN
295    s1='value of '//TRIM(str1)//' ='
296    s2=' read in starting file differs from parametrized '//TRIM(str2)//' ='
297    WRITE(msg,'(10x,a,i4,2x,a,i4)')TRIM(s1),n1,TRIM(s2),n2
298    CALL ABORT_gcm(TRIM(modname),TRIM(msg),1)
299  END IF
300END SUBROUTINE check_dim
301
302
303SUBROUTINE get_var1(var,v)
304  CHARACTER(LEN=*), INTENT(IN)  :: var
305  REAL,             INTENT(OUT) :: v(:)
306  CALL err(NF90_INQ_VARID(fID,var,vID),"inq",var)
307  CALL err(NF90_GET_VAR(fID,vID,v),"get",var)
308END SUBROUTINE get_var1
309
310
311SUBROUTINE get_var3(var,v) ! on U grid
312  CHARACTER(LEN=*), INTENT(IN)  :: var
313  REAL,             INTENT(OUT) :: v(:,:,:)
314  CALL err(NF90_INQ_VARID(fID,var,vID),"inq",var)
315  CALL err(NF90_GET_VAR(fID,vID,v),"get",var)
316
317END SUBROUTINE get_var3
318
319SUBROUTINE get_var4(var,v)
320  CHARACTER(LEN=*), INTENT(IN)  :: var
321  REAL,             INTENT(OUT) :: v(:,:,:,:)
322  CALL err(NF90_INQ_VARID(fID,var,vID),"inq",var)
323  CALL err(NF90_GET_VAR(fID,vID,v),"get",var)
324END SUBROUTINE get_var4
325
326SUBROUTINE err(ierr,typ,nam)
327  INTEGER,          INTENT(IN) :: ierr   !--- NetCDF ERROR CODE
328  CHARACTER(LEN=*), INTENT(IN) :: typ    !--- TYPE OF OPERATION
329  CHARACTER(LEN=*), INTENT(IN) :: nam    !--- FIELD/FILE NAME
330  IF(ierr==NF90_NoERR) RETURN
331  SELECT CASE(typ)
332    CASE('inq');   msg="Field <"//TRIM(nam)//"> is missing"
333    CASE('get');   msg="Reading failed for <"//TRIM(nam)//">"
334    CASE('open');  msg="File opening failed for <"//TRIM(nam)//">"
335    CASE('close'); msg="File closing failed for <"//TRIM(nam)//">"
336  END SELECT
337  CALL ABORT_gcm(TRIM(modname),TRIM(msg),ierr)
338END SUBROUTINE err
339
340END SUBROUTINE read_data_gcm
Note: See TracBrowser for help on using the repository browser.