source: trunk/LMDZ.MARS/libf/dynlonlat_phylonlat/phymars/start2archive.F @ 1422

Last change on this file since 1422 was 1422, checked in by milmd, 10 years ago

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