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

Last change on this file was 490, checked in by dumas, 4 months ago

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