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

Last change on this file since 3094 was 2999, checked in by llange, 18 months ago

Mars PCM
Include perenial_co2ice (equivalent of watercap) to distinguich between CO2 frost and perenial CO2 ice for paleoclimate studies.
When no frost is present and we dig into perenial ice, the surface albedo is changed. The albedo for seasonal ice is set to 0.65, and the perenial ice albedo can be fixed in the callphys.def. I recommand values between 0.8 and 0.9.
To use this, paleoclimate must be set to True and TESalbedo to false in the callphys.def. Else, it runs as usual with TES albedo
LL

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