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

Last change on this file since 2842 was 2842, checked in by romain.vande, 3 years ago

Mars PEM:
PEM is now adapted to run with XIOS diurnal averages (when they will work properly)
RV
File size: 9.9 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
7      use netcdf, only: nf90_open,NF90_NOWRITE,nf90_noerr,nf90_strerror, &
8                        nf90_get_var, nf90_inq_varid, nf90_inq_dimid, &
9                        nf90_inquire_dimension,nf90_close
10      use comsoil_h, only: nsoilmx
11      USE soil_evolution_mod, ONLY: soil_pem
12
13      IMPLICIT NONE
14
15!=======================================================================
16!
17! Read initial confitions file
18!
19!=======================================================================
20
21  include "dimensions.h"
22
23!===============================================================================
24! Arguments:
25  CHARACTER(LEN=*), INTENT(IN) :: fichnom          !--- FILE NAME
26  INTEGER, INTENT(IN) :: timelen                   ! number of times stored in the file
27
28  INTEGER :: iim_input,jjm_input,nlayer,nslope
29  REAL, ALLOCATABLE ::  h2o_ice_s(:,:,:)                       ! h2o_ice_s of the concatenated file
30  REAL, ALLOCATABLE ::  co2_ice_s(:,:,:)                       ! co2_ice_s of the concatenated file
31
32  REAL, ALLOCATABLE ::  h2o_ice_s_slope(:,:,:,:)                       ! co2_ice_s of the concatenated file
33
34  REAL, INTENT(OUT) ::  min_h2o_ice_s(iim_input+1,jjm_input+1) ! Minimum of h2o_ice_s of the year
35  REAL, INTENT(OUT) ::  min_co2_ice_s(iim_input+1,jjm_input+1) ! Minimum of co2_ice_s of the year
36  REAL, INTENT(OUT) ::  min_co2_ice_slope(iim_input+1,jjm_input+1,nslope) ! Minimum of co2_ice slope of the year
37  REAL, INTENT(OUT) ::  min_h2o_ice_slope(iim_input+1,jjm_input+1,nslope) ! Minimum of co2_ice slope of the year
38  REAL, INTENT(OUT) ::  vmr_co2_gcm(iim_input+1,jjm_input+1,timelen)      !!!!vmr_co2_phys_gcm(iim_input+1,jjm_input+1,timelen)
39  REAL, INTENT(OUT) ::  q_h2o_GCM(iim_input+1,jjm_input+1,timelen)
40  REAL, INTENT(OUT) ::  q_co2_GCM(iim_input+1,jjm_input+1,timelen)
41  REAL, ALLOCATABLE ::  q1_co2_GCM(:,:,:)
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
51  REAL ::  TI_gcm(iim_input+1,jjm_input+1,nsoilmx,nslope,timelen) ! Thermal Inertia  of the concatenated file
52  REAL, INTENT(OUT) ::  TI_ave(iim_input+1,jjm_input+1,nsoilmx,nslope) ! Average Thermal Inertia  of the concatenated file
53  REAL, INTENT(OUT) ::  co2_ice_slope(iim_input+1,jjm_input+1,nslope,timelen) ! Minimum of co2_ice slope of the year
54!===============================================================================
55!   Local Variables
56  CHARACTER(LEN=256) :: msg, var, modname
57  INTEGER,PARAMETER :: length=100
58  INTEGER :: iq, fID, vID, idecal
59  INTEGER :: ierr
60  CHARACTER(len=12) :: start_file_type="earth" ! default start file type
61
62  REAL,ALLOCATABLE :: time(:) ! times stored in start
63  INTEGER :: indextime ! index of selected time
64
65  INTEGER :: edges(4),corner(4)
66  INTEGER :: i,j,t
67  real,save :: m_co2, m_noco2, A , B, mmean
68
69  INTEGER :: islope
70  CHARACTER*2 :: num
71
72!-----------------------------------------------------------------------
73  modname="pemetat0"
74
75      m_co2 = 44.01E-3  ! CO2 molecular mass (kg/mol)   
76      m_noco2 = 33.37E-3  ! Non condensible mol mass (kg/mol)   
77      A =(1/m_co2 - 1/m_noco2)
78      B=1/m_noco2
79
80      allocate(co2_ice_s(iim+1,jjm+1,timelen))
81      allocate(q1_co2_GCM(iim+1,jjm+1,timelen))
82      allocate(h2o_ice_s_slope(iim+1,jjm+1,nslope,timelen))
83      allocate(h2o_ice_s(iim+1,jjm+1,timelen))
84
85  print *, "Opening ", fichnom, "..."
86
87!  Open initial state NetCDF file
88  var=fichnom
89  CALL err(NF90_OPEN(var,NF90_NOWRITE,fID),"open",var)
90
91     print *, "Downloading data for h2oice ..."
92
93! Get h2o_ice_s of the concatenated file
94  CALL get_var3("h2o_ice_s"   ,h2o_ice_s)
95
96     print *, "Downloading data for h2oice done"
97     print *, "Downloading data for co2ice ..."
98
99  CALL get_var3("co2ice"   ,co2_ice_s)
100
101     print *, "Downloading data for co2ice done"
102     print *, "Downloading data for vmr co2..."
103
104  CALL get_var3("co2_cropped"   ,q_co2_GCM)
105
106     print *, "Downloading data for vmr co2 done"
107     print *, "Downloading data for vmr h20..."
108
109  CALL get_var3("h2o_cropped"   ,q_h2o_GCM)
110
111     print *, "Downloading data for vmr h2o done"
112     print *, "Downloading data for surface pressure ..."
113
114  CALL get_var3("ps"   ,ps_GCM)
115
116     print *, "Downloading data for surface pressure done"
117     print *, "nslope=", nslope
118     print *, "Downloading data for co2ice_slope ..."
119
120if(nslope.gt.1) then
121
122DO islope=1,nslope
123  write(num,fmt='(i2.2)') islope
124  call get_var3("co2ice_slope"//num,co2_ice_slope(:,:,islope,:))
125ENDDO
126
127     print *, "Downloading data for co2ice_slope done"
128     print *, "Downloading data for h2o_ice_s_slope ..."
129
130DO islope=1,nslope
131  write(num,fmt='(i2.2)') islope
132  call get_var3("h2o_ice_s_slope"//num,h2o_ice_s_slope(:,:,islope,:))
133ENDDO
134
135     print *, "Downloading data for h2o_ice_s_slope done"
136     print *, "Downloading data for tsurf_slope ..."
137
138DO islope=1,nslope
139  write(num,fmt='(i2.2)') islope
140  call get_var3("tsurf_slope"//num,tsurf_gcm(:,:,islope,:))
141ENDDO
142
143     print *, "Downloading data for tsurf_slope done"
144
145     if(soil_pem) then
146
147     print *, "Downloading data for tsoil_slope ..."
148
149DO islope=1,nslope
150  write(num,fmt='(i2.2)') islope
151  call get_var4("tsoil_slope"//num,tsoil_gcm(:,:,:,islope,:))
152ENDDO
153
154     print *, "Downloading data for tsoil_slope done"
155     print *, "Downloading data for inertiesoil_slope ..."
156
157DO islope=1,nslope
158  write(num,fmt='(i2.2)') islope
159  call get_var4("inertiesoil_slope"//num,TI_gcm(:,:,:,islope,:))
160ENDDO
161
162     print *, "Downloading data for inertiesoil_slope done"
163
164  endif !soil_pem
165
166  else !nslope=1 no slope, we copy all the values
167    co2_ice_slope(:,:,1,:)=co2_ice_s(:,:,:)
168    h2o_ice_s_slope(:,:,1,:)=h2o_ice_s(:,:,:)
169    call get_var3("tsurf",tsurf_gcm(:,:,1,:))
170    if(soil_pem) then
171      call get_var4("tsoil",tsoil_gcm(:,:,:,1,:))
172      call get_var4("inertiesoil",TI_gcm(:,:,:,1,:))
173    endif !soil_pem
174  endif !nslope=1
175
176! Compute the minimum over the year for each point
177  print *, "Computing the min of h2o_ice"
178  min_h2o_ice_s(:,:)=minval(h2o_ice_s,3)
179  print *, "Computing the min of co2_ice"
180  min_co2_ice_s(:,:)=minval(co2_ice_s,3)
181
182  print *, "Computing the min of h2o_ice_slope"
183  min_h2o_ice_slope(:,:,:)=minval(h2o_ice_s_slope,4)
184  print *, "Computing the min of co2_ice_slope"
185  min_co2_ice_slope(:,:,:)=minval(co2_ice_slope,4)
186
187!Compute averages
188
189    print *, "Computing average of tsurf"
190    tsurf_ave(:,:,:)=SUM(tsurf_gcm(:,:,:,:),4)/timelen
191
192  if(soil_pem) then
193    print *, "Computing average of tsoil"
194    tsoil_ave(:,:,:,:)=SUM(tsoil_gcm(:,:,:,:,:),5)/timelen
195    print *, "Computing average of TI"
196    TI_ave(:,:,:,:)=SUM(TI_gcm(:,:,:,:,:),5)/timelen
197  endif
198
199! By definition, a density is positive, we get rid of the negative values
200  DO i=1,iim+1
201    DO j = 1, jjm+1
202       if (min_co2_ice_s(i,j).LT.0) then
203          min_h2o_ice_s(i,j)  = 0.
204          min_co2_ice_s(i,j)  = 0.
205       endif
206       DO islope=1,nslope
207          if (min_co2_ice_slope(i,j,islope).LT.0) then
208            min_co2_ice_slope(i,j,islope)  = 0.
209          endif
210          if (min_h2o_ice_slope(i,j,islope).LT.0) then
211            min_h2o_ice_slope(i,j,islope)  = 0.
212          endif
213       ENDDO
214    ENDDO
215  ENDDO
216
217  DO i=1,iim+1
218    DO j = 1, jjm+1
219      DO t = 1, timelen
220         if (q_co2_GCM(i,j,t).LT.0) then
221              q_co2_GCM(i,j,t)=1E-10
222         elseif (q_co2_GCM(i,j,t).GT.1) then
223              q_co2_GCM(i,j,t)=1.
224         endif
225         if (q_h2o_GCM(i,j,t).LT.0) then
226              q_h2o_GCM(i,j,t)=1E-30
227         elseif (q_h2o_GCM(i,j,t).GT.1) then
228              q_h2o_GCM(i,j,t)=1.
229         endif
230         mmean=1/(A*q_co2_GCM(i,j,t) +B)
231         vmr_co2_gcm(i,j,t) = q_co2_GCM(i,j,t)*mmean/m_co2
232      ENDDO
233    ENDDO
234  ENDDO
235
236  CONTAINS
237
238SUBROUTINE check_dim(n1,n2,str1,str2)
239  INTEGER,          INTENT(IN) :: n1, n2
240  CHARACTER(LEN=*), INTENT(IN) :: str1, str2
241  CHARACTER(LEN=256) :: s1, s2
242  IF(n1/=n2) THEN
243    s1='value of '//TRIM(str1)//' ='
244    s2=' read in starting file differs from parametrized '//TRIM(str2)//' ='
245    WRITE(msg,'(10x,a,i4,2x,a,i4)')TRIM(s1),n1,TRIM(s2),n2
246    CALL ABORT_gcm(TRIM(modname),TRIM(msg),1)
247  END IF
248END SUBROUTINE check_dim
249
250
251SUBROUTINE get_var1(var,v)
252  CHARACTER(LEN=*), INTENT(IN)  :: var
253  REAL,             INTENT(OUT) :: v(:)
254  CALL err(NF90_INQ_VARID(fID,var,vID),"inq",var)
255  CALL err(NF90_GET_VAR(fID,vID,v),"get",var)
256END SUBROUTINE get_var1
257
258
259SUBROUTINE get_var3(var,v) ! on U grid
260  CHARACTER(LEN=*), INTENT(IN)  :: var
261  REAL,             INTENT(OUT) :: v(:,:,:)
262  CALL err(NF90_INQ_VARID(fID,var,vID),"inq",var)
263  CALL err(NF90_GET_VAR(fID,vID,v),"get",var)
264
265END SUBROUTINE get_var3
266
267SUBROUTINE get_var4(var,v)
268  CHARACTER(LEN=*), INTENT(IN)  :: var
269  REAL,             INTENT(OUT) :: v(:,:,:,:)
270  CALL err(NF90_INQ_VARID(fID,var,vID),"inq",var)
271  CALL err(NF90_GET_VAR(fID,vID,v),"get",var)
272END SUBROUTINE get_var4
273
274SUBROUTINE err(ierr,typ,nam)
275  INTEGER,          INTENT(IN) :: ierr   !--- NetCDF ERROR CODE
276  CHARACTER(LEN=*), INTENT(IN) :: typ    !--- TYPE OF OPERATION
277  CHARACTER(LEN=*), INTENT(IN) :: nam    !--- FIELD/FILE NAME
278  IF(ierr==NF90_NoERR) RETURN
279  SELECT CASE(typ)
280    CASE('inq');   msg="Field <"//TRIM(nam)//"> is missing"
281    CASE('get');   msg="Reading failed for <"//TRIM(nam)//">"
282    CASE('open');  msg="File opening failed for <"//TRIM(nam)//">"
283    CASE('close'); msg="File closing failed for <"//TRIM(nam)//">"
284  END SELECT
285  CALL ABORT_gcm(TRIM(modname),TRIM(msg),ierr)
286END SUBROUTINE err
287
288END SUBROUTINE read_data_gcm
Note: See TracBrowser for help on using the repository browser.