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

Last change on this file since 1523 was 1478, checked in by bclmd, 9 years ago

bug for dimensions of tslab in start2archive

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','',-2,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.