source: branches/GRISLIv3/SOURCES/Alps_files/output_alps_mod.f90

Last change on this file was 484, checked in by aquiquet, 5 months ago

Cleaning branch: small corrections in lect alps/anteis/hemin40
File size: 11.5 KB
Line 
1!> \file output_alps.f90
2!! Module de
3!<
4
5!> \namespace output_alps_mod
6!! Module de
7!! \author ...
8!! \date ...
9!! @note Used module
10!! @note   - use module3D_phy
11!<
12
13module  output_alps_mod
14
15use module3D_phy,only: xlong,ylat,b,s,acc,h,t,tbdot,ux,uy,uzk,uzsdot, &
16              hdot,bdot,num_ritz,tafor,sealevel,ts
17use geography, only: nx,ny,nz,geoplace,dirnameinp
18use runparam, only: dirnameout,runname
19use bilan_eau_mod
20use netcdf
21use io_netcdf_grisli
22use bilan_flux_mod
23
24
25implicit none
26     
27real ::  vol
28integer :: nflot                      !< nbr de point flottant
29real ::  bmean                        !<
30real ::  accmean                      !< accumulation moyenne
31real ::  ablmean                      !< ablation moyenne
32real ::  calvmean                     !< moyenne calving
33real ::  ablbordmean                  !<
34real ::  bmeltmean                    !< moyenne fusion basale
35real ::  tbmean                       !< temperature basale moyenne
36real ::  tbdotmean                    !< moyenne variation / temps temperature basale
37real ::  vsmean                       !< vitesse de surface moyenne
38!real ::  vsdotmean                    !< moyenne variation / temps vitesse de surface
39!real ::  uzsmean   !!!! utilise ?     !< vitesse verticale de surface moyenne
40real ::  uzsdotmean                   !< moyenne variation / temps vitesse verticale de surface
41real ::  uzkmean                      !< moyenne vitesse verticale de surface
42real ::  hdotmean                     !< moyenne derivee / temps de h
43real ::  bdotmean                     !< moyenne bedrock derive / temps
44real ::  volf                         !< volume au dessus de la flottaison
45
46real :: lim                           !< Total ice mass
47real :: iareag                        !< surface posee
48real :: iareaf                        !< surface flottante
49real :: tendacabf                     !< Total SMB flux
50real :: tendlibmassbf                 !< Total Basal mass balance flux
51real :: tendlibmassbffl               !< Total Basal mass balance flux on floating ice
52real :: tendlicalvf                   !< Total calving flux
53real :: tendlifmassbf                 !< Total calving and ice front melting flux
54real :: tendligroundf                 !< Total grounding line flux
55
56real,dimension(nx,ny) :: corrsurf     !< facteur de correction de la surface
57!real,parameter :: ice_density=910.    !< densite de la glace pour conversion en masse
58
59! variables netcdf
60integer,parameter :: ncshortout=1 ! 1 sorties netcdf short initMIP
61integer,parameter :: nvar=10 ! nombre de variables dans le fichier de sortie temporel netcdf
62integer :: ncid
63integer :: status
64integer :: timeDimID
65integer :: timeVarID
66integer,dimension(nvar) :: varID
67integer :: nbtimeout  ! index time output
68
69real,dimension(nvar) :: var_shortoutput
70
71
72CONTAINS
73
74subroutine init_outshort
75
76double precision,dimension(:,:),pointer      :: tab               !< tableau 2d real pointer
77character(len=100),dimension(nvar) :: namevar ! name, standard_name, long_name, units
78character(len=100),dimension(nvar) :: standard_name
79character(len=100),dimension(nvar) :: long_name
80character(len=100),dimension(nvar) :: units
81
82integer :: k
83
84
85if (ncshortout.eq.1) then  ! ecriture netcdf
86
87! lecture du fichier avec les corrections de surface
88  if (geoplace.eq.'ant16km') then
89      call Read_Ncdf_var('z',trim(DIRNAMEINP)//'/corrsurf-initMIP-16km.grd',tab)
90      corrsurf(:,:)=tab(:,:)
91  else
92      corrsurf(:,:)= 1.
93  end if
94
95  open(568,file=trim(dirsource)//'/Fichiers-parametres/short-ISMIPnc.dat',status='old')
96! lecture en-tete
97  read(568,*)
98  read(568,*)
99  read(568,*)
100  read(568,*)
101  read(568,*)
102! lecture des infos sur les variables :
103  do k=1,nvar
104    read(568,'(a100)') namevar(k)
105    read(568,'(a100)') standard_name(k)
106    read(568,'(a100)') long_name(k)
107    read(568,'(a100)') units(k)
108    read(568,*)
109  enddo
110  close(568)
111! Fichier Netcdf initMIP
112! creation du fichier Netcdf :
113  status=nf90_create(path = trim(dirnameout)//'short'//runname//'.nc', cmode = nf90_clobber, ncid = ncid)
114  if (status /= nf90_noerr) call handle_err(status)
115
116! definition des dimension :
117  status=nf90_def_dim(ncid, name="time", len=NF90_UNLIMITED, dimid=timeDimID)
118  if (status /= nf90_noerr) call handle_err(status)
119  status=nf90_def_var(ncid, name="time", xtype=nf90_float, dimids=(/ timeDimID/), varid=timeVarID) 
120  if (status /= nf90_noerr) call handle_err(status)
121  status=nf90_put_att(ncid, timeVarID, "standard_name", "time")
122  if (status /= nf90_noerr) call handle_err(status)
123  status=nf90_put_att(ncid, timeVarID,"units", "years since 1995-01-01 00:00:00")
124  if (status /= nf90_noerr) call handle_err(status)
125  status=nf90_put_att(ncid, timeVarID,"calendar", "360_day")
126  if (status /= nf90_noerr) call handle_err(status)
127 
128! definition des variables de sortie :
129  do k=1,nvar  ! boucle sur le nbr de variable a definir
130    status=nf90_def_var(ncid, name=trim(namevar(k)), xtype=nf90_float, dimids= &
131          (/ timeDimID /), varid=varID(k))
132    if (status /= nf90_noerr) call handle_err(status)
133    status=nf90_put_att(ncid, varID(k), "standard_name", trim(standard_name(k)))
134    if (status /= nf90_noerr) call handle_err(status)
135    status=nf90_put_att(ncid, varID(k), "long_name", trim(long_name(k)))
136    if (status /= nf90_noerr) call handle_err(status)
137    status=nf90_put_att(ncid, varID(k), "units", trim(units(k)))
138    if (status /= nf90_noerr) call handle_err(status)
139  enddo
140
141! fin de la definition du fchier :
142  status=nf90_enddef(ncid)
143  if (status /= nf90_noerr) call handle_err(status)
144  nbtimeout = 0 ! initialisation compteur sorties axe time
145else ! pas de sortie netcdf et sans correction de surface
146  corrsurf(:,:)=1.
147endif
148
149end subroutine init_outshort
150
151
152
153!_________________________________________________________________________
154subroutine shortoutput
155
156! 1_initialization
157!------------------
158real ::  smax,hmax,hmean
159integer :: np,i,j,k
160
161
162      vol=0. 
163      np=0
164      nflot=0
165      hmax=0. 
166      smax=0.
167      bmean=0. 
168      accmean=0. 
169      ablmean=0. 
170      calvmean=0. 
171      ablbordmean=0.
172      bmeltmean=0.
173      tbmean=0.
174      tbdotmean=0.
175      vsmean=0.
176!      vsdotmean=0.
177!      uzsmean=0.
178      uzsdotmean=0.
179      uzkmean=0.
180      hdotmean=0.
181      bdotmean=0.
182      volf=0.
183      lim=0.
184      tendacabf=0.
185      iareag=0.
186      iareaf=0.
187      tendlicalvf=0.
188      tendligroundf=0.
189      tendlibmassbf=0.
190      tendlibmassbffl=0.
191      tendlifmassbf=0.
192! 2_preparing outputs
193!--------------------     
194    do j=1,ny
195      do i=1,nx
196        if (ice(i,j).eq.1) then ! point englace
197          if (.not.flot(i,j)) then ! point pose
198            np=np+1
199            vol=vol+h(i,j)
200            iareag=iareag+1.*corrsurf(i,j)                            ! surface englacee posee
201
202!         calcul de la hauteur au dessus de la flottaison
203            if (sealevel_2d(i,j)-B(i,j).le.0.) then    ! socle au dessus du niveau des mers
204              volf=volf+h(i,j)*corrsurf(i,j)                          ! volume au-dessus de la flottaison
205            else
206              volf=volf+(h(i,j)-row/ro*(sealevel_2d(i,j)-b(i,j)))*corrsurf(i,j) ! volume au-dessus de la flottaison
207            endif
208
209            if (h(i,j).gt.hmax) hmax=h(i,j) 
210            if (s(i,j).gt.smax) smax=s(i,j) 
211            bmean=bm(i,j)+bmean 
212            accmean=acc(i,j)+accmean
213            tbmean=tbmean+t(i,j,nz)
214            tbdotmean=tbdotmean+tbdot(i,j)
215            vsmean=vsmean+sqrt(ux(i,j,1)**2+uy(i,j,1)**2)
216!          vsdotmean=vsdotmean+vsdot(i,j)
217!          uzsmean=uzsmean+uz(i,j,1)
218            uzsdotmean=uzsdotmean+uzsdot(i,j)
219            uzkmean=uzkmean+uzk(i,j)
220            hdotmean=hdotmean+abs(hdot(i,j)) 
221            bdotmean=bdotmean+abs(bdot(i,j)) 
222            bmeltmean=bmeltmean+bmelt(i,j)
223          else ! point flottant
224            iareaf=iareaf+1.*corrsurf(i,j)                               ! surface flottante
225            tendlibmassbffl=tendlibmassbffl-bmelt_dtt(i,j)*corrsurf(i,j)/dtt_flux   ! fonte basale sur zone flottante
226          endif
227          lim=lim+h(i,j)*corrsurf(i,j)                                   ! volume total de glace
228          tendacabf=tendacabf+bm_dtt(i,j)*corrsurf(i,j)/dtt_flux              ! smb surface
229          tendlibmassbf=tendlibmassbf-bmelt_dtt(i,j)*corrsurf(i,j)/dtt_flux   ! fonte basale     
230        endif
231        tendlicalvf=tendlicalvf+calv_dtt(i,j)*corrsurf(i,j)/dtt_flux        ! calving
232        tendlifmassbf=tendlifmassbf+(calv_dtt(i,j)-ablbord_dtt(i,j))*corrsurf(i,j)/dtt_flux ! calving + ablbord (ablbord est positif)
233        tendligroundf=tendligroundf+grline_dtt(i,j)*corrsurf(i,j)/dtt_flux  ! flux a la grounding line
234      end do
235    end do
236
237
238      if (np.ne.0) then
239        hmean=vol/np 
240        vol=vol*dx*dy 
241        volf=volf*dx*dy*ice_density
242        bmean=bmean/np 
243        accmean=accmean/np 
244        ablmean=bmean-accmean 
245        calvmean=calvmean/np 
246        bmeltmean=bmeltmean/np
247        ablbordmean=ablbordmean/np
248        tbmean=tbmean/np
249        tbdotmean=tbdotmean/np
250        vsmean=vsmean/np
251!        vsdotmean=vsdotmean/np
252!        uzsmean=uzsmean/np
253        uzsdotmean=uzsdotmean/np
254        uzkmean=uzkmean/np
255        hdotmean=hdotmean/np
256      endif
257      lim=lim*dx*dy*ice_density
258      iareag=iareag*dx*dy 
259      iareaf=iareaf*dx*dy
260      tendacabf=tendacabf*dx*dy*ice_density/secyear
261      tendlibmassbf=tendlibmassbf*dx*dy*ice_density/secyear
262      tendlibmassbffl=tendlibmassbffl*dx*dy*ice_density/secyear
263      tendlicalvf=tendlicalvf*dx*dy*ice_density/secyear
264      tendlifmassbf=tendlifmassbf*dx*dy*ice_density/secyear
265      tendligroundf=tendligroundf*ice_density/secyear
266     
267     
268      bdotmean=bdotmean/nx/ny 
269
270
271! 2_writing outputs
272!------------------   
273!     **** short display ****
274
275        write(num_ritz,904) nt,time,tafor,sealevel,vol,volf,np, &
276          nint(hmean),nint(smax),                    &
277          bmean,tbmean,nint(vsmean),                 &
278                                  tbdotmean,hdotmean,dt,accmean, &
279                                  diff_H,water_bilan,sum(calv_dtt(:,:))/dtt_flux, &
280                                  sum(ablbord_dtt(:,:))/dtt_flux, &
281                                  sum(Bm_dtt(:,:))/dtt_flux,sum(bmelt_dtt(:,:))/dtt_flux
282
283
284903   format(i8,1x,f0.2,1x,f0.4,1x,f0.2,1x,2(e10.5,1x),i6,1x,i4,1x,i5,1x, &
285             f0.4,1x,f0.3,1x,i3,4(1x,e8.2),1x,f0.4,1x,f0.4) 
286904   format(i8,1x,f0.2,4(1x,e15.8),3(1x,i8),2(1x,e15.8),1x,i8,10(1x,e15.8))
287!940   format('%%%% ',a,'   time=',f8.0,' %%%%')
288
289
290if (ncshortout.eq.1) then  ! ecriture netcdf
291  ! Total ice mass
292  var_shortoutput(1)=lim
293  ! Mass above floatation
294  var_shortoutput(2)=volf
295  ! Grounded ice area
296  var_shortoutput(3)=iareag
297  ! Floating ice area
298  var_shortoutput(4)=iareaf
299  ! Total SMB flux
300  var_shortoutput(5)=tendacabf
301  ! Total Basal mass balance flux
302  var_shortoutput(6)=tendlibmassbf
303  ! Total Basal mass balance flux beneath floating ice
304  var_shortoutput(7)=tendlibmassbffl
305  ! Total calving flux
306  var_shortoutput(8)=tendlicalvf
307  ! Total calving and ice front melting flux
308  var_shortoutput(9)=tendlifmassbf
309  ! Total grounding line flux
310  var_shortoutput(10)=tendligroundf
311
312  nbtimeout=nbtimeout+1
313 
314  status=nf90_put_var(ncid, timeVarID, time, start=(/nbtimeout/))
315  if (status /= nf90_noerr) call handle_err(status)
316
317  do k=1,nvar  ! boucle sur le nbr de variable a ecrire
318    status=nf90_put_var(ncid, varID(k), var_shortoutput(k),start=(/nbtimeout/))
319    if (status /= nf90_noerr) call handle_err(status)
320  enddo
321  status=nf90_sync(ncid)
322  if (status /= nf90_noerr) call handle_err(status)
323endif
324
325
326end subroutine shortoutput
327
328subroutine handle_err(status)
329  integer, intent(in) :: status
330  if (status /= nf90_noerr) then
331     print*,trim(nf90_strerror(status))
332     stop "stopped"
333  end if
334end subroutine handle_err
335
336end module  output_alps_mod
Note: See TracBrowser for help on using the repository browser.