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

Last change on this file since 3316 was 3316, checked in by jbclement, 9 months ago

Mars PCM:
Reversion of r3305 and r3307.
JBC

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