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

Last change on this file since 2943 was 2943, checked in by llange, 20 months ago

Mars PCM
Following r-2942: Fix a bug in newstart when rewriting inertiesoil. Inertiesoil is now also managed in startarchive.
When using startarchive or newstart, inertiesoil is set to inertiedat.
LL

File size: 19.6 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
23      use surfdat_h, only: ini_surfdat_h, qsurf
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! AD: SSO parameters
39      USE surfdat_h, ONLY: phisfi, albedodat, z0, z0_default,
40     &    zmea, zstd, zsig, zgam, zthe, hmons, summit, base
41      implicit none
42
43      include "dimensions.h"
44      integer, parameter :: ngridmx = (2+(jjm-1)*iim - 1/jjm)
45      include "paramet.h"
46      include "comdissip.h"
47      include "comgeom.h"
48      include "netcdf.inc"
49
50c-----------------------------------------------------------------------
51c   Declarations
52c-----------------------------------------------------------------------
53
54c variables dynamiques du GCM
55c -----------------------------
56      REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm) ! vents covariants
57      REAL teta(ip1jmp1,llm)                    ! temperature potentielle
58      REAL,ALLOCATABLE :: q(:,:,:)   ! champs advectes
59      REAL pks(ip1jmp1)                      ! exner (f pour filtre)
60      REAL pk(ip1jmp1,llm)
61      REAL pkf(ip1jmp1,llm)
62      REAL phis(ip1jmp1)                     ! geopotentiel au sol
63      REAL masse(ip1jmp1,llm)                ! masse de l'atmosphere
64      REAL ps(ip1jmp1)                       ! pression au sol
65      REAL p3d(iip1, jjp1, llm+1)            ! pression aux interfaces
66     
67c Variable Physiques (grille physique)
68c ------------------------------------
69      REAL,ALLOCATABLE :: tsurf(:,:)        ! Surface temperature
70      REAL,ALLOCATABLE :: tsoil(:,:,:) ! Soil temperature
71      REAL,ALLOCATABLE :: inertiesoil(:,:,:) ! Soil thermal inertia (!= inertiedat which is for present day climate)
72      REAL,ALLOCATABLE :: watercap(:,:)        ! h2o ice layer
73      REAL tauscaling(ngridmx) ! dust conversion factor
74      REAL totcloudfrac(ngridmx) ! sub-grid cloud fraction
75      REAL q2(ngridmx,llm+1)
76      REAL,ALLOCATABLE :: emis(:,:)
77      REAL,ALLOCATABLE :: albedo(:,:,:)
78      REAL wstar(ngridmx)
79      DOUBLE PRECISION mem_Nccn_co2(ngridmx,llm)
80      DOUBLE PRECISION mem_Mccn_co2(ngridmx,llm)
81      DOUBLE PRECISION mem_Mh2o_co2(ngridmx,llm)
82      INTEGER start,length
83      PARAMETER (length = 100)
84      REAL tab_cntrl_fi(length) ! tableau des parametres de startfi
85      INTEGER*4 day_ini_fi
86
87c Variable naturelle / grille scalaire
88c ------------------------------------
89      REAL T(ip1jmp1,llm),us(ip1jmp1,llm),vs(ip1jmp1,llm)
90      REAL,ALLOCATABLE :: tsurfS(:,:)
91      REAL,ALLOCATABLE :: tsoilS(:,:,:)
92      REAL,ALLOCATABLE :: inertiesoilS(:,:,:)! Variable Soil Thermal Inertia  (obtained from PEM)
93      REAL ithS(ip1jmp1,nsoilmx) ! Soil Thermal Inertia for inertie dat (present day climate)
94      REAL,ALLOCATABLE :: watercapS(:,:)
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,
180     &                  1)
181
182      fichnom = 'startfi.nc'
183      Lmodif=0
184
185      allocate(tsurf(ngridmx,nslope))
186      allocate(tsoil(ngridmx,nsoilmx,nslope))
187      allocate(inertiesoil(ngridmx,nsoilmx,nslope))
188      allocate(watercap(ngridmx,nslope))
189      allocate(emis(ngridmx,nslope))
190      allocate(albedo(ngridmx,2,nslope))
191
192      allocate(qsurfS(ip1jmp1,nqtot,nslope))
193      allocate(tsurfS(ip1jmp1,nslope))
194      allocate(tsoilS(ip1jmp1,nsoilmx,nslope))
195      allocate(inertiesoilS(ip1jmp1,nsoilmx,nslope))
196      allocate(watercapS(ip1jmp1,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,def_slope,
204     &      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(nsoilmx,ngridmx,iip1,jjp1,tsoil(:,:,islope),
302     &    tsoilS(:,:,islope))
303      call gr_fi_dyn(nsoilmx,ngridmx,iip1,jjp1,inertiesoil(:,:,islope),
304     &    inertiesoil(:,:,islope))
305      ! Note: thermal inertia "inertiedat" is in comsoil.h
306      call gr_fi_dyn(1,ngridmx,iip1,jjp1,emis(:,islope),
307     &     emisS(:,islope))
308      call gr_fi_dyn(1,ngridmx,iip1,jjp1,albedo(:,1,islope),
309     &   albedoS(:,islope))
310      call gr_fi_dyn(nqtot,ngridmx,iip1,jjp1,qsurf(:,:,islope),
311     &   qsurfS(:,:,islope))
312      call gr_fi_dyn(1,ngridmx,iip1,jjp1,subslope_dist(:,islope),
313     &    subslope_distS(:,islope))
314      enddo
315      call gr_fi_dyn(nsoilmx,ngridmx,iip1,jjp1,inertiedat,ithS)
316      call gr_fi_dyn(llm+1,ngridmx,iip1,jjp1,q2,q2S)
317      call gr_fi_dyn(1,ngridmx,iip1,jjp1,tauscaling,tauscalingS)
318      call gr_fi_dyn(1,ngridmx,iip1,jjp1,totcloudfrac,totcloudfracS)
319
320! AD: for SSO parameters
321      call gr_fi_dyn(1,ngridmx,iip1,jjp1,zmea,zmeaS)
322      call gr_fi_dyn(1,ngridmx,iip1,jjp1,zstd,zstdS)
323      call gr_fi_dyn(1,ngridmx,iip1,jjp1,zsig,zsigS)
324      call gr_fi_dyn(1,ngridmx,iip1,jjp1,zthe,ztheS)
325      call gr_fi_dyn(1,ngridmx,iip1,jjp1,zgam,zgamS)
326      call gr_fi_dyn(1,ngridmx,iip1,jjp1,albedodat,albedodatS)
327      call gr_fi_dyn(1,ngridmx,iip1,jjp1,z0,z0S)
328      call gr_fi_dyn(1,ngridmx,iip1,jjp1,hmons,hmonsS)
329      call gr_fi_dyn(1,ngridmx,iip1,jjp1,summit,summitS)
330      call gr_fi_dyn(1,ngridmx,iip1,jjp1,base,baseS)
331
332c=======================================================================
333c Info pour controler
334c=======================================================================
335
336      DO iq=1,nqtot
337        if (trim(tname(iq)) .eq. "co2") then
338           igcm_co2=iq
339        endif
340      enddo
341
342      ptotal =  0.
343      co2icetotal = 0.
344      DO j=1,jjp1
345         DO i=1,iim
346           DO islope=1,nslope
347           ptotal=ptotal+aire(i+(iim+1)*(j-1))*ps(i+(iim+1)*(j-1))/g
348           co2icetotal = co2icetotal +
349     &        qsurfS(i+(iim+1)*(j-1),igcm_co2,islope)*
350     &        aire(i+(iim+1)*(j-1))*
351     &    subslope_distS(i+(iim+1)*(j-1),islope)/
352     &    cos(pi*def_slope_mean(islope))
353           ENDDO
354         ENDDO
355      ENDDO
356      write(*,*)'Ancienne grille : masse de l''atm :',ptotal
357      write(*,*)'Ancienne grille : masse de la glace CO2 :',co2icetotal
358
359c-----------------------------------------------------------------------
360c Passage de "ptotal" et "co2icetotal" par tab_cntrl_fi
361c-----------------------------------------------------------------------
362
363      tab_cntrl_fi(49) = ptotal
364      tab_cntrl_fi(50) = co2icetotal
365
366c=======================================================================
367c Ecriture dans le fichier  "start_archive"
368c=======================================================================
369
370c-----------------------------------------------------------------------
371c Ouverture de "start_archive"
372c-----------------------------------------------------------------------
373
374      ierr = NF_OPEN ('start_archive.nc', NF_WRITE,nid)
375 
376c-----------------------------------------------------------------------
377c  si "start_archive" n'existe pas:
378c    1_ ouverture
379c    2_ creation de l'entete dynamique ("ini_archive")
380c-----------------------------------------------------------------------
381c ini_archive:
382c On met dans l'entete le tab_cntrl dynamique (1 a 16)
383c  On y ajoute les valeurs du tab_cntrl_fi (a partir de 51)
384c  En plus les deux valeurs ptotal et co2icetotal (99 et 100)
385c-----------------------------------------------------------------------
386
387      if (ierr.ne.NF_NOERR) then
388         write(*,*)'OK, Could not open file "start_archive.nc"'
389         write(*,*)'So let s create a new "start_archive"'
390         ierr = NF_CREATE('start_archive.nc',
391     &  IOR(NF_CLOBBER,NF_64BIT_OFFSET), nid)
392         call ini_archive(nid,day_ini,phis,ithS,tab_cntrl_fi,
393     &         def_slope,subslope_distS)
394      endif
395
396c-----------------------------------------------------------------------
397c Ecriture de la coordonnee temps (date en jours)
398c-----------------------------------------------------------------------
399
400      date = day_ini + hour_ini
401      ierr= NF_INQ_VARID(nid,"Time",varid)
402      ierr= NF_INQ_DIMID(nid,"Time",dimid)
403      ierr= NF_INQ_DIMLEN(nid,dimid,timelen)
404      ntime=timelen+1
405
406      write(*,*) "******************"
407      write(*,*) "ntime",ntime
408      write(*,*) "******************"
409#ifdef NC_DOUBLE
410      ierr= NF_PUT_VARA_DOUBLE(nid,varid,ntime,1,date)
411#else
412      ierr= NF_PUT_VARA_REAL(nid,varid,ntime,1,date)
413#endif
414      if (ierr.ne.NF_NOERR) then
415         write(*,*) "time matter ",NF_STRERROR(ierr)
416         stop
417      endif
418
419c-----------------------------------------------------------------------
420c Ecriture des champs  (co2ice,emis,ps,Tsurf,T,u,v,q2,q,qsurf)
421c-----------------------------------------------------------------------
422c ATTENTION: q2 a une couche de plus!!!!
423c    Pour creer un fichier netcdf lisible par grads,
424c    On passe donc une des couches de q2 a part
425c    comme une variable 2D (la couche au sol: "q2surf")
426c    Les lmm autres couches sont nommees "q2atm" (3D)
427c-----------------------------------------------------------------------
428
429      call write_archive(nid,ntime,'watercap','couche de glace h2o',
430     &  'kg/m2',2,watercapS)
431      call write_archive(nid,ntime,'tauscaling',
432     &  'dust conversion factor',' ',2,tauscalingS)
433      call write_archive(nid,ntime,'totcloudfrac',
434     &  'sub grid cloud fraction',' ',2,totcloudfracS)
435      call write_archive(nid,ntime,'emis','grd emis',' ',2,emisS)
436      call write_archive(nid,ntime,'albedo','surface albedo',' ',
437     &                             2,albedoS)
438      call write_archive(nid,ntime,'ps','Psurf','Pa',2,ps)
439      call write_archive(nid,ntime,'tsurf','surf T','K',2,tsurfS)
440      call write_archive(nid,ntime,'temp','temperature','K',3,t)
441      call write_archive(nid,ntime,'u','Vent zonal','m.s-1',3,us)
442      call write_archive(nid,ntime,'v','Vent merid','m.s-1',3,vs)
443      call write_archive(nid,ntime,'q2surf','wind variance','m2.s-2',2,
444     .              q2S)
445      call write_archive(nid,ntime,'q2atm','wind variance','m2.s-2',3,
446     .              q2S(1,2))
447! AD: SSO parameters
448      call write_archive(nid,ntime,'ZMEA','zmea',' ',2,zmeaS)
449      call write_archive(nid,ntime,'ZSTD','zstd',' ',2,zstdS)
450      call write_archive(nid,ntime,'ZSIG','zsig',' ',2,zsigS)
451      call write_archive(nid,ntime,'ZTHE','zthe',' ',2,ztheS)
452      call write_archive(nid,ntime,'ZGAM','zgam',' ',2,zgamS)
453      call write_archive(nid,ntime,'albedodat','albedodat',
454     &                             ' ',2,albedodatS)
455      call write_archive(nid,ntime,'z0','z0',' ',2,z0S)
456      call write_archive(nid,ntime,'summit','summit',
457     &                             ' ',2,summitS)
458      call write_archive(nid,ntime,'hmons','hmons',' ',2,hmonsS)
459      call write_archive(nid,ntime,'base','base',' ',2,baseS)
460c-----------------------------------------------------------------------
461c Ecriture du champs  q  ( q[1,nqtot] )
462c-----------------------------------------------------------------------
463      do iq=1,nqtot
464c       write(str2,'(i2.2)') iq
465c        call write_archive(nid,ntime,'q'//str2,'tracer','kg/kg',
466c     .         3,q(1,1,iq))
467        call write_archive(nid,ntime,tname(iq),'tracer','kg/kg',
468     &         3,q(1,1,iq))
469      end do
470c-----------------------------------------------------------------------
471c Ecriture du champs  qsurf  ( qsurf[1,nqtot] )
472c-----------------------------------------------------------------------
473      do iq=1,nqtot
474c       write(str2,'(i2.2)') iq
475c       call write_archive(nid,ntime,'qsurf'//str2,'Tracer on surface',
476c     $  'kg.m-2',2,qsurfS(1,iq))
477        txt=trim(tname(iq))//"_surf"
478        call write_archive(nid,ntime,txt,'Tracer on surface',
479     &  'kg.m-2',2,qsurfS(:,iq,:))
480      enddo
481
482
483c-----------------------------------------------------------------------
484c Ecriture du champs  tsoil  ( Tg[1,10] )
485c-----------------------------------------------------------------------
486c "tsoil" Temperature au sol definie dans 10 couches dans le sol
487c   Les 10 couches sont lues comme 10 champs
488c  nommees Tg[1,10]
489
490c      do isoil=1,nsoilmx
491c       write(str2,'(i2.2)') isoil
492c       call write_archive(nid,ntime,'Tg'//str2,'Ground Temperature ',
493c     .   'K',2,tsoilS(1,isoil))
494c      enddo
495
496! Write soil temperatures tsoil
497      call write_archive(nid,ntime,'tsoil','Soil temperature',
498     &     'K',-3,tsoilS(:,:,:))
499! Write soil thermal inertia
500      call write_archive(nid,ntime,'inertiesoil','Soil TI',
501     &     'J.s-1/2.m-2.K-1',-3,inertiesoilS(:,:,:))
502! Write soil thermal inertia for current climate
503      call write_archive(nid,ntime,'inertiedat',
504     &     'Soil thermal inertia',
505     &     'J.s-1/2.m-2.K-1',-3,ithS)
506
507! Write (0D) volumetric heat capacity (stored in comsoil.h)
508!      call write_archive(nid,ntime,'volcapa',
509!     &     'Soil volumetric heat capacity',
510!     &     'J.m-3.K-1',0,volcapa)
511! Note: no need to write volcapa, it is stored in "controle" table
512
513      ierr=NF_CLOSE(nid)
514c-----------------------------------------------------------------------
515c Fin
516c-----------------------------------------------------------------------
517
518      write(*,*) "startarchive: all is well that ends well"
519     
520      end
Note: See TracBrowser for help on using the repository browser.