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

Last change on this file since 1647 was 1647, checked in by jvatant, 8 years ago

+ Major clean of the new LMDZ.TITAN from too-generic options and routines (water, co2, ocean, surface type ...)
+ From this revision LMDZ.TITAN begins to be really separated from LMDZ.GENERIC
+ Partial desactivation of aerosols, only the dummy case is still enabled to keep the code running ( new aerosol routines to come in followings commits )

JVO

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