source: trunk/LMDZ.TITAN/libf/dynphy_lonlat/phytitan/start2archive.F @ 1789

Last change on this file since 1789 was 1789, checked in by jvatant, 7 years ago

Added the surface methane tank and put it in start files
--JVO

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