source: trunk/LMDZ.GENERIC/libf/dyn3d/start2archive.F @ 832

Last change on this file since 832 was 787, checked in by aslmd, 12 years ago

LMDZ.GENERIC. (Sorry for long text but this is a quite major commit)

Paving the path for parallel computations. And moving towards a more flexible code.

Automatic allocation is used within all routines in phystd. No further mention to ngridmx and nqmx.

  1. ngridmx and nqmx are still used in LMDZ.GENERIC in the dyn3d part
  2. if the LMDZ4/LMDZ5 dynamical core is used, there is no more fixed dimensions ngridmx and nqmx --> a fully flexible parallel implementation is now possible (e.g. no need to recompile when changing numbers of processors)

The important stuff :

  • Compilation checked with ifort. OK with and without debug mode. No errors. Checked for: gcm, newstart, rcm1d, kcm1d
  • RUN GCM: Running an Earth test case. Comparison with previous revision --> debug mode : perfect match. bit by bit (diff command). checked with plots --> O1 mode : close match (checked with plots) --> O2 mode : sometimes up to 0.5 K departure.... BUT in this new version O2 and O1 are quite close while in previous version O1 and O2 differed by about, well, typically 0.5 K (pictures available on request)
  • RUN NEWSTART : perfect match (bit-by-bit) in either debug or normal mode.
  • RUN RCM1D : perfect match in normal mode.
  • RUN KCM1D : not tested (I don't know what is the use of kcm1d)

List of main changes :

  • Additional arguments to some subroutines (ngrid and nq)
  • F77 include strategy is obsolete and replaced by F90 module strategy In this new strategy arrays are allocatable and allocated once at first use This has to be done for all common featuring arrays defined with ngridmx or nqmx

surfdat.h >> surfdat_h.F90
tracer.h >> tracer_h.F90
comsaison.h >> comsaison_h.F90
comgeomfi.h >> comgeomfi_h.F90
comsoil.h >> comsoil_h.F90
comdiurn.h >> comdiurn_h.F90
fisice.h >> DELETED. was not used. probably a fossil.
watercap.h >> DELETED. variable put in surfdat_h.F90

  • F77 'save' strategy is obsolete and replaced by F90 'allocatable save' strategy (see previous point and e.g. new version of physiq.F90)
  • Suppressing any mention to advtrac.h which is a common in the dynamics and needs nqmx This was easily solved by adding an argument with tracer names, coming from the dynamics This is probably not a definitive solution, ... but this allows for generic physics to work easily with either LMDZ.GENERIC or LMDZ dynamical cores
  • Removing consistency tests between nq and nqmx ; and ngrid and ngridmx. No use now!
  • Adaptation of rcm1d, kcm1d, newstart given above-mentioned changes

A note on phyetat0 and soil_setting:

  • Now written so that a slice of horizontal size 'ngrid' starting at grid point 'cursor' is read in startfi.nc 'cursor' is defined in dimphys.h and initialized by inifis (or in newstart) this is useful for parallel computations. default behavior is the usual one : sequential runs, cursor is 1, size ngrid is the whole global domain

A note on an additional change :

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