source: trunk/LMDZ.MARS/libf/dyn3d/start2archive.F @ 664

Last change on this file since 664 was 410, checked in by acolaitis, 13 years ago

Modified NCDF norm of our files from classic to 64bit offset to support variables indices of more than integer*4 maximum length. This is a priori retrocompatible with classic format. I have tested it, it works for all file outputs, newstart (on classic or 64-bit offset files). One can check the format of his .nc with ncdump -k file.nc

File size: 13.8 KB
Line 
1c=======================================================================
2      PROGRAM start2archive
3c=======================================================================
4c
5c
6c   Date:    01/1997
7c   ----
8c
9c
10c   Objet:   Passage des  fichiers netcdf d'etat initial "start" et
11c   -----    "startfi" a un fichier netcdf unique "start_archive"
12c
13c  "start_archive" est une banque d'etats initiaux:
14c  On peut stocker plusieurs etats initiaux dans un meme fichier "start_archive"
15c    (Veiller dans ce cas avoir un day_ini different pour chacun des start)
16c
17c
18c
19c=======================================================================
20
21      implicit none
22
23#include "dimensions.h"
24#include "paramet.h"
25#include "comconst.h"
26#include "comdissip.h"
27#include "comvert.h"
28#include "comgeom.h"
29#include "logic.h"
30#include "temps.h"
31#include "control.h"
32#include "ener.h"
33#include "description.h"
34
35#include "dimphys.h"
36#include "comsoil.h"
37#include"advtrac.h"
38#include "netcdf.inc"
39
40c-----------------------------------------------------------------------
41c   Declarations
42c-----------------------------------------------------------------------
43
44c variables dynamiques du GCM
45c -----------------------------
46      REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm) ! vents covariants
47      REAL teta(ip1jmp1,llm)                    ! temperature potentielle
48      REAL q(ip1jmp1,llm,nqmx)               ! champs advectes
49      REAL pks(ip1jmp1)                      ! exner (f pour filtre)
50      REAL pk(ip1jmp1,llm)
51      REAL pkf(ip1jmp1,llm)
52      REAL beta(iip1,jjp1,llm)
53      REAL phis(ip1jmp1)                     ! geopotentiel au sol
54      REAL masse(ip1jmp1,llm)                ! masse de l'atmosphere
55      REAL ps(ip1jmp1)                       ! pression au sol
56      REAL p3d(iip1, jjp1, llm+1)            ! pression aux interfaces
57     
58c Variable Physiques (grille physique)
59c ------------------------------------
60      REAL tsurf(ngridmx)       ! Surface temperature
61      REAL tsoil(ngridmx,nsoilmx) ! Soil temperature
62      REAL co2ice(ngridmx)      ! CO2 ice layer
63      REAL q2(ngridmx,nlayermx+1),qsurf(ngridmx,nqmx)
64      REAL emis(ngridmx)
65      INTEGER start,length
66      PARAMETER (length = 100)
67      REAL tab_cntrl_fi(length) ! tableau des parametres de startfi
68      INTEGER*4 day_ini_fi
69
70c Variable naturelle / grille scalaire
71c ------------------------------------
72      REAL T(ip1jmp1,llm),us(ip1jmp1,llm),vs(ip1jmp1,llm)
73      REAL tsurfS(ip1jmp1)
74      REAL tsoilS(ip1jmp1,nsoilmx)
75      REAL ithS(ip1jmp1,nsoilmx) ! Soil Thermal Inertia
76      REAL co2iceS(ip1jmp1)
77      REAL q2S(ip1jmp1,llm+1),qsurfS(ip1jmp1,nqmx)
78      REAL emisS(ip1jmp1)
79
80c Variables intermediaires : vent naturel, mais pas coord scalaire
81c----------------------------------------------------------------
82      REAL vn(ip1jm,llm),un(ip1jmp1,llm)
83
84c Autres  variables
85c -----------------
86      LOGICAL startdrs
87      INTEGER Lmodif
88
89      REAL ptotal, co2icetotal
90      REAL timedyn,timefi !fraction du jour dans start, startfi
91      REAL date
92
93      CHARACTER*2 str2
94      CHARACTER*80 fichier
95      data  fichier /'startfi'/
96
97      INTEGER ij, l,i,j,isoil,iq
98      character*80      fichnom
99      integer :: ierr,ntime
100      integer :: nq,numvanle
101      character(len=30) :: txt ! to store some text
102
103c Netcdf
104c-------
105      integer varid,dimid,timelen
106      INTEGER nid,nid1
107
108c-----------------------------------------------------------------------
109c   Initialisations
110c-----------------------------------------------------------------------
111
112      grireg   = .TRUE.
113
114c=======================================================================
115c Lecture des donnees
116c=======================================================================
117! Load tracer names:
118      call iniadvtrac(nq,numvanle)
119
120      fichnom = 'start.nc'
121      CALL dynetat0(fichnom,nqmx,vcov,ucov,teta,q,masse,
122     .       ps,phis,timedyn)
123
124
125      fichnom = 'startfi.nc'
126      Lmodif=0
127
128      CALL phyetat0 (fichnom,0,Lmodif,nsoilmx,nqmx,day_ini_fi,timefi,
129     .      tsurf,tsoil,emis,q2,qsurf,co2ice)
130
131       ierr = NF_OPEN (fichnom, NF_NOWRITE,nid1)
132       IF (ierr.NE.NF_NOERR) THEN
133         write(6,*)' Pb d''ouverture du fichier'//fichnom
134        CALL ABORT
135       ENDIF
136                                               
137      ierr = NF_INQ_VARID (nid1, "controle", varid)
138      IF (ierr .NE. NF_NOERR) THEN
139       PRINT*, "start2archive: Le champ <controle> est absent"
140       CALL abort
141      ENDIF
142#ifdef NC_DOUBLE
143       ierr = NF_GET_VAR_DOUBLE(nid1, varid, tab_cntrl_fi)
144#else
145      ierr = NF_GET_VAR_REAL(nid1, varid, tab_cntrl_fi)
146#endif
147       IF (ierr .NE. NF_NOERR) THEN
148          PRINT*, "start2archive: Lecture echoue pour <controle>"
149          CALL abort
150       ENDIF
151
152      ierr = NF_CLOSE(nid1)
153
154c-----------------------------------------------------------------------
155c Controle de la synchro
156c-----------------------------------------------------------------------
157!mars a voir      if ((day_ini_fi.ne.day_ini).or.(abs(timefi-timedyn).gt.1.e-10))
158      if ((day_ini_fi.ne.day_ini))
159     &  stop ' Probleme de Synchro entre start et startfi !!!'
160
161
162c *****************************************************************
163c    Option : Reinitialisation des dates dans la premieres annees :
164       do while (day_ini.ge.669)
165          day_ini=day_ini-669
166       enddo
167c *****************************************************************
168
169c-----------------------------------------------------------------------
170c   Initialisations
171c-----------------------------------------------------------------------
172
173      CALL defrun_new(99, .FALSE. )
174      call iniconst
175      call inigeom
176      call inifilr
177      CALL pression(ip1jmp1, ap, bp, ps, p3d)
178      call exner_hyb(ip1jmp1, ps, p3d, beta, pks, pk, pkf)
179
180c=======================================================================
181c Transformation EN VARIABLE NATURELLE / GRILLE SCALAIRE si necessaire
182c=======================================================================
183c  Les variables modeles dependent de la resolution. Il faut donc
184c  eliminer les facteurs responsables de cette dependance
185c  (pour utiliser newstart)
186c=======================================================================
187
188c-----------------------------------------------------------------------
189c Vent   (depend de la resolution horizontale)
190c-----------------------------------------------------------------------
191c
192c ucov --> un  et  vcov --> vn
193c un --> us  et   vn --> vs
194c
195c-----------------------------------------------------------------------
196
197      call covnat(llm,ucov, vcov, un, vn)
198      call wind_scal(un,vn,us,vs)
199
200c-----------------------------------------------------------------------
201c Temperature  (depend de la resolution verticale => de "sigma.def")
202c-----------------------------------------------------------------------
203c
204c h --> T
205c
206c-----------------------------------------------------------------------
207
208      DO l=1,llm
209         DO ij=1,ip1jmp1
210            T(ij,l)=teta(ij,l)*pk(ij,l)/cpp !mars deduit de l'equation dans newstart
211         ENDDO
212      ENDDO
213
214c-----------------------------------------------------------------------
215c Variable physique
216c-----------------------------------------------------------------------
217c
218c tsurf --> tsurfS
219c co2ice --> co2iceS
220c tsoil --> tsoilS
221c emis --> emisS
222c q2 --> q2S
223c qsurf --> qsurfS
224c
225c-----------------------------------------------------------------------
226
227      call gr_fi_dyn(1,ngridmx,iip1,jjp1,tsurf,tsurfS)
228      call gr_fi_dyn(1,ngridmx,iip1,jjp1,co2ice,co2iceS)
229      call gr_fi_dyn(nsoilmx,ngridmx,iip1,jjp1,tsoil,tsoilS)
230      ! Note: thermal inertia "inertiedat" is in comsoil.h
231      call gr_fi_dyn(nsoilmx,ngridmx,iip1,jjp1,inertiedat,ithS)
232      call gr_fi_dyn(1,ngridmx,iip1,jjp1,emis,emisS)
233      call gr_fi_dyn(llm+1,ngridmx,iip1,jjp1,q2,q2S)
234      call gr_fi_dyn(nqmx,ngridmx,iip1,jjp1,qsurf,qsurfS)
235
236c=======================================================================
237c Info pour controler
238c=======================================================================
239
240      ptotal =  0.
241      co2icetotal = 0.
242      DO j=1,jjp1
243         DO i=1,iim
244           ptotal=ptotal+aire(i+(iim+1)*(j-1))*ps(i+(iim+1)*(j-1))/g
245           co2icetotal = co2icetotal +
246     &            co2iceS(i+(iim+1)*(j-1))*aire(i+(iim+1)*(j-1))
247         ENDDO
248      ENDDO
249      write(*,*)'Ancienne grille : masse de l''atm :',ptotal
250      write(*,*)'Ancienne grille : masse de la glace CO2 :',co2icetotal
251
252c-----------------------------------------------------------------------
253c Passage de "ptotal" et "co2icetotal" par tab_cntrl_fi
254c-----------------------------------------------------------------------
255
256      tab_cntrl_fi(49) = ptotal
257      tab_cntrl_fi(50) = co2icetotal
258
259c=======================================================================
260c Ecriture dans le fichier  "start_archive"
261c=======================================================================
262
263c-----------------------------------------------------------------------
264c Ouverture de "start_archive"
265c-----------------------------------------------------------------------
266
267      ierr = NF_OPEN ('start_archive.nc', NF_WRITE,nid)
268 
269c-----------------------------------------------------------------------
270c  si "start_archive" n'existe pas:
271c    1_ ouverture
272c    2_ creation de l'entete dynamique ("ini_archive")
273c-----------------------------------------------------------------------
274c ini_archive:
275c On met dans l'entete le tab_cntrl dynamique (1 a 16)
276c  On y ajoute les valeurs du tab_cntrl_fi (a partir de 51)
277c  En plus les deux valeurs ptotal et co2icetotal (99 et 100)
278c-----------------------------------------------------------------------
279
280      if (ierr.ne.NF_NOERR) then
281         write(*,*)'OK, Could not open file "start_archive.nc"'
282         write(*,*)'So let s create a new "start_archive"'
283         ierr = NF_CREATE('start_archive.nc',
284     &  IOR(NF_CLOBBER,NF_64BIT_OFFSET), nid)
285         call ini_archive(nid,day_ini,phis,ithS,tab_cntrl_fi)
286      endif
287
288c-----------------------------------------------------------------------
289c Ecriture de la coordonnee temps (date en jours)
290c-----------------------------------------------------------------------
291
292      date = day_ini
293      ierr= NF_INQ_VARID(nid,"Time",varid)
294      ierr= NF_INQ_DIMID(nid,"Time",dimid)
295      ierr= NF_INQ_DIMLEN(nid,dimid,timelen)
296      ntime=timelen+1
297
298      write(*,*) "******************"
299      write(*,*) "ntime",ntime
300      write(*,*) "******************"
301#ifdef NC_DOUBLE
302      ierr= NF_PUT_VARA_DOUBLE(nid,varid,ntime,1,date)
303#else
304      ierr= NF_PUT_VARA_REAL(nid,varid,ntime,1,date)
305#endif
306      if (ierr.ne.NF_NOERR) then
307         write(*,*) "time matter ",NF_STRERROR(ierr)
308         stop
309      endif
310
311c-----------------------------------------------------------------------
312c Ecriture des champs  (co2ice,emis,ps,Tsurf,T,u,v,q2,q,qsurf)
313c-----------------------------------------------------------------------
314c ATTENTION: q2 a une couche de plus!!!!
315c    Pour creer un fichier netcdf lisible par grads,
316c    On passe donc une des couches de q2 a part
317c    comme une variable 2D (la couche au sol: "q2surf")
318c    Les lmm autres couches sont nommees "q2atm" (3D)
319c-----------------------------------------------------------------------
320
321      call write_archive(nid,ntime,'co2ice','couche de glace co2',
322     &  'kg/m2',2,co2iceS)
323      call write_archive(nid,ntime,'emis','grd emis',' ',2,emisS)
324      call write_archive(nid,ntime,'ps','Psurf','Pa',2,ps)
325      call write_archive(nid,ntime,'tsurf','surf T','K',2,tsurfS)
326      call write_archive(nid,ntime,'temp','temperature','K',3,t)
327      call write_archive(nid,ntime,'u','Vent zonal','m.s-1',3,us)
328      call write_archive(nid,ntime,'v','Vent merid','m.s-1',3,vs)
329      call write_archive(nid,ntime,'q2surf','wind variance','m2.s-2',2,
330     .              q2S)
331      call write_archive(nid,ntime,'q2atm','wind variance','m2.s-2',3,
332     .              q2S(1,2))
333
334c-----------------------------------------------------------------------
335c Ecriture du champs  q  ( q[1,nqmx] )
336c-----------------------------------------------------------------------
337      do iq=1,nqmx
338c       write(str2,'(i2.2)') iq
339c        call write_archive(nid,ntime,'q'//str2,'tracer','kg/kg',
340c     .         3,q(1,1,iq))
341        call write_archive(nid,ntime,tnom(iq),'tracer','kg/kg',
342     &         3,q(1,1,iq))
343      end do
344c-----------------------------------------------------------------------
345c Ecriture du champs  qsurf  ( qsurf[1,nqmx] )
346c-----------------------------------------------------------------------
347      do iq=1,nqmx
348c       write(str2,'(i2.2)') iq
349c       call write_archive(nid,ntime,'qsurf'//str2,'Tracer on surface',
350c     $  'kg.m-2',2,qsurfS(1,iq))
351        txt=trim(tnom(iq))//"_surf"
352        call write_archive(nid,ntime,txt,'Tracer on surface',
353     &  'kg.m-2',2,qsurfS(1,iq))
354      enddo
355
356
357c-----------------------------------------------------------------------
358c Ecriture du champs  tsoil  ( Tg[1,10] )
359c-----------------------------------------------------------------------
360c "tsoil" Temperature au sol definie dans 10 couches dans le sol
361c   Les 10 couches sont lues comme 10 champs
362c  nommees Tg[1,10]
363
364c      do isoil=1,nsoilmx
365c       write(str2,'(i2.2)') isoil
366c       call write_archive(nid,ntime,'Tg'//str2,'Ground Temperature ',
367c     .   'K',2,tsoilS(1,isoil))
368c      enddo
369
370! Write soil temperatures tsoil
371      call write_archive(nid,ntime,'tsoil','Soil temperature',
372     &     'K',-3,tsoilS)
373
374! Write soil thermal inertia
375      call write_archive(nid,ntime,'inertiedat',
376     &     'Soil thermal inertia',
377     &     'J.s-1/2.m-2.K-1',-3,ithS)
378
379! Write (0D) volumetric heat capacity (stored in comsoil.h)
380!      call write_archive(nid,ntime,'volcapa',
381!     &     'Soil volumetric heat capacity',
382!     &     'J.m-3.K-1',0,volcapa)
383! Note: no need to write volcapa, it is stored in "controle" table
384
385      ierr=NF_CLOSE(nid)
386c-----------------------------------------------------------------------
387c Fin
388c-----------------------------------------------------------------------
389
390      end
Note: See TracBrowser for help on using the repository browser.