!> \file bmelt-ismip6-param.f90 !! bmelt computed from the non-local formulation suggested in ismip6 !< !> \namespace bmelt-ismip6-param !! Module for sub-shelf basal melting (grounded or ice shelves) !! \author aquiquet !! \date April 2019 !! @note from ismip6 suggested parametrisation !! Should be chosen in the module_choix !! @note Used modules !! @note - use module3D_phy !! @note - use netcdf !! @note - use io_netcdf_grisli !< !ncwa -a time TO_file_mean.nc TO_file_mean2.nc !cdo yearmean TO_file.nc TO_file_mean.nc module bmelt_beckmann_gcm_mod !$ USE OMP_LIB 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 ! note: the geom. (nx,ny,dx,dy) come from module_geoplace ! note: the densities come from param_phy_mod use netcdf use io_netcdf_grisli implicit none integer, parameter :: nbassins = 10 !< number of sectors in the ocean integer :: nzoc !< number of vertical levels in the ocean (read in netcdf T/S file) real*8 :: coef_OM !combined coefficient by DeConto et Pollard (m/yr/C2) real*8 :: K_t real*8 :: n_tour real*8, dimension (:,:,:), allocatable :: temp_ocean , mask_oce2, temp_ocean_2 !< thermal forcing , input real*8, dimension (:,:,:), allocatable :: salinity_ocean !< thermal forcing , input real :: bmelt_empty ! bmelt value for bassins without any GCM data ! integer, dimension (nx,ny) :: profondeur integer, dimension(nx,ny) :: bassin real*8, dimension (:), allocatable :: zoc !< depth of oceanic levels , input real*8, dimension (nx,ny) :: mesh_area !< grid cell area real*8, dimension (nx,ny) :: ice_draft !< ice draft (S-H) character(len=200) :: TO_file, SO_file, Bassin_file ! fichiers de forcage logical, dimension(nbassins) :: mask_bassin_vide ! mask true if bassin without any GCM data contains subroutine init_bmelt ! this routine is used to initialise the sub-shelf basal melting rate real*8, dimension(:), pointer :: tab1d => null() !< 2d array real pointer, needed for netcdf readings real*8, dimension(:,:,:), pointer :: tab3d => null() !< 3d array real pointer, needed for netcdf readings real*8, dimension(:,:), pointer :: tab2d => null() !< 3d array real pointer, needed for netcdf readings character(len=100) :: file_inputs !< files read character(len=100) :: TO_file !< files read character(len=100) :: SO_file !< files read character(len=100) :: Bassin_file !< files read integer :: i,j,k integer :: n,imin,imax,jmin,jmax logical :: test ! pour stoper boucle integer :: nbr_pts integer :: depth_max_grid ! max depth level were T S are defined (all grid) integer, dimension(nx,ny) :: depth_max ! max depth were T S are defined for every point integer, dimension(:,:,:), allocatable :: mask_gcm ! mask were T S are defined integer, dimension(:,:,:), allocatable :: mask_gcm_init ! mask were T S are defined in GCM before interpolation integer, dimension(nbassins) :: max_depth_bassin ! max level were T S are defined for every bassins integer :: nmax ! to avoid inifite loop integer :: err namelist/bmelt_beckmann_gcm_mod/TO_file,SO_file,Bassin_file,bmelt_empty rewind(num_param) ! loop back at the beginning of the param_list.dat file read(num_param,bmelt_beckmann_gcm_mod) write(num_rep_42,'(A)')'!___________________________________________________________' write(num_rep_42,*) write(num_rep_42,bmelt_beckmann_gcm_mod) K_t = 15.77 file_inputs=TRIM(DIRNAMEINP)//'Snapshots-GCM/'//TRIM(TO_file) ! read the number of levels in T/S files : nzoc = read_nzoc(file_inputs) ! print*, 'nzoc = ',nzoc ! allocate tabs : if (.not.allocated(zoc)) then allocate(zoc(nzoc),stat=err) if (err/=0) then print *,"erreur a l'allocation du tableau zoc",err stop 4 end if end if if (.not.allocated(temp_ocean)) then allocate(temp_ocean(nx,ny,nzoc),stat=err) if (err/=0) then print *,"erreur a l'allocation du tableau temp_ocean",err stop 4 end if end if if (.not.allocated(mask_oce2)) then allocate(mask_oce2(nx,ny,nzoc),stat=err) if (err/=0) then print *,"erreur a l'allocation du tableau mask_oce2",err stop 4 end if end if if (.not.allocated(temp_ocean_2)) then allocate(temp_ocean_2(nx,ny,nzoc),stat=err) if (err/=0) then print *,"erreur a l'allocation du tableau temp_ocean_2",err stop 4 end if end if if (.not.allocated(salinity_ocean)) then allocate(salinity_ocean(nx,ny,nzoc),stat=err) if (err/=0) then print *,"erreur a l'allocation du tableau salinity_ocean",err stop 4 end if end if if (.not.allocated(mask_gcm)) then allocate(mask_gcm(nx,ny,nzoc),stat=err) if (err/=0) then print *,"erreur a l'allocation du tableau mask_gcm",err stop 4 end if end if if (.not.allocated(mask_gcm_init)) then allocate(mask_gcm_init(nx,ny,nzoc),stat=err) if (err/=0) then print *,"erreur a l'allocation du tableau mask_gcm_init",err stop 4 end if end if call Read_Ncdf_var('lev',file_inputs,tab1d) zoc(:) = -tab1d(:) call Read_Ncdf_var('thetao',file_inputs,tab3d) temp_ocean(:,:,:) = tab3d(:,:,:) file_inputs=TRIM(DIRNAMEINP)//'Snapshots-GCM/'//TRIM(SO_file) call Read_Ncdf_var('so',file_inputs,tab3d) salinity_ocean(:,:,:) = tab3d(:,:,:) file_inputs=TRIM(DIRNAMEINP)//'Snapshots-GCM/'//TRIM(Bassin_file) call Read_Ncdf_var('mask',file_inputs,tab2d) bassin(:,:) = tab2d(:,:) ! calcul de la profondeur max pour chaque point : !$OMP PARALLEL !$OMP WORKSHARE depth_max(:,:) = 0 !$OMP END WORKSHARE !$OMP END PARALLEL !$OMP PARALLEL PRIVATE(k) !$OMP DO do j=1,ny do i=1,nx if (temp_ocean(i,j,1).lt.1.e10) then k=2 do while (temp_ocean(i,j,k).lt.1.e10 .and. (k.le.nzoc)) k=k+1 enddo depth_max(i,j) = k-1 ! max level were T and S are defined for this point endif enddo enddo !$OMP END DO !$OMP END PARALLEL depth_max_grid=maxval(depth_max(:,:)) ! max depth were T and S are defined in the grid !$OMP PARALLEL if (depth_max_grid.lt.nzoc) then !$OMP DO do j=1,ny do i=1,nx 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 ! extension of the values downwards do k=depth_max(i,j)+1,nzoc temp_ocean(i,j,k)=temp_ocean(i,j,depth_max_grid) salinity_ocean(i,j,k) = salinity_ocean(i,j,depth_max_grid) enddo endif enddo enddo !$OMP END DO endif !$OMP DO do n=1,nbassins max_depth_bassin(n)=maxval(depth_max(:,:),mask=bassin(:,:).eq.n) enddo !$OMP END DO ! print*,'max_depth_bassin bassin ', max_depth_bassin ! mask_gcm_init = 1 si valeur dans gcm, 0 si pas de valeur => tableau avec uniquement les valeurs initiales du GCM !$OMP WORKSHARE mask_gcm_init(:,:,:) = 1 where (temp_ocean(:,:,:).gt.1.e10) mask_gcm_init(:,:,:) = 0 endwhere ! mask_gcm = 1 si valeur dans gcm, 0 si pas de valeur => tableau update lorsqu'on rempli la grille mask_gcm(:,:,:) = 1 where (temp_ocean(:,:,:).gt.1.e10) mask_gcm(:,:,:) = 0 endwhere !$OMP END WORKSHARE ! mask_bassin_vide : identifie les bassins oceaniques sans aucune valeur du GCM en surface !$OMP BARRIER !$OMP DO do n=1,nbassins if (sum(mask_gcm(:,:,1),mask=bassin(:,:).eq.n).eq.0) then mask_bassin_vide(n)=.true. else mask_bassin_vide(n)=.false. endif enddo !$OMP END DO !$OMP END PARALLEL nmax = max(nx,ny) !$OMP PARALLEL PRIVATE(n,test) !$OMP DO do k=1,nzoc do j=1,ny do i=1,nx 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 test=.true. n=1 do while (test.and.(n.lt.nmax)) imin=max(1,i-n) imax=min(NX,i+n) jmin=max(1,j-n) jmax=min(NY,j+n) nbr_pts = count((bassin(imin:imax,jmin:jmax).eq.bassin(i,j)) .and. (mask_gcm_init(imin:imax,jmin:jmax,k).eq.1)) if (nbr_pts.ge.1) then temp_ocean(i,j,k)=sum(temp_ocean(imin:imax,jmin:jmax,k),mask = (bassin(imin:imax,jmin:jmax).eq.bassin(i,j)) & .and. (mask_gcm_init(imin:imax,jmin:jmax,k).eq.1)) / nbr_pts ! calcul la moyenne de temp_ocean sur les pts non masques salinity_ocean(i,j,k)=sum(salinity_ocean(imin:imax,jmin:jmax,k),mask = (bassin(imin:imax,jmin:jmax).eq.bassin(i,j)) & .and. (mask_gcm_init(imin:imax,jmin:jmax,k).eq.1)) / nbr_pts ! calcul la moyenne de salinity_ocean sur les pts non masques mask_gcm(i,j,k)=1 test=.false. endif n=n+1 enddo 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) endif enddo enddo enddo !$OMP END DO ! extension des valeurs de temp salinité sur les points des secteurs vides !$OMP DO do k=1,nzoc do j=1,ny do i=1,nx !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 if (mask_gcm(i,j,k).eq.0) then ! point sans valeur dans bassin sans aucun points GCM test=.true. n=1 do while (test.and.(n.lt.nmax)) imin=max(1,i-n) imax=min(NX,i+n) jmin=max(1,j-n) jmax=min(NY,j+n) nbr_pts = count(mask_gcm_init(imin:imax,jmin:jmax,k).eq.1) if (nbr_pts.ge.1) then 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 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 mask_gcm(i,j,k)=1 test=.false. endif n=n+1 enddo endif enddo enddo enddo !$OMP END DO !$OMP WORKSHARE debug_3d(:,:,129) = temp_ocean(:,:,19) debug_3d(:,:,130) = salinity_ocean(:,:,19) mesh_area(:,:) = dx*dy ! this could be corrected to account for projection deformation !$OMP END WORKSHARE !$OMP END PARALLEL ! print*,'zoc profondeur : ',zoc ! print*, 'temp_ocean', temp_ocean(1,1,:) ! print*, 'salinity_ocean', salinity_ocean(1,1,:) !temp_ocean(:,:,:) = temp_ocean(:,:,:) - 273.15 coef_OM = K_t * 0.01420418516 if ( ubound(zoc,dim=1).ne.nzoc) then write (*,*) "bmelt_beckmann_gcm: pb with the number of oceanic layers! abort..." STOP endif end subroutine init_bmelt subroutine bmeltshelf ! this routine is used to compute the sub-shelf basal melting rate real*8, dimension(nx,ny) :: TO_draft real*8, dimension(nx,ny) :: SO_draft real*8, dimension(nx,ny) :: T_freez integer :: i,j,k,kinf,ksup,ngr real*8 :: bmloc !$OMP PARALLEL PRIVATE(ksup,kinf,ngr,bmloc) !$OMP WORKSHARE ice_draft(:,:) = S(:,:)-H(:,:)-sealevel_2d(:,:) !$OMP END WORKSHARE !$OMP DO do j=1,ny do i=1,nx if (flot(i,j).and.(temp_ocean(i,j,1).lt.1.e10)) then if (H(i,j).gt.0.d0) then !limit on the critical thickness of ice to define the ice shelf mask ! we should use Hcalv ! 1 - Linear interpolation of the thermal forcing on the ice draft depth : ksup=nzoc do k=nzoc-1,2,-1 if ( zoc(k) .le. ice_draft(i,j) ) ksup = k enddo kinf = ksup - 1 if ( ice_draft(i,j) .gt. zoc(1) ) then TO_draft(i,j) = temp_ocean(i,j,1) SO_draft(i,j) = salinity_ocean(i,j,1) T_freez(i,j) = 0.0939 - 0.057 * salinity_ocean(i,j,1) + 0.000764 * zoc(1) elseif ( ice_draft(i,j) .lt. zoc(nzoc) ) then TO_draft(i,j) = temp_ocean(i,j,nzoc) SO_draft(i,j) = salinity_ocean(i,j,nzoc) T_freez(i,j) = 0.0939 - 0.057 * salinity_ocean(i,j,nzoc) + 0.000764 * zoc(nzoc) else !TO_draft(i,j) = ( (zoc(ksup)-ice_draft(i,j)) * temp_ocean(i,j,kinf) & ! & + (ice_draft(i,j)-zoc(kinf)) * temp_ocean(i,j,ksup) ) / (zoc(ksup)-zoc(kinf)) !SO_draft(i,j) = ( (zoc(ksup)-ice_draft(i,j)) * salinity_ocean(i,j,kinf) & ! & + (ice_draft(i,j)-zoc(kinf)) * salinity_ocean(i,j,ksup) ) / (zoc(ksup)-zoc(kinf)) 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))) 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))) T_freez(i,j) = 0.0939 - 0.057 * SO_draft(i,j) + 0.000764 * ice_draft(i,j) endif else TO_draft(i,j) = temp_ocean(i,j,1) SO_draft(i,j) = salinity_ocean(i,j,1) T_freez(i,j) = 0.0939 - 0.057 * salinity_ocean(i,j,1) + 0.000764 * zoc(1) endif else TO_draft(i,j) = -9999.9d0 SO_draft(i,j) = -9999.9d0 T_freez(i,j) = -9999.9d0 endif enddo enddo !$OMP END DO ! 3 - Calculation of melting rate : ! melt rate in m/yr (meters of pure water per year) : ! [ * rhofw_SI / rhoi_SI to get it in meters of ice per year ] !$OMP DO do j=1,ny do i=1,nx if ( TO_draft(i,j) .gt. -9000.d0 ) then ! floating points if (mask_bassin_vide(bassin(i,j)) .or. (bassin(i,j).eq.1)) then ! bassins without any GCM data or bassin 1 bmelt(i,j) = bmelt_empty else ! Bassins with GCM data bmelt(i,j) = coef_OM * ( TO_draft(i,j) - T_freez(i,j) )* abs( TO_draft(i,j) - T_freez(i,j) ) endif else bmelt(i,j)=0. endif enddo enddo !$OMP END DO !$OMP DO do j=2,ny-1 do i=2,nx-1 if ( TO_draft(i,j) .le. -9000.d0 ) then !grounded points ngr=0 bmloc=0.d0 if (TO_draft(i+1,j) .le. -9000.d0 ) then ngr=ngr+1 bmloc=bmloc+bmelt(i+1,j) endif if (TO_draft(i-1,j) .le. -9000.d0 ) then ngr=ngr+1 bmloc=bmloc+bmelt(i-1,j) endif if (TO_draft(i,j+1) .le. -9000.d0 ) then ngr=ngr+1 bmloc=bmloc+bmelt(i,j+1) endif if (TO_draft(i,j-1) .le. -9000.d0 ) then ngr=ngr+1 bmloc=bmloc+bmelt(i,j-1) endif bmelt(i,j)=(ngr/4.)*bmloc+(1-ngr/4.)*bmelt(i,j) endif enddo enddo !$OMP END DO !$OMP END PARALLEL end subroutine bmeltshelf function read_nzoc(filename) ! Reads number of vertical levels in ocean forcing NetCDF File. ! Returns the levels number ! Reads use netcdf implicit none ! input character (len = *), intent(in) :: filename integer :: read_nzoc ! This will be the netCDF ID for the file and data variable. integer :: ncid, varid, status integer, dimension(3) :: count integer, dimension(3) :: start ! Open the file. status = nf90_open(filename, nf90_nowrite, ncid) if (status/=nf90_noerr) then write(*,*)"unable to open netcdf file : ",filename stop endif ! Get the varids of the netCDF variable. status = nf90_inq_dimid(ncid, "lev", varid) status = nf90_inquire_dimension(ncid, varid, len = read_nzoc) ! Close the file. This frees up any internal netCDF resources ! associated with the file. status = nf90_close(ncid) end function read_nzoc end module bmelt_beckmann_gcm_mod