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

Last change on this file since 1758 was 1711, checked in by mvals, 8 years ago

Mars GCM

  • implementation of a sub-grid water cloud fraction scheme (by A. Pottier): If CLFvarying is set to true in callphys.def, then the sub-grid cloud fraction routine is applied in watercloud.F and aeropacity.F.
  • accordingly modified files are: watercloud.F, aeropacity.F, callradite.F, conf_phys.F, phyetat0.F, phyredem.F90, physiq_mod.F, testphys1d.F, callkeys.h, newtsart.F, start2archive.F, lect_start_archive.F
  • added file: tcondwater.F90, used by watercloud.F to calculate the condensation temperature of water
  • watercloud.F, aeropacity.F, callradite.F are converted to module files as watercloud_mod.F, aeropacity_mod.F, callradite_mod.F

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