| 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 | module bmelt_ismip6_param_mod |
|---|
| 18 | |
|---|
| 19 | 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 |
|---|
| 20 | ! note: the geom. (nx,ny,dx,dy) come from module_geoplace |
|---|
| 21 | ! note: the densities come from param_phy_mod |
|---|
| 22 | use netcdf |
|---|
| 23 | use io_netcdf_grisli |
|---|
| 24 | |
|---|
| 25 | implicit none |
|---|
| 26 | |
|---|
| 27 | integer, parameter :: nzoc = 30 !< number of vertical levels in the ocean |
|---|
| 28 | real*8,parameter :: cpw = 3974.d0 !Specific heat of sea water (J/kg/K) |
|---|
| 29 | |
|---|
| 30 | real*8, dimension (nx,ny,nzoc) :: thermal_forcing !< thermal forcing , input |
|---|
| 31 | real*8, dimension (nzoc) :: zoc !< depth of oceanic levels , input |
|---|
| 32 | real*8, dimension (nx,ny) :: basinNumber !< IMBIE2 oc. basin identifier, input |
|---|
| 33 | real*8, dimension (nx,ny) :: deltaT_basin !< deltaT , input |
|---|
| 34 | real*8 :: gamma0 !< gamma0 , input |
|---|
| 35 | |
|---|
| 36 | real*8, dimension (nx,ny) :: mesh_area !< grid cell area |
|---|
| 37 | real*8, dimension (nx,ny) :: ice_draft !< ice draft (S-H) |
|---|
| 38 | |
|---|
| 39 | real*8 :: cste !< a pre-computed constant |
|---|
| 40 | integer :: Nbasin !< number of oceanic basin |
|---|
| 41 | |
|---|
| 42 | ! For initMIP, we need to read an anomaly of bmelt...: |
|---|
| 43 | integer :: bmelt_time !< pour selectionner le type de calcul de bmelt |
|---|
| 44 | double precision,dimension(nx,ny) :: bmelt_anom !< anomalie de bmelt pour exp initMIP abmb |
|---|
| 45 | ! double precision,dimension(nx,ny,nzoc) :: TF_0 !> thermal_forcing initial |
|---|
| 46 | double precision,dimension(:,:,:,:),allocatable :: TF_time !< snapshots thermal_forcing pour exp ISMIP |
|---|
| 47 | real,dimension(:),allocatable :: time_snap !> date des snapshots indice : nb de nb_snap |
|---|
| 48 | real :: time_depart_snaps !> temps du debut premier snapshot |
|---|
| 49 | integer :: nb_snap !> nombre de snapshots |
|---|
| 50 | contains |
|---|
| 51 | |
|---|
| 52 | |
|---|
| 53 | |
|---|
| 54 | subroutine init_bmelt |
|---|
| 55 | |
|---|
| 56 | ! this routine is used to initialise the sub-shelf basal melting rate |
|---|
| 57 | |
|---|
| 58 | |
|---|
| 59 | real*8, dimension(:), pointer :: tab1d => null() !< 2d array real pointer, needed for netcdf readings |
|---|
| 60 | real*8, dimension(:,:), pointer :: tab2d => null() !< 2d array real pointer, needed for netcdf readings |
|---|
| 61 | real*8, dimension(:,:,:), pointer :: tab3d => null() !< 3d array real pointer, needed for netcdf readings |
|---|
| 62 | real*8, dimension(:,:,:,:), pointer :: tab4d => null() !< 4d array real pointer, needed for netcdf readings |
|---|
| 63 | |
|---|
| 64 | character(len=100) :: file_inputs !< files read |
|---|
| 65 | character(len=100) :: file_TF !< files read |
|---|
| 66 | character(len=100) :: file_basinNumbers !< files read |
|---|
| 67 | character(len=100) :: file_coef !< files read |
|---|
| 68 | character(len=100) :: file_bmelt_anom !< files read |
|---|
| 69 | |
|---|
| 70 | integer :: i,j !juste pour les verifs |
|---|
| 71 | integer :: err ! recuperation d'erreur |
|---|
| 72 | |
|---|
| 73 | namelist/bmelt_ismip6_param/file_TF,file_basinNumbers,file_coef,gamma0 |
|---|
| 74 | namelist/bmelt_anom_initMIP/file_bmelt_anom,nb_snap,time_depart_snaps,bmelt_time |
|---|
| 75 | |
|---|
| 76 | rewind(num_param) ! loop back at the beginning of the param_list.dat file |
|---|
| 77 | read(num_param,bmelt_ismip6_param) |
|---|
| 78 | write(num_rep_42,'(A)') '! module bmelt_ismip6_param' |
|---|
| 79 | write(num_rep_42,bmelt_ismip6_param) |
|---|
| 80 | write(num_rep_42,'(A)') '! file_TF: 3d thermal forcing' |
|---|
| 81 | write(num_rep_42,'(A)') '! file_basinNumbers: identifier for the Imbie2 basins' |
|---|
| 82 | write(num_rep_42,'(A)') '! file_coef: deltaT' |
|---|
| 83 | write(num_rep_42,'(A)') '! gamma0: value associated with the deltaT file' |
|---|
| 84 | |
|---|
| 85 | |
|---|
| 86 | file_inputs=TRIM(DIRNAMEINP)//trim(file_TF) |
|---|
| 87 | call Read_Ncdf_var('z',file_inputs,tab1d) |
|---|
| 88 | zoc(:) = tab1d(:) |
|---|
| 89 | call Read_Ncdf_var('thermal_forcing',file_inputs,tab3d) |
|---|
| 90 | thermal_forcing(:,:,:) = tab3d(:,:,:) |
|---|
| 91 | |
|---|
| 92 | file_inputs=TRIM(DIRNAMEINP)//trim(file_basinNumbers) |
|---|
| 93 | call Read_Ncdf_var('basinNumber',file_inputs,tab2d) |
|---|
| 94 | basinNumber(:,:) = tab2d(:,:) |
|---|
| 95 | |
|---|
| 96 | file_inputs=TRIM(DIRNAMEINP)//trim(file_coef) |
|---|
| 97 | call Read_Ncdf_var('deltaT_basin',file_inputs,tab2d) |
|---|
| 98 | deltaT_basin(:,:) = tab2d(:,:) |
|---|
| 99 | |
|---|
| 100 | if ( ubound(zoc,dim=1).ne.nzoc) then |
|---|
| 101 | write (*,*) "bmelt-ismip6-param: pb with the number of oceanic layers! abort..." |
|---|
| 102 | STOP |
|---|
| 103 | endif |
|---|
| 104 | |
|---|
| 105 | mesh_area(:,:) = dx*dy ! this could be corrected to account for projection deformation |
|---|
| 106 | |
|---|
| 107 | Nbasin=maxval(basinNumber)+1 ! number of basins |
|---|
| 108 | |
|---|
| 109 | cste = (row*cpw/(ro*cl))**2 ! in K^(-2) |
|---|
| 110 | |
|---|
| 111 | where (thermal_forcing(:,:,:).lt.0.d0) |
|---|
| 112 | thermal_forcing(:,:,:) = 3.5 |
|---|
| 113 | endwhere |
|---|
| 114 | ! TF_0(:,:,:) = thermal_forcing(:,:,:) |
|---|
| 115 | |
|---|
| 116 | ! ______ Anomalies... |
|---|
| 117 | |
|---|
| 118 | rewind(num_param) ! pour revenir au debut du fichier param_list.dat |
|---|
| 119 | read(num_param,bmelt_anom_initMIP) |
|---|
| 120 | |
|---|
| 121 | write(num_rep_42,'(A)')'! module bmelt_ismip6_param ' |
|---|
| 122 | write(num_rep_42,bmelt_anom_initMIP) |
|---|
| 123 | write(num_rep_42,'(A)')'! file_bmelt_anom = fichier anomalie bmelt ' |
|---|
| 124 | write(num_rep_42,'(A)')'! nb_snap = nombre de snapshots ' |
|---|
| 125 | write(num_rep_42,'(A)')'! time_depart_snaps = debut du forçage ' |
|---|
| 126 | write(num_rep_42,'(A)')'! bmelt_time = 0:fixe, 1:anomalies, 2:interp snapshots ' |
|---|
| 127 | write(num_rep_42,'(A)')'!_______________________________________________________________________' |
|---|
| 128 | |
|---|
| 129 | file_bmelt_anom = trim(dirnameinp)//trim(file_bmelt_anom) |
|---|
| 130 | |
|---|
| 131 | if ( bmelt_time == 1 ) then |
|---|
| 132 | call Read_Ncdf_var('abmb',file_bmelt_anom,tab2d) |
|---|
| 133 | bmelt_anom (:,:) = tab2d(:,:) |
|---|
| 134 | elseif (bmelt_time == 2) then ! fichier 3D de TF_time : lecture des snapshots |
|---|
| 135 | ! allocation dynamique de time_snap |
|---|
| 136 | if (allocated(time_snap)) then |
|---|
| 137 | deallocate(time_snap,stat=err) |
|---|
| 138 | if (err/=0) then |
|---|
| 139 | print *,"Erreur à la desallocation de time_snap",err |
|---|
| 140 | stop |
|---|
| 141 | end if |
|---|
| 142 | end if |
|---|
| 143 | |
|---|
| 144 | allocate(time_snap(nb_snap),stat=err) |
|---|
| 145 | if (err/=0) then |
|---|
| 146 | print *,"erreur a l'allocation du tableau time_snap ",err |
|---|
| 147 | print *,"nb_snap = ",nb_snap |
|---|
| 148 | stop |
|---|
| 149 | end if |
|---|
| 150 | ! allocation dynamique de TF_time |
|---|
| 151 | if (allocated(TF_time)) then |
|---|
| 152 | deallocate(TF_time,stat=err) |
|---|
| 153 | if (err/=0) then |
|---|
| 154 | print *,"Erreur à la desallocation de TF_time",err |
|---|
| 155 | stop |
|---|
| 156 | end if |
|---|
| 157 | end if |
|---|
| 158 | |
|---|
| 159 | allocate(TF_time(nx,ny,nzoc,nb_snap),stat=err) |
|---|
| 160 | if (err/=0) then |
|---|
| 161 | print *,"erreur a l'allocation du tableau TF_time ",err |
|---|
| 162 | print *,"nx,ny,nb_snap = ",nx,',',ny,',',nb_snap |
|---|
| 163 | stop |
|---|
| 164 | end if |
|---|
| 165 | ! lecture de tann_snap |
|---|
| 166 | call Read_Ncdf_var('thermal_forcing',file_bmelt_anom,tab4D) |
|---|
| 167 | TF_time(:,:,:,:) = tab4D(:,:,:,:) |
|---|
| 168 | |
|---|
| 169 | ! lecture de time_snap |
|---|
| 170 | call Read_Ncdf_var('time',file_bmelt_anom,tab1D) |
|---|
| 171 | time_snap(:) = tab1D(:) |
|---|
| 172 | time_snap(:) = time_snap(:) + time_depart_snaps |
|---|
| 173 | endif |
|---|
| 174 | end subroutine init_bmelt |
|---|
| 175 | |
|---|
| 176 | |
|---|
| 177 | |
|---|
| 178 | subroutine bmeltshelf |
|---|
| 179 | |
|---|
| 180 | ! this routine is used to compute the sub-shelf basal melting rate |
|---|
| 181 | |
|---|
| 182 | real*8, dimension(nx,ny) :: TF_draft |
|---|
| 183 | real*8,allocatable,dimension(:) :: mean_TF, IS_area |
|---|
| 184 | |
|---|
| 185 | integer :: i,j,k,kinf,ksup,ngr |
|---|
| 186 | real*8 :: bmloc |
|---|
| 187 | |
|---|
| 188 | real :: coefanomtime ! pour initMIP abmb |
|---|
| 189 | |
|---|
| 190 | allocate( mean_TF(Nbasin), IS_area(Nbasin) ) |
|---|
| 191 | |
|---|
| 192 | if (bmelt_time == 0) then |
|---|
| 193 | coefanomtime = 0. |
|---|
| 194 | else if (bmelt_time == 1) then |
|---|
| 195 | coefanomtime = min ( real(time/40.) , 1. ) |
|---|
| 196 | else if (bmelt_time == 2) then |
|---|
| 197 | call TF_ISMIP_RCM |
|---|
| 198 | end if |
|---|
| 199 | |
|---|
| 200 | mean_TF(:) = 0.d0 |
|---|
| 201 | IS_area(:) = 0.d0 |
|---|
| 202 | |
|---|
| 203 | ice_draft(:,:) = S(:,:)-H(:,:)-sealevel_2d(:,:) |
|---|
| 204 | |
|---|
| 205 | do j=1,ny |
|---|
| 206 | do i=1,nx |
|---|
| 207 | |
|---|
| 208 | if (flot(i,j)) then |
|---|
| 209 | |
|---|
| 210 | if (H(i,j).gt.0.d0) then !limit on the critical thickness of ice to define the ice shelf mask |
|---|
| 211 | ! we should use Hcalv |
|---|
| 212 | |
|---|
| 213 | ! 1 - Linear interpolation of the thermal forcing on the ice draft depth : |
|---|
| 214 | ksup=nzoc |
|---|
| 215 | do k=nzoc-1,2,-1 |
|---|
| 216 | if ( zoc(k) .le. ice_draft(i,j) ) ksup = k |
|---|
| 217 | enddo |
|---|
| 218 | kinf = ksup - 1 |
|---|
| 219 | if ( ice_draft(i,j) .gt. zoc(1) ) then |
|---|
| 220 | TF_draft(i,j) = thermal_forcing(i,j,1) |
|---|
| 221 | elseif ( ice_draft(i,j) .lt. zoc(nzoc) ) then |
|---|
| 222 | TF_draft(i,j) = thermal_forcing(i,j,nzoc) |
|---|
| 223 | else |
|---|
| 224 | TF_draft(i,j) = ( (zoc(ksup)-ice_draft(i,j)) * thermal_forcing(i,j,kinf) & |
|---|
| 225 | & + (ice_draft(i,j)-zoc(kinf)) * thermal_forcing(i,j,ksup) ) / (zoc(ksup)-zoc(kinf)) |
|---|
| 226 | endif |
|---|
| 227 | |
|---|
| 228 | ! 2 - Mean Thermal forcing in individual basins (NB: fortran norm while basins start at zero) : |
|---|
| 229 | mean_TF( basinNumber(i,j)+1 ) = mean_TF( basinNumber(i,j)+1 ) + mesh_area(i,j) * TF_draft(i,j) |
|---|
| 230 | IS_area( basinNumber(i,j)+1 ) = IS_area( basinNumber(i,j)+1 ) + mesh_area(i,j) |
|---|
| 231 | |
|---|
| 232 | else |
|---|
| 233 | TF_draft(i,j) = thermal_forcing(i,j,1) |
|---|
| 234 | endif |
|---|
| 235 | |
|---|
| 236 | else |
|---|
| 237 | |
|---|
| 238 | TF_draft(i,j) = -9999.9d0 |
|---|
| 239 | |
|---|
| 240 | endif |
|---|
| 241 | |
|---|
| 242 | enddo |
|---|
| 243 | enddo |
|---|
| 244 | |
|---|
| 245 | where ( IS_area(:).gt.0.d0) ! for all basins that have ice shelves |
|---|
| 246 | mean_TF(:) = mean_TF(:) / IS_area(:) |
|---|
| 247 | elsewhere ! we have no floating points over the whole basin |
|---|
| 248 | mean_TF(:) = -9999.d0 |
|---|
| 249 | endwhere |
|---|
| 250 | |
|---|
| 251 | ! 3 - Calculation of melting rate : |
|---|
| 252 | |
|---|
| 253 | ! melt rate in m/yr (meters of pure water per year) : |
|---|
| 254 | ! [ * rhofw_SI / rhoi_SI to get it in meters of ice per year ] |
|---|
| 255 | do j=1,ny |
|---|
| 256 | do i=1,nx |
|---|
| 257 | if ( TF_draft(i,j) .gt. -5.d0 ) then !floating points |
|---|
| 258 | if ((mean_TF(basinNumber(i,j)+1).gt.-9999.d0).and.(H(i,j).gt.0.)) then !should use Hcoup |
|---|
| 259 | bmelt(i,j) = gamma0 * cste * ( TF_draft(i,j) + deltaT_basin(i,j) )* abs( mean_TF(basinNumber(i,j)+1) + deltaT_basin(i,j) ) |
|---|
| 260 | !if (bmelt(i,j).lt.-0.2) then |
|---|
| 261 | ! write(*,*) "Strong sub-shelf refreezing: ", bmelt(i,j), i,j,flot(i,j),ice_draft(i,j),TF_draft(i,j), mean_TF(basinNumber(i,j)+1), deltaT_basin(i,j) |
|---|
| 262 | !endif |
|---|
| 263 | else ! in case we have no floating points over a whole basin, we take the local expression |
|---|
| 264 | bmelt(i,j) = gamma0 * cste * max( TF_draft(i,j) + deltaT_basin(i,j) , 0.d0 ) * max( TF_draft(i,j) + deltaT_basin(i,j) , 0.d0 ) |
|---|
| 265 | endif |
|---|
| 266 | if (bmelt_time.eq.1) bmelt(i,j) = bmelt(i,j) + (coefanomtime * bmelt_anom(i,j)) |
|---|
| 267 | endif |
|---|
| 268 | enddo |
|---|
| 269 | enddo |
|---|
| 270 | |
|---|
| 271 | do j=1,ny |
|---|
| 272 | do i=1,nx |
|---|
| 273 | if ( TF_draft(i,j) .le. -5.d0 ) then !grounded points |
|---|
| 274 | ngr=0 |
|---|
| 275 | bmloc=0.d0 |
|---|
| 276 | if (flot(i+1,j)) then |
|---|
| 277 | ngr=ngr+1 |
|---|
| 278 | bmloc=bmloc+bmelt(i+1,j) |
|---|
| 279 | endif |
|---|
| 280 | if (flot(i-1,j)) then |
|---|
| 281 | ngr=ngr+1 |
|---|
| 282 | bmloc=bmloc+bmelt(i-1,j) |
|---|
| 283 | endif |
|---|
| 284 | if (flot(i,j+1)) then |
|---|
| 285 | ngr=ngr+1 |
|---|
| 286 | bmloc=bmloc+bmelt(i,j+1) |
|---|
| 287 | endif |
|---|
| 288 | if (flot(i,j-1)) then |
|---|
| 289 | ngr=ngr+1 |
|---|
| 290 | bmloc=bmloc+bmelt(i,j-1) |
|---|
| 291 | endif |
|---|
| 292 | bmelt(i,j)=(ngr/4.)*bmloc+(1-ngr/4.)*bmelt(i,j) |
|---|
| 293 | if (bmelt_time.eq.1) bmelt(i,j) = bmelt(i,j) + ((ngr/4.) * coefanomtime * bmelt_anom(i,j)) |
|---|
| 294 | endif |
|---|
| 295 | enddo |
|---|
| 296 | enddo |
|---|
| 297 | |
|---|
| 298 | |
|---|
| 299 | end subroutine bmeltshelf |
|---|
| 300 | |
|---|
| 301 | subroutine TF_ISMIP_RCM ! calcule le thermal forcing a partir des snapshots |
|---|
| 302 | |
|---|
| 303 | implicit none |
|---|
| 304 | integer :: k ! pour calculer les indices de temps |
|---|
| 305 | real,dimension(nx,ny,nzoc) :: thermal_forcing_time ! pour calcul Thermal forcing |
|---|
| 306 | |
|---|
| 307 | ! calcul TF a partir fichier snapshots TF_ISMIP_RCM |
|---|
| 308 | ! Calcule le mass balance d'apres un fichier de snapshots |
|---|
| 309 | |
|---|
| 310 | ! calcule thermal_forcing_time par interpolation entre deux snapshots |
|---|
| 311 | ! avant prend la valeur de reference |
|---|
| 312 | ! apres prend la derniere valeur |
|---|
| 313 | ! en general les snapshots vont de 1995 a 2100 (time_snap de 0 Ã 105) |
|---|
| 314 | if(time.lt.time_snap(1)) then ! time avant le forcage |
|---|
| 315 | thermal_forcing_time(:,:,:) = TF_time(:,:,:,1) |
|---|
| 316 | else if (time.ge.time_snap(nb_snap)) then ! time apres le forcage |
|---|
| 317 | thermal_forcing_time(:,:,:) = TF_time(:,:,:,nb_snap) |
|---|
| 318 | else ! cas general |
|---|
| 319 | do k = 1 , nb_snap-1 |
|---|
| 320 | if((time.ge.time_snap(k)).and.(time.lt.time_snap(k+1))) then ! entre k et k+1 |
|---|
| 321 | !cdc version avec interpolation lineaire entre 2 annees |
|---|
| 322 | !~ thermal_forcing_time(:,:,:) = TF_time(:,:,:,k) + (TF_time(:,:,:,k+1)-TF_time(:,:,:,k)) * & |
|---|
| 323 | !~ (time-time_snap(k))/(time_snap(k+1)-time_snap(k)) |
|---|
| 324 | !cdc version IMSIP avec annee n => TF(n) pendant un an |
|---|
| 325 | thermal_forcing_time(:,:,:) = TF_time(:,:,:,k) |
|---|
| 326 | endif |
|---|
| 327 | end do |
|---|
| 328 | endif |
|---|
| 329 | |
|---|
| 330 | thermal_forcing(:,:,:) = thermal_forcing_time(:,:,:) |
|---|
| 331 | !~ print*,'TF_ISMIP_RCM time_snap TF', k, time_snap(k), time_snap(k+1), thermal_forcing_time(190,220,1), thermal_forcing_time(190,220,30),TF_time(190,220,1,1),TF_time(190,220,30,1) |
|---|
| 332 | |
|---|
| 333 | end subroutine TF_ISMIP_RCM |
|---|
| 334 | |
|---|
| 335 | |
|---|
| 336 | |
|---|
| 337 | end module bmelt_ismip6_param_mod |
|---|
![(please configure the [header_logo] section in trac.ini)](/grisli/chrome/site/Laurentide-16km-3D-Umag-water-bordure-GRISLI-small.png)
