source: trunk/LMDZ.PLUTO/libf/dynphy_lonlat/phypluto/start2archive.F @ 3411

Last change on this file since 3411 was 3380, checked in by tbertrand, 7 months ago

LMDZ.PLUTO
Removing useless stuff for Pluto in start2archive.F
TB

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