source: codes/icosagcm/branches/SATURN_DYNAMICO/LMDZ.COMMON/libf/phystd/start2archive.F @ 224

Last change on this file since 224 was 222, checked in by ymipsl, 10 years ago

Creating temporary dynamico/lmdz/saturn branche

YM

File size: 18.1 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, tname
22      USE comsoil_h
23      USE comgeomfi_h, ONLY: lati, long, area
24!      use control_mod
25      use comgeomphy, only: initcomgeomphy
26      use slab_ice_h, only: noceanmx
27! to use  'getin'
28      USE ioipsl_getincom
29
30      implicit none
31
32#include "dimensions.h"
33#include "paramet.h"
34#include "comconst.h"
35#include "comdissip.h"
36#include "comvert.h"
37#include "comgeom.h"
38#include "logic.h"
39#include "temps.h"
40!#include "control.h"
41#include "ener.h"
42
43#include "dimphys.h"
44#include "planete.h"
45!#include"advtrac.h"
46#include "netcdf.inc"
47#include "callkeys.h"
48c-----------------------------------------------------------------------
49c   Declarations
50c-----------------------------------------------------------------------
51
52c variables dynamiques du GCM
53c -----------------------------
54      REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm) ! vents covariants
55      REAL teta(ip1jmp1,llm)                    ! temperature potentielle
56      REAL,ALLOCATABLE :: q(:,:,:)   ! champs advectes
57      REAL pks(ip1jmp1)                      ! exner (f pour filtre)
58      REAL pk(ip1jmp1,llm)
59      REAL pkf(ip1jmp1,llm)
60      REAL beta(iip1,jjp1,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 tsurf(ngridmx)       ! Surface temperature
69      REAL tsoil(ngridmx,nsoilmx) ! Soil temperature
70      REAL co2ice(ngridmx)      ! CO2 ice layer
71      REAL q2(ngridmx,nlayermx+1)
72      REAL,ALLOCATABLE :: qsurf(:,:)
73      REAL emis(ngridmx)
74      INTEGER start,length
75      PARAMETER (length = 100)
76      REAL tab_cntrl_fi(length) ! tableau des parametres de startfi
77      REAL tab_cntrl_dyn(length) ! tableau des parametres de start
78      INTEGER*4 day_ini_fi
79
80!     added by FF for cloud fraction setup
81      REAL hice(ngridmx)
82      REAL cloudfrac(ngridmx,nlayermx),totalcloudfrac(ngridmx)
83
84!     added by BC for slab ocean
85      REAL rnat(ngridmx),pctsrf_sic(ngridmx),sea_ice(ngridmx)
86      REAL tslab(ngridmx,noceanmx),tsea_ice(ngridmx)
87
88
89c Variable naturelle / grille scalaire
90c ------------------------------------
91      REAL T(ip1jmp1,llm),us(ip1jmp1,llm),vs(ip1jmp1,llm)
92      REAL tsurfS(ip1jmp1)
93      REAL tsoilS(ip1jmp1,nsoilmx)
94      REAL ithS(ip1jmp1,nsoilmx) ! Soil Thermal Inertia
95      REAL co2iceS(ip1jmp1)
96      REAL q2S(ip1jmp1,llm+1)
97      REAL,ALLOCATABLE :: qsurfS(:,:)
98      REAL emisS(ip1jmp1)
99
100!     added by FF for cloud fraction setup
101      REAL hiceS(ip1jmp1)
102      REAL cloudfracS(ip1jmp1,llm),totalcloudfracS(ip1jmp1)
103
104!     added by BC for slab ocean
105      REAL rnatS(ip1jmp1),pctsrf_sicS(ip1jmp1),sea_iceS(ip1jmp1)
106      REAL tslabS(ip1jmp1,noceanmx),tsea_iceS(ip1jmp1)
107
108
109c Variables intermediaires : vent naturel, mais pas coord scalaire
110c----------------------------------------------------------------
111      REAL vn(ip1jm,llm),un(ip1jmp1,llm)
112
113c Autres  variables
114c -----------------
115      LOGICAL startdrs
116      INTEGER Lmodif
117
118      REAL ptotal, co2icetotal
119      REAL timedyn,timefi !fraction du jour dans start, startfi
120      REAL date
121
122      CHARACTER*2 str2
123      CHARACTER*80 fichier 
124      data  fichier /'startfi'/
125
126      INTEGER ij, l,i,j,isoil,iq
127      character*80      fichnom
128      integer :: ierr,ntime
129      integer :: nq,numvanle
130      character(len=30) :: txt ! to store some text
131
132c Netcdf
133c-------
134      integer varid,dimid,timelen 
135      INTEGER nid,nid1
136
137c-----------------------------------------------------------------------
138c   Initialisations 
139c-----------------------------------------------------------------------
140
141      CALL defrun_new(99, .TRUE. )
142      grireg   = .TRUE.
143
144! initialize "serial/parallel" related stuff
145      CALL init_phys_lmdz(iim,jjp1,llm,1,(/(jjm-1)*iim+2/))
146      call initcomgeomphy
147
148      ! ALLOCATE ARRAYS IN comgeomfi_h (usually done in inifis)
149      ! this must be here for start2archive to work
150      IF (.not. ALLOCATED(lati)) ALLOCATE(lati(ngridmx))
151      IF (.not. ALLOCATED(long)) ALLOCATE(long(ngridmx))
152      IF (.not. ALLOCATED(area)) ALLOCATE(area(ngridmx))
153
154c=======================================================================
155c Lecture des donnees
156c=======================================================================
157! Load tracer number and names:
158      call iniadvtrac(nqtot,numvanle)
159
160! allocate arrays:
161      allocate(q(ip1jmp1,llm,nqtot))
162      allocate(qsurf(ngridmx,nqtot))
163      allocate(qsurfS(ip1jmp1,nqtot))
164! other array allocations:
165      call ini_comsoil_h(ngridmx)
166
167      fichnom = 'start.nc'
168      CALL dynetat0(fichnom,nqtot,vcov,ucov,teta,q,masse,
169     .       ps,phis,timedyn)
170
171! load 'controle' array from dynamics start file
172
173       ierr = NF_OPEN (fichnom, NF_NOWRITE,nid1)
174       IF (ierr.NE.NF_NOERR) THEN
175         write(6,*)' Pb d''ouverture du fichier'//trim(fichnom)
176        CALL ABORT
177       ENDIF
178                                               
179      ierr = NF_INQ_VARID (nid1, "controle", varid)
180      IF (ierr .NE. NF_NOERR) THEN
181       PRINT*, "start2archive: Le champ <controle> est absent"
182       CALL abort
183      ENDIF
184#ifdef NC_DOUBLE
185       ierr = NF_GET_VAR_DOUBLE(nid1, varid, tab_cntrl_dyn)
186#else
187      ierr = NF_GET_VAR_REAL(nid1, varid, tab_cntrl_dyn)
188#endif
189       IF (ierr .NE. NF_NOERR) THEN
190          PRINT*, "start2archive: Lecture echoue pour <controle>"
191          CALL abort
192       ENDIF
193
194      ierr = NF_CLOSE(nid1)
195     
196
197      fichnom = 'startfi.nc'
198      Lmodif=0
199
200
201      CALL phyetat0 (ngridmx,fichnom,0,Lmodif,nsoilmx,nqtot,day_ini_fi,
202     .      timefi,
203     .      tsurf,tsoil,emis,q2,qsurf,
204!       change FF 05/2011
205     .       cloudfrac,totalcloudfrac,hice,
206!       change BC 05/2014
207     .       rnat,pctsrf_sic,tslab,tsea_ice,sea_ice)
208
209
210
211
212! load 'controle' array from physics start file
213
214       ierr = NF_OPEN (fichnom, NF_NOWRITE,nid1)
215       IF (ierr.NE.NF_NOERR) THEN
216         write(6,*)' Pb d''ouverture du fichier'//trim(fichnom)
217        CALL ABORT
218       ENDIF
219                                               
220      ierr = NF_INQ_VARID (nid1, "controle", varid)
221      IF (ierr .NE. NF_NOERR) THEN
222       PRINT*, "start2archive: Le champ <controle> est absent"
223       CALL abort
224      ENDIF
225#ifdef NC_DOUBLE
226       ierr = NF_GET_VAR_DOUBLE(nid1, varid, tab_cntrl_fi)
227#else
228      ierr = NF_GET_VAR_REAL(nid1, varid, tab_cntrl_fi)
229#endif
230       IF (ierr .NE. NF_NOERR) THEN
231          PRINT*, "start2archive: Lecture echoue pour <controle>"
232          CALL abort
233       ENDIF
234
235      ierr = NF_CLOSE(nid1)
236
237
238c-----------------------------------------------------------------------
239c Controle de la synchro
240c-----------------------------------------------------------------------
241!mars a voir      if ((day_ini_fi.ne.day_ini).or.(abs(timefi-timedyn).gt.1.e-10))
242      if ((day_ini_fi.ne.day_ini)) 
243     &  stop ' Probleme de Synchro entre start et startfi !!!'
244
245
246c *****************************************************************
247c    Option : Reinitialisation des dates dans la premieres annees :
248       do while (day_ini.ge.year_day)
249          day_ini=day_ini-year_day
250       enddo
251c *****************************************************************
252
253c-----------------------------------------------------------------------
254c   Initialisations 
255c-----------------------------------------------------------------------
256
257      CALL defrun_new(99, .FALSE. )
258      call iniconst
259      call inigeom
260      call inifilr
261      CALL pression(ip1jmp1, ap, bp, ps, p3d)
262      call exner_hyb(ip1jmp1, ps, p3d, beta, pks, pk, pkf)
263
264c=======================================================================
265c Transformation EN VARIABLE NATURELLE / GRILLE SCALAIRE si necessaire
266c=======================================================================
267c  Les variables modeles dependent de la resolution. Il faut donc
268c  eliminer les facteurs responsables de cette dependance
269c  (pour utiliser newstart)
270c=======================================================================
271
272c-----------------------------------------------------------------------
273c Vent   (depend de la resolution horizontale) 
274c-----------------------------------------------------------------------
275c
276c ucov --> un  et  vcov --> vn
277c un --> us  et   vn --> vs
278c
279c-----------------------------------------------------------------------
280
281      call covnat(llm,ucov, vcov, un, vn) 
282      call wind_scal(un,vn,us,vs) 
283
284c-----------------------------------------------------------------------
285c Temperature  (depend de la resolution verticale => de "sigma.def")
286c-----------------------------------------------------------------------
287c
288c h --> T
289c
290c-----------------------------------------------------------------------
291
292      DO l=1,llm
293         DO ij=1,ip1jmp1
294            T(ij,l)=teta(ij,l)*pk(ij,l)/cpp !mars deduit de l'equation dans newstart
295         ENDDO
296      ENDDO
297
298c-----------------------------------------------------------------------
299c Variable physique 
300c-----------------------------------------------------------------------
301c
302c tsurf --> tsurfS
303c co2ice --> co2iceS
304c tsoil --> tsoilS
305c emis --> emisS
306c q2 --> q2S
307c qsurf --> qsurfS
308c
309c-----------------------------------------------------------------------
310
311      call gr_fi_dyn(1,ngridmx,iip1,jjp1,tsurf,tsurfS)
312!      call gr_fi_dyn(1,ngridmx,iip1,jjp1,co2ice,co2iceS)
313      call gr_fi_dyn(nsoilmx,ngridmx,iip1,jjp1,tsoil,tsoilS)
314      ! Note: thermal inertia "inertiedat" is in comsoil.h
315      call gr_fi_dyn(nsoilmx,ngridmx,iip1,jjp1,inertiedat,ithS)
316      call gr_fi_dyn(1,ngridmx,iip1,jjp1,emis,emisS)
317      call gr_fi_dyn(llm+1,ngridmx,iip1,jjp1,q2,q2S)
318      call gr_fi_dyn(nqtot,ngridmx,iip1,jjp1,qsurf,qsurfS)
319      call gr_fi_dyn(llm,ngridmx,iip1,jjp1,cloudfrac,cloudfracS)
320      call gr_fi_dyn(1,ngridmx,iip1,jjp1,hice,hiceS)
321      call gr_fi_dyn(1,ngridmx,iip1,jjp1,totalcloudfrac,totalcloudfracS)
322
323      call gr_fi_dyn(1,ngridmx,iip1,jjp1,rnat,rnatS)
324      call gr_fi_dyn(1,ngridmx,iip1,jjp1,pctsrf_sic,pctsrf_sicS)
325      call gr_fi_dyn(1,ngridmx,iip1,jjp1,tsea_ice,tsea_iceS)
326      call gr_fi_dyn(1,ngridmx,iip1,jjp1,sea_ice,sea_iceS)
327      call gr_fi_dyn(noceanmx,ngridmx,iip1,jjp1,tslab,tslabS)
328
329c=======================================================================
330c Info pour controler
331c=======================================================================
332
333      ptotal =  0.
334      co2icetotal = 0.
335      DO j=1,jjp1
336         DO i=1,iim
337           ptotal=ptotal+aire(i+(iim+1)*(j-1))*ps(i+(iim+1)*(j-1))/g
338!           co2icetotal = co2icetotal +
339!     &            co2iceS(i+(iim+1)*(j-1))*aire(i+(iim+1)*(j-1))
340         ENDDO
341      ENDDO
342      write(*,*)'Ancienne grille : masse de l''atm :',ptotal
343!      write(*,*)'Ancienne grille : masse de la glace CO2 :',co2icetotal
344
345c-----------------------------------------------------------------------
346c Passage de "ptotal" et "co2icetotal" par tab_cntrl_fi
347c-----------------------------------------------------------------------
348
349      tab_cntrl_fi(49) = ptotal
350      tab_cntrl_fi(50) = co2icetotal
351
352c=======================================================================
353c Ecriture dans le fichier  "start_archive"
354c=======================================================================
355
356c-----------------------------------------------------------------------
357c Ouverture de "start_archive" 
358c-----------------------------------------------------------------------
359
360      ierr = NF_OPEN ('start_archive.nc', NF_WRITE,nid)
361 
362c-----------------------------------------------------------------------
363c  si "start_archive" n'existe pas:
364c    1_ ouverture
365c    2_ creation de l'entete dynamique ("ini_archive")
366c-----------------------------------------------------------------------
367c ini_archive:
368c On met dans l'entete le tab_cntrl dynamique (1 a 16)
369c  On y ajoute les valeurs du tab_cntrl_fi (a partir de 51)
370c  En plus les deux valeurs ptotal et co2icetotal (99 et 100)
371c-----------------------------------------------------------------------
372
373      if (ierr.ne.NF_NOERR) then
374         write(*,*)'OK, Could not open file "start_archive.nc"'
375         write(*,*)'So let s create a new "start_archive"'
376         ierr = NF_CREATE('start_archive.nc', NF_CLOBBER, nid)
377         call ini_archive(nid,day_ini,phis,ithS,tab_cntrl_fi,
378     &                                          tab_cntrl_dyn)
379      endif
380
381c-----------------------------------------------------------------------
382c Ecriture de la coordonnee temps (date en jours)
383c-----------------------------------------------------------------------
384
385      date = day_ini
386      ierr= NF_INQ_VARID(nid,"Time",varid)
387      ierr= NF_INQ_DIMID(nid,"Time",dimid)
388      ierr= NF_INQ_DIMLEN(nid,dimid,timelen)
389      ntime=timelen+1
390
391      write(*,*) "******************"
392      write(*,*) "ntime",ntime
393      write(*,*) "******************"
394#ifdef NC_DOUBLE
395      ierr= NF_PUT_VARA_DOUBLE(nid,varid,ntime,1,date)
396#else
397      ierr= NF_PUT_VARA_REAL(nid,varid,ntime,1,date)
398#endif
399      if (ierr.ne.NF_NOERR) then
400         write(*,*) "time matter ",NF_STRERROR(ierr)
401         stop
402      endif
403
404c-----------------------------------------------------------------------
405c Ecriture des champs  (co2ice,emis,ps,Tsurf,T,u,v,q2,q,qsurf)
406c-----------------------------------------------------------------------
407c ATTENTION: q2 a une couche de plus!!!!
408c    Pour creer un fichier netcdf lisible par grads,
409c    On passe donc une des couches de q2 a part
410c    comme une variable 2D (la couche au sol: "q2surf")
411c    Les lmm autres couches sont nommees "q2atm" (3D)
412c-----------------------------------------------------------------------
413
414!      call write_archive(nid,ntime,'co2ice','couche de glace co2',
415!     &  'kg/m2',2,co2iceS)
416      call write_archive(nid,ntime,'emis','grd emis',' ',2,emisS)
417      call write_archive(nid,ntime,'ps','Psurf','Pa',2,ps)
418      call write_archive(nid,ntime,'tsurf','surf T','K',2,tsurfS)
419      call write_archive(nid,ntime,'temp','temperature','K',3,t)
420      call write_archive(nid,ntime,'u','Vent zonal','m.s-1',3,us)
421      call write_archive(nid,ntime,'v','Vent merid','m.s-1',3,vs)
422      call write_archive(nid,ntime,'q2surf','wind variance','m2.s-2',2,
423     .              q2S)
424      call write_archive(nid,ntime,'q2atm','wind variance','m2.s-2',3,
425     .              q2S(1,2))
426
427c-----------------------------------------------------------------------
428c Ecriture du champs  q  ( q[1,nqtot] )
429c-----------------------------------------------------------------------
430      do iq=1,nqtot
431        call write_archive(nid,ntime,tname(iq),'tracer','kg/kg',
432     &         3,q(1,1,iq))
433      end do
434c-----------------------------------------------------------------------
435c Ecriture du champs  qsurf  ( qsurf[1,nqtot] )
436c-----------------------------------------------------------------------
437      do iq=1,nqtot
438        txt=trim(tname(iq))//"_surf"
439        call write_archive(nid,ntime,txt,'Tracer on surface',
440     &  'kg.m-2',2,qsurfS(1,iq))
441      enddo
442
443
444c-----------------------------------------------------------------------
445c Ecriture du champs  tsoil  ( Tg[1,10] )
446c-----------------------------------------------------------------------
447c "tsoil" Temperature au sol definie dans 10 couches dans le sol
448c   Les 10 couches sont lues comme 10 champs
449c  nommees Tg[1,10]
450
451c      do isoil=1,nsoilmx
452c       write(str2,'(i2.2)') isoil
453c       call write_archive(nid,ntime,'Tg'//str2,'Ground Temperature ',
454c     .   'K',2,tsoilS(1,isoil))
455c      enddo
456
457! Write soil temperatures tsoil
458      call write_archive(nid,ntime,'tsoil','Soil temperature',
459     &     'K',-3,tsoilS)
460
461! Write soil thermal inertia
462      call write_archive(nid,ntime,'inertiedat',
463     &     'Soil thermal inertia',
464     &     'J.s-1/2.m-2.K-1',-3,ithS)
465
466! Write (0D) volumetric heat capacity (stored in comsoil.h)
467!      call write_archive(nid,ntime,'volcapa',
468!     &     'Soil volumetric heat capacity',
469!     &     'J.m-3.K-1',0,volcapa)
470! Note: no need to write volcapa, it is stored in "controle" table
471
472c-----------------------------------------------------------------------
473c Ecriture du champs  cloudfrac,hice,totalcloudfrac
474c-----------------------------------------------------------------------
475      call write_archive(nid,ntime,'hice',
476     &         'Height of oceanic ice','m',2,hiceS)
477      call write_archive(nid,ntime,'totalcloudfrac',
478     &        'Total cloud Fraction','',2,totalcloudfracS)
479      call write_archive(nid,ntime,'cloudfrac'
480     &        ,'Cloud fraction','',3,cloudfracS)
481
482c-----------------------------------------------------------------------
483c Slab ocean
484c-----------------------------------------------------------------------
485      OPEN(99,file='callphys.def',status='old',form='formatted'
486     &     ,iostat=ierr)
487      CLOSE(99)
488
489      IF(ierr.EQ.0) THEN
490
491
492         write(*,*) "Use slab-ocean ?"
493         ok_slab_ocean=.false.         ! default value
494         call getin("ok_slab_ocean",ok_slab_ocean)
495         write(*,*) "ok_slab_ocean = ",ok_slab_ocean
496
497      if(ok_slab_ocean) then
498      call write_archive(nid,ntime,'rnat'
499     &        ,'rnat','',2,rnatS)
500      call write_archive(nid,ntime,'pctsrf_sic'
501     &        ,'pctsrf_sic','',2,pctsrf_sicS)
502      call write_archive(nid,ntime,'sea_ice'
503     &        ,'sea_ice','',2,sea_iceS)
504      call write_archive(nid,ntime,'tslab'
505     &        ,'tslab','',3,tslabS)
506      call write_archive(nid,ntime,'tsea_ice'
507     &        ,'tsea_ice','',2,tsea_iceS)
508      endif !ok_slab_ocean
509      ENDIF
510c-----------------------------------------------------------------------
511c Fin 
512c-----------------------------------------------------------------------
513      ierr=NF_CLOSE(nid)
514
515      write(*,*) "start2archive: All is well that ends well."
516
517      end 
Note: See TracBrowser for help on using the repository browser.