source: trunk/LMDZ.MARS/libf/phymars/inistats.F @ 1532

Last change on this file since 1532 was 1532, checked in by emillour, 9 years ago

Mars GCM:

  • Some fixes for buggy outputs in 1D introduced by previous code modifications.
  • made statto.h a module.
  • ecritphy in dyn3d/control_mod.F90 should be an integer, not a real.

EM

File size: 4.4 KB
Line 
1      subroutine inistats(ierr)
2
3      use statto_mod, only: istats,istime
4      use mod_phys_lmdz_para, only : is_master
5      USE comvert_mod, ONLY: ap,bp,aps,bps,preff,pseudoalt,presnivs
6      USE comcstfi_h, ONLY: pi
7      USE time_phylmdz_mod, ONLY: daysec,dtphys
8      USE regular_lonlat_mod, ONLY: lon_reg, lat_reg
9      USE mod_grid_phy_lmdz, ONLY: nbp_lon, nbp_lat, nbp_lev
10      implicit none
11
12      include "netcdf.inc"
13
14      integer,intent(out) :: ierr
15      integer :: nid
16      integer :: l,nsteppd
17      real, dimension(nbp_lev) ::  sig_s
18      real,allocatable :: lon_reg_ext(:) ! extended longitudes
19      integer :: idim_lat,idim_lon,idim_llm,idim_llmp1,idim_time
20      real, dimension(istime) :: lt
21      integer :: nvarid
22
23
24      IF (nbp_lon*nbp_lat==1) THEN
25        ! 1D model
26        ALLOCATE(lon_reg_ext(1))
27      ELSE
28        ! 3D model
29        ALLOCATE(lon_reg_ext(nbp_lon+1))
30      ENDIF
31     
32      write (*,*)
33      write (*,*) '                        || STATS ||'
34      write (*,*)
35      write (*,*) 'daysec',daysec
36      write (*,*) 'dtphys',dtphys
37      nsteppd=nint(daysec/dtphys)
38      write (*,*) 'nsteppd=',nsteppd
39      if (abs(float(nsteppd)-daysec/dtphys).gt.1.e-8*daysec)
40     &   stop'Dans Instat:  1jour .ne. n pas physiques'
41
42      if(mod(nsteppd,istime).ne.0)
43     &   stop'Dans Instat:  1jour .ne. n*istime pas physiques'
44
45      istats=nsteppd/istime
46      write (*,*) 'istats=',istats
47      write (*,*) 'Storing ',istime,'times per day'
48      write (*,*) 'thus every ',istats,'physical timestep '
49      write (*,*)
50
51      do l= 1, nbp_lev
52         sig_s(l)=((ap(l)+ap(l+1))/preff+bp(l)+bp(l+1))/2.
53         pseudoalt(l)=-10.*log(presnivs(l)/preff)   
54      enddo
55     
56      lon_reg_ext(1:nbp_lon)=lon_reg(1:nbp_lon)
57      IF (nbp_lon*nbp_lat/=1) THEN
58        ! In 3D, add extra redundant point (180 degrees,
59        ! since lon_reg starts at -180)
60        lon_reg_ext(nbp_lon+1)=-lon_reg_ext(1)
61      ENDIF
62
63      if (is_master) then
64      ! only the master needs do this
65
66      ierr = NF_CREATE("stats.nc",IOR(NF_CLOBBER,NF_64BIT_OFFSET),nid)
67      if (ierr.ne.NF_NOERR) then
68         write (*,*) NF_STRERROR(ierr)
69         stop ""
70      endif
71
72      ierr = NF_DEF_DIM (nid, "latitude", nbp_lat, idim_lat)
73      IF (nbp_lon*nbp_lat==1) THEN
74        ierr = NF_DEF_DIM (nid, "longitude", 1, idim_lon)
75      ELSE
76        ierr = NF_DEF_DIM (nid, "longitude", nbp_lon+1, idim_lon)
77      ENDIF
78      ierr = NF_DEF_DIM (nid, "altitude", nbp_lev, idim_llm)
79      ierr = NF_DEF_DIM (nid, "llmp1", nbp_lev+1, idim_llmp1)
80      ierr = NF_DEF_DIM (nid, "Time", NF_UNLIMITED, idim_time)
81
82      ierr = NF_ENDDEF(nid)
83      call def_var_stats(nid,"Time","Time",
84     &            "days since 0000-00-0 00:00:00",1,
85     &            idim_time,nvarid,ierr)
86! Time is initialised later by mkstats subroutine
87
88      call def_var_stats(nid,"latitude","latitude",
89     &            "degrees_north",1,idim_lat,nvarid,ierr)
90#ifdef NC_DOUBLE
91      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,lat_reg/pi*180)
92#else
93      ierr = NF_PUT_VAR_REAL (nid,nvarid,lat_reg/pi*180)
94#endif
95      call def_var_stats(nid,"longitude","East longitude",
96     &            "degrees_east",1,idim_lon,nvarid,ierr)
97#ifdef NC_DOUBLE
98      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,lon_reg_ext/pi*180)
99#else
100      ierr = NF_PUT_VAR_REAL (nid,nvarid,lon_reg_ext/pi*180)
101#endif
102
103! Niveaux verticaux, aps et bps
104      ierr = NF_REDEF (nid)
105#ifdef NC_DOUBLE
106      ierr = NF_DEF_VAR (nid,"altitude", NF_DOUBLE, 1,idim_llm,nvarid)
107#else
108      ierr = NF_DEF_VAR (nid,"altitude", NF_FLOAT, 1,idim_llm,nvarid)
109#endif
110      ierr = NF_PUT_ATT_TEXT (nid,nvarid,"long_name",8,"altitude")
111      ierr = NF_PUT_ATT_TEXT (nid,nvarid,'units',2,"km")
112      ierr = NF_PUT_ATT_TEXT (nid,nvarid,'positive',2,"up")
113      ierr = NF_ENDDEF(nid)
114#ifdef NC_DOUBLE
115      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,pseudoalt)
116#else
117      ierr = NF_PUT_VAR_REAL (nid,nvarid,pseudoalt)
118#endif
119      call def_var_stats(nid,"aps","hybrid pressure at midlayers"
120     & ," ",1,idim_llm,nvarid,ierr)
121#ifdef NC_DOUBLE
122      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,aps)
123#else
124      ierr = NF_PUT_VAR_REAL (nid,nvarid,aps)
125#endif
126
127      call def_var_stats(nid,"bps","hybrid sigma at midlayers"
128     & ," ",1,idim_llm,nvarid,ierr)
129#ifdef NC_DOUBLE
130      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,bps)
131#else
132      ierr = NF_PUT_VAR_REAL (nid,nvarid,bps)
133#endif
134
135      ierr=NF_CLOSE(nid)
136
137      endif ! of if (is_master)
138      end subroutine inistats
139
Note: See TracBrowser for help on using the repository browser.