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

Last change on this file since 465 was 465, checked in by aquiquet, 4 months ago

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