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

Last change on this file since 2322 was 2263, checked in by jnaar, 5 years ago

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