[335] | 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 | |
---|
| 20 | module bmelt_beckmann_gcm_mod |
---|
| 21 | |
---|
[339] | 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 |
---|
[335] | 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 | |
---|
[346] | 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) |
---|
[335] | 33 | real*8 :: coef_OM !combined coefficient by DeConto et Pollard (m/yr/C2) |
---|
| 34 | real*8 :: K_t |
---|
[339] | 35 | real*8 :: n_tour |
---|
[346] | 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 | |
---|
[339] | 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 |
---|
[346] | 43 | real*8, dimension (:), allocatable :: zoc !< depth of oceanic levels , input |
---|
[335] | 44 | real*8, dimension (nx,ny) :: mesh_area !< grid cell area |
---|
| 45 | real*8, dimension (nx,ny) :: ice_draft !< ice draft (S-H) |
---|
[339] | 46 | character(len=200) :: TO_file, SO_file, Bassin_file ! fichiers de forcage |
---|
[346] | 47 | logical, dimension(nbassins) :: mask_bassin_vide ! mask true if bassin without any GCM data |
---|
[339] | 48 | |
---|
[335] | 49 | contains |
---|
| 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 |
---|
[339] | 59 | real*8, dimension(:,:), pointer :: tab2d => null() !< 3d array real pointer, needed for netcdf readings |
---|
[335] | 60 | character(len=100) :: file_inputs !< files read |
---|
| 61 | character(len=100) :: TO_file !< files read |
---|
| 62 | character(len=100) :: SO_file !< files read |
---|
[339] | 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 |
---|
[346] | 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 |
---|
[339] | 72 | integer, dimension(nbassins) :: max_depth_bassin ! max level were T S are defined for every bassins |
---|
| 73 | integer :: nmax ! to avoid inifite loop |
---|
[346] | 74 | integer :: err |
---|
[339] | 75 | namelist/bmelt_beckmann_gcm_mod/TO_file,SO_file,Bassin_file,bmelt_empty |
---|
[335] | 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) |
---|
[339] | 82 | |
---|
[335] | 83 | K_t = 15.77 |
---|
| 84 | |
---|
[339] | 85 | file_inputs=TRIM(DIRNAMEINP)//'Snapshots-GCM/'//TRIM(TO_file) |
---|
[346] | 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 | |
---|
[335] | 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 | |
---|
[339] | 147 | file_inputs=TRIM(DIRNAMEINP)//'Snapshots-GCM/'//TRIM(SO_file) |
---|
[335] | 148 | call Read_Ncdf_var('so',file_inputs,tab3d) |
---|
| 149 | salinity_ocean(:,:,:) = tab3d(:,:,:) |
---|
[339] | 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 | |
---|
[335] | 198 | |
---|
[339] | 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 |
---|
[335] | 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 | |
---|
[339] | 334 | |
---|
| 335 | !$OMP PARALLEL PRIVATE(ksup,kinf,ngr,bmloc) |
---|
| 336 | !$OMP WORKSHARE |
---|
[335] | 337 | ice_draft(:,:) = S(:,:)-H(:,:)-sealevel_2d(:,:) |
---|
[339] | 338 | !$OMP END WORKSHARE |
---|
[335] | 339 | |
---|
[339] | 340 | !$OMP DO |
---|
[335] | 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 |
---|
[339] | 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) |
---|
[335] | 359 | elseif ( ice_draft(i,j) .lt. zoc(nzoc) ) then |
---|
[339] | 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) |
---|
[335] | 363 | else |
---|
[339] | 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) |
---|
[335] | 371 | endif |
---|
[339] | 372 | |
---|
[335] | 373 | else |
---|
| 374 | TO_draft(i,j) = temp_ocean(i,j,1) |
---|
[339] | 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) |
---|
[335] | 377 | endif |
---|
| 378 | |
---|
| 379 | else |
---|
[339] | 380 | TO_draft(i,j) = -9999.9d0 |
---|
| 381 | SO_draft(i,j) = -9999.9d0 |
---|
| 382 | T_freez(i,j) = -9999.9d0 |
---|
[335] | 383 | endif |
---|
| 384 | enddo |
---|
| 385 | enddo |
---|
[339] | 386 | !$OMP END DO |
---|
| 387 | |
---|
[335] | 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 ] |
---|
[339] | 392 | !$OMP DO |
---|
[335] | 393 | do j=1,ny |
---|
| 394 | do i=1,nx |
---|
[339] | 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. |
---|
[335] | 403 | endif |
---|
| 404 | enddo |
---|
| 405 | enddo |
---|
[339] | 406 | !$OMP END DO |
---|
[335] | 407 | |
---|
[339] | 408 | !$OMP DO |
---|
[335] | 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 |
---|
[339] | 414 | if (TO_draft(i+1,j) .le. -9000.d0 ) then |
---|
[335] | 415 | ngr=ngr+1 |
---|
| 416 | bmloc=bmloc+bmelt(i+1,j) |
---|
| 417 | endif |
---|
[339] | 418 | if (TO_draft(i-1,j) .le. -9000.d0 ) then |
---|
[335] | 419 | ngr=ngr+1 |
---|
| 420 | bmloc=bmloc+bmelt(i-1,j) |
---|
| 421 | endif |
---|
[339] | 422 | if (TO_draft(i,j+1) .le. -9000.d0 ) then |
---|
[335] | 423 | ngr=ngr+1 |
---|
| 424 | bmloc=bmloc+bmelt(i,j+1) |
---|
| 425 | endif |
---|
[339] | 426 | if (TO_draft(i,j-1) .le. -9000.d0 ) then |
---|
[335] | 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 |
---|
[339] | 434 | !$OMP END DO |
---|
| 435 | !$OMP END PARALLEL |
---|
[335] | 436 | |
---|
| 437 | end subroutine bmeltshelf |
---|
[346] | 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 | |
---|
[335] | 473 | end module bmelt_beckmann_gcm_mod |
---|