source: trunk/LMDZ.MARS/libf/dyn3d/start2archive.F @ 1047

Last change on this file since 1047 was 1047, checked in by emillour, 11 years ago

Mars GCM:

  • IMPORTANT CHANGE: Removed all reference/use of ngridmx (dimphys.h) in routines (necessary prerequisite to using parallel dynamics); in most cases this just means adding 'ngrid' as routine argument, and making local saved variables allocatable (and allocated at first call). In the process, had to convert many *.h files to equivalent modules: yomaer.h => yomaer_h.F90 , surfdat.h => surfdat_h.F90 , comsaison.h => comsaison_h.F90 , yomlw.h => yomlw_h.F90 , comdiurn.h => comdiurn_h.F90 , dimradmars.h => dimradmars_mod.F90 , comgeomfi.h => comgeomfi_h.F90, comsoil.h => comsoil_h.F90 , slope.h => slope_mod.F90
  • Also updated EOF routines, everything is now in eofdump_mod.F90
  • Removed unused routine lectfux.F (in dyn3d)

EM
File size: 14.3 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: iniadvtrac, nqtot, tnom
22      use comsoil_h, only: nsoilmx, inertiedat
23      use surfdat_h, only: ini_surfdat_h
24      use comsoil_h, only: ini_comsoil_h
25      implicit none
26
27#include "dimensions.h"
28      integer, parameter :: ngridmx = (2+(jjm-1)*iim - 1/jjm)
29#include "paramet.h"
30#include "comconst.h"
31#include "comdissip.h"
32#include "comvert.h"
33#include "comgeom.h"
34#include "logic.h"
35#include "temps.h"
36#include "control.h"
37#include "ener.h"
38#include "description.h"
39
40!#include "dimphys.h"
41!#include "comsoil.h"
42!#include"advtrac.h"
43#include "netcdf.inc"
44
45c-----------------------------------------------------------------------
46c   Declarations
47c-----------------------------------------------------------------------
48
49c variables dynamiques du GCM
50c -----------------------------
51      REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm) ! vents covariants
52      REAL teta(ip1jmp1,llm)                    ! temperature potentielle
53      REAL,ALLOCATABLE :: q(:,:,:)   ! champs advectes
54      REAL pks(ip1jmp1)                      ! exner (f pour filtre)
55      REAL pk(ip1jmp1,llm)
56      REAL pkf(ip1jmp1,llm)
57      REAL beta(iip1,jjp1,llm)
58      REAL phis(ip1jmp1)                     ! geopotentiel au sol
59      REAL masse(ip1jmp1,llm)                ! masse de l'atmosphere
60      REAL ps(ip1jmp1)                       ! pression au sol
61      REAL p3d(iip1, jjp1, llm+1)            ! pression aux interfaces
62     
63c Variable Physiques (grille physique)
64c ------------------------------------
65      REAL tsurf(ngridmx)       ! Surface temperature
66      REAL tsoil(ngridmx,nsoilmx) ! Soil temperature
67      REAL co2ice(ngridmx)      ! CO2 ice layer
68      REAL q2(ngridmx,llm+1)
69      REAL,ALLOCATABLE :: qsurf(:,:)
70      REAL emis(ngridmx)
71      INTEGER start,length
72      PARAMETER (length = 100)
73      REAL tab_cntrl_fi(length) ! tableau des parametres de startfi
74      INTEGER*4 day_ini_fi
75
76c Variable naturelle / grille scalaire
77c ------------------------------------
78      REAL T(ip1jmp1,llm),us(ip1jmp1,llm),vs(ip1jmp1,llm)
79      REAL tsurfS(ip1jmp1)
80      REAL tsoilS(ip1jmp1,nsoilmx)
81      REAL ithS(ip1jmp1,nsoilmx) ! Soil Thermal Inertia
82      REAL co2iceS(ip1jmp1)
83      REAL q2S(ip1jmp1,llm+1)
84      REAL,ALLOCATABLE :: qsurfS(:,:)
85      REAL emisS(ip1jmp1)
86
87c Variables intermediaires : vent naturel, mais pas coord scalaire
88c----------------------------------------------------------------
89      REAL vn(ip1jm,llm),un(ip1jmp1,llm)
90
91c Autres  variables
92c -----------------
93      LOGICAL startdrs
94      INTEGER Lmodif
95
96      REAL ptotal, co2icetotal
97      REAL timedyn,timefi !fraction du jour dans start, startfi
98      REAL date
99
100      CHARACTER*2 str2
101      CHARACTER*80 fichier
102      data  fichier /'startfi'/
103
104      INTEGER ij, l,i,j,isoil,iq
105      character*80      fichnom
106      integer :: ierr,ntime
107      integer :: nq,numvanle
108      character(len=30) :: txt ! to store some text
109
110c Netcdf
111c-------
112      integer varid,dimid,timelen
113      INTEGER nid,nid1
114
115c-----------------------------------------------------------------------
116c   Initialisations
117c-----------------------------------------------------------------------
118
119      CALL defrun_new(99, .TRUE. )
120      grireg   = .TRUE.
121
122c=======================================================================
123c Lecture des donnees
124c=======================================================================
125! Load tracer number and names:
126      call iniadvtrac(nqtot,numvanle)
127
128! allocate arrays:
129      allocate(q(ip1jmp1,llm,nqtot))
130      allocate(qsurf(ngridmx,nqtot))
131      allocate(qsurfS(ip1jmp1,nqtot))
132      call ini_surfdat_h(ngridmx)
133      call ini_comsoil_h(ngridmx)
134     
135
136      fichnom = 'start.nc'
137      CALL dynetat0(fichnom,nqtot,vcov,ucov,teta,q,masse,
138     .       ps,phis,timedyn)
139
140
141      fichnom = 'startfi.nc'
142      Lmodif=0
143
144      CALL phyetat0 (fichnom,0,Lmodif,nsoilmx,ngridmx,llm,nqtot,
145     &      day_ini_fi,timefi,tsurf,tsoil,emis,q2,qsurf,co2ice)
146
147       ierr = NF_OPEN (fichnom, NF_NOWRITE,nid1)
148       IF (ierr.NE.NF_NOERR) THEN
149         write(6,*)' Pb d''ouverture du fichier'//fichnom
150        CALL ABORT
151       ENDIF
152                                               
153      ierr = NF_INQ_VARID (nid1, "controle", varid)
154      IF (ierr .NE. NF_NOERR) THEN
155       PRINT*, "start2archive: Le champ <controle> est absent"
156       CALL abort
157      ENDIF
158#ifdef NC_DOUBLE
159       ierr = NF_GET_VAR_DOUBLE(nid1, varid, tab_cntrl_fi)
160#else
161      ierr = NF_GET_VAR_REAL(nid1, varid, tab_cntrl_fi)
162#endif
163       IF (ierr .NE. NF_NOERR) THEN
164          PRINT*, "start2archive: Lecture echoue pour <controle>"
165          CALL abort
166       ENDIF
167
168      ierr = NF_CLOSE(nid1)
169
170c-----------------------------------------------------------------------
171c Controle de la synchro
172c-----------------------------------------------------------------------
173!mars a voir      if ((day_ini_fi.ne.day_ini).or.(abs(timefi-timedyn).gt.1.e-10))
174      if ((day_ini_fi.ne.day_ini))
175     &  stop ' Probleme de Synchro entre start et startfi !!!'
176
177
178c *****************************************************************
179c    Option : Reinitialisation des dates dans la premieres annees :
180       do while (day_ini.ge.669)
181          day_ini=day_ini-669
182       enddo
183c *****************************************************************
184
185c-----------------------------------------------------------------------
186c   Initialisations
187c-----------------------------------------------------------------------
188
189      CALL defrun_new(99, .FALSE. )
190      call iniconst
191      call inigeom
192      call inifilr
193      CALL pression(ip1jmp1, ap, bp, ps, p3d)
194      call exner_hyb(ip1jmp1, ps, p3d, beta, pks, pk, pkf)
195
196c=======================================================================
197c Transformation EN VARIABLE NATURELLE / GRILLE SCALAIRE si necessaire
198c=======================================================================
199c  Les variables modeles dependent de la resolution. Il faut donc
200c  eliminer les facteurs responsables de cette dependance
201c  (pour utiliser newstart)
202c=======================================================================
203
204c-----------------------------------------------------------------------
205c Vent   (depend de la resolution horizontale)
206c-----------------------------------------------------------------------
207c
208c ucov --> un  et  vcov --> vn
209c un --> us  et   vn --> vs
210c
211c-----------------------------------------------------------------------
212
213      call covnat(llm,ucov, vcov, un, vn)
214      call wind_scal(un,vn,us,vs)
215
216c-----------------------------------------------------------------------
217c Temperature  (depend de la resolution verticale => de "sigma.def")
218c-----------------------------------------------------------------------
219c
220c h --> T
221c
222c-----------------------------------------------------------------------
223
224      DO l=1,llm
225         DO ij=1,ip1jmp1
226            T(ij,l)=teta(ij,l)*pk(ij,l)/cpp !mars deduit de l'equation dans newstart
227         ENDDO
228      ENDDO
229
230c-----------------------------------------------------------------------
231c Variable physique
232c-----------------------------------------------------------------------
233c
234c tsurf --> tsurfS
235c co2ice --> co2iceS
236c tsoil --> tsoilS
237c emis --> emisS
238c q2 --> q2S
239c qsurf --> qsurfS
240c
241c-----------------------------------------------------------------------
242
243      call gr_fi_dyn(1,ngridmx,iip1,jjp1,tsurf,tsurfS)
244      call gr_fi_dyn(1,ngridmx,iip1,jjp1,co2ice,co2iceS)
245      call gr_fi_dyn(nsoilmx,ngridmx,iip1,jjp1,tsoil,tsoilS)
246      ! Note: thermal inertia "inertiedat" is in comsoil.h
247      call gr_fi_dyn(nsoilmx,ngridmx,iip1,jjp1,inertiedat,ithS)
248      call gr_fi_dyn(1,ngridmx,iip1,jjp1,emis,emisS)
249      call gr_fi_dyn(llm+1,ngridmx,iip1,jjp1,q2,q2S)
250      call gr_fi_dyn(nqtot,ngridmx,iip1,jjp1,qsurf,qsurfS)
251
252c=======================================================================
253c Info pour controler
254c=======================================================================
255
256      ptotal =  0.
257      co2icetotal = 0.
258      DO j=1,jjp1
259         DO i=1,iim
260           ptotal=ptotal+aire(i+(iim+1)*(j-1))*ps(i+(iim+1)*(j-1))/g
261           co2icetotal = co2icetotal +
262     &            co2iceS(i+(iim+1)*(j-1))*aire(i+(iim+1)*(j-1))
263         ENDDO
264      ENDDO
265      write(*,*)'Ancienne grille : masse de l''atm :',ptotal
266      write(*,*)'Ancienne grille : masse de la glace CO2 :',co2icetotal
267
268c-----------------------------------------------------------------------
269c Passage de "ptotal" et "co2icetotal" par tab_cntrl_fi
270c-----------------------------------------------------------------------
271
272      tab_cntrl_fi(49) = ptotal
273      tab_cntrl_fi(50) = co2icetotal
274
275c=======================================================================
276c Ecriture dans le fichier  "start_archive"
277c=======================================================================
278
279c-----------------------------------------------------------------------
280c Ouverture de "start_archive"
281c-----------------------------------------------------------------------
282
283      ierr = NF_OPEN ('start_archive.nc', NF_WRITE,nid)
284 
285c-----------------------------------------------------------------------
286c  si "start_archive" n'existe pas:
287c    1_ ouverture
288c    2_ creation de l'entete dynamique ("ini_archive")
289c-----------------------------------------------------------------------
290c ini_archive:
291c On met dans l'entete le tab_cntrl dynamique (1 a 16)
292c  On y ajoute les valeurs du tab_cntrl_fi (a partir de 51)
293c  En plus les deux valeurs ptotal et co2icetotal (99 et 100)
294c-----------------------------------------------------------------------
295
296      if (ierr.ne.NF_NOERR) then
297         write(*,*)'OK, Could not open file "start_archive.nc"'
298         write(*,*)'So let s create a new "start_archive"'
299         ierr = NF_CREATE('start_archive.nc',
300     &  IOR(NF_CLOBBER,NF_64BIT_OFFSET), nid)
301         call ini_archive(nid,day_ini,phis,ithS,tab_cntrl_fi)
302      endif
303
304c-----------------------------------------------------------------------
305c Ecriture de la coordonnee temps (date en jours)
306c-----------------------------------------------------------------------
307
308      date = day_ini + hour_ini
309      ierr= NF_INQ_VARID(nid,"Time",varid)
310      ierr= NF_INQ_DIMID(nid,"Time",dimid)
311      ierr= NF_INQ_DIMLEN(nid,dimid,timelen)
312      ntime=timelen+1
313
314      write(*,*) "******************"
315      write(*,*) "ntime",ntime
316      write(*,*) "******************"
317#ifdef NC_DOUBLE
318      ierr= NF_PUT_VARA_DOUBLE(nid,varid,ntime,1,date)
319#else
320      ierr= NF_PUT_VARA_REAL(nid,varid,ntime,1,date)
321#endif
322      if (ierr.ne.NF_NOERR) then
323         write(*,*) "time matter ",NF_STRERROR(ierr)
324         stop
325      endif
326
327c-----------------------------------------------------------------------
328c Ecriture des champs  (co2ice,emis,ps,Tsurf,T,u,v,q2,q,qsurf)
329c-----------------------------------------------------------------------
330c ATTENTION: q2 a une couche de plus!!!!
331c    Pour creer un fichier netcdf lisible par grads,
332c    On passe donc une des couches de q2 a part
333c    comme une variable 2D (la couche au sol: "q2surf")
334c    Les lmm autres couches sont nommees "q2atm" (3D)
335c-----------------------------------------------------------------------
336
337      call write_archive(nid,ntime,'co2ice','couche de glace co2',
338     &  'kg/m2',2,co2iceS)
339      call write_archive(nid,ntime,'emis','grd emis',' ',2,emisS)
340      call write_archive(nid,ntime,'ps','Psurf','Pa',2,ps)
341      call write_archive(nid,ntime,'tsurf','surf T','K',2,tsurfS)
342      call write_archive(nid,ntime,'temp','temperature','K',3,t)
343      call write_archive(nid,ntime,'u','Vent zonal','m.s-1',3,us)
344      call write_archive(nid,ntime,'v','Vent merid','m.s-1',3,vs)
345      call write_archive(nid,ntime,'q2surf','wind variance','m2.s-2',2,
346     .              q2S)
347      call write_archive(nid,ntime,'q2atm','wind variance','m2.s-2',3,
348     .              q2S(1,2))
349
350c-----------------------------------------------------------------------
351c Ecriture du champs  q  ( q[1,nqtot] )
352c-----------------------------------------------------------------------
353      do iq=1,nqtot
354c       write(str2,'(i2.2)') iq
355c        call write_archive(nid,ntime,'q'//str2,'tracer','kg/kg',
356c     .         3,q(1,1,iq))
357        call write_archive(nid,ntime,tnom(iq),'tracer','kg/kg',
358     &         3,q(1,1,iq))
359      end do
360c-----------------------------------------------------------------------
361c Ecriture du champs  qsurf  ( qsurf[1,nqtot] )
362c-----------------------------------------------------------------------
363      do iq=1,nqtot
364c       write(str2,'(i2.2)') iq
365c       call write_archive(nid,ntime,'qsurf'//str2,'Tracer on surface',
366c     $  'kg.m-2',2,qsurfS(1,iq))
367        txt=trim(tnom(iq))//"_surf"
368        call write_archive(nid,ntime,txt,'Tracer on surface',
369     &  'kg.m-2',2,qsurfS(1,iq))
370      enddo
371
372
373c-----------------------------------------------------------------------
374c Ecriture du champs  tsoil  ( Tg[1,10] )
375c-----------------------------------------------------------------------
376c "tsoil" Temperature au sol definie dans 10 couches dans le sol
377c   Les 10 couches sont lues comme 10 champs
378c  nommees Tg[1,10]
379
380c      do isoil=1,nsoilmx
381c       write(str2,'(i2.2)') isoil
382c       call write_archive(nid,ntime,'Tg'//str2,'Ground Temperature ',
383c     .   'K',2,tsoilS(1,isoil))
384c      enddo
385
386! Write soil temperatures tsoil
387      call write_archive(nid,ntime,'tsoil','Soil temperature',
388     &     'K',-3,tsoilS)
389
390! Write soil thermal inertia
391      call write_archive(nid,ntime,'inertiedat',
392     &     'Soil thermal inertia',
393     &     'J.s-1/2.m-2.K-1',-3,ithS)
394
395! Write (0D) volumetric heat capacity (stored in comsoil.h)
396!      call write_archive(nid,ntime,'volcapa',
397!     &     'Soil volumetric heat capacity',
398!     &     'J.m-3.K-1',0,volcapa)
399! Note: no need to write volcapa, it is stored in "controle" table
400
401      ierr=NF_CLOSE(nid)
402c-----------------------------------------------------------------------
403c Fin
404c-----------------------------------------------------------------------
405
406      end
Note: See TracBrowser for help on using the repository browser.