Index: LMDZ5/branches/LMDZ5V2.0-dev/libf/dyn3d/etat0_netcdf.F90 =================================================================== --- LMDZ5/branches/LMDZ5V2.0-dev/libf/dyn3d/etat0_netcdf.F90 (revision 1485) +++ LMDZ5/branches/LMDZ5V2.0-dev/libf/dyn3d/etat0_netcdf.F90 (revision 1486) @@ -98,5 +98,5 @@ REAL :: dummy LOGICAL :: ok_newmicro, ok_journe, ok_mensuel, ok_instan, ok_hf - LOGICAL :: ok_LES, ok_ade, ok_aie, aerosol_couple, new_aod + LOGICAL :: ok_LES, ok_ade, ok_aie, aerosol_couple, new_aod, callstats INTEGER :: iflag_radia, flag_aerosol REAL :: bl95_b0, bl95_b1, fact_cldcon, facttemps, ratqsbas, ratqshaut @@ -130,4 +130,5 @@ !--- CONSTRUCT A GRID CALL conf_phys( ok_journe, ok_mensuel, ok_instan, ok_hf, ok_LES, & + callstats, & solarlong0,seuil_inversion, & fact_cldcon, facttemps,ok_newmicro,iflag_radia, & Index: LMDZ5/branches/LMDZ5V2.0-dev/libf/dyn3dpar/etat0_netcdf.F90 =================================================================== --- LMDZ5/branches/LMDZ5V2.0-dev/libf/dyn3dpar/etat0_netcdf.F90 (revision 1485) +++ LMDZ5/branches/LMDZ5V2.0-dev/libf/dyn3dpar/etat0_netcdf.F90 (revision 1486) @@ -98,5 +98,5 @@ REAL :: dummy LOGICAL :: ok_newmicro, ok_journe, ok_mensuel, ok_instan, ok_hf - LOGICAL :: ok_LES, ok_ade, ok_aie, aerosol_couple, new_aod + LOGICAL :: ok_LES, ok_ade, ok_aie, aerosol_couple, new_aod, callstats INTEGER :: iflag_radia, flag_aerosol REAL :: bl95_b0, bl95_b1, fact_cldcon, facttemps, ratqsbas, ratqshaut @@ -130,4 +130,5 @@ !--- CONSTRUCT A GRID CALL conf_phys( ok_journe, ok_mensuel, ok_instan, ok_hf, ok_LES, & + callstats, & solarlong0,seuil_inversion, & fact_cldcon, facttemps,ok_newmicro,iflag_radia, & Index: LMDZ5/branches/LMDZ5V2.0-dev/libf/phylmd/conf_phys.F90 =================================================================== --- LMDZ5/branches/LMDZ5V2.0-dev/libf/phylmd/conf_phys.F90 (revision 1485) +++ LMDZ5/branches/LMDZ5V2.0-dev/libf/phylmd/conf_phys.F90 (revision 1486) @@ -13,4 +13,5 @@ subroutine conf_phys(ok_journe, ok_mensuel, ok_instan, ok_hf, & ok_LES,& + callstats,& solarlong0,seuil_inversion, & fact_cldcon, facttemps,ok_newmicro,iflag_radia,& @@ -66,4 +67,5 @@ logical :: ok_journe, ok_mensuel, ok_instan, ok_hf logical :: ok_LES + LOGICAL :: callstats LOGICAL :: ok_ade, ok_aie, aerosol_couple INTEGER :: flag_aerosol @@ -79,4 +81,5 @@ logical,SAVE :: ok_journe_omp, ok_mensuel_omp, ok_instan_omp, ok_hf_omp logical,SAVE :: ok_LES_omp + LOGICAL,SAVE :: callstats_omp LOGICAL,SAVE :: ok_ade_omp, ok_aie_omp, aerosol_couple_omp INTEGER, SAVE :: flag_aerosol_omp @@ -1418,4 +1421,13 @@ ok_LES_omp = .false. call getin('OK_LES', ok_LES_omp) + +!Config Key = callstats +!Config Desc = Pour des sorties callstats +!Config Def = .false. +!Config Help = Pour creer le fichier stats contenant les sorties +! stats +! + callstats_omp = .false. + call getin('callstats', callstats_omp) ! !Config Key = ecrit_LES @@ -1581,4 +1593,5 @@ ok_hines = ok_hines_omp ok_LES = ok_LES_omp + callstats = callstats_omp ecrit_LES = ecrit_LES_omp carbon_cycle_tr = carbon_cycle_tr_omp Index: LMDZ5/branches/LMDZ5V2.0-dev/libf/phylmd/inistats.F =================================================================== --- LMDZ5/branches/LMDZ5V2.0-dev/libf/phylmd/inistats.F (revision 1486) +++ LMDZ5/branches/LMDZ5V2.0-dev/libf/phylmd/inistats.F (revision 1486) @@ -0,0 +1,133 @@ + subroutine inistats(ierr) + + implicit none + +#include "dimensions.h" +#include "paramet.h" +#include "comgeom.h" +#include "comvert.h" +#include "comconst.h" +#include "statto.h" +#include "netcdf.inc" + + integer,intent(out) :: ierr + integer :: nid + integer :: l,nsteppd + real, dimension(llm) :: sig_s + integer :: idim_lat,idim_lon,idim_llm,idim_llmp1,idim_time + real, dimension(istime) :: lt + integer :: nvarid + real, dimension(llm) :: pseudoalt + + write (*,*) + write (*,*) ' || STATS ||' + write (*,*) + write (*,*) 'daysec',daysec + write (*,*) 'dtphys',dtphys + nsteppd=nint(daysec/dtphys) + write (*,*) 'nsteppd=',nsteppd + if (abs(float(nsteppd)-daysec/dtphys).gt.1.e-8*daysec) + & stop'Dans Instat: 1jour .ne. n pas physiques' + + if(mod(nsteppd,istime).ne.0) + & stop'Dans Instat: 1jour .ne. n*istime pas physiques' + + istats=nsteppd/istime + write (*,*) 'istats=',istats + write (*,*) 'Storing ',istime,'times per day' + write (*,*) 'thus every ',istats,'physical timestep ' + write (*,*) + + do l= 1, llm + sig_s(l)=((ap(l)+ap(l+1))/preff+bp(l)+bp(l+1))/2. + pseudoalt(l)=log(preff/presnivs(l))*8. + enddo + + ierr = NF_CREATE("stats.nc",NF_CLOBBER,nid) + if (ierr.ne.NF_NOERR) then + write (*,*) NF_STRERROR(ierr) + stop "" + endif + + ierr = NF_DEF_DIM (nid, "latitude", jjp1, idim_lat) + ierr = NF_DEF_DIM (nid, "longitude", iip1, idim_lon) + ierr = NF_DEF_DIM (nid, "altitude", llm, idim_llm) + ierr = NF_DEF_DIM (nid, "llmp1", llm+1, idim_llmp1) + ierr = NF_DEF_DIM (nid, "Time", NF_UNLIMITED, idim_time) + + ierr = NF_ENDDEF(nid) + call def_var_stats(nid,"Time","Time", + & "hours since 0000-00-0 00:00:00",1, + & idim_time,nvarid,ierr) +! Time is initialised later by mkstats subroutine + + call def_var_stats(nid,"latitude","latitude", + & "degrees_north",1,idim_lat,nvarid,ierr) +#ifdef NC_DOUBLE + ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,rlatu/pi*180) +#else + ierr = NF_PUT_VAR_REAL (nid,nvarid,rlatu/pi*180) +#endif + call def_var_stats(nid,"longitude","East longitude", + & "degrees_east",1,idim_lon,nvarid,ierr) +#ifdef NC_DOUBLE + ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,rlonv/pi*180) +#else + ierr = NF_PUT_VAR_REAL (nid,nvarid,rlonv/pi*180) +#endif + +! Niveaux verticaux, aps et bps + ierr = NF_REDEF (nid) +! presnivs +#ifdef NC_DOUBLE + ierr = NF_DEF_VAR (nid,"presnivs", NF_DOUBLE, 1,idim_llm,nvarid) +#else + ierr = NF_DEF_VAR (nid,"presnivs", NF_FLOAT, 1,idim_llm,nvarid) +#endif + ierr = NF_PUT_ATT_TEXT (nid,nvarid,"long_name",15, + & "Vertical levels") + ierr = NF_PUT_ATT_TEXT (nid,nvarid,'units',2,"Pa") + ierr = NF_PUT_ATT_TEXT (nid,nvarid,'positive',4,"down") + ierr = NF_ENDDEF(nid) +#ifdef NC_DOUBLE + ierr = NF_PUT_VAR_DOUBLE (nid,nvarid, + & presnivs(1:llm)) +#else + ierr = NF_PUT_VAR_REAL (nid,nvarid, + & presnivs(1:llm)) +#endif +! Pseudo alts +#ifdef NC_DOUBLE + ierr = NF_DEF_VAR (nid,"altitude", NF_DOUBLE, 1,idim_llm,nvarid) +#else + ierr = NF_DEF_VAR (nid,"altitude", NF_FLOAT, 1,idim_llm,nvarid) +#endif + ierr = NF_PUT_ATT_TEXT (nid,nvarid,"long_name",8,"altitude") + ierr = NF_PUT_ATT_TEXT (nid,nvarid,'units',2,"km") + ierr = NF_PUT_ATT_TEXT (nid,nvarid,'positive',2,"up") + ierr = NF_ENDDEF(nid) +#ifdef NC_DOUBLE + ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,pseudoalt) +#else + ierr = NF_PUT_VAR_REAL (nid,nvarid,pseudoalt) +#endif +! call def_var_stats(nid,"aps","hybrid pressure at midlayers"," ", +! & 1,idim_llm,nvarid,ierr) +!#ifdef NC_DOUBLE +! ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,aps) +!#else +! ierr = NF_PUT_VAR_REAL (nid,nvarid,aps) +!#endif + +! call def_var_stats(nid,"bps","hybrid sigma at midlayers"," ", +! & 1,idim_llm,nvarid,ierr) +!#ifdef NC_DOUBLE +! ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,bps) +!#else +! ierr = NF_PUT_VAR_REAL (nid,nvarid,bps) +!#endif + + ierr=NF_CLOSE(nid) + + end subroutine inistats + Index: LMDZ5/branches/LMDZ5V2.0-dev/libf/phylmd/mkstat.F90 =================================================================== --- LMDZ5/branches/LMDZ5V2.0-dev/libf/phylmd/mkstat.F90 (revision 1486) +++ LMDZ5/branches/LMDZ5V2.0-dev/libf/phylmd/mkstat.F90 (revision 1486) @@ -0,0 +1,164 @@ +subroutine mkstats(ierr) + + +! +! This program writes a stats.nc file from sums and sums of squares +! to means and standard deviations and also writes netcdf style +! file so that the data can be viewed easily. The data file is +! overwritten in place. +! SRL 21 May 1996 +! Yann W. july 2003 + + +implicit none + +#include "dimensions.h" +#include "statto.h" +#include "netcdf.inc" + +integer,parameter :: iip1=iim+1 +integer,parameter :: jjp1=jjm+1 +integer :: ierr,nid,nbvar,i,ndims,lt,nvarid +integer, dimension(4) :: id,varid,start,size +integer, dimension(5) :: dimids +character (len=50) :: name,nameout,units,title +real, dimension(iip1,jjp1,llm) :: sum3d,square3d,mean3d,sd3d +real, dimension(iip1,jjp1) :: sum2d,square2d,mean2d,sd2d +real, dimension(istime) :: time +real, dimension(jjp1) :: lat +real, dimension(iip1) :: lon +real, dimension(llm) :: alt +logical :: lcopy=.true. +!integer :: latid,lonid,altid,timeid +integer :: meanid,sdid +!integer, dimension(4) :: dimout + +! Incrementation of count for the last step, which is not done in wstats +count(istime)=count(istime)+1 + +ierr = NF_OPEN("stats.nc",NF_WRITE,nid) + +! We catch the id of dimensions of the stats file + +ierr= NF_INQ_DIMID(nid,"latitude",id(1)) +ierr= NF_INQ_DIMID(nid,"longitude",id(2)) +ierr= NF_INQ_DIMID(nid,"altitude",id(3)) +ierr= NF_INQ_DIMID(nid,"Time",id(4)) + +ierr= NF_INQ_VARID(nid,"latitude",varid(1)) +ierr= NF_INQ_VARID(nid,"longitude",varid(2)) +ierr= NF_INQ_VARID(nid,"altitude",varid(3)) +ierr= NF_INQ_VARID(nid,"Time",varid(4)) + +! Time initialisation + +do i=1,istime + time(i)=i*24./istime +#ifdef NC_DOUBLE + ierr= NF_PUT_VARA_DOUBLE(nid,varid(4),i,1,time(i)) +#else + ierr= NF_PUT_VARA_REAL(nid,varid(4),i,1,time(i)) +#endif +enddo + +! We catche the values of the variables + +#ifdef NC_DOUBLE + ierr = NF_GET_VAR_DOUBLE(nid,varid(1),lat) + ierr = NF_GET_VAR_DOUBLE(nid,varid(2),lon) + ierr = NF_GET_VAR_DOUBLE(nid,varid(3),alt) +#else + ierr = NF_GET_VAR_REAL(nid,varid(1),lat) + ierr = NF_GET_VAR_REAL(nid,varid(2),lon) + ierr = NF_GET_VAR_REAL(nid,varid(3),alt) +#endif + +! We catch the number of variables in the stats file +ierr = NF_INQ_NVARS(nid,nbvar) + +! to catche the "real" number of variables (without the "additionnal variables") +nbvar=(nbvar-4)/2 + +do i=1,nbvar + varid=(i-1)*2+5 + + ! What's the variable's name? + ierr=NF_INQ_VARNAME(nid,varid,name) + write(*,*) "OK variable ",name + ! Its units? + units=" " + ierr=NF_GET_ATT_TEXT(nid,varid,"units",units) + ! Its title? + title=" " + ierr=NF_GET_ATT_TEXT(nid,varid,"title",title) + ! Its number of dimensions? + ierr=NF_INQ_VARNDIMS(nid,varid,ndims) + ! Its values? + + if(ndims==4) then ! lat, lon, alt & time + +! dimout(1)=lonid +! dimout(2)=latid +! dimout(3)=altid +! dimout(4)=timeid + + size=(/iip1,jjp1,llm,1/) + do lt=1,istime + start=(/1,1,1,lt/) + ! Extraction of the "source" variables +#ifdef NC_DOUBLE + ierr = NF_GET_VARA_DOUBLE(nid,varid,start,size,sum3d) + ierr = NF_GET_VARA_DOUBLE(nid,varid+1,start,size,square3d) +#else + ierr = NF_GET_VARA_REAL(nid,varid,start,size,sum3d) + ierr = NF_GET_VARA_REAL(nid,varid+1,start,size,square3d) +#endif + ! Calculation of these variables + mean3d=sum3d/count(lt) + sd3d=sqrt(max(0.,square3d/count(lt)-mean3d**2)) + ! Writing of the variables +#ifdef NC_DOUBLE + ierr = NF_PUT_VARA_DOUBLE(nid,varid,start,size,mean3d) + ierr = NF_PUT_VARA_DOUBLE(nid,varid+1,start,size,sd3d) +#else + ierr = NF_PUT_VARA_REAL(nid,varid,start,size,mean3d) + ierr = NF_PUT_VARA_REAL(nid,varid+1,start,size,sd3d) +#endif + enddo + + else if (ndims.eq.3) then + +! dimout(1)=lonid +! dimout(2)=latid +! dimout(3)=timeid + + size=(/iip1,jjp1,1,0/) + do lt=1,istime + start=(/1,1,lt,0/) + ! Extraction of the "source" variables +#ifdef NC_DOUBLE + ierr = NF_GET_VARA_DOUBLE(nid,varid,start,size,sum2d) + ierr = NF_GET_VARA_DOUBLE(nid,varid+1,start,size,square2d) +#else + ierr = NF_GET_VARA_REAL(nid,varid,start,size,sum2d) + ierr = NF_GET_VARA_REAL(nid,varid+1,start,size,square2d) +#endif + ! Calculation of these variables + mean2d=sum2d/count(lt) + sd2d=sqrt(max(0.,square2d/count(lt)-mean2d**2)) + ! Writing of the variables +#ifdef NC_DOUBLE + ierr = NF_PUT_VARA_DOUBLE(nid,varid,start,size,mean2d) + ierr = NF_PUT_VARA_DOUBLE(nid,varid+1,start,size,sd2d) +#else + ierr = NF_PUT_VARA_REAL(nid,varid,start,size,mean2d) + ierr = NF_PUT_VARA_REAL(nid,varid+1,start,size,sd2d) +#endif + enddo + + endif +enddo + +ierr= NF_CLOSE(nid) + +end Index: LMDZ5/branches/LMDZ5V2.0-dev/libf/phylmd/physiq.F =================================================================== --- LMDZ5/branches/LMDZ5V2.0-dev/libf/phylmd/physiq.F (revision 1485) +++ LMDZ5/branches/LMDZ5V2.0-dev/libf/phylmd/physiq.F (revision 1486) @@ -158,4 +158,8 @@ save ok_LES c$OMP THREADPRIVATE(ok_LES) +c + LOGICAL callstats ! sortir le fichier stats + save callstats +c$OMP THREADPRIVATE(callstats) c LOGICAL ok_region ! sortir le fichier regional @@ -1150,4 +1154,5 @@ ! and 360 + INTEGER ierr #include "YOMCST.h" #include "YOETHF.h" @@ -1222,4 +1227,5 @@ . ok_instan, ok_hf, . ok_LES, + . callstats, . solarlong0,seuil_inversion, . fact_cldcon, facttemps,ok_newmicro,iflag_radia, @@ -3695,5 +3701,58 @@ c==================================================================== c - + +c ----------------------------------------------------------------- +c WSTATS: Saving statistics +c ----------------------------------------------------------------- +c ("stats" stores and accumulates 8 key variables in file "stats.nc" +c which can later be used to make the statistic files of the run: +c "stats") only possible in 3D runs ! + + + IF (callstats) THEN + + call wstats(klon,o_psol%name,"Surface pressure","Pa" + & ,2,paprs(:,1)) + call wstats(klon,o_tsol%name,"Surface temperature","K", + & 2,zxtsol) + zx_tmp_fi2d(:) = rain_fall(:) + snow_fall(:) + call wstats(klon,o_precip%name,"Precip Totale liq+sol", + & "kg/(s*m2)",2,zx_tmp_fi2d) + zx_tmp_fi2d(:) = rain_lsc(:) + snow_lsc(:) + call wstats(klon,o_plul%name,"Large-scale Precip", + & "kg/(s*m2)",2,zx_tmp_fi2d) + zx_tmp_fi2d(:) = rain_con(:) + snow_con(:) + call wstats(klon,o_pluc%name,"Convective Precip", + & "kg/(s*m2)",2,zx_tmp_fi2d) + call wstats(klon,o_sols%name,"Solar rad. at surf.", + & "W/m2",2,solsw) + call wstats(klon,o_soll%name,"IR rad. at surf.", + & "W/m2",2,sollw) + zx_tmp_fi2d(:) = topsw(:)-toplw(:) + call wstats(klon,o_nettop%name,"Net dn radiatif flux at TOA", + & "W/m2",2,zx_tmp_fi2d) + + + + call wstats(klon,o_temp%name,"Air temperature","K", + & 3,t_seri) + call wstats(klon,o_vitu%name,"Zonal wind","m.s-1", + & 3,u_seri) + call wstats(klon,o_vitv%name,"Meridional wind", + & "m.s-1",3,v_seri) + call wstats(klon,o_vitw%name,"Vertical wind", + & "m.s-1",3,omega) + call wstats(klon,o_ovap%name,"Specific humidity", "kg/kg", + & 3,q_seri) + + + + IF(lafin) THEN + write (*,*) "Writing stats..." + call mkstats(ierr) + ENDIF + + ENDIF !if callstats + IF (lafin) THEN Index: LMDZ5/branches/LMDZ5V2.0-dev/libf/phylmd/statto.h =================================================================== --- LMDZ5/branches/LMDZ5V2.0-dev/libf/phylmd/statto.h (revision 1486) +++ LMDZ5/branches/LMDZ5V2.0-dev/libf/phylmd/statto.h (revision 1486) @@ -0,0 +1,34 @@ + +! statto: +! This include file controls the production of statistics. +! Some variables could be set in a namelist, but it is easier to +! do it here since arrays can then be dimensioned using parameters +! and values shouldn't have to change too often. SRL + +! Calculate stats every istats physics timesteps, starting at first +! call. If istats=0 then don't do statistics at all. Check value +! if number of physics timesteps changes. + integer istats + +! Calculate itime independent sums and sums of squares, +! example, istat=1,istime=1 gives a single time mean + integer, parameter :: istime=12 + +! Number of 2D and 3D variables on which to do statistics. + integer n2dvar, n3dvar + parameter (n2dvar = 8, n3dvar = 5) + +! Units for writing stats header and data + integer usdata + +! count tab to know the variable record + integer count(istime) + +! Record of the number of stores made for each time. + integer nstore(istime) + +! Size of the "controle" array + integer, parameter :: cntrlsize=15 + +! common /sttcom/ dummy,nstore,istats,usdata + common /sttcom/ nstore,istats,usdata,count Index: LMDZ5/branches/LMDZ5V2.0-dev/libf/phylmd/wstats.F90 =================================================================== --- LMDZ5/branches/LMDZ5V2.0-dev/libf/phylmd/wstats.F90 (revision 1486) +++ LMDZ5/branches/LMDZ5V2.0-dev/libf/phylmd/wstats.F90 (revision 1486) @@ -0,0 +1,351 @@ +subroutine wstats(ngrid,nom,titre,unite,dim,px) + +implicit none + +#include "dimensions.h" +#include "statto.h" +#include "netcdf.inc" + +integer,intent(in) :: ngrid +character (len=*),intent(in) :: nom,titre,unite +integer,intent(in) :: dim +real, dimension(ngrid,llm),intent(in) :: px +integer,parameter :: iip1=iim+1 +integer,parameter :: jjp1=jjm+1 +real, dimension(iip1,jjp1,llm) :: mean3d,sd3d,dx3 +real, dimension(iip1,jjp1) :: mean2d,sd2d,dx2 +character (len=50) :: namebis +character (len=50), save :: firstvar +integer :: ierr,varid,nbdim,nid +integer :: meanid,sdid +integer, dimension(4) :: id,start,size +logical, save :: firstcall=.TRUE. +integer :: l,i,j,ig0 +integer,save :: index + +integer, save :: step=0 + + +if (firstcall) then + firstcall=.false. + firstvar=trim((nom)) + call inistats(ierr) +endif + +if (firstvar==nom) then ! If we're back to the first variable + step = step + 1 +endif + +if (mod(step,istats).ne.0) then + RETURN +endif + +ierr = NF_OPEN("stats.nc",NF_WRITE,nid) + +namebis=trim(nom) +ierr= NF_INQ_VARID(nid,namebis,meanid) + +if (ierr.ne.NF_NOERR) then + + if (firstvar==nom) then + index=1 + count=0 + endif + + +!declaration de la variable + +! choix du nom des coordonnees + ierr= NF_INQ_DIMID(nid,"longitude",id(1)) + ierr= NF_INQ_DIMID(nid,"latitude",id(2)) + if (dim.eq.3) then + ierr= NF_INQ_DIMID(nid,"altitude",id(3)) + ierr= NF_INQ_DIMID(nid,"Time",id(4)) + nbdim=4 + else if (dim==2) then + ierr= NF_INQ_DIMID(nid,"Time",id(3)) + nbdim=3 + endif + + write (*,*) "=====================" + write (*,*) "STATS: creation de ",nom + namebis=trim(nom) + call def_var_stats(nid,namebis,titre,unite,nbdim,id,meanid,ierr) + call inivar(nid,meanid,ngrid,dim,index,px,ierr) + namebis=trim(nom)//"_sd" + call def_var_stats(nid,namebis,trim(titre)//" total standard deviation over the season",unite,nbdim,id,sdid,ierr) + call inivar(nid,sdid,ngrid,dim,index,px,ierr) + + ierr= NF_CLOSE(nid) + return + +else + namebis=trim(nom)//"_sd" + ierr= NF_INQ_VARID(nid,namebis,sdid) + +endif + +if (firstvar==nom) then + count(index)=count(int(index))+1 + index=index+1 + if (index>istime) then + index=1 + endif +endif + +if (count(index)==0) then + if (dim.eq.3) then + start=(/1,1,1,index/) + size=(/iip1,jjp1,llm,1/) + mean3d=0 + sd3d=0 + else if (dim.eq.2) then + start=(/1,1,index,0/) + size=(/iip1,jjp1,1,0/) + mean2d=0 + sd2d=0 + endif +else + if (dim.eq.3) then + start=(/1,1,1,index/) + size=(/iip1,jjp1,llm,1/) +#ifdef NC_DOUBLE + ierr = NF_GET_VARA_DOUBLE(nid,meanid,start,size,mean3d) + ierr = NF_GET_VARA_DOUBLE(nid,sdid,start,size,sd3d) +#else + ierr = NF_GET_VARA_REAL(nid,meanid,start,size,mean3d) + ierr = NF_GET_VARA_REAL(nid,sdid,start,size,sd3d) +#endif + if (ierr.ne.NF_NOERR) then + write (*,*) NF_STRERROR(ierr) + stop "" + endif + + else if (dim.eq.2) then + start=(/1,1,index,0/) + size=(/iip1,jjp1,1,0/) +#ifdef NC_DOUBLE + ierr = NF_GET_VARA_DOUBLE(nid,meanid,start,size,mean2d) + ierr = NF_GET_VARA_DOUBLE(nid,sdid,start,size,sd2d) +#else + ierr = NF_GET_VARA_REAL(nid,meanid,start,size,mean2d) + ierr = NF_GET_VARA_REAL(nid,sdid,start,size,sd2d) +#endif + if (ierr.ne.NF_NOERR) then + write (*,*) NF_STRERROR(ierr) + stop "" + endif + endif +endif + +if (dim.eq.3) then + +! Passage variable physique --> variable dynamique + + DO l=1,llm + DO i=1,iip1 + dx3(i,1,l)=px(1,l) + dx3(i,jjp1,l)=px(ngrid,l) + ENDDO + DO j=2,jjm + ig0= 1+(j-2)*iim + DO i=1,iim + dx3(i,j,l)=px(ig0+i,l) + ENDDO + dx3(iip1,j,l)=dx3(1,j,l) + ENDDO + ENDDO + + mean3d(:,:,:)= mean3d(:,:,:)+dx3(:,:,:) + sd3d(:,:,:)= sd3d(:,:,:)+dx3(:,:,:)**2 + +#ifdef NC_DOUBLE + ierr = NF_PUT_VARA_DOUBLE(nid,meanid,start,size,mean3d) + ierr = NF_PUT_VARA_DOUBLE(nid,sdid,start,size,sd3d) +#else + ierr = NF_PUT_VARA_REAL(nid,meanid,start,size,mean3d) + ierr = NF_PUT_VARA_REAL(nid,sdid,start,size,sd3d) +#endif + if (ierr.ne.NF_NOERR) then + write (*,*) NF_STRERROR(ierr) + stop "" + endif + +else if (dim.eq.2) then + +! Passage variable physique --> physique dynamique + + DO i=1,iip1 + dx2(i,1)=px(1,1) + dx2(i,jjp1)=px(ngrid,1) + ENDDO + DO j=2,jjm + ig0= 1+(j-2)*iim + DO i=1,iim + dx2(i,j)=px(ig0+i,1) + ENDDO + dx2(iip1,j)=dx2(1,j) + ENDDO + + mean2d(:,:)= mean2d(:,:)+dx2(:,:) + sd2d(:,:)= sd2d(:,:)+dx2(:,:)**2 + +#ifdef NC_DOUBLE + ierr = NF_PUT_VARA_DOUBLE(nid,meanid,start,size,mean2d) + ierr = NF_PUT_VARA_DOUBLE(nid,sdid,start,size,sd2d) +#else + ierr = NF_PUT_VARA_REAL(nid,meanid,start,size,mean2d) + ierr = NF_PUT_VARA_REAL(nid,sdid,start,size,sd2d) +#endif + if (ierr.ne.NF_NOERR) then + write (*,*) NF_STRERROR(ierr) + stop "" + endif + +endif + +ierr= NF_CLOSE(nid) + +end subroutine wstats + +!====================================================== +subroutine inivar(nid,varid,ngrid,dim,index,px,ierr) + +implicit none + +include "dimensions.h" +!!!!!!include "dimphys.h" +include "netcdf.inc" + +integer, intent(in) :: nid,varid,dim,index,ngrid +real, dimension(ngrid,llm), intent(in) :: px +integer, intent(out) :: ierr + +integer,parameter :: iip1=iim+1 +integer,parameter :: jjp1=jjm+1 + +integer :: l,i,j,ig0 +integer, dimension(4) :: start,size +real, dimension(iip1,jjp1,llm) :: dx3 +real, dimension(iip1,jjp1) :: dx2 + +if (dim.eq.3) then + + start=(/1,1,1,index/) + size=(/iip1,jjp1,llm,1/) + +! Passage variable physique --> variable dynamique + + DO l=1,llm + DO i=1,iip1 + dx3(i,1,l)=px(1,l) + dx3(i,jjp1,l)=px(ngrid,l) + ENDDO + DO j=2,jjm + ig0= 1+(j-2)*iim + DO i=1,iim + dx3(i,j,l)=px(ig0+i,l) + ENDDO + dx3(iip1,j,l)=dx3(1,j,l) + ENDDO + ENDDO + +#ifdef NC_DOUBLE + ierr = NF_PUT_VARA_DOUBLE(nid,varid,start,size,dx3) +#else + ierr = NF_PUT_VARA_REAL(nid,varid,start,size,dx3) +#endif + +else if (dim.eq.2) then + + start=(/1,1,index,0/) + size=(/iip1,jjp1,1,0/) + +! Passage variable physique --> physique dynamique + + DO i=1,iip1 + dx2(i,1)=px(1,1) + dx2(i,jjp1)=px(ngrid,1) + ENDDO + DO j=2,jjm + ig0= 1+(j-2)*iim + DO i=1,iim + dx2(i,j)=px(ig0+i,1) + ENDDO + dx2(iip1,j)=dx2(1,j) + ENDDO + +#ifdef NC_DOUBLE + ierr = NF_PUT_VARA_DOUBLE(nid,varid,start,size,dx2) +#else + ierr = NF_PUT_VARA_REAL(nid,varid,start,size,dx2) +#endif + +endif + +end subroutine inivar + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + +subroutine def_var_stats(nid,name,title,units,nbdim,dimids,nvarid,ierr) + +! This subroutine defines variable 'name' in a (pre-existing and opened) +! NetCDF file (known from its NetCDF ID 'nid'). +! The number of dimensions 'nbdim' of the variable, as well as the IDs of +! corresponding dimensions must be set (in array 'dimids'). +! Upon successfull definition of the variable, 'nvarid' contains the +! NetCDF ID of the variable. +! The variables' attributes 'title' (Note that 'long_name' would be more +! appropriate) and 'units' are also set. + +implicit none + +#include "netcdf.inc" + +integer,intent(in) :: nid ! NetCDF file ID +character(len=*),intent(in) :: name ! the variable's name +character(len=*),intent(in) :: title ! 'title' attribute of variable +character(len=*),intent(in) :: units ! 'units' attribute of variable +integer,intent(in) :: nbdim ! number of dimensions of the variable +integer,dimension(nbdim),intent(in) :: dimids ! NetCDF IDs of the dimensions + ! the variable is defined along +integer,intent(out) :: nvarid ! NetCDF ID of the variable +integer,intent(out) :: ierr ! returned NetCDF staus code + +! 1. Switch to NetCDF define mode +ierr=NF_REDEF(nid) + +! 2. Define the variable +#ifdef NC_DOUBLE +ierr = NF_DEF_VAR (nid,adjustl(name),NF_DOUBLE,nbdim,dimids,nvarid) +#else +ierr = NF_DEF_VAR (nid,adjustl(name),NF_FLOAT,nbdim,dimids,nvarid) +#endif +if(ierr/=NF_NOERR) then + write(*,*) "def_var_stats: Failed defining variable "//trim(name) + write(*,*) NF_STRERROR(ierr) + stop "" +endif + +! 3. Write attributes +ierr=NF_PUT_ATT_TEXT(nid,nvarid,"title",& + len_trim(adjustl(title)),adjustl(title)) +if(ierr/=NF_NOERR) then + write(*,*) "def_var_stats: Failed writing title attribute for "//trim(name) + write(*,*) NF_STRERROR(ierr) + stop "" +endif + +ierr=NF_PUT_ATT_TEXT(nid,nvarid,"units",& + len_trim(adjustl(units)),adjustl(units)) +if(ierr/=NF_NOERR) then + write(*,*) "def_var_stats: Failed writing units attribute for "//trim(name) + write(*,*) NF_STRERROR(ierr) + stop "" +endif + +! 4. Switch out of NetCDF define mode +ierr = NF_ENDDEF(nid) + +end subroutine def_var_stats +