source: trunk/SOURCES/Ant40_files/output_anta40_mod-0.4.f90 @ 237

Last change on this file since 237 was 237, checked in by aquiquet, 6 years ago

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