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

Last change on this file since 2849 was 2849, checked in by llange, 2 years ago

PEM
Update the codes for subsurface evolution to run with XIOS
1) Water density at the surface and in the soil is now read in the XIOS file
2) Reshape routine created as XIOS output have one element less in the longitude grid
3) The ice table is now computed using these water densities
4) Minor fixs in the main, adsorption module, and tendencies evolutions

LL
File size: 10.8 KB
Line 
1!
2! $Id $
3!
4SUBROUTINE read_data_GCM(fichnom,min_h2o_ice_s,min_co2_ice_s,iim_input,jjm_input,nlayer,vmr_co2_gcm,ps_GCM,timelen, &
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 soil_evolution_mod, ONLY: soil_pem
13
14      IMPLICIT NONE
15
16!=======================================================================
17!
18! Read initial confitions file
19!
20!=======================================================================
21
22  include "dimensions.h"
23
24!===============================================================================
25! Arguments:
26  CHARACTER(LEN=*), INTENT(IN) :: fichnom          !--- FILE NAME
27  INTEGER, INTENT(IN) :: timelen                   ! number of times stored in the file
28
29  INTEGER :: iim_input,jjm_input,nlayer,nslope
30  REAL, ALLOCATABLE ::  h2o_ice_s(:,:,:)                       ! h2o_ice_s of the concatenated file
31  REAL, ALLOCATABLE ::  co2_ice_s(:,:,:)                       ! co2_ice_s of the concatenated file
32
33  REAL, ALLOCATABLE ::  h2o_ice_s_slope(:,:,:,:)                       ! co2_ice_s of the concatenated file
34
35  REAL, INTENT(OUT) ::  min_h2o_ice_s(iim_input+1,jjm_input+1) ! Minimum of h2o_ice_s of the year
36  REAL, INTENT(OUT) ::  min_co2_ice_s(iim_input+1,jjm_input+1) ! Minimum of co2_ice_s of the year
37  REAL, INTENT(OUT) ::  min_co2_ice_slope(iim_input+1,jjm_input+1,nslope) ! Minimum of co2_ice slope of the year
38  REAL, INTENT(OUT) ::  min_h2o_ice_slope(iim_input+1,jjm_input+1,nslope) ! Minimum of co2_ice slope of the year
39  REAL, INTENT(OUT) ::  vmr_co2_gcm(iim_input+1,jjm_input+1,timelen)      !!!!vmr_co2_phys_gcm(iim_input+1,jjm_input+1,timelen)
40  REAL, INTENT(OUT) ::  q_h2o_GCM(iim_input+1,jjm_input+1,timelen)
41  REAL, INTENT(OUT) ::  q_co2_GCM(iim_input+1,jjm_input+1,timelen)
42  REAL,  INTENT(OUT) ::  ps_GCM(iim_input+1,jjm_input+1,timelen)
43
44!SOIL
45  REAL, INTENT(OUT) ::  tsurf_ave(iim_input+1,jjm_input+1,nslope) ! Average surface temperature of the concatenated file
46  REAL, INTENT(OUT) ::  tsoil_ave(iim_input+1,jjm_input+1,nsoilmx,nslope) ! Average soil temperature of the concatenated file
47
48  REAL ,INTENT(OUT) ::  tsurf_gcm(iim_input+1,jjm_input+1,nslope,timelen) ! Surface temperature of the concatenated file
49  REAL , INTENT(OUT) ::  tsoil_gcm(iim_input+1,jjm_input+1,nsoilmx,nslope,timelen) ! Soil temperature of the concatenated file
50  REAL , INTENT(OUT) ::  watersurf_density(iim_input+1,jjm_input+1,nslope,timelen) ! Soil temperature of the concatenated file
51  REAL , INTENT(OUT) ::  watersoil_density(iim_input+1,jjm_input+1,nsoilmx,nslope,timelen) ! Soil temperature of the concatenated file
52
53  REAL ::  TI_gcm(iim_input+1,jjm_input+1,nsoilmx,nslope,timelen) ! Thermal Inertia  of the concatenated file
54  REAL, INTENT(OUT) ::  TI_ave(iim_input+1,jjm_input+1,nsoilmx,nslope) ! Average Thermal Inertia  of the concatenated file
55  REAL, INTENT(OUT) ::  co2_ice_slope(iim_input+1,jjm_input+1,nslope,timelen) ! Minimum of co2_ice slope of the year
56!===============================================================================
57!   Local Variables
58  CHARACTER(LEN=256) :: msg, var, modname
59  INTEGER,PARAMETER :: length=100
60  INTEGER :: iq, fID, vID, idecal
61  INTEGER :: ierr
62  CHARACTER(len=12) :: start_file_type="earth" ! default start file type
63
64  REAL,ALLOCATABLE :: time(:) ! times stored in start
65  INTEGER :: indextime ! index of selected time
66
67  INTEGER :: edges(4),corner(4)
68  INTEGER :: i,j,t
69  real,save :: m_co2, m_noco2, A , B, mmean
70
71  INTEGER :: islope
72  CHARACTER*2 :: num
73
74!-----------------------------------------------------------------------
75  modname="pemetat0"
76
77      m_co2 = 44.01E-3  ! CO2 molecular mass (kg/mol)   
78      m_noco2 = 33.37E-3  ! Non condensible mol mass (kg/mol)   
79      A =(1/m_co2 - 1/m_noco2)
80      B=1/m_noco2
81
82      allocate(co2_ice_s(iim+1,jjm+1,timelen))
83      allocate(h2o_ice_s_slope(iim+1,jjm+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 tsurf_slope ..."
138
139DO islope=1,nslope
140  write(num,fmt='(i2.2)') islope
141  call get_var3("tsurf_slope"//num,tsurf_gcm(:,:,islope,:))
142ENDDO
143
144     print *, "Downloading data for tsurf_slope done"
145
146     if(soil_pem) then
147
148     print *, "Downloading data for tsoil_slope ..."
149
150DO islope=1,nslope
151  write(num,fmt='(i2.2)') islope
152  call get_var4("tsoil_slope"//num,tsoil_gcm(:,:,:,islope,:))
153ENDDO
154
155     print *, "Downloading data for tsoil_slope done"
156     print *, "Downloading data for inertiesoil_slope ..."
157
158DO islope=1,nslope
159  write(num,fmt='(i2.2)') islope
160  call get_var4("inertiesoil_slope"//num,TI_gcm(:,:,:,islope,:))
161ENDDO
162
163     print *, "Downloading data for inertiesoil_slope done"
164
165     print *, "Downloading data for watersoil_density ..."
166
167DO islope=1,nslope
168  write(num,fmt='(i2.2)') islope
169  call get_var4("Waterdensity_soil_slope"//num,watersoil_density(:,:,:,islope,:))
170ENDDO
171
172     print *, "Downloading data for  watersoil_density  done"
173
174     print *, "Downloading data for  watersurf_density  ..."
175
176DO islope=1,nslope
177  write(num,fmt='(i2.2)') islope
178  call get_var3("Waterdensity_surface"//num,watersurf_density(:,:,islope,:))
179ENDDO
180
181     print *, "Downloading data for  watersurf_density  done"
182
183  endif !soil_pem
184
185  else !nslope=1 no slope, we copy all the values
186    co2_ice_slope(:,:,1,:)=co2_ice_s(:,:,:)
187    h2o_ice_s_slope(:,:,1,:)=h2o_ice_s(:,:,:)
188    call get_var3("tsurf",tsurf_gcm(:,:,1,:))
189    if(soil_pem) then
190      call get_var4("tsoil",tsoil_gcm(:,:,:,1,:))
191      call get_var4("inertiesoil",TI_gcm(:,:,:,1,:))
192    endif !soil_pem
193  endif !nslope=1
194
195! Compute the minimum over the year for each point
196  print *, "Computing the min of h2o_ice"
197  min_h2o_ice_s(:,:)=minval(h2o_ice_s,3)
198  print *, "Computing the min of co2_ice"
199  min_co2_ice_s(:,:)=minval(co2_ice_s,3)
200
201  print *, "Computing the min of h2o_ice_slope"
202  min_h2o_ice_slope(:,:,:)=minval(h2o_ice_s_slope,4)
203  print *, "Computing the min of co2_ice_slope"
204  min_co2_ice_slope(:,:,:)=minval(co2_ice_slope,4)
205
206!Compute averages
207
208    print *, "Computing average of tsurf"
209    tsurf_ave(:,:,:)=SUM(tsurf_gcm(:,:,:,:),4)/timelen
210
211  if(soil_pem) then
212    print *, "Computing average of tsoil"
213    tsoil_ave(:,:,:,:)=SUM(tsoil_gcm(:,:,:,:,:),5)/timelen
214    print *, "Computing average of TI"
215    TI_ave(:,:,:,:)=SUM(TI_gcm(:,:,:,:,:),5)/timelen
216  endif
217
218! By definition, a density is positive, we get rid of the negative values
219  DO i=1,iim+1
220    DO j = 1, jjm+1
221       if (min_co2_ice_s(i,j).LT.0) then
222          min_h2o_ice_s(i,j)  = 0.
223          min_co2_ice_s(i,j)  = 0.
224       endif
225       DO islope=1,nslope
226          if (min_co2_ice_slope(i,j,islope).LT.0) then
227            min_co2_ice_slope(i,j,islope)  = 0.
228          endif
229          if (min_h2o_ice_slope(i,j,islope).LT.0) then
230            min_h2o_ice_slope(i,j,islope)  = 0.
231          endif
232       ENDDO
233    ENDDO
234  ENDDO
235
236  DO i=1,iim+1
237    DO j = 1, jjm+1
238      DO t = 1, timelen
239         if (q_co2_GCM(i,j,t).LT.0) then
240              q_co2_GCM(i,j,t)=1E-10
241         elseif (q_co2_GCM(i,j,t).GT.1) then
242              q_co2_GCM(i,j,t)=1.
243         endif
244         if (q_h2o_GCM(i,j,t).LT.0) then
245              q_h2o_GCM(i,j,t)=1E-30
246         elseif (q_h2o_GCM(i,j,t).GT.1) then
247              q_h2o_GCM(i,j,t)=1.
248         endif
249         mmean=1/(A*q_co2_GCM(i,j,t) +B)
250         vmr_co2_gcm(i,j,t) = q_co2_GCM(i,j,t)*mmean/m_co2
251      ENDDO
252    ENDDO
253  ENDDO
254
255
256      deallocate(co2_ice_s)
257      deallocate(h2o_ice_s_slope)
258      deallocate(h2o_ice_s)
259
260
261  CONTAINS
262
263SUBROUTINE check_dim(n1,n2,str1,str2)
264  INTEGER,          INTENT(IN) :: n1, n2
265  CHARACTER(LEN=*), INTENT(IN) :: str1, str2
266  CHARACTER(LEN=256) :: s1, s2
267  IF(n1/=n2) THEN
268    s1='value of '//TRIM(str1)//' ='
269    s2=' read in starting file differs from parametrized '//TRIM(str2)//' ='
270    WRITE(msg,'(10x,a,i4,2x,a,i4)')TRIM(s1),n1,TRIM(s2),n2
271    CALL ABORT_gcm(TRIM(modname),TRIM(msg),1)
272  END IF
273END SUBROUTINE check_dim
274
275
276SUBROUTINE get_var1(var,v)
277  CHARACTER(LEN=*), INTENT(IN)  :: var
278  REAL,             INTENT(OUT) :: v(:)
279  CALL err(NF90_INQ_VARID(fID,var,vID),"inq",var)
280  CALL err(NF90_GET_VAR(fID,vID,v),"get",var)
281END SUBROUTINE get_var1
282
283
284SUBROUTINE get_var3(var,v) ! on U grid
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)
289
290END SUBROUTINE get_var3
291
292SUBROUTINE get_var4(var,v)
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)
297END SUBROUTINE get_var4
298
299SUBROUTINE err(ierr,typ,nam)
300  INTEGER,          INTENT(IN) :: ierr   !--- NetCDF ERROR CODE
301  CHARACTER(LEN=*), INTENT(IN) :: typ    !--- TYPE OF OPERATION
302  CHARACTER(LEN=*), INTENT(IN) :: nam    !--- FIELD/FILE NAME
303  IF(ierr==NF90_NoERR) RETURN
304  SELECT CASE(typ)
305    CASE('inq');   msg="Field <"//TRIM(nam)//"> is missing"
306    CASE('get');   msg="Reading failed for <"//TRIM(nam)//">"
307    CASE('open');  msg="File opening failed for <"//TRIM(nam)//">"
308    CASE('close'); msg="File closing failed for <"//TRIM(nam)//">"
309  END SELECT
310  CALL ABORT_gcm(TRIM(modname),TRIM(msg),ierr)
311END SUBROUTINE err
312
313END SUBROUTINE read_data_gcm
Note: See TracBrowser for help on using the repository browser.