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

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

MARS PEM:

  • Add a PEMETAT0 that read "startfi_pem.nc"
  • Add the soil in the model: soil temperature, thermal properties, ice table
  • Add a routine that compute CO2 + H2O adsorption
  • Minor corrections in PEM.F90

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