source: trunk/LMDZ.GENERIC/libf/dynlonlat_phylonlat/phystd/start2archive.F @ 1428

Last change on this file since 1428 was 1422, checked in by milmd, 10 years ago

In GENERIC, MARS and COMMON models replace some include files by modules (usefull for decoupling physics with dynamics).

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