source: trunk/LMDZ.MARS/libf/phymars/start2archive.F @ 1242

Last change on this file since 1242 was 1230, checked in by emillour, 11 years ago

Mars GCM:
Fix the bug fix in start2archive.
EM

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