source: trunk/LMDZ.COMMON/libf/dyn3d_common/dynredem.F90 @ 3509

Last change on this file since 3509 was 3509, checked in by jbclement, 13 days ago

Dynamic + Mars PCM:
Addition of the description for the 'controle' array in the "start.nc" and "startfi.nc" files. It is given by the variable 'controle_descriptor' whose the element 'controle_descriptor(i)' explains 'controle(i)'.
JBC

File size: 16.5 KB
Line 
1!
2! $Id: dynredem.F 1635 2012-07-12 11:37:16Z lguez $
3!
4!
5SUBROUTINE dynredem0(fichnom,iday_end,phis)
6#ifdef CPP_IOIPSL
7  USE IOIPSL
8#endif
9  USE infotrac, ONLY: nqtot, tname, ttext
10  USE netcdf, ONLY: NF90_CREATE, NF90_DEF_DIM, NF90_INQ_VARID, NF90_GLOBAL, &
11                    NF90_CLOSE,  NF90_PUT_ATT, NF90_UNLIMITED, NF90_CLOBBER
12  USE dynredem_mod, ONLY: cre_var, put_var1, put_var2, put_char1, err, modname, fil
13  use netcdf95, only: NF95_PUT_VAR
14  use control_mod, only : planet_type
15  USE comvert_mod, ONLY: ap,bp,aps,bps,presnivs,pseudoalt,pa,preff, &
16                        nivsig,nivsigs
17  USE comconst_mod, ONLY: daysec,dtvr,rad,omeg,g,cpp,kappa,pi
18  USE logic_mod, ONLY: fxyhypb,ysinus
19  USE serre_mod, ONLY: clon,clat,grossismx,grossismy,dzoomx,dzoomy, &
20                        taux,tauy
21  USE temps_mod, ONLY: annee_ref,day_ref,itau_dyn,itaufin, &
22                        start_time,hour_ini
23  USE ener_mod, ONLY: etot0,ptot0,ztot0,stot0,ang0
24
25  IMPLICIT NONE
26!=======================================================================
27! Writting the NetCDF restart file (initialisation)
28!=======================================================================
29!   Declarations:
30!   -------------
31  include "dimensions.h"
32  include "paramet.h"
33  include "comgeom2.h"
34  include "netcdf.inc"
35  include "iniprint.h"
36
37!===============================================================================
38! Arguments:
39  CHARACTER(LEN=*), INTENT(IN) :: fichnom          !--- FILE NAME
40  INTEGER,          INTENT(IN) :: iday_end         !---
41  REAL,             INTENT(IN) :: phis(iip1, jjp1) !--- GROUND GEOPOTENTIAL
42!===============================================================================
43! Local variables:
44  INTEGER :: iq,l
45  INTEGER, PARAMETER :: length=100
46  INTEGER, PARAMETER :: ldscrpt = 28, ndscrpt = 42
47  REAL :: tab_cntrl(length) ! run parameters
48  character(ndscrpt), dimension(ldscrpt), parameter :: dscrpt_tab_cntrl = (/ &
49      "(1)  Number of nodes along longitude      ", &
50      "(2)  Number of nodes along latitude       ", &
51      "(3)  Number of atmospheric layers         ", &
52      "(4)  Initial day                          ", &
53      "(5)  Radius of the planet                 ", &
54      "(6)  Rotation of the planet (rad/s)       ", &
55      "(7)  Gravity (m/s2)                       ", &
56      "(8)  Specific heat Cp (J.kg-1.K-1)        ", &
57      "(9)  = r/Cp (=kappa)                      ", &
58      "(10) Lenght of a sol (s)                  ", &
59      "(11) Dynamical time step (s)              ", &
60      "(12) Total energy                         ", &
61      "(13) Total pressure                       ", &
62      "(14) Total enstrophy                      ", &
63      "(15) Total enthalpy                       ", &
64      "(16) Total angular momentum               ", &
65      "(17) Reference pressure (Pa)              ", &
66      "(18) Reference surface pressure (Pa)      ", &
67      "(19) Longitude of center of zoom          ", &
68      "(20) Latitude of center of zoom           ", &
69      "(21) Zooming factor, along longitude      ", &
70      "(22) Zooming factor, along latitude       ", &
71      "(23) -                                    ", &
72      "(24) Extention (in longitude) of zoom     ", &
73      "(25) Extention (in latitude) of zoom      ", &
74      "(26) -                                    ", &
75      "(27) Stiffness factor of zoom in longitude", &
76      "(28) Stiffness factor of zoom in latitude "/)
77  INTEGER :: ierr
78  CHARACTER(LEN=80) :: abort_message
79
80!   For NetCDF:
81  CHARACTER(LEN=30) :: unites
82  INTEGER :: indexID, descrptID, dscrpt_sID
83  INTEGER :: rlonuID, rlonvID, rlatuID, rlatvID
84  INTEGER :: sID, sigID, nID, vID, timID
85  INTEGER :: yyears0, jjour0, mmois0
86  REAL    :: zan0, zjulian, hours
87
88  CHARACTER(len=12) :: start_file_type="earth" ! default start file type
89  INTEGER :: idecal
90
91!===============================================================================
92  ! fill dynredem_mod module variables
93  modname='dynredem0'; fil=fichnom
94
95#ifdef CPP_IOIPSL
96  call ymds2ju(annee_ref, 1, iday_end, 0.0, zjulian)
97  call ju2ymds(zjulian, yyears0, mmois0, jjour0, hours)
98#else
99! set yyears0, mmois0, jjour0 to 0,1,1 (hours is not used)
100  yyears0=0
101  mmois0=1
102  jjour0=1
103#endif       
104
105  !!! AS: idecal is a hack to be able to read planeto starts...
106  !!!     .... while keeping everything OK for LMDZ EARTH
107  if ((planet_type.eq."generic").or.(planet_type.eq."mars")) then
108    write(lunout,*) trim(modname),' : Planeto-like start file'
109    start_file_type="planeto"
110    idecal = 4
111  else
112    if ( planet_type.eq."titan" ) then
113      ! Titan inherited Earth-like start files with idecal=5
114      write(lunout,*) trim(modname),' : Titan start file'
115    else
116      write(lunout,*) trim(modname),' : Earth-like start file'
117    endif
118    idecal = 5
119  endif
120
121  tab_cntrl(:)  = 0.
122  tab_cntrl(1)  = REAL(iim)
123  tab_cntrl(2)  = REAL(jjm)
124  tab_cntrl(3)  = REAL(llm)
125  if (start_file_type.eq."earth") then
126    tab_cntrl(4)=REAL(day_ref)
127  else
128    !tab_cntrl(4)=REAL(day_end)
129    tab_cntrl(4)=REAL(iday_end)
130  endif
131  tab_cntrl(5)  = REAL(annee_ref)
132  tab_cntrl(idecal+1)  = rad
133  tab_cntrl(idecal+2)  = omeg
134  tab_cntrl(idecal+3)  = g
135  tab_cntrl(idecal+4)  = cpp
136  tab_cntrl(idecal+5) = kappa
137  tab_cntrl(idecal+6) = daysec
138  tab_cntrl(idecal+7) = dtvr
139  tab_cntrl(idecal+8) = etot0
140  tab_cntrl(idecal+9) = ptot0
141  tab_cntrl(idecal+10) = ztot0
142  tab_cntrl(idecal+11) = stot0
143  tab_cntrl(idecal+12) = ang0
144  tab_cntrl(idecal+13) = pa
145  tab_cntrl(idecal+14) = preff
146
147!    .....    parameters for the zoom      ......   
148  tab_cntrl(idecal+15)  = clon
149  tab_cntrl(idecal+16)  = clat
150  tab_cntrl(idecal+17)  = grossismx
151  tab_cntrl(idecal+18)  = grossismy
152!
153  IF ( fxyhypb )   THEN
154    tab_cntrl(idecal+19) = 1.
155    tab_cntrl(idecal+20) = dzoomx
156    tab_cntrl(idecal+21) = dzoomy
157    tab_cntrl(idecal+22) = 0.
158    tab_cntrl(idecal+23) = taux
159    tab_cntrl(idecal+24) = tauy
160  ELSE
161    tab_cntrl(idecal+19) = 0.
162    tab_cntrl(idecal+20) = dzoomx
163    tab_cntrl(idecal+21) = dzoomy
164    tab_cntrl(idecal+22) = 0.
165    tab_cntrl(idecal+23) = 0.
166    tab_cntrl(idecal+24) = 0.
167    IF( ysinus )  tab_cntrl(idecal+22) = 1.
168  ENDIF
169
170  if (start_file_type.eq."earth") then
171    tab_cntrl(idecal+25) = REAL(iday_end)
172    tab_cntrl(idecal+26) = REAL(itau_dyn + itaufin)
173! start_time: start_time of simulation (not necessarily 0.)
174    tab_cntrl(idecal+27) = start_time
175  endif
176
177  if (planet_type=="mars") then ! For Mars only
178    !tab_cntrl(29)=hour_ini
179    tab_cntrl(29)=0
180  endif
181
182!--- File creation
183  CALL err(NF90_CREATE(fichnom,NF90_CLOBBER,nid))
184
185!--- Some global attributes
186  CALL err(NF90_PUT_ATT(nid,NF90_GLOBAL,"title","Fichier demarrage dynamique"))
187
188!--- Dimensions
189  if (start_file_type.eq."earth") then
190    CALL err(NF90_DEF_DIM(nid,"index", length, indexID))
191    CALL err(NF90_DEF_DIM(nid,"rlonu", iip1,   rlonuID))
192    CALL err(NF90_DEF_DIM(nid,"rlatu", jjp1,   rlatuID))
193    CALL err(NF90_DEF_DIM(nid,"rlonv", iip1,   rlonvID))
194    CALL err(NF90_DEF_DIM(nid,"rlatv", jjm,    rlatvID))
195    CALL err(NF90_DEF_DIM(nid,"sigs",  llm,        sID))
196    CALL err(NF90_DEF_DIM(nid,"sig",   llmp1,    sigID))
197    CALL err(NF90_DEF_DIM(nid,"temps", NF90_UNLIMITED, timID))
198    CALL err(NF90_DEF_DIM(nid,"descriptor", ldscrpt, descrptID))
199    CALL err(NF90_DEF_DIM(nid,"description_size", ndscrpt, dscrpt_sID))
200  else
201    CALL err(NF90_DEF_DIM(nid,"index", length, indexID))
202    CALL err(NF90_DEF_DIM(nid,"rlonu", iip1,   rlonuID))
203    CALL err(NF90_DEF_DIM(nid,"latitude", jjp1,   rlatuID))
204    CALL err(NF90_DEF_DIM(nid,"longitude", iip1,   rlonvID))
205    CALL err(NF90_DEF_DIM(nid,"rlatv", jjm,    rlatvID))
206    CALL err(NF90_DEF_DIM(nid,"altitude",  llm,        sID))
207    CALL err(NF90_DEF_DIM(nid,"interlayer",   llmp1,    sigID))
208    CALL err(NF90_DEF_DIM(nid,"Time", NF90_UNLIMITED, timID))
209    CALL err(NF90_DEF_DIM(nid,"descriptor", ldscrpt, descrptID))
210    CALL err(NF90_DEF_DIM(nid,"description_size", ndscrpt, dscrpt_sID))
211  endif
212
213!--- Define and save invariant fields
214  CALL put_var1(nid,"controle","Parametres de controle",[indexID],tab_cntrl)
215  CALL put_char1(nid,"controle_descriptor","Description of control parameters",[dscrpt_sID,descrptID],dscrpt_tab_cntrl)
216  CALL put_var1(nid,"rlonu"   ,"Longitudes des points U",[rlonuID],rlonu)
217  CALL put_var1(nid,"rlatu"   ,"Latitudes des points U" ,[rlatuID],rlatu)
218  CALL put_var1(nid,"rlonv"   ,"Longitudes des points V",[rlonvID],rlonv)
219  CALL put_var1(nid,"rlatv"   ,"Latitudes des points V" ,[rlatvID],rlatv)
220  if (start_file_type.eq."earth") then
221    CALL put_var1(nid,"nivsigs" ,"Numero naturel des couches s"    ,[sID]  ,nivsigs)
222    CALL put_var1(nid,"nivsig"  ,"Numero naturel des couches sigma",[sigID],nivsig)
223  endif ! of if (start_file_type.eq."earth")
224  CALL put_var1(nid,"ap"      ,"Coefficient A pour hybride"      ,[sigID],ap)
225  CALL put_var1(nid,"bp"      ,"Coefficient B pour hybride"      ,[sigID],bp)
226  if (start_file_type.ne."earth") then
227    CALL put_var1(nid,"aps","Coef AS: hybrid pressure at midlayers",[sID],aps)
228    CALL put_var1(nid,"bps","Coef BS: hybrid sigma at midlayers",[sID],bps)
229  endif ! of if (start_file_type.eq."earth")
230  CALL put_var1(nid,"presnivs",""                                ,[sID]  ,presnivs)
231  if (start_file_type.ne."earth") then
232        ierr = NF_REDEF (nid)
233#ifdef NC_DOUBLE
234        ierr = NF_DEF_VAR(nid,"latitude",NF_DOUBLE,1,rlatuID,vID)
235#else
236        ierr = NF_DEF_VAR(nid,"latitude",NF_FLOAT,1,rlatuID,vID)
237#endif
238        ierr =NF_PUT_ATT_TEXT(nid,vID,'units',13,"degrees_north")
239        ierr = NF_PUT_ATT_TEXT (nid,vID,"long_name", 14, &
240              "North latitude")
241        ierr = NF_ENDDEF(nid)
242        call NF95_PUT_VAR(nid,vID,rlatu*180/pi)
243!
244        ierr = NF_REDEF (nid)
245#ifdef NC_DOUBLE
246        ierr=NF_DEF_VAR(nid,"longitude",NF_DOUBLE,1,rlonvID,vID)
247#else
248        ierr=NF_DEF_VAR(nid,"longitude",NF_FLOAT,1,rlonvID,vID)
249#endif
250        ierr = NF_PUT_ATT_TEXT (nid,vID,"long_name", 14, &
251              "East longitude")
252        ierr = NF_PUT_ATT_TEXT(nid,vID,'units',12,"degrees_east")
253        ierr = NF_ENDDEF(nid)
254        call NF95_PUT_VAR(nid,vID,rlonv*180/pi)
255!
256        ierr = NF_REDEF (nid)
257#ifdef NC_DOUBLE
258        ierr = NF_DEF_VAR (nid, "altitude", NF_DOUBLE, 1, &
259             sID,vID)
260#else
261        ierr = NF_DEF_VAR (nid, "altitude", NF_FLOAT, 1, &
262             sID,vID)
263#endif
264        ierr = NF_PUT_ATT_TEXT(nid,vID,"long_name",10,"pseudo-alt")
265        ierr = NF_PUT_ATT_TEXT (nid,vID,'units',2,"km")
266        ierr = NF_PUT_ATT_TEXT (nid,vID,'positive',2,"up")
267        ierr = NF_ENDDEF(nid)
268        call NF95_PUT_VAR(nid,vID,pseudoalt)
269        CALL err(NF_REDEF(nid))
270  endif ! of if (start_file_type.ne."earth")
271
272! covariant <-> contravariant <-> natural conversion coefficients
273  CALL put_var2(nid,"cu","Coefficient de passage pour U",[rlonuID,rlatuID],cu)
274  CALL put_var2(nid,"cv","Coefficient de passage pour V",[rlonvID,rlatvID],cv)
275  CALL put_var2(nid,"aire","Aires de chaque maille"     ,[rlonvID,rlatuID],aire)
276  CALL put_var2(nid,"phisinit","Geopotentiel au sol"    ,[rlonvID,rlatuID],phis)
277
278
279! Define variables that will be stored later:
280  WRITE(unites,"('days since ',i4,'-',i2.2,'-',i2.2,' 00:00:00')") &
281               yyears0,mmois0,jjour0
282  IF (start_file_type.eq."earth") THEN
283    CALL cre_var(nid,"temps","Temps de simulation",[timID],unites)
284  ELSE
285    CALL cre_var(nid,"Time","Temps de simulation",[timID],unites)
286  ENDIF
287
288  CALL cre_var(nid,"ucov" ,"Vitesse U"  ,[rlonuID,rlatuID,sID,timID])
289  CALL cre_var(nid,"vcov" ,"Vitesse V"  ,[rlonvID,rlatvID,sID,timID])
290  CALL cre_var(nid,"teta" ,"Temperature",[rlonvID,rlatuID,sID,timID])
291
292  IF(nqtot.GE.1) THEN
293    DO iq=1,nqtot
294      CALL cre_var(nid,tname(iq),ttext(iq),[rlonvID,rlatuID,sID,timID])
295    END DO
296  ENDIF
297
298  CALL cre_var(nid,"masse","Masse d air"    ,[rlonvID,rlatuID,sID,timID])
299  CALL cre_var(nid,"ps"   ,"Pression au sol",[rlonvID,rlatuID    ,timID])
300
301  CALL err(NF90_CLOSE (nid)) ! close file
302
303  WRITE(lunout,*)TRIM(modname)//': iim,jjm,llm,iday_end',iim,jjm,llm,iday_end
304  WRITE(lunout,*)TRIM(modname)//': rad,omeg,g,cpp,kappa',rad,omeg,g,cpp,kappa
305
306END SUBROUTINE dynredem0
307
308!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
309
310SUBROUTINE dynredem1(fichnom,time,vcov,ucov,teta,q,masse,ps)
311!
312!-------------------------------------------------------------------------------
313! Purpose: Write the NetCDF restart file (append).
314!-------------------------------------------------------------------------------
315  USE infotrac, ONLY: nqtot, tname, type_trac
316  USE control_mod, only : planet_type
317  USE netcdf,   ONLY: NF90_OPEN,  NF90_NOWRITE, NF90_GET_VAR, NF90_INQ_VARID,  &
318                      NF90_CLOSE, NF90_WRITE,   NF90_PUT_VAR, NF90_NoErr
319  use netcdf95, only: NF95_PUT_VAR
320  USE temps_mod, ONLY: itaufin,itau_dyn
321  USE dynredem_mod, ONLY: dynredem_write_u, dynredem_write_v, dynredem_read_u, &
322                          err, modname, fil, msg
323 
324  IMPLICIT NONE
325  include "dimensions.h"
326  include "paramet.h"
327  include "netcdf.inc"
328  include "comgeom.h"
329  include "iniprint.h"
330!===============================================================================
331! Arguments:
332  CHARACTER(LEN=*), INTENT(IN) :: fichnom              !-- FILE NAME
333  REAL, INTENT(IN)    ::  time                         !-- TIME
334  REAL, INTENT(IN)    ::  vcov(iip1,jjm, llm)          !-- V COVARIANT WIND
335  REAL, INTENT(IN)    ::  ucov(iip1,jjp1,llm)          !-- U COVARIANT WIND
336  REAL, INTENT(IN)    ::  teta(iip1,jjp1,llm)          !-- POTENTIAL TEMPERATURE
337  REAL, INTENT(INOUT) ::     q(iip1,jjp1,llm,nqtot)    !-- TRACERS
338  REAL, INTENT(IN)    :: masse(iip1,jjp1,llm)          !-- MASS PER CELL
339  REAL, INTENT(IN)    ::    ps(iip1,jjp1)              !-- GROUND PRESSURE
340!===============================================================================
341! Local variables:
342  INTEGER :: l, iq, nid, vID, ierr, nid_trac, vID_trac
343  INTEGER,SAVE :: nb=0
344  INTEGER, PARAMETER :: length=100
345  REAL               :: tab_cntrl(length) ! tableau des parametres du run
346  CHARACTER(LEN=256) :: var, dum
347  LOGICAL            :: lread_inca
348  CHARACTER(LEN=80) :: abort_message
349  CHARACTER(len=12) :: start_file_type="earth" ! default start file type
350
351  ! fill dynredem_mod module variables
352  modname='dynredem1'; fil=fichnom
353
354  if ((planet_type.eq."generic").or.(planet_type.eq."mars")) then
355      write(lunout,*) trim(modname),' : Planeto-like start file'
356      start_file_type="planeto"
357  else
358      write(lunout,*) trim(modname),' : Earth-like start file'
359  endif
360
361  CALL err(NF90_OPEN(fil,NF90_WRITE,nid),"open",fil)
362
363!--- Write/extend time coordinate
364  nb = nb + 1
365  if (start_file_type.eq."earth") then
366        ierr = NF_INQ_VARID(nid, "temps", vID)
367        IF (ierr .NE. NF_NOERR) THEN
368          write(lunout,*) NF_STRERROR(ierr)
369          abort_message='Variable temps n est pas definie'
370          CALL abort_gcm(modname,abort_message,ierr)
371        ENDIF
372 else
373        ierr = NF_INQ_VARID(nid,"Time", vID)
374        IF (ierr .NE. NF_NOERR) THEN
375          write(lunout,*) NF_STRERROR(ierr)
376          abort_message='Variable Time not found'
377          CALL abort_gcm(modname,abort_message,ierr)
378        ENDIF
379  endif ! of if (start_file_type.eq."earth")
380  call NF95_PUT_VAR(nid,vID,time,start=(/nb/))
381  WRITE(lunout,*)TRIM(modname)//": Saving for ", nb, time
382
383!--- Rewrite control table (itaufin undefined in dynredem0)
384  var="controle"
385  CALL err(NF90_INQ_VARID(nid,var,vID),"inq",var)
386  CALL err(NF90_GET_VAR(nid,vID,tab_cntrl),"get",var)
387  if (start_file_type=="earth") then
388    tab_cntrl(31) = REAL(itau_dyn + itaufin)
389  else
390    tab_cntrl(31) = 0
391  endif
392  CALL err(NF90_INQ_VARID(nid,var,vID),"inq",var)
393  CALL err(NF90_PUT_VAR(nid,vID,tab_cntrl),"put",var)
394
395!--- Save fields
396  CALL dynredem_write_u(nid,"ucov" ,ucov ,llm, nb)
397  CALL dynredem_write_v(nid,"vcov" ,vcov ,llm, nb)
398  CALL dynredem_write_u(nid,"teta" ,teta ,llm, nb)
399  CALL dynredem_write_u(nid,"masse",masse,llm, nb)
400  CALL dynredem_write_u(nid,"ps"   ,ps   ,1, nb)
401
402!--- Tracers in file "start_trac.nc" (added by Anne)
403  lread_inca=.FALSE.; fil="start_trac.nc"
404  IF(type_trac=='inca') INQUIRE(FILE=fil,EXIST=lread_inca)
405  IF(lread_inca) CALL err(NF90_OPEN(fil,NF90_NOWRITE,nid_trac),"open")
406
407!--- Save tracers
408  IF(nqtot.GE.1) THEN
409    DO iq=1,nqtot
410      var=tname(iq); ierr=-1
411      IF(lread_inca) THEN                  !--- Possibly read from "start_trac.nc"
412        fil="start_trac.nc"
413        ierr=NF90_INQ_VARID(nid_trac,var,vID_trac)
414        dum='inq'; IF(ierr==NF90_NoErr) dum='fnd'
415        WRITE(lunout,*)msg(dum,var)
416
417
418        IF(ierr==NF90_NoErr) CALL dynredem_read_u(nid_trac,var,q(:,:,:,iq),llm)
419      END IF ! of IF(lread_inca)
420      fil=fichnom
421      CALL dynredem_write_u(nid,var,q(:,:,:,iq),llm,nb)
422    END DO ! of DO iq=1,nqtot
423  ENDIF ! of IF(nqtot.GE.1)
424
425  CALL err(NF90_CLOSE(nid),"close")
426  fil="start_trac.nc"
427  IF(lread_inca) CALL err(NF90_CLOSE(nid_trac),"close")
428
429END SUBROUTINE dynredem1
430
Note: See TracBrowser for help on using the repository browser.