source: trunk/LMDZ.COMMON/libf/dyn3dpar/dynredem_p.F90 @ 2236

Last change on this file since 2236 was 1850, checked in by jvatant, 7 years ago

Follow-up of r1849 : understood the black magic !
--JVO

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