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

Last change on this file since 2885 was 2885, checked in by romain.vande, 22 months ago

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