source: trunk/LMDZ.MARS/libf/dynphy_lonlat/phymars/start2archive.F @ 1711

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