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

Last change on this file since 2997 was 2959, checked in by romain.vande, 20 months ago

Mars PCM :
Correct start2archive to write watercaptag correctly.
Watercaptag will be set to false and correctly handle by the PCM in the case where we change resolution.
+ Correct inertiesoil writting in start2archive
RV

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