source: trunk/MESOSCALE/LMDZ.MARS/libf_gcm/dyn3d/start2archive.F @ 3584

Last change on this file since 3584 was 57, checked in by aslmd, 14 years ago

mineur LMD_MM_MARS: ajout du GCM ancienne physique, systeme maintenant complet sur SVN (ne manque que la base de donnees d'etats initiaux)

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