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

Last change on this file since 3026 was 2366, checked in by jvatant, 5 years ago

Titan GCM : Major maintenance catching up commits from the generic including :

  • r2356 and 2354 removing obsolete old dynamical core
  • various minor addition to physics and gestion of phys_state_var_mode, especially in dyn1d
  • adding MESOSCALE CPP keys around chemistry and microphysics (disabled in mesoscale for now)
File size: 18.9 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      USE comchem_h, only : cnames, nkim, nlaykim_up, preskim, ykim_up
24      USE ioipsl_getincom, only: getin
25      USE planete_mod, only: year_day
26      USE mod_const_mpi, ONLY: COMM_LMDZ
27      USE control_mod, only: planet_type
28      use filtreg_mod, only: inifilr
29      USE comvert_mod, ONLY: ap,bp
30      USE comconst_mod, ONLY: daysec,dtphys,rad,g,r,cpp
31      USE temps_mod, ONLY: day_ini
32      USE iniphysiq_mod, ONLY: iniphysiq
33      use phys_state_var_mod, only: phys_state_var_init
34      use phyetat0_mod, only: phyetat0
35      use tracer_h
36      use exner_hyb_m, only: exner_hyb
37      implicit none
38
39      include "dimensions.h"
40      integer, parameter :: ngridmx = (2+(jjm-1)*iim - 1/jjm)
41      include "paramet.h"
42      include "comdissip.h"
43      include "comgeom.h"
44
45      include "netcdf.inc"
46c-----------------------------------------------------------------------
47c   Declarations
48c-----------------------------------------------------------------------
49
50c variables dynamiques du GCM
51c -----------------------------
52      REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm) ! vents covariants
53      REAL teta(ip1jmp1,llm)                    ! temperature potentielle
54      REAL,ALLOCATABLE :: q(:,:,:)   ! champs advectes
55      REAL pks(ip1jmp1)                      ! exner (f pour filtre)
56      REAL pk(ip1jmp1,llm)
57      REAL pkf(ip1jmp1,llm)
58      REAL phis(ip1jmp1)                     ! geopotentiel au sol
59      REAL masse(ip1jmp1,llm)                ! masse de l'atmosphere
60      REAL ps(ip1jmp1)                       ! pression au sol
61      REAL p3d(iip1, jjp1, llm+1)            ! pression aux interfaces
62     
63c Variable Physiques (grille physique)
64c ------------------------------------
65      REAL tsurf(ngridmx)        ! Surface temperature
66      REAL,ALLOCATABLE :: tsoil(:,:) ! Soil temperature
67      REAL q2(ngridmx,llm+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
76c     Added by JVO for Titan specifities
77      REAL tankCH4(ngridmx) ! Depth of surface methane tank
78
79c Variable naturelle / grille scalaire
80c ------------------------------------
81      REAL T(ip1jmp1,llm),us(ip1jmp1,llm),vs(ip1jmp1,llm)
82      REAL tsurfS(ip1jmp1)
83      REAL,ALLOCATABLE :: tsoilS(:,:)
84      REAL,ALLOCATABLE :: ithS(:,:) ! Soil Thermal Inertia
85      REAL q2S(ip1jmp1,llm+1)
86      REAL,ALLOCATABLE :: qsurfS(:,:)
87      REAL emisS(ip1jmp1)
88
89c     Added by JVO for Titan specifities
90      REAL tankCH4S(ip1jmp1)  ! Depth of surface methane tank
91
92      REAL,ALLOCATABLE :: ykim_upS(:,:,:)  ! Upper chemistry fields
93 
94c Variables intermediaires : vent naturel, mais pas coord scalaire
95c----------------------------------------------------------------
96      REAL vn(ip1jm,llm),un(ip1jmp1,llm)
97
98c Autres  variables
99c -----------------
100      LOGICAL startdrs
101      INTEGER Lmodif
102
103      LOGICAL kim
104
105      REAL ptotal
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
130      planet_type="titan"
131
132c=======================================================================
133c Lecture des donnees
134c=======================================================================
135! Load tracer number and names:
136      call infotrac_init
137
138! allocate arrays:
139      allocate(q(ip1jmp1,llm,nqtot))
140      allocate(qsurf(ngridmx,nqtot))
141      allocate(qsurfS(ip1jmp1,nqtot))
142! other array allocations:
143!      call ini_comsoil_h(ngridmx) ! done via iniphysiq
144
145      fichnom = 'start.nc'
146      CALL dynetat0(fichnom,vcov,ucov,teta,q,masse,
147     .       ps,phis,timedyn)
148
149! load 'controle' array from dynamics start file
150
151       ierr = NF_OPEN (fichnom, NF_NOWRITE,nid1)
152       IF (ierr.NE.NF_NOERR) THEN
153         write(6,*)' Pb d''ouverture du fichier'//trim(fichnom)
154        CALL ABORT
155       ENDIF
156                                               
157      ierr = NF_INQ_VARID (nid1, "controle", varid)
158      IF (ierr .NE. NF_NOERR) THEN
159       PRINT*, "start2archive: Le champ <controle> est absent"
160       CALL abort
161      ENDIF
162#ifdef NC_DOUBLE
163       ierr = NF_GET_VAR_DOUBLE(nid1, varid, tab_cntrl_dyn)
164#else
165      ierr = NF_GET_VAR_REAL(nid1, varid, tab_cntrl_dyn)
166#endif
167       IF (ierr .NE. NF_NOERR) THEN
168          PRINT*, "start2archive: Lecture echoue pour <controle>"
169          CALL abort
170       ENDIF
171
172      ierr = NF_CLOSE(nid1)
173
174! Get value of the "subsurface_layers" dimension from physics start file
175      fichnom = 'startfi.nc'
176      ierr = NF_OPEN (fichnom, NF_NOWRITE,nid1)
177       IF (ierr.NE.NF_NOERR) THEN
178         write(6,*)' Pb d''ouverture du fichier'//trim(fichnom)
179        CALL ABORT
180       ENDIF
181      ierr = NF_INQ_DIMID(nid1,"subsurface_layers",varid)
182      IF (ierr .NE. NF_NOERR) THEN
183       PRINT*, "start2archive: No subsurface_layers dimension!!"
184       CALL abort
185      ENDIF
186      ierr = NF_INQ_DIMLEN(nid1,varid,nsoilmx)
187      IF (ierr .NE. NF_NOERR) THEN
188       PRINT*, "start2archive: Failed reading subsurface_layers value!!"
189       CALL abort
190      ENDIF
191     
192      ! allocate arrays of nsoilmx size
193      allocate(tsoil(ngridmx,nsoilmx))
194      allocate(tsoilS(ip1jmp1,nsoilmx))
195      allocate(ithS(ip1jmp1,nsoilmx))
196
197! Get value of the "upper_chemistry_layers" dimension from physics start file
198
199      ! NB : nlaykim_up is recalculated in phyetat0 ... this is redundant for newstart
200      ! If you want to do it properly put ykim_upS in a specific separate header
201      ! Anyway it should be safe. -- JVO 18
202
203      ierr = NF_INQ_DIMID(nid1,"upper_chemistry_layers",varid)
204      IF (ierr .NE. NF_NOERR) THEN
205       PRINT*, "start2archive: No upper_chemistry_layers dimension!!"
206       CALL abort
207      ENDIF
208      ierr = NF_INQ_DIMLEN(nid1,varid,nlaykim_up)
209      IF (ierr .NE. NF_NOERR) THEN
210       PRINT*, "start2archive: Failed reading
211     . upper_chemistry_layers value!!"
212       CALL abort
213      ENDIF
214
215      ALLOCATE(preskim(nlaykim_up))
216     
217      ! Allocate other arrays of nlaykim_up size, only if they're present
218      ! The test is on H but could be on any as we assume we can't do incomplete chemistry
219
220      ierr = NF_INQ_VARID(nid1,'H_up',varid)
221      IF (ierr .NE. NF_NOERR) THEN
222        PRINT*, "start2archive: Missing field(s) for upper chemistry ...
223     . I presume they're all absent !"
224        kim=.FALSE.
225      ELSE
226        PRINT*,"start2archive: I found a field for upper chemistry ...
227     . I presume they're all here as you can't do uncomplete chemistry!"
228     
229        ALLOCATE(ykim_up(nkim,ngridmx,nlaykim_up))
230        ALLOCATE(ykim_upS(nkim,ip1jmp1,nlaykim_up))
231       
232        kim=.TRUE.
233      ENDIF
234
235      ierr = NF_CLOSE(nid1)
236
237c-----------------------------------------------------------------------
238c   Initialisations
239c-----------------------------------------------------------------------
240
241      CALL defrun_new(99, .FALSE. )
242      call iniconst
243      call inigeom
244      call inifilr
245
246! Initialize the physics
247         CALL iniphysiq(iim,jjm,llm,
248     &                  (jjm-1)*iim+2,comm_lmdz,
249     &                  daysec,day_ini,dtphys,
250     &                  rlatu,rlatv,rlonu,rlonv,
251     &                  aire,cu,cv,rad,g,r,cpp,
252     &                  1)
253
254      fichnom = 'startfi.nc'
255      Lmodif=0
256
257! Allocate saved arrays (as in firstcall of physiq)
258      call phys_state_var_init(nqtot)
259     
260! Initialize tracer names, indexes and properties
261      CALL initracer2(nqtot,tname)
262
263      CALL phyetat0(.true.,ngridmx,llm,fichnom,0,Lmodif,nsoilmx,nqtot,
264     .      day_ini_fi,timefi,
265     .      tsurf,tsoil,emis,q2,qsurf,tankCH4)
266
267
268! load 'controle' array from physics start file
269
270       ierr = NF_OPEN (fichnom, NF_NOWRITE,nid1)
271       IF (ierr.NE.NF_NOERR) THEN
272         write(6,*)' Pb d''ouverture du fichier'//trim(fichnom)
273        CALL ABORT
274       ENDIF
275                                               
276      ierr = NF_INQ_VARID (nid1, "controle", varid)
277      IF (ierr .NE. NF_NOERR) THEN
278       PRINT*, "start2archive: Le champ <controle> est absent"
279       CALL abort
280      ENDIF
281#ifdef NC_DOUBLE
282       ierr = NF_GET_VAR_DOUBLE(nid1, varid, tab_cntrl_fi)
283#else
284      ierr = NF_GET_VAR_REAL(nid1, varid, tab_cntrl_fi)
285#endif
286       IF (ierr .NE. NF_NOERR) THEN
287          PRINT*, "start2archive: Lecture echoue pour <controle>"
288          CALL abort
289       ENDIF
290
291      ierr = NF_CLOSE(nid1)
292
293
294c-----------------------------------------------------------------------
295c Controle de la synchro
296c-----------------------------------------------------------------------
297!mars a voir      if ((day_ini_fi.ne.day_ini).or.(abs(timefi-timedyn).gt.1.e-10))
298      if ((day_ini_fi.ne.day_ini))
299     &  stop ' Probleme de Synchro entre start et startfi !!!'
300
301
302c *****************************************************************
303c    Option : Reinitialisation des dates dans la premieres annees :
304       do while (day_ini.ge.year_day)
305          day_ini=day_ini-year_day
306       enddo
307c *****************************************************************
308
309      CALL pression(ip1jmp1, ap, bp, ps, p3d)
310      call exner_hyb(ip1jmp1, ps, p3d, pks, pk, pkf)
311
312c=======================================================================
313c Transformation EN VARIABLE NATURELLE / GRILLE SCALAIRE si necessaire
314c=======================================================================
315c  Les variables modeles dependent de la resolution. Il faut donc
316c  eliminer les facteurs responsables de cette dependance
317c  (pour utiliser newstart)
318c=======================================================================
319
320c-----------------------------------------------------------------------
321c Vent   (depend de la resolution horizontale)
322c-----------------------------------------------------------------------
323c
324c ucov --> un  et  vcov --> vn
325c un --> us  et   vn --> vs
326c
327c-----------------------------------------------------------------------
328
329      call covnat(llm,ucov, vcov, un, vn)
330      call wind_scal(un,vn,us,vs)
331
332c-----------------------------------------------------------------------
333c Temperature  (depend de la resolution verticale => de "sigma.def")
334c-----------------------------------------------------------------------
335c
336c h --> T
337c
338c-----------------------------------------------------------------------
339
340      DO l=1,llm
341         DO ij=1,ip1jmp1
342            T(ij,l)=teta(ij,l)*pk(ij,l)/cpp !mars deduit de l'equation dans newstart
343         ENDDO
344      ENDDO
345
346c-----------------------------------------------------------------------
347c Variable physique
348c-----------------------------------------------------------------------
349c
350c tsurf --> tsurfS
351c tsoil --> tsoilS
352c emis --> emisS
353c q2 --> q2S
354c qsurf --> qsurfS
355c tankCH4 --> tankCH4S
356c ykim_up --> ykim_upS ( nkim=44 upper chemistry fields )
357c
358c-----------------------------------------------------------------------
359
360      call gr_fi_dyn(1,ngridmx,iip1,jjp1,tsurf,tsurfS)
361      call gr_fi_dyn(nsoilmx,ngridmx,iip1,jjp1,tsoil,tsoilS)
362      ! Note: thermal inertia "inertiedat" is in comsoil.h
363      call gr_fi_dyn(nsoilmx,ngridmx,iip1,jjp1,inertiedat,ithS)
364      call gr_fi_dyn(1,ngridmx,iip1,jjp1,emis,emisS)
365      call gr_fi_dyn(llm+1,ngridmx,iip1,jjp1,q2,q2S)
366      call gr_fi_dyn(nqtot,ngridmx,iip1,jjp1,qsurf,qsurfS)
367      call gr_fi_dyn(1,ngridmx,iip1,jjp1,tankCH4,tankCH4S)
368
369      IF (kim) THEN ! NB : fields are in comchem_startarch_h
370         DO iq=1,nkim
371           call gr_fi_dyn(nlaykim_up,ngridmx,iip1,jjp1,ykim_up(iq,:,:),
372     &                                                 ykim_upS(iq,:,:))
373         ENDDO
374      ENDIF
375
376c=======================================================================
377c Info pour controler
378c=======================================================================
379
380      ptotal =  0.
381      DO j=1,jjp1
382         DO i=1,iim
383           ptotal=ptotal+aire(i+(iim+1)*(j-1))*ps(i+(iim+1)*(j-1))/g
384         ENDDO
385      ENDDO
386      write(*,*)'Ancienne grille : masse de l''atm :',ptotal
387
388c-----------------------------------------------------------------------
389c Passage de "ptotal"  par tab_cntrl_fi
390c-----------------------------------------------------------------------
391
392      tab_cntrl_fi(49) = ptotal
393      tab_cntrl_fi(50) = 0.
394
395c=======================================================================
396c Ecriture dans le fichier  "start_archive"
397c=======================================================================
398
399c-----------------------------------------------------------------------
400c Ouverture de "start_archive"
401c-----------------------------------------------------------------------
402
403      ierr = NF_OPEN ('start_archive.nc', NF_WRITE,nid)
404 
405c-----------------------------------------------------------------------
406c  si "start_archive" n'existe pas:
407c    1_ ouverture
408c    2_ creation de l'entete dynamique ("ini_archive")
409c-----------------------------------------------------------------------
410c ini_archive:
411c On met dans l'entete le tab_cntrl dynamique (1 a 16)
412c  On y ajoute les valeurs du tab_cntrl_fi (a partir de 51)
413c  En plus les deux valeurs ptotal et co2icetotal (99 et 100)
414c-----------------------------------------------------------------------
415
416      if (ierr.ne.NF_NOERR) then
417         write(*,*)'OK, Could not open file "start_archive.nc"'
418         write(*,*)'So let s create a new "start_archive"'
419         ierr = NF_CREATE('start_archive.nc', NF_CLOBBER, nid)
420         call ini_archive(nid,day_ini,phis,ithS,tab_cntrl_fi,
421     &                                          tab_cntrl_dyn)
422      endif
423
424c-----------------------------------------------------------------------
425c Ecriture de la coordonnee temps (date en jours)
426c-----------------------------------------------------------------------
427
428      date = day_ini
429      ierr= NF_INQ_VARID(nid,"Time",varid)
430      ierr= NF_INQ_DIMID(nid,"Time",dimid)
431      ierr= NF_INQ_DIMLEN(nid,dimid,timelen)
432      ntime=timelen+1
433
434      write(*,*) "******************"
435      write(*,*) "ntime",ntime
436      write(*,*) "******************"
437#ifdef NC_DOUBLE
438      ierr= NF_PUT_VARA_DOUBLE(nid,varid,ntime,1,date)
439#else
440      ierr= NF_PUT_VARA_REAL(nid,varid,ntime,1,date)
441#endif
442      if (ierr.ne.NF_NOERR) then
443         write(*,*) "time matter ",NF_STRERROR(ierr)
444         stop
445      endif
446
447c-----------------------------------------------------------------------
448c Ecriture des champs  (emis,ps,Tsurf,T,u,v,q2,q,qsurf,tankCH4)
449c-----------------------------------------------------------------------
450c ATTENTION: q2 a une couche de plus!!!!
451c    Pour creer un fichier netcdf lisible par grads,
452c    On passe donc une des couches de q2 a part
453c    comme une variable 2D (la couche au sol: "q2surf")
454c    Les lmm autres couches sont nommees "q2atm" (3D)
455c-----------------------------------------------------------------------
456
457      call write_archive(nid,ntime,'emis','grd emis',' ',2,emisS)
458      call write_archive(nid,ntime,'ps','Psurf','Pa',2,ps)
459      call write_archive(nid,ntime,'tsurf','surf T','K',2,tsurfS)
460      call write_archive(nid,ntime,'temp','temperature','K',3,t)
461      call write_archive(nid,ntime,'u','Vent zonal','m.s-1',3,us)
462      call write_archive(nid,ntime,'v','Vent merid','m.s-1',3,vs)
463      call write_archive(nid,ntime,'q2surf','wind variance','m2.s-2',2,
464     .              q2S)
465      call write_archive(nid,ntime,'q2atm','wind variance','m2.s-2',3,
466     .              q2S(1,2))
467
468c-----------------------------------------------------------------------
469c Ecriture du champs  q  ( q[1,nqtot] )
470c-----------------------------------------------------------------------
471      do iq=1,nqtot
472        call write_archive(nid,ntime,tname(iq),'tracer','kg/kg',
473     &         3,q(1,1,iq))
474      end do
475c-----------------------------------------------------------------------
476c Ecriture du champs  qsurf  ( qsurf[1,nqtot] )
477c-----------------------------------------------------------------------
478      do iq=1,nqtot
479        txt=trim(tname(iq))//"_surf"
480        call write_archive(nid,ntime,txt,'Tracer on surface',
481     &  'kg.m-2',2,qsurfS(1,iq))
482      enddo
483
484
485c-----------------------------------------------------------------------
486c Ecriture du champs  tsoil  ( Tg[1,10] )
487c-----------------------------------------------------------------------
488c "tsoil" Temperature au sol definie dans 10 couches dans le sol
489c   Les 10 couches sont lues comme 10 champs
490c  nommees Tg[1,10]
491
492c      do isoil=1,nsoilmx
493c       write(str2,'(i2.2)') isoil
494c       call write_archive(nid,ntime,'Tg'//str2,'Ground Temperature ',
495c     .   'K',2,tsoilS(1,isoil))
496c      enddo
497
498! Write soil temperatures tsoil
499      call write_archive(nid,ntime,'tsoil','Soil temperature',
500     &     'K',-3,tsoilS)
501
502! Write soil thermal inertia
503      call write_archive(nid,ntime,'inertiedat',
504     &     'Soil thermal inertia',
505     &     'J.s-1/2.m-2.K-1',-3,ithS)
506
507! Write (0D) volumetric heat capacity (stored in comsoil.h)
508!      call write_archive(nid,ntime,'volcapa',
509!     &     'Soil volumetric heat capacity',
510!     &     'J.m-3.K-1',0,volcapa)
511! Note: no need to write volcapa, it is stored in "controle" table
512
513c-----------------------------------------------------------------
514c Ecriture du champs  tankCH4
515c-----------------------------------------------------------------
516      call write_archive(nid,ntime,'tankCH4',
517     &         'Depth of surface methane tank','m',2,tankCH4S)
518
519c-----------------------------------------------------------------
520c Ecriture des champs upper_chemistry
521c NB : We assume a given order of the nkim=44 chemistry tracers !
522c ( H=1, H2=2 ..., C4N2=nkim=44) -> hardcoded in comchem_h.
523c-----------------------------------------------------------------
524
525      IF (kim) THEN
526        DO iq=1,nkim
527          call write_archive(nid,ntime,trim(cnames(iq))//'_up',
528     .                trim(cnames(iq))//' in upper atmosphere',
529     .                'mol/mol',4,ykim_upS(iq,:,:))
530        ENDDO
531      ENDIF
532
533c Fin
534c-----------------------------------------------------------------------
535      ierr=NF_CLOSE(nid)
536
537      write(*,*) "start2archive: All is well that ends well."
538
539      end
Note: See TracBrowser for help on using the repository browser.