source: trunk/LMDZ.GENERIC/libf/phystd/start2archive.F @ 993

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

Generic GCM:

  • Some more cleanup in dynamics:
    • Moved "start2archive" (and auxilliary routines) to phystd
    • removed unused (obsolete) testharm.F , para_netcdf.h , readhead_NC.F , angtot.h from dyn3d
    • removed obsolete addit.F (and change corresponding lines in gcm)
    • remove unused "description.h" (and many places where it was "included")

EM

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