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

Last change on this file since 847 was 837, checked in by aslmd, 13 years ago

LMDZ.GENERIC. Corrected problems with allocated arrays in start2archive and newstart. Applied a workaround to make those work without tracers (-cpp NOTRAC -- perhaps there is a better solution). Checked that everything works in debug mode.

File size: 15.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 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
140      fichnom = 'start.nc'
141      CALL dynetat0(fichnom,nqmx,vcov,ucov,teta,q,masse,
142     .       ps,phis,timedyn)
143
144! load 'controle' array from dynamics start file
145
146       ierr = NF_OPEN (fichnom, NF_NOWRITE,nid1)
147       IF (ierr.NE.NF_NOERR) THEN
148         write(6,*)' Pb d''ouverture du fichier'//trim(fichnom)
149        CALL ABORT
150       ENDIF
151                                               
152      ierr = NF_INQ_VARID (nid1, "controle", varid)
153      IF (ierr .NE. NF_NOERR) THEN
154       PRINT*, "start2archive: Le champ <controle> est absent"
155       CALL abort
156      ENDIF
157#ifdef NC_DOUBLE
158       ierr = NF_GET_VAR_DOUBLE(nid1, varid, tab_cntrl_dyn)
159#else
160      ierr = NF_GET_VAR_REAL(nid1, varid, tab_cntrl_dyn)
161#endif
162       IF (ierr .NE. NF_NOERR) THEN
163          PRINT*, "start2archive: Lecture echoue pour <controle>"
164          CALL abort
165       ENDIF
166
167      ierr = NF_CLOSE(nid1)
168     
169
170      fichnom = 'startfi.nc'
171      Lmodif=0
172
173      CALL phyetat0 (ngridmx,fichnom,0,Lmodif,nsoilmx,nqmx,day_ini_fi,
174     .      timefi,
175     .      tsurf,tsoil,emis,q2,qsurf,
176!       change FF 05/2011
177     .       cloudfrac,totalcloudfrac,hice)
178
179
180! load 'controle' array from physics start file
181
182       ierr = NF_OPEN (fichnom, NF_NOWRITE,nid1)
183       IF (ierr.NE.NF_NOERR) THEN
184         write(6,*)' Pb d''ouverture du fichier'//trim(fichnom)
185        CALL ABORT
186       ENDIF
187                                               
188      ierr = NF_INQ_VARID (nid1, "controle", varid)
189      IF (ierr .NE. NF_NOERR) THEN
190       PRINT*, "start2archive: Le champ <controle> est absent"
191       CALL abort
192      ENDIF
193#ifdef NC_DOUBLE
194       ierr = NF_GET_VAR_DOUBLE(nid1, varid, tab_cntrl_fi)
195#else
196      ierr = NF_GET_VAR_REAL(nid1, varid, tab_cntrl_fi)
197#endif
198       IF (ierr .NE. NF_NOERR) THEN
199          PRINT*, "start2archive: Lecture echoue pour <controle>"
200          CALL abort
201       ENDIF
202
203      ierr = NF_CLOSE(nid1)
204
205
206c-----------------------------------------------------------------------
207c Controle de la synchro
208c-----------------------------------------------------------------------
209!mars a voir      if ((day_ini_fi.ne.day_ini).or.(abs(timefi-timedyn).gt.1.e-10))
210      if ((day_ini_fi.ne.day_ini))
211     &  stop ' Probleme de Synchro entre start et startfi !!!'
212
213
214c *****************************************************************
215c    Option : Reinitialisation des dates dans la premieres annees :
216       do while (day_ini.ge.year_day)
217          day_ini=day_ini-year_day
218       enddo
219c *****************************************************************
220
221c-----------------------------------------------------------------------
222c   Initialisations
223c-----------------------------------------------------------------------
224
225      CALL defrun_new(99, .FALSE. )
226      call iniconst
227      call inigeom
228      call inifilr
229      CALL pression(ip1jmp1, ap, bp, ps, p3d)
230      call exner_hyb(ip1jmp1, ps, p3d, beta, pks, pk, pkf)
231
232c=======================================================================
233c Transformation EN VARIABLE NATURELLE / GRILLE SCALAIRE si necessaire
234c=======================================================================
235c  Les variables modeles dependent de la resolution. Il faut donc
236c  eliminer les facteurs responsables de cette dependance
237c  (pour utiliser newstart)
238c=======================================================================
239
240c-----------------------------------------------------------------------
241c Vent   (depend de la resolution horizontale)
242c-----------------------------------------------------------------------
243c
244c ucov --> un  et  vcov --> vn
245c un --> us  et   vn --> vs
246c
247c-----------------------------------------------------------------------
248
249      call covnat(llm,ucov, vcov, un, vn)
250      call wind_scal(un,vn,us,vs)
251
252c-----------------------------------------------------------------------
253c Temperature  (depend de la resolution verticale => de "sigma.def")
254c-----------------------------------------------------------------------
255c
256c h --> T
257c
258c-----------------------------------------------------------------------
259
260      DO l=1,llm
261         DO ij=1,ip1jmp1
262            T(ij,l)=teta(ij,l)*pk(ij,l)/cpp !mars deduit de l'equation dans newstart
263         ENDDO
264      ENDDO
265
266c-----------------------------------------------------------------------
267c Variable physique
268c-----------------------------------------------------------------------
269c
270c tsurf --> tsurfS
271c co2ice --> co2iceS
272c tsoil --> tsoilS
273c emis --> emisS
274c q2 --> q2S
275c qsurf --> qsurfS
276c
277c-----------------------------------------------------------------------
278
279      call gr_fi_dyn(1,ngridmx,iip1,jjp1,tsurf,tsurfS)
280!      call gr_fi_dyn(1,ngridmx,iip1,jjp1,co2ice,co2iceS)
281      call gr_fi_dyn(nsoilmx,ngridmx,iip1,jjp1,tsoil,tsoilS)
282      ! Note: thermal inertia "inertiedat" is in comsoil.h
283      call gr_fi_dyn(nsoilmx,ngridmx,iip1,jjp1,inertiedat,ithS)
284      call gr_fi_dyn(1,ngridmx,iip1,jjp1,emis,emisS)
285      call gr_fi_dyn(llm+1,ngridmx,iip1,jjp1,q2,q2S)
286      call gr_fi_dyn(nqmx,ngridmx,iip1,jjp1,qsurf,qsurfS)
287      call gr_fi_dyn(llm,ngridmx,iip1,jjp1,cloudfrac,cloudfracS)
288      call gr_fi_dyn(1,ngridmx,iip1,jjp1,hice,hiceS)
289      call gr_fi_dyn(1,ngridmx,iip1,jjp1,totalcloudfrac,totalcloudfracS)
290
291c=======================================================================
292c Info pour controler
293c=======================================================================
294
295      ptotal =  0.
296      co2icetotal = 0.
297      DO j=1,jjp1
298         DO i=1,iim
299           ptotal=ptotal+aire(i+(iim+1)*(j-1))*ps(i+(iim+1)*(j-1))/g
300!           co2icetotal = co2icetotal +
301!     &            co2iceS(i+(iim+1)*(j-1))*aire(i+(iim+1)*(j-1))
302         ENDDO
303      ENDDO
304      write(*,*)'Ancienne grille : masse de l''atm :',ptotal
305!      write(*,*)'Ancienne grille : masse de la glace CO2 :',co2icetotal
306
307c-----------------------------------------------------------------------
308c Passage de "ptotal" et "co2icetotal" par tab_cntrl_fi
309c-----------------------------------------------------------------------
310
311      tab_cntrl_fi(49) = ptotal
312      tab_cntrl_fi(50) = co2icetotal
313
314c=======================================================================
315c Ecriture dans le fichier  "start_archive"
316c=======================================================================
317
318c-----------------------------------------------------------------------
319c Ouverture de "start_archive"
320c-----------------------------------------------------------------------
321
322      ierr = NF_OPEN ('start_archive.nc', NF_WRITE,nid)
323 
324c-----------------------------------------------------------------------
325c  si "start_archive" n'existe pas:
326c    1_ ouverture
327c    2_ creation de l'entete dynamique ("ini_archive")
328c-----------------------------------------------------------------------
329c ini_archive:
330c On met dans l'entete le tab_cntrl dynamique (1 a 16)
331c  On y ajoute les valeurs du tab_cntrl_fi (a partir de 51)
332c  En plus les deux valeurs ptotal et co2icetotal (99 et 100)
333c-----------------------------------------------------------------------
334
335      if (ierr.ne.NF_NOERR) then
336         write(*,*)'OK, Could not open file "start_archive.nc"'
337         write(*,*)'So let s create a new "start_archive"'
338         ierr = NF_CREATE('start_archive.nc', NF_CLOBBER, nid)
339         call ini_archive(nid,day_ini,phis,ithS,tab_cntrl_fi,
340     &                                          tab_cntrl_dyn)
341      endif
342
343c-----------------------------------------------------------------------
344c Ecriture de la coordonnee temps (date en jours)
345c-----------------------------------------------------------------------
346
347      date = day_ini
348      ierr= NF_INQ_VARID(nid,"Time",varid)
349      ierr= NF_INQ_DIMID(nid,"Time",dimid)
350      ierr= NF_INQ_DIMLEN(nid,dimid,timelen)
351      ntime=timelen+1
352
353      write(*,*) "******************"
354      write(*,*) "ntime",ntime
355      write(*,*) "******************"
356#ifdef NC_DOUBLE
357      ierr= NF_PUT_VARA_DOUBLE(nid,varid,ntime,1,date)
358#else
359      ierr= NF_PUT_VARA_REAL(nid,varid,ntime,1,date)
360#endif
361      if (ierr.ne.NF_NOERR) then
362         write(*,*) "time matter ",NF_STRERROR(ierr)
363         stop
364      endif
365
366c-----------------------------------------------------------------------
367c Ecriture des champs  (co2ice,emis,ps,Tsurf,T,u,v,q2,q,qsurf)
368c-----------------------------------------------------------------------
369c ATTENTION: q2 a une couche de plus!!!!
370c    Pour creer un fichier netcdf lisible par grads,
371c    On passe donc une des couches de q2 a part
372c    comme une variable 2D (la couche au sol: "q2surf")
373c    Les lmm autres couches sont nommees "q2atm" (3D)
374c-----------------------------------------------------------------------
375
376!      call write_archive(nid,ntime,'co2ice','couche de glace co2',
377!     &  'kg/m2',2,co2iceS)
378      call write_archive(nid,ntime,'emis','grd emis',' ',2,emisS)
379      call write_archive(nid,ntime,'ps','Psurf','Pa',2,ps)
380      call write_archive(nid,ntime,'tsurf','surf T','K',2,tsurfS)
381      call write_archive(nid,ntime,'temp','temperature','K',3,t)
382      call write_archive(nid,ntime,'u','Vent zonal','m.s-1',3,us)
383      call write_archive(nid,ntime,'v','Vent merid','m.s-1',3,vs)
384      call write_archive(nid,ntime,'q2surf','wind variance','m2.s-2',2,
385     .              q2S)
386      call write_archive(nid,ntime,'q2atm','wind variance','m2.s-2',3,
387     .              q2S(1,2))
388
389c-----------------------------------------------------------------------
390c Ecriture du champs  q  ( q[1,nqmx] )
391c-----------------------------------------------------------------------
392      do iq=1,nqmx
393        call write_archive(nid,ntime,tnom(iq),'tracer','kg/kg',
394     &         3,q(1,1,iq))
395      end do
396c-----------------------------------------------------------------------
397c Ecriture du champs  qsurf  ( qsurf[1,nqmx] )
398c-----------------------------------------------------------------------
399      do iq=1,nqmx
400        txt=trim(tnom(iq))//"_surf"
401        call write_archive(nid,ntime,txt,'Tracer on surface',
402     &  'kg.m-2',2,qsurfS(1,iq))
403      enddo
404
405
406c-----------------------------------------------------------------------
407c Ecriture du champs  tsoil  ( Tg[1,10] )
408c-----------------------------------------------------------------------
409c "tsoil" Temperature au sol definie dans 10 couches dans le sol
410c   Les 10 couches sont lues comme 10 champs
411c  nommees Tg[1,10]
412
413c      do isoil=1,nsoilmx
414c       write(str2,'(i2.2)') isoil
415c       call write_archive(nid,ntime,'Tg'//str2,'Ground Temperature ',
416c     .   'K',2,tsoilS(1,isoil))
417c      enddo
418
419! Write soil temperatures tsoil
420      call write_archive(nid,ntime,'tsoil','Soil temperature',
421     &     'K',-3,tsoilS)
422
423! Write soil thermal inertia
424      call write_archive(nid,ntime,'inertiedat',
425     &     'Soil thermal inertia',
426     &     'J.s-1/2.m-2.K-1',-3,ithS)
427
428! Write (0D) volumetric heat capacity (stored in comsoil.h)
429!      call write_archive(nid,ntime,'volcapa',
430!     &     'Soil volumetric heat capacity',
431!     &     'J.m-3.K-1',0,volcapa)
432! Note: no need to write volcapa, it is stored in "controle" table
433
434c-----------------------------------------------------------------------
435c Ecriture du champs  cloudfrac,hice,totalcloudfrac
436c-----------------------------------------------------------------------
437      call write_archive(nid,ntime,'hice',
438     &         'Height of oceanic ice','m',2,hiceS)
439      call write_archive(nid,ntime,'totalcloudfrac',
440     &        'Total cloud Fraction','',2,totalcloudfracS)
441      call write_archive(nid,ntime,'cloudfrac'
442     &        ,'Cloud fraction','',3,cloudfracS)
443c-----------------------------------------------------------------------
444c Fin
445c-----------------------------------------------------------------------
446      ierr=NF_CLOSE(nid)
447
448      end
Note: See TracBrowser for help on using the repository browser.