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

Last change on this file since 3302 was 3216, checked in by emillour, 11 months ago

Mars PCM:
Fix start2archive; routine surfini is now in module surfini_mod.
EM

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