source: branches/GRISLIv3/SOURCES/bmelt-beckmann-gcm_mod.f90 @ 364

Last change on this file since 364 was 346, checked in by dumas, 3 years ago

nzoc, number of vertical levels in T/S files is read in netcdf file
  • Property svn:executable set to *
File size: 16.9 KB
Line 
1!> \file bmelt-ismip6-param.f90
2!! bmelt computed from the non-local formulation suggested in ismip6
3!<
4
5!> \namespace bmelt-ismip6-param
6!! Module for sub-shelf basal melting (grounded or ice shelves)
7!! \author aquiquet
8!! \date April 2019
9!! @note from ismip6 suggested parametrisation
10!! Should be chosen in the module_choix
11!! @note Used modules
12!! @note   - use module3D_phy
13!! @note   - use netcdf
14!! @note   - use io_netcdf_grisli
15!<
16
17!ncwa -a time TO_file_mean.nc TO_file_mean2.nc
18!cdo yearmean TO_file.nc TO_file_mean.nc
19
20module  bmelt_beckmann_gcm_mod
21
22  !$ USE OMP_LIB
23  use module3D_phy,only: nx,ny,dx,dy,ro,rofresh,row,cl,S,H,sealevel_2d,flot,bmelt,dirnameinp,num_param,num_rep_42,time,dt,debug_3d
24  ! note: the geom. (nx,ny,dx,dy) come from module_geoplace
25  ! note: the densities come from param_phy_mod
26  use netcdf
27  use io_netcdf_grisli
28
29  implicit none
30
31  integer, parameter :: nbassins = 10 !< number of sectors in the ocean
32  integer :: nzoc !< number of vertical levels in the ocean (read in netcdf T/S file)
33  real*8 :: coef_OM !combined coefficient by DeConto et Pollard (m/yr/C2)
34  real*8 :: K_t
35  real*8 :: n_tour
36  real*8, dimension (:,:,:), allocatable :: temp_ocean , mask_oce2, temp_ocean_2 !< thermal forcing            , input
37  real*8, dimension (:,:,:), allocatable :: salinity_ocean !< thermal forcing            , input
38 
39  real :: bmelt_empty ! bmelt value for bassins without any GCM data
40     
41!  integer, dimension (nx,ny) :: profondeur
42  integer, dimension(nx,ny) :: bassin
43  real*8, dimension (:), allocatable :: zoc             !< depth of oceanic levels    , input
44  real*8, dimension (nx,ny)      :: mesh_area       !< grid cell area
45  real*8, dimension (nx,ny)      :: ice_draft       !< ice draft (S-H)
46  character(len=200) :: TO_file, SO_file, Bassin_file   ! fichiers de forcage
47  logical, dimension(nbassins) :: mask_bassin_vide  ! mask true if bassin without any GCM data
48 
49contains
50
51 
52  subroutine init_bmelt
53   
54    ! this routine is used to initialise the sub-shelf basal melting rate
55
56
57    real*8, dimension(:),       pointer :: tab1d => null()   !< 2d array real pointer, needed for netcdf readings
58    real*8, dimension(:,:,:),   pointer :: tab3d => null()   !< 3d array real pointer, needed for netcdf readings
59    real*8, dimension(:,:),   pointer :: tab2d => null()   !< 3d array real pointer, needed for netcdf readings
60    character(len=100) :: file_inputs               !< files read
61    character(len=100) :: TO_file                   !< files read
62    character(len=100) :: SO_file                   !< files read
63    character(len=100) :: Bassin_file               !< files read
64    integer            :: i,j,k               
65    integer :: n,imin,imax,jmin,jmax
66    logical :: test ! pour stoper boucle
67    integer :: nbr_pts
68    integer :: depth_max_grid ! max depth level were T S are defined (all grid)
69    integer, dimension(nx,ny) :: depth_max ! max depth were T S are defined for every point
70    integer, dimension(:,:,:), allocatable :: mask_gcm  ! mask were T S are defined
71    integer, dimension(:,:,:), allocatable :: mask_gcm_init ! mask were T S are defined in GCM before interpolation
72    integer, dimension(nbassins) :: max_depth_bassin ! max level were T S are defined for every bassins
73    integer :: nmax ! to avoid inifite loop
74    integer :: err
75    namelist/bmelt_beckmann_gcm_mod/TO_file,SO_file,Bassin_file,bmelt_empty
76   
77    rewind(num_param)        ! loop back at the beginning of the param_list.dat file
78    read(num_param,bmelt_beckmann_gcm_mod)
79    write(num_rep_42,'(A)')'!___________________________________________________________'
80    write(num_rep_42,*)
81    write(num_rep_42,bmelt_beckmann_gcm_mod)
82
83    K_t = 15.77
84
85    file_inputs=TRIM(DIRNAMEINP)//'Snapshots-GCM/'//TRIM(TO_file)
86   
87    ! read the number of levels in T/S files :
88    nzoc = read_nzoc(file_inputs)
89    ! print*, 'nzoc = ',nzoc
90   
91    ! allocate tabs :
92    if (.not.allocated(zoc)) then
93      allocate(zoc(nzoc),stat=err)
94      if (err/=0) then
95          print *,"erreur a l'allocation du tableau zoc",err
96          stop 4
97       end if
98    end if
99    if (.not.allocated(temp_ocean)) then
100      allocate(temp_ocean(nx,ny,nzoc),stat=err)
101      if (err/=0) then
102          print *,"erreur a l'allocation du tableau temp_ocean",err
103          stop 4
104       end if
105    end if
106    if (.not.allocated(mask_oce2)) then
107      allocate(mask_oce2(nx,ny,nzoc),stat=err)
108      if (err/=0) then
109          print *,"erreur a l'allocation du tableau mask_oce2",err
110          stop 4
111       end if
112    end if
113    if (.not.allocated(temp_ocean_2)) then
114      allocate(temp_ocean_2(nx,ny,nzoc),stat=err)
115      if (err/=0) then
116          print *,"erreur a l'allocation du tableau temp_ocean_2",err
117          stop 4
118       end if
119    end if
120    if (.not.allocated(salinity_ocean)) then
121      allocate(salinity_ocean(nx,ny,nzoc),stat=err)
122      if (err/=0) then
123          print *,"erreur a l'allocation du tableau salinity_ocean",err
124          stop 4
125       end if
126    end if
127    if (.not.allocated(mask_gcm)) then
128      allocate(mask_gcm(nx,ny,nzoc),stat=err)
129      if (err/=0) then
130          print *,"erreur a l'allocation du tableau mask_gcm",err
131          stop 4
132       end if
133    end if
134    if (.not.allocated(mask_gcm_init)) then
135      allocate(mask_gcm_init(nx,ny,nzoc),stat=err)
136      if (err/=0) then
137          print *,"erreur a l'allocation du tableau mask_gcm_init",err
138          stop 4
139       end if
140    end if
141     
142    call Read_Ncdf_var('lev',file_inputs,tab1d)
143    zoc(:)  = -tab1d(:)
144    call Read_Ncdf_var('thetao',file_inputs,tab3d)
145    temp_ocean(:,:,:)  = tab3d(:,:,:) 
146       
147    file_inputs=TRIM(DIRNAMEINP)//'Snapshots-GCM/'//TRIM(SO_file)
148    call Read_Ncdf_var('so',file_inputs,tab3d)
149    salinity_ocean(:,:,:)  = tab3d(:,:,:)
150   
151    file_inputs=TRIM(DIRNAMEINP)//'Snapshots-GCM/'//TRIM(Bassin_file)
152    call Read_Ncdf_var('mask',file_inputs,tab2d)
153    bassin(:,:)  = tab2d(:,:)
154       
155   ! calcul de la profondeur max pour chaque point :
156   
157   !$OMP PARALLEL
158   !$OMP WORKSHARE
159   depth_max(:,:) = 0
160   !$OMP END WORKSHARE
161   !$OMP END PARALLEL
162   
163   !$OMP PARALLEL PRIVATE(k)
164   !$OMP DO
165   do j=1,ny
166      do i=1,nx
167         if (temp_ocean(i,j,1).lt.1.e10) then
168            k=2
169            do while (temp_ocean(i,j,k).lt.1.e10 .and. (k.le.nzoc))
170               k=k+1
171            enddo
172            depth_max(i,j) = k-1   ! max level were T and S are defined for this point
173          endif
174      enddo
175   enddo
176   !$OMP END DO
177   !$OMP END PARALLEL
178       
179   depth_max_grid=maxval(depth_max(:,:)) ! max depth were T and S are defined in the grid
180   
181   !$OMP PARALLEL
182   if (depth_max_grid.lt.nzoc) then
183      !$OMP DO
184      do j=1,ny
185         do i=1,nx
186            if (temp_ocean(i,j,1).lt.1.e10 .and. depth_max(i,j).eq.depth_max_grid) then ! T and S are defined up to depth_max_grid
187               ! extension of the values downwards
188               do k=depth_max(i,j)+1,nzoc
189                  temp_ocean(i,j,k)=temp_ocean(i,j,depth_max_grid)
190                  salinity_ocean(i,j,k) = salinity_ocean(i,j,depth_max_grid)
191               enddo
192            endif
193         enddo
194      enddo
195      !$OMP END DO
196   endif
197   
198
199   !$OMP DO
200   do n=1,nbassins
201      max_depth_bassin(n)=maxval(depth_max(:,:),mask=bassin(:,:).eq.n)
202   enddo
203   !$OMP END DO
204   
205!   print*,'max_depth_bassin bassin ', max_depth_bassin
206   
207   
208! mask_gcm_init = 1 si valeur dans gcm, 0 si pas de valeur => tableau avec uniquement les valeurs initiales du GCM
209   !$OMP WORKSHARE 
210   mask_gcm_init(:,:,:) = 1
211   where (temp_ocean(:,:,:).gt.1.e10)
212      mask_gcm_init(:,:,:) = 0
213        endwhere
214   
215! mask_gcm = 1 si valeur dans gcm, 0 si pas de valeur  => tableau update lorsqu'on rempli la grille   
216   mask_gcm(:,:,:) = 1
217   where (temp_ocean(:,:,:).gt.1.e10)
218      mask_gcm(:,:,:) = 0
219        endwhere
220   !$OMP END WORKSHARE 
221   
222   
223! mask_bassin_vide : identifie les bassins oceaniques sans aucune valeur du GCM en surface
224   !$OMP BARRIER
225   !$OMP DO
226   do n=1,nbassins
227      if (sum(mask_gcm(:,:,1),mask=bassin(:,:).eq.n).eq.0) then
228         mask_bassin_vide(n)=.true.
229      else
230         mask_bassin_vide(n)=.false.
231      endif
232   enddo
233   !$OMP END DO
234   !$OMP END PARALLEL
235   nmax = max(nx,ny)   
236
237   
238   !$OMP PARALLEL PRIVATE(n,test)
239   !$OMP DO
240   do k=1,nzoc
241      do j=1,ny
242         do i=1,nx
243            if ((mask_gcm(i,j,k).eq.0) .and. (.not.mask_bassin_vide(bassin(i,j))) .and. (k.le.max_depth_bassin(bassin(i,j)))) then ! point sans valeur dans bassin avec des points GCM
244               test=.true.
245               n=1
246               do while (test.and.(n.lt.nmax))
247                  imin=max(1,i-n)
248                  imax=min(NX,i+n)
249                  jmin=max(1,j-n)
250                  jmax=min(NY,j+n)
251                  nbr_pts = count((bassin(imin:imax,jmin:jmax).eq.bassin(i,j)) .and. (mask_gcm_init(imin:imax,jmin:jmax,k).eq.1))
252                  if (nbr_pts.ge.1) then
253                     temp_ocean(i,j,k)=sum(temp_ocean(imin:imax,jmin:jmax,k),mask = (bassin(imin:imax,jmin:jmax).eq.bassin(i,j)) &
254                           .and. (mask_gcm_init(imin:imax,jmin:jmax,k).eq.1)) / nbr_pts ! calcul la moyenne de temp_ocean sur les pts non masques
255                     salinity_ocean(i,j,k)=sum(salinity_ocean(imin:imax,jmin:jmax,k),mask = (bassin(imin:imax,jmin:jmax).eq.bassin(i,j)) &
256                           .and. (mask_gcm_init(imin:imax,jmin:jmax,k).eq.1)) / nbr_pts ! calcul la moyenne de salinity_ocean sur les pts non masques
257                     mask_gcm(i,j,k)=1
258                     test=.false.
259                  endif
260                  n=n+1
261               enddo
262               if (n.eq.nmax) print*,'bug boucle infinie',n,i,j,k, mask_gcm(i,j,k), mask_bassin_vide(bassin(i,j)),max_depth_bassin(bassin(i,j)), bassin(i,j)
263            endif
264         enddo
265      enddo
266   enddo
267   !$OMP END DO
268 
269! extension des valeurs de temp salinité sur les points des secteurs vides
270   !$OMP DO
271   do k=1,nzoc
272      do j=1,ny
273         do i=1,nx
274            !if ((mask_gcm(i,j,k).eq.0) .and. (mask_bassin_vide(bassin(i,j)))) then ! point sans valeur dans bassin sans aucun points GCM
275            if (mask_gcm(i,j,k).eq.0) then ! point sans valeur dans bassin sans aucun points GCM
276               test=.true.
277               n=1
278               do while (test.and.(n.lt.nmax))
279                  imin=max(1,i-n)
280                  imax=min(NX,i+n)
281                  jmin=max(1,j-n)
282                  jmax=min(NY,j+n)
283                  nbr_pts = count(mask_gcm_init(imin:imax,jmin:jmax,k).eq.1)
284                  if (nbr_pts.ge.1) then
285                     temp_ocean(i,j,k)=sum(temp_ocean(imin:imax,jmin:jmax,k),mask = mask_gcm_init(imin:imax,jmin:jmax,k).eq.1) / nbr_pts ! calcul la moyenne de temp_ocean sur les pts non masques
286                     salinity_ocean(i,j,k)=sum(salinity_ocean(imin:imax,jmin:jmax,k),mask = mask_gcm_init(imin:imax,jmin:jmax,k).eq.1) / nbr_pts! calcul la moyenne de salinity_ocean sur les pts non masques
287                     mask_gcm(i,j,k)=1
288                     test=.false.
289                  endif
290                  n=n+1
291               enddo
292            endif
293         enddo
294      enddo
295   enddo
296   !$OMP END DO
297   
298   !$OMP WORKSHARE
299   debug_3d(:,:,129) =  temp_ocean(:,:,19)
300   debug_3d(:,:,130) =  salinity_ocean(:,:,19)
301   mesh_area(:,:) = dx*dy ! this could be corrected to account for projection deformation
302   !$OMP END WORKSHARE
303   !$OMP END PARALLEL
304   
305!   print*,'zoc profondeur : ',zoc
306   
307!   print*, 'temp_ocean', temp_ocean(1,1,:)
308!   print*, 'salinity_ocean', salinity_ocean(1,1,:)
309
310        !temp_ocean(:,:,:) = temp_ocean(:,:,:) - 273.15
311       
312    coef_OM = K_t * 0.01420418516
313
314    if ( ubound(zoc,dim=1).ne.nzoc) then
315       write (*,*) "bmelt_beckmann_gcm: pb with the number of oceanic layers! abort..."
316       STOP
317    endif
318
319  end subroutine init_bmelt
320
321
322
323  subroutine bmeltshelf
324
325    ! this routine is used to compute the sub-shelf basal melting rate
326
327    real*8, dimension(nx,ny) :: TO_draft
328    real*8, dimension(nx,ny) :: SO_draft
329    real*8, dimension(nx,ny) :: T_freez
330
331    integer :: i,j,k,kinf,ksup,ngr
332    real*8  :: bmloc
333
334
335    !$OMP PARALLEL PRIVATE(ksup,kinf,ngr,bmloc)
336    !$OMP WORKSHARE
337    ice_draft(:,:) = S(:,:)-H(:,:)-sealevel_2d(:,:)
338    !$OMP END WORKSHARE
339
340    !$OMP DO
341    do j=1,ny
342       do i=1,nx
343
344          if (flot(i,j).and.(temp_ocean(i,j,1).lt.1.e10)) then
345
346             if (H(i,j).gt.0.d0) then !limit on the critical thickness of ice to define the ice shelf mask
347                                        ! we should use Hcalv
348
349                ! 1 -  Linear interpolation of the thermal forcing on the ice draft depth :
350                ksup=nzoc
351                do k=nzoc-1,2,-1
352                   if ( zoc(k) .le. ice_draft(i,j) ) ksup = k
353                enddo
354                kinf = ksup - 1
355                if ( ice_draft(i,j) .gt. zoc(1) ) then
356                  TO_draft(i,j) = temp_ocean(i,j,1)
357                  SO_draft(i,j) = salinity_ocean(i,j,1)
358                  T_freez(i,j) = 0.0939 - 0.057 * salinity_ocean(i,j,1) + 0.000764 * zoc(1) 
359                elseif ( ice_draft(i,j) .lt. zoc(nzoc) ) then
360                  TO_draft(i,j) = temp_ocean(i,j,nzoc)
361                  SO_draft(i,j) = salinity_ocean(i,j,nzoc)
362                  T_freez(i,j) = 0.0939 - 0.057 * salinity_ocean(i,j,nzoc) + 0.000764 * zoc(nzoc) 
363                else
364                  !TO_draft(i,j) = (   (zoc(ksup)-ice_draft(i,j)) * temp_ocean(i,j,kinf)   &
365                  !    &                   + (ice_draft(i,j)-zoc(kinf)) * temp_ocean(i,j,ksup) ) / (zoc(ksup)-zoc(kinf))
366                  !SO_draft(i,j) = (   (zoc(ksup)-ice_draft(i,j)) * salinity_ocean(i,j,kinf)   &
367                  !    &                   + (ice_draft(i,j)-zoc(kinf)) * salinity_ocean(i,j,ksup) ) / (zoc(ksup)-zoc(kinf))
368                  TO_draft(i,j) = temp_ocean(i,j,ksup) + ((ice_draft(i,j)-zoc(ksup))*(temp_ocean(i,j,kinf)-temp_ocean(i,j,ksup))/(zoc(kinf)-zoc(ksup))) 
369                  SO_draft(i,j) = salinity_ocean(i,j,ksup) + ((ice_draft(i,j)-zoc(ksup))*(salinity_ocean(i,j,kinf)-salinity_ocean(i,j,ksup))/(zoc(kinf)-zoc(ksup)))
370                  T_freez(i,j) = 0.0939 - 0.057 * SO_draft(i,j) + 0.000764 * ice_draft(i,j)
371                endif
372               
373             else
374                TO_draft(i,j) = temp_ocean(i,j,1)
375                SO_draft(i,j) = salinity_ocean(i,j,1)
376                T_freez(i,j) = 0.0939 - 0.057 * salinity_ocean(i,j,1) + 0.000764 * zoc(1) 
377             endif
378
379          else
380            TO_draft(i,j) = -9999.9d0
381            SO_draft(i,j) = -9999.9d0
382            T_freez(i,j) = -9999.9d0
383          endif
384       enddo
385    enddo
386    !$OMP END DO
387   
388    ! 3  - Calculation of melting rate :
389
390    ! melt rate in m/yr (meters of pure water per year) :
391    ! [ * rhofw_SI / rhoi_SI to get it in meters of ice per year ]
392    !$OMP DO
393    do j=1,ny
394       do i=1,nx
395          if ( TO_draft(i,j) .gt. -9000.d0 ) then ! floating points
396            if (mask_bassin_vide(bassin(i,j)) .or. (bassin(i,j).eq.1)) then ! bassins without any GCM data or bassin 1
397               bmelt(i,j) = bmelt_empty
398            else  ! Bassins with GCM data
399               bmelt(i,j) = coef_OM * ( TO_draft(i,j) - T_freez(i,j) )* abs( TO_draft(i,j) - T_freez(i,j) )
400            endif
401          else
402            bmelt(i,j)=0.
403          endif
404       enddo
405    enddo
406    !$OMP END DO
407   
408    !$OMP DO
409    do j=2,ny-1
410       do i=2,nx-1
411          if ( TO_draft(i,j) .le. -9000.d0 ) then !grounded points
412             ngr=0
413             bmloc=0.d0
414             if (TO_draft(i+1,j) .le. -9000.d0 ) then
415                ngr=ngr+1
416                bmloc=bmloc+bmelt(i+1,j)
417             endif
418             if (TO_draft(i-1,j) .le. -9000.d0 ) then
419                ngr=ngr+1
420                bmloc=bmloc+bmelt(i-1,j)
421             endif
422             if (TO_draft(i,j+1) .le. -9000.d0 ) then
423                ngr=ngr+1
424                bmloc=bmloc+bmelt(i,j+1)
425             endif
426             if (TO_draft(i,j-1) .le. -9000.d0 ) then
427                ngr=ngr+1
428                bmloc=bmloc+bmelt(i,j-1)
429             endif
430             bmelt(i,j)=(ngr/4.)*bmloc+(1-ngr/4.)*bmelt(i,j) 
431          endif
432       enddo
433    enddo
434    !$OMP END DO
435    !$OMP END PARALLEL
436
437  end subroutine bmeltshelf
438
439  function read_nzoc(filename)
440    ! Reads number of vertical levels in ocean forcing NetCDF File.
441    ! Returns the levels number
442
443    ! Reads
444    use netcdf
445
446    implicit none
447    ! input
448    character (len = *), intent(in) :: filename
449    integer :: read_nzoc
450
451    ! This will be the netCDF ID for the file and data variable.
452    integer :: ncid, varid, status
453    integer, dimension(3) :: count
454    integer, dimension(3) :: start
455
456    ! Open the file.
457    status = nf90_open(filename, nf90_nowrite, ncid)
458    if (status/=nf90_noerr) then
459       write(*,*)"unable to open netcdf file : ",filename
460       stop
461    endif
462
463    ! Get the varids of the netCDF variable.
464    status = nf90_inq_dimid(ncid, "lev", varid)
465    status = nf90_inquire_dimension(ncid, varid, len = read_nzoc)
466
467    ! Close the file. This frees up any internal netCDF resources
468    ! associated with the file.
469    status = nf90_close(ncid)
470
471  end function read_nzoc
472
473end module bmelt_beckmann_gcm_mod
Note: See TracBrowser for help on using the repository browser.