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

Last change on this file since 1255 was 1252, checked in by aslmd, 11 years ago

LMDZ.GENERIC LMDZ.COMMON LMDZ.UNIVERSAL. Bye Bye LMDZ.UNIVERSAL. Go to LMDZ.COMMON!

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