source: trunk/LMDZ.MARS/libf/dynphy_lonlat/phymars/start2archive_SSO.F @ 2959

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