Documentation/Forcings: mk_grid.f90

File mk_grid.f90, 9.2 KB (added by maignan, 12 years ago)
Line 
1PROGRAM ANA_FORCING
2
3IMPLICIT NONE
4INCLUDE "/usr/local/install/netcdf/include/netcdf.inc" !@FM sur obelix
5                                         
6INTEGER, PARAMETER                    :: IGRILLE= 8602    ! Numbre de points de grille
7                                                 !(nb de ligne du fichier coord_xxx correspondant)
8INTEGER, PARAMETER                    :: NBLON=143
9INTEGER, PARAMETER                    :: NBLAT=134
10REAL, DIMENSION(NBLON)                :: LON
11REAL, DIMENSION(NBLAT)                :: LAT
12REAL(kind=4), DIMENSION(IGRILLE)      :: indice_ligne,indice_col
13integer,DIMENSION(NBLON,NBLAT)        :: mask
14REAL,DIMENSION(NBLON,NBLAT)           :: mask0
15INTEGER                               :: JPOINT   ! points à extraire (fichier coord_xxx)
16INTEGER                               :: JPARA    ! varaible à extraire
17INTEGER, DIMENSION (IGRILLE)          :: I,IMAILLE
18REAL, DIMENSION(IGRILLE)              :: XLON, XLAT
19integer                               :: nlongi,nlati            !@FM
20integer                               :: day_ref,l,c
21integer                               :: month
22INTEGER                               :: II,dataid 
23INTEGER                               :: IANEXT   ! annee  à extraire
24CHARACTER(len=*),PARAMETER            :: pathout='/home/orchidee01/zhao/carboFrance/Safran/'
25CHARACTER(len=200)                    :: wname
26CHARACTER(len=200)                    :: str,str2 
27REAL, DIMENSION(IGRILLE,365*24)       :: Tair,Psurf,Qair,Rainf,Snowf,Wind,Swdown,Lwdown 
28integer,       dimension(8)           :: nuvar      !@F
29integer,dimension(4)                  :: corner,edge
30integer                               :: dims(4),ndim
31real(kind=4), parameter               :: missing_value= 1.E20
32real(kind=4), parameter               :: Pconst= 101400.        ! Surface Pressure (Pa)
33integer                               :: mois,jour,nval
34integer                               :: id, day, hour, step,y,nbyear
35integer                               :: hh 
36integer                               :: ncid,nxi,nyi,nti,nri,nlon,nlat,nmask,nindligne,nindcol,nbi,ndlon,ndlat,nuvar_time,nuvar_cell    !@FM
37integer,DIMENSION(IGRILLE)            :: ncell,cell
38integer                               :: ntime1,ntime2 
39integer                               :: status,ier
40real(kind=4), dimension(365*24)       :: timestp,time
41real(kind=4), parameter               :: secs_step=3600.
42CHARACTER(len=50), DIMENSION(8)       :: nomvar,nomunit,nomvar_long  !@FM
43
44DATA nomvar /'Tair','PSurf','Qair','Wind','Rainf','Snowf','SWdown','LWdown'/
45DATA nomvar_long /'Near surface air temperature','Surface pressure',& 
46               'near surface specific humidity','Near surface wind speed at 10m',&
47               'Rainfall rate','Snowfall rate',&
48               'Surface incident shortwave radiation',&
49               'Surface incident longwave radiation'/
50DATA nomunit /'K','Pa','kg/kg','m/s','kg/m^2/s','kg/m^2/s','W/m^2','W/m^2'/
51!------------------------------------------------------------------------------------------
52! Read in mask:
53     !status=nf_open('/home/orchidee01/viovy/cerfacs/arp_A1B_swdown_1950_2099.nc', 0, ncid)
54     !See mk_mask_new.R ???=> mask_Safran.nc)
55      status=nf_open('/home/orchidee01/zhao/carboFrance/Safran/mask_Safran.nc', 0, ncid)
56      corner(1)=1
57      edge(1)=NBLON
58         status=nf_inq_varid(ncid,"longitude",dataid)
59         status=nf_get_vara_real(ncid, dataid, corner, edge, LON)
60      corner(1)=1
61      edge(1)=NBLAT
62         status=nf_inq_varid(ncid,"latitude",dataid)
63         status=nf_get_vara_real(ncid, dataid, corner, edge, LAT)
64
65      corner(1)=1
66      edge(1)=IGRILLE
67         status=nf_inq_varid(ncid,"cell",dataid)
68          status=nf_get_vara_int(ncid, dataid, corner, edge, cell)
69         !print*,cell(1:10)
70         !print*,cell(8593:8602)
71         !pause
72
73      corner(1)=1
74      edge(1)=IGRILLE
75         status=nf_inq_varid(ncid,"indice_ligne",dataid)
76         status=nf_get_vara_real(ncid, dataid, corner, edge, indice_ligne)
77         status=nf_inq_varid(ncid,"indice_col",dataid)
78         status=nf_get_vara_real(ncid, dataid, corner, edge, indice_col)
79
80        !print*,"ligne:"
81        !print*,indice_ligne(1:10)
82        !print*,"col:"
83        !print*,indice_col(1:10)
84
85
86      corner(1)=1
87      corner(2)=1
88      edge(1)=NBLON
89      edge(2)=NBLAT
90          status=nf_inq_varid(ncid,"mask",dataid)
91          status=nf_get_vara_real(ncid, dataid, corner, edge, mask0)
92          mask(:,:)=int(mask0(:,:))
93          print*,(mask(:,11:12))
94        !open(20,file="test_mask.txt")
95        !write(20,*) mask
96        !close(20)
97
98! LECTURE DES NUMEROS DE MAILLE A EXTRAIRE
99OPEN (UNIT=50, FILE='/home/orchidee01/viovy/safran/coord_france', FORM='FORMATTED', STATUS='OLD')
100      DO JPOINT=1,IGRILLE
101      READ(50,*) I(JPOINT), IMAILLE(JPOINT), xlon(JPOINT), xlat(JPOINT)
102      ENDDO
103
104
105!Do IANEXT=1994,2007
106 Do IANEXT=2001,2007
107 ! Process Safran data year by year
108
109   !IANEXT=1994
110!  Extract the variables of the selected grid from Safran dataset
111     call SAFRAN_FORCING(IANEXT,1,Tair)
112    !print*,(Tair(IGRILLE,ii),ii=1,10)
113     call SAFRAN_FORCING(IANEXT,2,Qair)
114     call SAFRAN_FORCING(IANEXT,3,Wind)
115     call SAFRAN_FORCING(IANEXT,4,Rainf)
116     call SAFRAN_FORCING(IANEXT,5,Snowf)
117     call SAFRAN_FORCING(IANEXT,6,SWdown)
118     call SAFRAN_FORCING(IANEXT,7,LWdown)
119! No Surface pressure in Safran dataset, so set it as a constant
120     PSurf(:,:)=Pconst
121
122     WRITE(WNAME,1001) PATHOUT,IANEXT
123     print*, (wname) 
124
125! Save in NetCDF format as required for running orchidee
126 ! Set up the netCDF file:
127
128   ! Define dimensions:
129     ier=NF_CREATE(wname, NF_SHARE, ncid)
130     ier=NF_DEF_DIM(ncid,'longitude',NBLON,nxi)
131     ier=NF_DEF_DIM(ncid,'latitude',NBLAT,nyi)
132     ier=NF_DEF_DIM(ncid,'cell',IGRILLE,nri)
133     ier=NF_DEF_DIM(ncid,'tstep',NF_UNLIMITED,nti)
134
135   ! The longitude
136     ndim=1
137     dims(1)=nxi
138     ier=NF_DEF_VAR(ncid,'longitude',NF_FLOAT,ndim,dims,nlongi)
139     status=NF_PUT_ATT_TEXT (ncid,nlongi,'units',12,'degrees_east')
140     status=NF_PUT_ATT_TEXT (ncid,nlongi,'long_name',9,'Longitude')
141
142   ! The latitude
143     ndim=1
144     dims(1)=nyi
145     ier=NF_DEF_VAR(ncid,'latitude',NF_FLOAT,ndim,dims,nlati)
146     status=NF_PUT_ATT_TEXT (ncid,nlati,'units',13,'degrees_north')
147     status=NF_PUT_ATT_TEXT (ncid,nlati,'long_name',8,'Latitude')
148
149   ! The mask
150     ndim=2
151     dims(1)=nxi
152     dims(2)=nyi
153     ier=NF_DEF_VAR(ncid,'mask',NF_INT,ndim,dims,nmask)
154
155   ! The cell, indice_ligne, indice_col
156     ndim=1
157     dims(1)=nri
158     ier=NF_DEF_VAR(ncid,'cell',NF_INT,ndim,dims,ncell)
159     status=NF_PUT_ATT_TEXT (ncid,ncell,'compress',18,'latitude longitude')
160     ier=NF_DEF_VAR(ncid,'indice_ligne',NF_FLOAT,ndim,dims,nindligne)
161     ier=NF_DEF_VAR(ncid,'indice_col',NF_FLOAT,ndim,dims,nindcol)
162
163   ! The step:
164     day=1
165     month=1
166     ndim=1
167     dims(1)=nti
168     ier=NF_DEF_VAR(ncid,'time',NF_FLOAT,ndim,dims,ntime1)
169     write(str,1002) IANEXT,month,day
170     status=NF_PUT_ATT_TEXT (ncid,ntime1,'units',LEN_TRIM(str),str)
171     status=NF_PUT_ATT_TEXT (ncid,ntime1,'calendar',6,'noleap')
172     status=NF_PUT_ATT_TEXT (ncid,ntime1,'long_name',9,'Time axis')
173     status=NF_PUT_ATT_REAL (ncid,ntime1,'tstep_sec',NF_FLOAT,1,secs_step)
174
175   ! The meterology variables:
176     ndim=2
177     dims(1)=nri
178     dims(2)=nti
179     DO ii=1,8
180            ier=NF_DEF_VAR(ncid,nomvar(ii),NF_FLOAT,ndim,dims,nuvar(ii))
181            status=NF_PUT_ATT_TEXT (ncid,nuvar(ii),'long_name',LEN_TRIM(nomvar_long(ii)),nomvar_long(ii))
182            status=NF_PUT_ATT_TEXT (ncid,nuvar(ii),'units',LEN_TRIM(nomunit(ii)),nomunit(ii))
183            status=NF_PUT_ATT_TEXT (ncid,nuvar(ii),'coordinates',10,'cell tstep')
184            nval=1
185            status=NF_PUT_ATT_REAL (ncid,nuvar(ii),'missing_value',NF_FLOAT,nval,missing_value)
186     ENDDO
187
188   ! End definition
189     ier=NF_ENDDEF(ncid)
190
191     Do ii=1,24*365
192        timestp(ii)=ii-1
193        time(ii)=timestp(ii)*secs_step
194     Enddo
195     corner(1)=1
196     edge(1)=24*365
197     ier =NF_PUT_VARA_REAL(ncid, ntime1, corner, edge, time)
198
199
200     corner(1)=1
201     edge(1)=NBLAT
202     ier =NF_PUT_VARA_REAL(ncid, nlati, corner, edge, LAT)
203
204     corner(1)=1
205     edge(1)=NBLON
206     ier =NF_PUT_VARA_REAL(ncid, nlongi, corner, edge,LON)
207
208     corner(1)=1
209     corner(2)=1
210     edge(1)=NBLON
211     edge(2)=NBLAT
212     ier =NF_PUT_VARA_INT(ncid, nmask, corner, edge, mask)
213
214     corner(1)=1
215     edge(1)=IGRILLE
216     ier =NF_PUT_VARA_INT(ncid, ncell, corner, edge, cell)
217     ier =NF_PUT_VARA_REAL(ncid, nindligne, corner, edge, indice_ligne)
218     ier =NF_PUT_VARA_REAL(ncid, nindcol, corner, edge, indice_col)
219         
220     corner(1)=1
221     corner(2)=1
222     edge(1)=IGRILLE
223     edge(2)=24*365
224
225     ier =NF_PUT_VARA_REAL(ncid,nuvar(1),corner, edge,Tair)
226     ier =NF_PUT_VARA_REAL(ncid,nuvar(2),corner, edge,Psurf)
227     ier =NF_PUT_VARA_REAL(ncid,nuvar(3),corner, edge,Qair)
228     ier =NF_PUT_VARA_REAL(ncid,nuvar(4),corner, edge,Wind)
229     ier =NF_PUT_VARA_REAL(ncid,nuvar(5),corner, edge,Rainf)
230     ier =NF_PUT_VARA_REAL(ncid,nuvar(6),corner, edge,Snowf)
231     ier =NF_PUT_VARA_REAL(ncid,nuvar(7),corner, edge,SWdown)
232     ier =NF_PUT_VARA_REAL(ncid,nuvar(8),corner, edge,LWdown)
233     ier=NF_CLOSE(ncid)
234
235     print*,"End of ", IANEXT
236 Enddo  ! IANEXT
237
238
2391001 FORMAT(A,'FR.',I4.4,'_Safran.nc')
2401002 FORMAT('seconds since ',I4.4,'-',I2.2,'-',I2.2,' 00:00:00')
2411003 FORMAT('timesteps since ',I4.4,'-',I2.2,'-',I2.2,' 00:00:00')
242   
243
244END PROGRAM ANA_FORCING
245
246