source: trunk/LMDZ.GENERIC/libf/dyn3d/start2archive.F @ 937

Last change on this file since 937 was 850, checked in by aslmd, 12 years ago

LMDZ.GENERIC. bug fix for start2archive (there was a problem because initracer was not called). also corrected iniadvtrac so that it could be compiled with -t 0. compiler was complaining about lines with iadv(nqmx). those lines are actually not necessary (this could be included in the loop above). added an if statement in iniadvtrac so that traceur.def does not even need to be here when the model is compiled with -t 0

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