source: trunk/SOURCES/Ant40_files/bmelt-ant-regions-oce_mod.f90 @ 159

Last change on this file since 159 was 34, checked in by dumas, 8 years ago

Ajout du module Ant40_files/bmelt-ant-regions-oce_mod.f90 pour forcer la fonte basale avec des fichiers de temperature oceanique (M2 Ning) | lecture fichier netcdf pour Antarctique | Ajout dans les sortie de ghf : flux geothermique | climat_GrIce2sea_years_mod.f90 : lecture de smb et Tann au lieu de z | correction bug sorties Netdf dans initial-0.3.f90

File size: 12.1 KB
Line 
1!> \file bmelt-ant-regions-oce_mod.f90
2!! Module qui calcule la fusion basale en fonction de la temp oceanique (grounded ou ice shelves).
3!<
4
5!> \namespace bmelt_ant_region_oce
6!! Module qui calcule la fusion basale (grounded ou ice shelves)
7!! \author C. Dumas
8!! \date juillet
9!! @note pour les ice shelves Antarctique, tient compte des différentes régions
10!! A choisir dans le module_choix
11!! @note Used module
12!! @note   - use module3D_phy
13!<
14MODULE  bmelt_ant_regions_oce         ! cat juillet 2005
15
16       use module3D_phy
17       use interface_input
18       use io_netcdf_grisli
19
20implicit none
21
22REAL,dimension(nx,ny) ::  bmgrz       !< fusion basale a la grounding zone
23real,dimension(nx,ny) ::  bmshelf     !< fusion basale sous shelf
24
25REAL,dimension(nx,ny) :: dist_talu    !< distance du point au talu continental
26REAL,dimension(nx,ny) :: typeshelf    !< Type de shelf Ronne->1 Ross ->2 ....
27                                      !< utilises pour moduler la fusion sous le shelf
28integer, dimension(10) :: region      !< pour écrire dans le fichier param
29character(len=30),dimension(10) :: regname !< nom des régions
30real :: bsupshelf
31integer :: grd                        !< pour une sortie
32
33real :: bmelt_Ross,bmgrz_Ross         !< fusion basale Ross
34real :: bmelt_FRis,bmgrz_FRis         !< fusion basale filchner-Ronne ice shelf
35real :: bmelt_Amery,bmgrz_Amery       !< Amery ice shelves
36real :: bmelt_PIG,bmgrz_PIG           !< Pine Island glacier
37real :: bmelt_Pen,bmgrz_Pen           !< Petits ice shelves Peninsule
38real :: bmelt_other,bmgrz_other       !< Autre ice shelves
39real :: bmelt_talus,bmgrz_talus       !< Au dela du talus continental
40
41real,dimension(nx,ny) :: temp_oce1, temp_oce2 ! temerature oceanique fichier 1 et 2
42real,dimension(nx,ny) :: diff_temp_oce ! variation temperature ocean => pour calcul variation fusion basale
43real :: beta_coefbmelt                ! coef variation fusion basale fonction de variation de temp oceanique (m.yr-1.K-1)
44
45
46CONTAINS
47!-------------------------------------------------------------------------------
48subroutine init_bmelt
49
50  character(len=100) :: file_temp_oce1      ! fichier de ref (ctrl) contenant la carte des temperatures de l'ocean de subsurface
51  character(len=100) :: file_temp_oce2      ! 2eme fichier : c'est la difference qui permet de calculer la variation de fusion basale
52  character(len=100) :: file_ncdf      !< fichier netcdf issue des fichiers .dat
53  integer :: nbr_pts ! nbr de points utilises pour calcul temperature oceanique sur points sans valeurs
54  logical,dimension(nx,ny) :: mask_oce1,mask_oce2 ! masque des zones avec temp ocean (true si existe valeur)
55  real :: flpath ! distance de recherche (rayon pour les temp oceanique non presentes dans le fichier)
56  integer :: n,imin,imax,jmin,jmax
57  logical :: test ! pour stoper boucle
58
59! Cette routine fait l'initialisation pour la fusion basale.
60! Elle est appelée par inputfile-vec-0.5.f90
61
62
63namelist/bmelt_ant_reg/bmelt_Ross,bmgrz_Ross,bmelt_FRis,bmgrz_FRis,  & 
64                       bmelt_Amery,bmgrz_Amery,bmelt_PIG,bmgrz_PIG,  &
65                       bmelt_Pen,bmgrz_Pen,bmelt_other,bmgrz_other,  &       
66                       bmelt_talus,bmgrz_talus     
67
68namelist/bmelt_ocean_temp/file_temp_oce1,file_temp_oce2,beta_coefbmelt
69
70
71rewind(num_param)        ! pour revenir au debut du fichier param_list.dat
72read(num_param,bmelt_ant_reg)
73rewind(num_param)
74read(num_param,bmelt_ocean_temp)
75
76
77!    ecriture dans le fichier parametres
78     write(num_rep_42,*)'fusion basale sous les ice shelves : bmelt-ant-regions_mod'
79     write(num_rep_42,*)'----------------------------------------------------------'
80
81! lecture du fichier contenant les distances au talu continental (m)
82      open(88,file=TRIM(DIRNAMEINP)//'distance_talu-40km.xy')
83
84      do j=1,ny
85         do i=1,nx
86            read(88,'(i3,1x,i3,1x,f10.2)') k,k,dist_talu(i,j)
87         enddo
88      enddo
89      close(88)
90
91! les bords de la grille sont mis a dist_talus=0 de force
92
93
94      dist_talu(1,:)=0.   !a gauche
95      dist_talu(2,:)=0.
96     
97      dist_talu(nx-1,:)=0. !  a droite
98      dist_talu(nx,:)=0. 
99     
100      dist_talu(:,1)=0.    ! en bas   
101      dist_talu(:,2)=0.
102
103      dist_talu(:,ny-1)=0.    ! en haut
104      dist_talu(:,ny)=0.
105
106
107      debug_3D(:,:,32)=dist_talu(:,:)
108
109! lecture du fichier contenant les types de shelves
110!  Ronne-Flichner ->1, Ross -> 2 , Amery -> 3,
111!  PIG-> 4, les petits shelves au dessus de PIG -> 5
112
113      open(88,file=TRIM(DIRNAMEINP)//'numer-ice-shelves-juil07.dat',status='OLD')
114      typeshelf(:,:)=100
115      do k=1,nx*ny
116        read(88,*,end=36) i,j,typeshelf(i,j)
117      end do
11836      close(88)
119
120
121
122  region(:)=0
123  regname(1)='Ronne-Fichner'
124  regname(2)='Ross'
125  regname(3)='Amery'
126  regname(4)='Pig'
127  regname(5)='Petits ice shelves peninsule'
128  regname(6)='autres ice shelves'
129  regname(7)='Au dela du talus continental'
130 
131bms_init:    do j=1,ny
132         do i=1,nx
133
134talus:    if (dist_talu(i,j).GT.5000.) then
135            if (nint(typeshelf(i,j)).eq.1) then ! Ronne-Filchner FRis
136                    bmshelf(i,j)=bmelt_FRis         
137                    bmgrz(i,j)=bmgrz_FRis         
138
139
140               if (region(1).eq.0) then          ! pour l'impression dans 42
141                  region(1)=1
142                  write(num_rep_42,80)regname(1),bmshelf(i,j),bmgrz(i,j)
143               endif
144
14580    format(a32,' bmshelf=',f8.4,'  bmgrz=',f8.4)
146
147            else  if (nint(typeshelf(i,j)).eq.2) then ! Ross
148               bmshelf(i,j)=bmelt_Ross 
149               bmgrz(i,j)=bmgrz_Ross     
150
151               if (region(2).eq.0) then
152                  region(2)=1
153                  write(num_rep_42,80)regname(2),bmshelf(i,j),bmgrz(i,j)
154               endif
155
156
157            else  if (nint(typeshelf(i,j)).eq.3) then ! Amery
158                    bmshelf(i,j)=bmelt_Amery   
159                    bmgrz(i,j)=bmgrz_Amery         
160               if (region(3).eq.0) then
161                  region(3)=1
162                  write(num_rep_42,80)regname(3),bmshelf(i,j),bmgrz(i,j)
163               endif
164
165
166            else  if (nint(typeshelf(i,j)).eq.4) then ! Pig
167                    bmshelf(i,j)=bmelt_PIG 
168                    bmgrz(i,j)=bmgrz_PIG       
169
170               if (region(4).eq.0) then
171                  region(4)=1
172                  write(num_rep_42,80)regname(4),bmshelf(i,j),bmgrz(i,j)
173               endif
174
175            else  if (nint(typeshelf(i,j)).eq.5) then ! petits shelves peninsule
176                    bmshelf(i,j)=bmelt_Pen     
177                    bmgrz(i,j)=bmgrz_Pen       
178
179               if (region(5).eq.0) then
180                  region(5)=1
181                  write(num_rep_42,80)regname(5),bmshelf(i,j),bmgrz(i,j)
182               endif
183
184
185            else !
186                    typeshelf(i,j)=6
187                    bmshelf(i,j)=bmelt_other     
188                    bmgrz(i,j)=bmgrz_other       
189
190               if (region(6).eq.0) then
191                  region(6)=1
192                  write(num_rep_42,80)regname(6),bmshelf(i,j),bmgrz(i,j)
193               endif
194
195
196            endif
197
198
199       else   ! au dela du talus continental
200
201                    bmshelf(i,j)=bmelt_talus 
202                    bmgrz(i,j)=bmgrz_talus   
203 
204               if (region(7).eq.0) then
205                  region(7)=1
206                  write(num_rep_42,80)regname(7),bmshelf(i,j),bmgrz(i,j)
207               endif
208
209
210            endif talus   ! fin du test sur dist_talus
211
212         enddo
213      enddo bms_init
214
215
216    call lect_input(1,'VOTEMP',1,temp_oce1,TRIM(DIRNAMEINP)//trim(file_temp_oce1),file_ncdf)    ! temp ocean ctrl
217    call lect_input(1,'VOTEMP',1,temp_oce2,TRIM(DIRNAMEINP)//trim(file_temp_oce2),file_ncdf)    ! temp ocean
218
219
220
221
222    ! recherche de la valeur du plus proche voisin
223! initialisation
224    mask_oce1(:,:)=.true.
225    where (temp_oce1(:,:).lt.-9999.)
226       mask_oce1(:,:)=.false.
227    endwhere
228
229!    print*, 'temp_oce1(141,141)',temp_oce1(141,141)
230!    print*, 'mask_oce1(141,141)',mask_oce1(141,141)
231
232!    print*, 'temp_oce1(5,5)',temp_oce1(5,5)
233!    print*, 'mask_oce1(5,5)',mask_oce1(5,5)
234
235
236    do j=1,ny
237       do i=1,nx
238          if (.not.mask_oce1(i,j)) then
239!             print*,'boucle mask i j ',i,j,temp_oce1(i,j)
240             test=.true.
241             n=1
242             do while (test)
243                imin=max(1,i-n)
244                imax=min(NX,i+n)
245                jmin=max(1,j-n)
246                jmax=min(NY,j+n)
247                nbr_pts = count(mask_oce1(imin:imax,jmin:jmax))
248!                print*,'nbr_pts', nbr_pts
249                if (nbr_pts.GT.1) then
250                   temp_oce1(i,j)=sum(temp_oce1(imin:imax,jmin:jmax),mask=mask_oce1(imin:imax,jmin:jmax))/nbr_pts ! calcul la moyenne de temp_oce1 sur les pts non masques
251!                   print*,'fin boucle mask i j n ',i,j,n,temp_oce1(i,j)
252                   test=.false.
253                endif
254                n=n+1
255             enddo
256          endif
257       enddo
258    enddo
259
260   
261!    print*, 'temp_oce1 :',temp_oce1(5,5),TRIM(DIRNAMEINP)//trim(file_temp_oce1)
262
263!    stop
264! initialisation
265    mask_oce2(:,:)=.true.
266    where (temp_oce2(:,:).lt.-9999.)
267       mask_oce2(:,:)=.false.
268    endwhere
269
270    do j=1,ny
271       do i=1,nx
272          if (.not.mask_oce2(i,j)) then
273             test=.true.
274             n=1
275             do while (test)
276                imin=max(1,i-n)
277                imax=min(NX,i+n)
278                jmin=max(1,j-n)
279                jmax=min(NY,j+n)
280                nbr_pts = count(mask_oce2(imin:imax,jmin:jmax))
281                if (nbr_pts.GT.1) then
282                   temp_oce2(i,j)=sum(temp_oce2(imin:imax,jmin:jmax),mask=mask_oce2(imin:imax,jmin:jmax))/nbr_pts ! calcul la moyenne de temp_oce1 sur les pts non masques
283                   test=.false.
284                endif
285                n=n+1
286             enddo
287          endif
288       enddo
289    enddo
290
291
292    diff_temp_oce(:,:)=temp_oce2(:,:)-temp_oce1(:,:)
293
294    return
295  end subroutine init_bmelt
296
297
298!________________________________________________________________________________
299subroutine bmeltshelf
300
301
302! cette routine calcule la fusion basale proprement dite
303
304integer :: ngr           ! nombre de voisins flottants
305
306    do j=2,ny-1
307      do i=2,nx-1
308 
309
310         shelf:   if (flot(i,j)) then    ! partie flottante
311           
312            bmelt(i,j)=diff_temp_oce(i,j)*beta_coefbmelt + bmshelf(i,j)
313
314            if (fbm(i,j)) then
315               bmelt(i,j)=diff_temp_oce(i,j)*beta_coefbmelt + bmgrz(i,j)
316            endif
317
318
319! ATTENTION le bloc suivant sert a determiner la fusion basale d'equilibre
320! pour les shelves stationnaires
321
322            if (igrdline.eq.1) then
323               corrbmelt(i,j)=hdot(i,j)+bmelt(i,j)   ! le bmelt d'equilibre
324               debug_3D(i,j,28)=corrbmelt(i,j)
325            endif
326
327
328         else                   ! point posé, on compte le nombre de voisins flottants
329            ngr=0
330            if (flot(i+1,j)) ngr=ngr+1
331            if (flot(i-1,j)) ngr=ngr+1
332            if (flot(i,j+1)) ngr=ngr+1
333            if (flot(i,j-1)) ngr=ngr+1
334
335!   la fusion des points limites est une combinaison entre valeur posée et valeur flottante
336!   en fonction du nombre de points flottants
337
338            bmelt(i,j)= ngr/4.* (bmgrz(i,j)+diff_temp_oce(i,j)*beta_coefbmelt) + (1.-ngr/4.)*bmelt(i,j)
339         endif shelf
340      end do
341   end do
342
343! les bords de la grille sont mis a dist_talus=0 de force
344! attention, pourrait poser des problemes avec AGRIF
345
346   bmelt(1,:)=bmelt_talus   !a gauche
347   bmelt(2,:)=bmelt_talus
348     
349   bmelt(nx-1,:)=bmelt_talus !  a droite
350   bmelt(nx,:)=bmelt_talus 
351     
352   bmelt(:,1)=bmelt_talus    ! en bas   
353   bmelt(:,2)=bmelt_talus
354
355   bmelt(:,ny-1)=bmelt_talus    ! en haut
356   bmelt(:,ny)=bmelt_talus
357
358!cdc test
359!      bmelt(:,:)=diff_temp_oce(:,:)
360!      bmelt(:,:)=bmelt(:,:) + beta_coefbmelt*diff_temp_oce(:,:)
361
362
363if (igrdline.eq.1) then
364    bmelt(:,:)=0.      ! hdot donne alors le -bmelt
365endif
366
367! ecriture des bmelt pour en faire des statistiques
368!--------------------------------------------------
369!open(144,file='bmelt-test')
370!do j=1,ny
371!   do i=1,nx
372!      if (fbm(i,j)) then
373!         grd=1
374!      else
375!         grd=0
376!      endif
377!      if ((flot(i,j)).and.(H(i,j).gt.1.)) then
378!         write(144,999) i,j,grd,nint(typeshelf(i,j)), &
379!              bmelt(i,j),bmgrz(i,j),bmshelf(i,j),hdot(i,j),dist_talu(i,j)
380!      endif
381!   end do
382!end do
383!999 format(4(i0,2x),5(e12.3,2x))
384!close(144)
385
386    return
387    end subroutine  bmeltshelf
388
389END MODULE bmelt_ant_regions_oce
Note: See TracBrowser for help on using the repository browser.