source: trunk/SOURCES/GrIce2sea_files/output_Grice2sea_mod.f90 @ 55

Last change on this file since 55 was 11, checked in by dumas, 9 years ago

Grice2sea compilé et validé avec le module climat_Grice2sea_years_mod. climat_GrIce2sea_years_mod.f90 inclus massb_Ice2sea_RCM et massb_Ice2sea_fixe. pdd_declar_mod.f90 supprimé, les déclarations de variables concernant le pdd sont maintenant dans le module ablation_mod.f90.
File size: 12.2 KB
Line 
1!> \file GrIce2sea_files/output_Grice2sea-regions_mod.f90
2!> \namespace output_grice2sea_mod
3
4!! Module de sortie temporelle Grice2sea avec regions
5!! \author Cat
6!! \date ...
7!! @note Used module
8!! @note   - use module3D_phy
9!! @note   - use param_phy_mod
10!<
11
12
13module output_grice2sea_mod
14
15  USE module3D_phy
16  USE param_phy_mod
17  use io_netcdf_grisli
18
19  implicit none
20     
21!real ::  vol ; !integer :: np
22  real ::  bmean                        !<
23  real ::  accmean                      !< accumulation moyenne
24  real ::  ablmean                      !< ablation moyenne
25  real ::  calvmean                     !< moyenne calving
26  real ::  ablbordmean                  !<
27  real ::  bmeltmean                    !< moyenne fusion basale
28  real ::  tbmean                       !< temperature basale moyenne
29  real ::  tbdotmean                    !< moyenne variation / temps temperature basale
30  real ::  vsmean                       !< vitesse de surface moyenne
31  real ::  uzsdotmean                   !< moyenne variation / temps vitesse verticale de surface
32  real ::  uzkmean                      !< moyenne vitesse verticale de surface
33  real ::  hdotmean                     !< moyenne derivee / temps de h
34  real ::  bdotmean                     !< moyenne bedrock derive / temps
35  real ::  volf                         !< volume au dessus de la flottaison
36 
37 
38! dimensionnement pour les bassins ou regions
39 
40  integer,parameter              :: nbreg_max=20      !< nombre maxi de regions
41  integer                        :: nbregions         !< nombre de region
42  integer                        :: i_region          !< index des regions
43  integer                        :: scal_np           !< nombre de points grounded par region
44  integer, dimension(nbreg_max)  :: np_reg            !< nombre de points grounded par region
45  integer, dimension(nx,ny)      :: mask_region       !< carte des region
46  integer, dimension(nbreg_max)  :: nbpoint_reg       !< nb de points par region
47 
48  real, dimension(nbreg_max)     :: vol_reg           !< volumes par region (points grounded)
49  real, dimension(nbreg_max)     :: volbuoy_reg       !< volumes au dessus flottaison par region
50  real, dimension(nbreg_max)     :: meanhdot_reg      !< hdot moyens par region
51  real, dimension(nbreg_max)     :: sigma_hdot_reg    !< sigma hdot  par region
52  real, dimension(nbreg_max)     :: smb_reg           !< bilan de masse en surface
53 
54  real                           :: scal_vol           !< volumes par region (points grounded)
55  real                           :: scal_volbuoy       !< volumes au dessus flottaison par region
56  real                           :: scal_meanhdot      !< hdot moyens par region
57  real                           :: scal_sigma_hdot    !< sigma hdot  par region
58  real                           :: scal_smb           ! bilan de mass en surface par region
59 
60 
61  integer                        :: nvol        = 250 !< numero du fichier volume     
62  integer                        :: nvolbuoy    = 251 !< numero du fichier volume au dessus flottaison   
63  integer                        :: nhdot       = 252 !< numero du fichier hdot mean 
64  integer                        :: nsigma_hdot = 253 !< numero du fichier hdot sigma
65  integer                        :: npoints     = 254 !< numero du fichier nb de points
66  integer                        :: nsmb        = 255 !< numero du fichier smb
67 
68 
69contains
70 
71!_______________________________________________________________________________
72!< initialise les sorties
73 
74  subroutine init_outshort
75   
76    implicit none
77   
78    real*8, dimension(:,:),   pointer  :: work_tab    !< tableau 2d real ecrit dans le fichier   
79    character(len=100)    :: file_ncdf      !< fichier netcdf issue des fichiers .dat
80   
81   
82   
83   
84    character(len=100) :: region_file                  !  Fichier dans lequel sont les regions
85    character(len=100) :: filout                       !  Pour fichier ecriture
86    namelist/output_regions/nbregions,region_file
87   
88    if (itracebug.eq.1)  call tracebug('output_grice2sea_mod   subroutine init_out_short')
89   
90   
91    ndisp = 1                                          ! ndisp sortie courte tous les ndisp
92   
93   
94!    lecture des regions     -----------------------------------------------------
95   
96    rewind(num_param)        ! pour revenir au debut du fichier param_list.dat
97428 format(A)
98    read(num_param,output_regions)
99   
100    write(num_rep_42,428) '!___________________________________________________________' 
101    write(num_rep_42,428) '!  sorties temporelles par regions'
102    write(num_rep_42,output_regions)                   
103    write(num_rep_42,428) '! nbregions = nombre de regions'
104    write(num_rep_42,428) '! region_file = carte des regions '
105    write(num_rep_42,428) '!___________________________________________________________' 
106   
107   
108! lecture du masque region
109    file_ncdf   = trim(dirnameinp)//'test_map_region'//'.nc'
110    region_file = trim(dirnameinp)//region_file
111    write(6,*) file_ncdf
112    write(6,*) region_file
113!    call lect_input(1,'mask_region',1,work_tab,region_file,file_ncdf)
114    call Read_Ncdf_var('z',region_file,work_tab)
115
116
117    mask_region(:,:) = work_tab(:,:)
118    if (mask_region((nx-1)/2,(ny-1)/2).eq.0) then
119       write(6,*) 'mauvaise lecture du fichier region'
120       STOP
121    end if
122   
123   
124! ouverture des fichiers par type de variable, toutes les regions dans le meme fichier
125   
126    filout   = runname//'_vol_regions.dat'
127    filout    = TRIM(DIRNAMEOUT)//TRIM(filout)
128    open(nvol,file=filout)
129   
130    filout   = runname// '_volbuoy_regions.dat'
131    filout    = TRIM(DIRNAMEOUT)//TRIM(filout)
132    open(nvolbuoy,file=filout)
133   
134    filout   = runname//'_meanhdot_regions.dat'
135    filout    = TRIM(DIRNAMEOUT)//TRIM(filout)
136    open(nhdot,file=filout)
137   
138    filout   = runname//'_sigmahdot_regions.dat'
139    filout    = TRIM(DIRNAMEOUT)//TRIM(filout)
140    open(nsigma_hdot,file=filout)
141   
142    filout   = runname//'_nbpoints_regions.dat'
143    filout    = TRIM(DIRNAMEOUT)//TRIM(filout)
144    open(npoints,file=filout)
145   
146    filout   = runname//'_smb_regions.dat'
147    filout    = TRIM(DIRNAMEOUT)//TRIM(filout)
148    open(nsmb,file=filout)
149   
150  end subroutine init_outshort
151 
152 
153 
154!_________________________________________________________________________
155  subroutine shortoutput
156
157    implicit none
158    real ::  smax
159
160
161    if (itracebug.eq.1)  call tracebug('output_grice2sea_mod   subroutine shortoutput')
162
163   
164! Sorties pour toute la calotte (ancienne version)
165   
166! 1_initialization
167!------------------
168
169    vol=0. 
170    np=0
171    hmax=0. 
172    smax=0.
173    bmean=0. 
174    accmean=0. 
175    ablmean=0. 
176    calvmean=0. 
177    ablbordmean=0.
178    bmeltmean=0.
179    tbmean=0.
180    tbdotmean=0.
181    vsmean=0.
182    uzsdotmean=0.
183    uzkmean=0.
184    hdotmean=0.
185    bdotmean=0.
186    volf=0.
187   
188! 2_preparing outputs
189!--------------------     
190    do i=1,nx 
191       do j=1,ny
192          if (.not.flot(i,j)) then 
193             !       if (h(i,j).gt.1.) then
194             np=np+1
195             vol=vol+h(i,j) 
196             
197             !        calcul de la hauteur au dessus de la flottaison
198             if (sealevel-B(i,j).le.0.) then    ! socle au dessus du niveau des mers
199                volf=volf+h(i,j)
200             else
201                volf=volf+h(i,j)-row/ro*(sealevel-b(i,j))
202             endif
203             
204             
205             if (h(i,j).gt.hmax) hmax=h(i,j) 
206             if (s(i,j).gt.smax) smax=s(i,j) 
207             bmean=bm(i,j)+bmean 
208             accmean=acc(i,j)+accmean
209             tbmean=tbmean+t(i,j,nz)
210             tbdotmean=tbdotmean+tbdot(i,j)
211             vsmean=vsmean+sqrt(ux(i,j,1)**2+uy(i,j,1)**2)
212             uzsdotmean=uzsdotmean+uzsdot(i,j)
213             uzkmean=uzkmean+uzk(i,j)
214             hdotmean=hdotmean+abs(hdot(i,j)) 
215             bdotmean=bdotmean+abs(bdot(i,j)) 
216             bmeltmean=bmeltmean+bmelt(i,j)
217          endif
218          calvmean=calvmean+calv(i,j) 
219          ablbordmean=ablbordmean+ablbord(i,j)
220       end do
221    end do
222   
223   
224    if (np.ne.0) then
225       hmean=vol/np 
226       vol=vol*dx*dy 
227       volf=volf*dx*dy
228       bmean=bmean/np 
229       accmean=accmean/np 
230       ablmean=bmean-accmean 
231       calvmean=calvmean/np 
232       bmeltmean=bmeltmean/np
233       ablbordmean=ablbordmean/np
234       tbmean=tbmean/np
235       tbdotmean=tbdotmean/np
236       vsmean=vsmean/np
237       !        vsdotmean=vsdotmean/np
238       !        uzsmean=uzsmean/np
239       uzsdotmean=uzsdotmean/np
240       uzkmean=uzkmean/np
241       hdotmean=hdotmean/np 
242    endif
243   
244    bdotmean=bdotmean/nx/ny 
245   
246   
247! 2_writing outputs
248!------------------   
249!     **** short display ****
250   
251    write(num_ritz,903) nt,time,tafor,sealevel,vol,volf,np, &
252         nint(hmean),nint(smax),                    &
253         bmean,tbmean,nint(vsmean),                 &
254         !         tbdotmean,vsdotmean,hdotmean,bdotmean,    &
255         !          tbdotmean,hdotmean,dt,bmeltmean,accmean
256         
257         tbdotmean,hdotmean,dt,bmelt(3,3),accmean 
258   
259   
260903 format(i8,1x,f0.2,1x,f0.4,1x,f0.2,1x,2(e15.8,1x),i0,1x,i0,1x,i0,1x, &
261         f0.4,1x,f0.3,1x,i0,4(1x,e11.4),1x,f0.4,1x,f0.4) 
262    !940   format('%%%% ',a,'   time=',f8.0,' %%%%')
263   
264   
265    call output_regions
266   
267   
268  end subroutine shortoutput
269 
270!----------------------------------------------------------------------------
271  subroutine output_regions
272   
273    implicit none
274    character(len=100)    :: fmt1,fmt2,fmt3
275    character(len=2)      :: string_region
276    if (itracebug.eq.1)  call tracebug('output_grice2sea_mod   subroutine output_regions')
277
278
279    region:  do  i_region = 1, nbregions
280       
281! initialisation
282       scal_vol        = 0.
283       scal_volbuoy    = 0.
284       scal_meanhdot   = 0.
285       scal_sigma_hdot = 0.
286       scal_np         = 0
287       scal_smb        = 0.
288       
289       do i=1,nx 
290          do j=1,ny
291!             write(165,*) i_region,i,j,mask_region(i,j),flot(i,j)
292
293! Attention, pour le Groenland le test porte aussi sur l'epaisseur
294! pour le smb, on perd sans doute la demi maille des points ou H=0
295
296             if ((mask_region(i,j).eq.i_region).and.(.not.flot(i,j)).and.H(i,j).gt.1.) then
297 
298                scal_np         =  scal_np + 1
299                scal_vol        =  scal_vol + H(i,j)
300                scal_volbuoy    =  scal_volbuoy + H(i,j) - row/ro * max(0.,sealevel-b(i,j))
301                scal_meanhdot   =  scal_meanhdot + Hdot(i,j)
302                scal_sigma_hdot =  scal_sigma_hdot + Hdot(i,j)**2
303                scal_smb        =  scal_smb + bm(i,j)
304             end if
305          end do
306       end do
307       
308       np_reg (i_region)         = scal_np
309       vol_reg(i_region)         = scal_vol*dx*dx
310       volbuoy_reg(i_region)     = scal_volbuoy*dx*dx
311       smb_reg(i_region)         = scal_smb*dx*dx
312
313       if (scal_np.ne.0) then
314          meanhdot_reg(i_region)    = scal_meanhdot/scal_np
315          sigma_hdot_reg(i_region)  = (scal_sigma_hdot/scal_np)**0.5
316       else
317          write(6,*) 'region',i_region,'  n a aucun point grounded'
318       end if
319    end do region
320   
321   
322! ecriture dans les fichiers
323    write(string_region,"(i0)")  nbregions+1
324    fmt1   = '(f10.1,4x,'//string_region//'(es15.8,1x))'
325    fmt2   = '(f10.1,4x,'//string_region//'(i8,1x))'
326
327    write(string_region,"(i0)")  nbregions
328    fmt3   = '(f10.1,4x,'//string_region//'(es15.8,1x))'
329!    write(6,*) fmt1
330!    write(6,*) fmt2
331
332
333    write(nvol,fmt1)        time, sum(vol_reg(1:nbregions)),         &
334                                     (vol_reg(i_region)       ,  i_region = 1,nbregions) 
335    write(nvolbuoy,fmt1)    time, sum(volbuoy_reg(1:nbregions)),     &
336                                     (volbuoy_reg(i_region)   ,  i_region = 1,nbregions) 
337    write(nsmb,fmt1)        time, sum(smb_reg(1:nbregions)),     &
338                                     (smb_reg(i_region)   ,  i_region = 1,nbregions) 
339    write(nhdot,fmt1)       time, sum(meanhdot_reg(1:nbregions)),    &
340                                     (meanhdot_reg(i_region)  ,  i_region = 1,nbregions) 
341    write(npoints,fmt2)     time, sum(np_reg(1:nbregions)),          &
342                                     (np_reg(i_region)        ,  i_region = 1,nbregions) 
343
344    write(nsigma_hdot,fmt3) time, (sigma_hdot_reg(i_region),  i_region = 1,nbregions) 
345   
346   
347  end subroutine output_regions
348 
349end module output_grice2sea_mod
Note: See TracBrowser for help on using the repository browser.