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

Last change on this file since 3050 was 3050, checked in by jbclement, 14 months ago

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