source: trunk/SOURCES/Alps_files/output_alps_mod.f90 @ 356

Last change on this file since 356 was 356, checked in by aquiquet, 2 years ago

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