![(please configure the [header_logo] section in trac.ini)](/grisli/chrome/site/Laurentide-16km-3D-Umag-water-bordure-GRISLI-small.png)
Changeset 460
- Timestamp:
- 12/15/23 17:14:02 (6 months ago)
- File:
-
- 1 edited
-
trunk/SOURCES/bmelt-beckmann-gcm_mod.f90 (modified) (16 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/SOURCES/bmelt-beckmann-gcm_mod.f90
r346 r460 21 21 22 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 23 use module3D_phy,only: cl,S,H,sealevel_2d,flot,bmelt,num_param,num_rep_42,time,dt,debug_3d 24 use geography, only: nx,ny,dx,dy,dirnameinp 24 25 ! 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 26 use param_phy_mod,only:ro,rofresh,row 27 use io_netcdf_grisli,only: read_ncdf_var 28 28 29 29 implicit none 30 30 31 integer, parameter :: nbassins = 10 !< number of sectors in the ocean 31 real*8,parameter :: cpw = 4182.d0 32 integer :: nbassins !< number of sectors in the ocean 32 33 integer :: nzoc !< number of vertical levels in the ocean (read in netcdf T/S file) 33 34 real*8 :: coef_OM !combined coefficient by DeConto et Pollard (m/yr/C2) 34 35 real*8 :: K_t 35 36 real*8 :: n_tour 36 real*8, dimension (:,:,:), allocatable :: temp_ocean , mask_oce2, temp_ocean_2!< thermal forcing , input37 real*8, dimension (:,:,:), allocatable :: temp_ocean !< thermal forcing , input 37 38 real*8, dimension (:,:,:), allocatable :: salinity_ocean !< thermal forcing , input 38 39 … … 40 41 41 42 ! integer, dimension (nx,ny) :: profondeur 43 real*8, dimension(nx,ny) :: deltaT_basin 42 44 integer, dimension(nx,ny) :: bassin 43 45 real*8, dimension (:), allocatable :: zoc !< depth of oceanic levels , input … … 45 47 real*8, dimension (nx,ny) :: ice_draft !< ice draft (S-H) 46 48 character(len=200) :: TO_file, SO_file, Bassin_file ! fichiers de forcage 47 logical, dimension(nbassins) :: mask_bassin_vide ! mask true if bassin without any GCM data 49 logical, dimension(:), allocatable :: mask_bassin_vide ! mask true if bassin without any GCM data 50 integer :: grlmelt_bool !< do we apply ocean melt at the grounded grounding line points (1=yes) 48 51 49 52 contains … … 62 65 character(len=100) :: SO_file !< files read 63 66 character(len=100) :: Bassin_file !< files read 67 character(len=100) :: deltaT_file !< files read, spatial correction for Toc 68 integer :: deltaT_bool !< Do we apply a spatial correction for Toc (1=yes) 64 69 integer :: i,j,k 65 70 integer :: n,imin,imax,jmin,jmax … … 70 75 integer, dimension(:,:,:), allocatable :: mask_gcm ! mask were T S are defined 71 76 integer, dimension(:,:,:), allocatable :: mask_gcm_init ! mask were T S are defined in GCM before interpolation 72 integer, dimension( nbassins):: max_depth_bassin ! max level were T S are defined for every bassins77 integer, dimension(:), allocatable :: max_depth_bassin ! max level were T S are defined for every bassins 73 78 integer :: nmax ! to avoid inifite loop 74 79 integer :: err 75 namelist/bmelt_beckmann_gcm_mod/TO_file,SO_file,Bassin_file,bmelt_empty 80 81 namelist/bmelt_beckmann_gcm_mod/K_t,TO_file,SO_file,Bassin_file,nbassins,deltaT_bool,deltaT_file,bmelt_empty,grlmelt_bool 76 82 77 83 rewind(num_param) ! loop back at the beginning of the param_list.dat file … … 81 87 write(num_rep_42,bmelt_beckmann_gcm_mod) 82 88 83 K_t = 15.7789 !K_t = 15.77 84 90 85 91 file_inputs=TRIM(DIRNAMEINP)//'Snapshots-GCM/'//TRIM(TO_file) … … 97 103 end if 98 104 end if 105 if (.not.allocated(mask_bassin_vide)) then 106 allocate(mask_bassin_vide(nbassins),stat=err) 107 if (err/=0) then 108 print *,"erreur a l'allocation du tableau mask_bassin_vide",err 109 stop 4 110 end if 111 end if 99 112 if (.not.allocated(temp_ocean)) then 100 113 allocate(temp_ocean(nx,ny,nzoc),stat=err) … … 104 117 end if 105 118 end if 106 if (.not.allocated(mask_oce2)) then107 allocate(mask_oce2(nx,ny,nzoc),stat=err)108 if (err/=0) then109 print *,"erreur a l'allocation du tableau mask_oce2",err110 stop 4111 end if112 end if113 if (.not.allocated(temp_ocean_2)) then114 allocate(temp_ocean_2(nx,ny,nzoc),stat=err)115 if (err/=0) then116 print *,"erreur a l'allocation du tableau temp_ocean_2",err117 stop 4118 end if119 end if120 119 if (.not.allocated(salinity_ocean)) then 121 120 allocate(salinity_ocean(nx,ny,nzoc),stat=err) … … 136 135 if (err/=0) then 137 136 print *,"erreur a l'allocation du tableau mask_gcm_init",err 137 stop 4 138 end if 139 end if 140 if (.not.allocated(max_depth_bassin)) then 141 allocate(max_depth_bassin(nbassins),stat=err) 142 if (err/=0) then 143 print *,"erreur a l'allocation du tableau max_depth_bassin",err 138 144 stop 4 139 145 end if … … 148 154 call Read_Ncdf_var('so',file_inputs,tab3d) 149 155 salinity_ocean(:,:,:) = tab3d(:,:,:) 150 156 157 if (deltaT_bool.eq.1) then 158 file_inputs=TRIM(DIRNAMEINP)//'Snapshots-GCM/'//TRIM(deltaT_file) 159 call Read_Ncdf_var('deltaT_basin',file_inputs,tab2d) 160 deltaT_basin(:,:) = tab2d(:,:) 161 else 162 deltaT_basin(:,:) = 0. 163 endif 164 151 165 file_inputs=TRIM(DIRNAMEINP)//'Snapshots-GCM/'//TRIM(Bassin_file) 152 166 call Read_Ncdf_var('mask',file_inputs,tab2d) … … 178 192 179 193 depth_max_grid=maxval(depth_max(:,:)) ! max depth were T and S are defined in the grid 180 194 181 195 !$OMP PARALLEL 182 196 if (depth_max_grid.lt.nzoc) then … … 211 225 where (temp_ocean(:,:,:).gt.1.e10) 212 226 mask_gcm_init(:,:,:) = 0 213 227 endwhere 214 228 215 229 ! mask_gcm = 1 si valeur dans gcm, 0 si pas de valeur => tableau update lorsqu'on rempli la grille … … 217 231 where (temp_ocean(:,:,:).gt.1.e10) 218 232 mask_gcm(:,:,:) = 0 219 233 endwhere 220 234 !$OMP END WORKSHARE 221 235 … … 310 324 !temp_ocean(:,:,:) = temp_ocean(:,:,:) - 273.15 311 325 312 coef_OM = K_t * 0.01420418516 326 !coef_OM = K_t * 0.01420418516 327 coef_OM = K_t * row * cpw / (cl * ro) 313 328 314 329 if ( ubound(zoc,dim=1).ne.nzoc) then … … 393 408 do j=1,ny 394 409 do i=1,nx 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 410 !if ( TO_draft(i,j) .gt. -9000.d0 ) then ! floating points 411 if ( flot(i,j) ) then ! floating points 412 !if (mask_bassin_vide(bassin(i,j)) .or. (bassin(i,j).eq.1)) then ! bassins without any GCM data or bassin 1 413 !if (mask_bassin_vide(bassin(i,j)) .or. (bassin(i,j).eq.1)) then !bassins without any GCM data 414 if (mask_bassin_vide(bassin(i,j))) then !bassins without any GCM data 397 415 bmelt(i,j) = bmelt_empty 398 416 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) )417 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) ) 400 418 endif 401 419 else … … 409 427 do j=2,ny-1 410 428 do i=2,nx-1 411 if ( TO_draft(i,j) .le. -9000.d0 ) then !grounded points 412 ngr=0 413 bmloc=0.d0 414 if (TO_draft(i+1,j) .le. -9000.d0 ) then 415 ngr=ngr+1 416 bmloc=bmloc+bmelt(i+1,j) 417 endif 418 if (TO_draft(i-1,j) .le. -9000.d0 ) then 419 ngr=ngr+1 420 bmloc=bmloc+bmelt(i-1,j) 421 endif 422 if (TO_draft(i,j+1) .le. -9000.d0 ) then 423 ngr=ngr+1 424 bmloc=bmloc+bmelt(i,j+1) 425 endif 426 if (TO_draft(i,j-1) .le. -9000.d0 ) then 427 ngr=ngr+1 428 bmloc=bmloc+bmelt(i,j-1) 429 !if ( TO_draft(i,j) .le. -9000.d0 ) then !grounded points 430 if ( .not. flot(i,j) ) then 431 if (grlmelt_bool.eq.1) then 432 ngr=0 433 bmloc=0.d0 434 !if (TO_draft(i+1,j) .le. -9000.d0 ) then 435 if (flot(i+1,j)) then 436 ngr=ngr+1 437 bmloc=bmloc+bmelt(i+1,j) 438 endif 439 !if (TO_draft(i-1,j) .le. -9000.d0 ) then 440 if (flot(i-1,j)) then 441 ngr=ngr+1 442 bmloc=bmloc+bmelt(i-1,j) 443 endif 444 !if (TO_draft(i,j+1) .le. -9000.d0 ) then 445 if (flot(i,j+1)) then 446 ngr=ngr+1 447 bmloc=bmloc+bmelt(i,j+1) 448 endif 449 !if (TO_draft(i,j-1) .le. -9000.d0 ) then 450 if (flot(i,j-1)) then 451 ngr=ngr+1 452 bmloc=bmloc+bmelt(i,j-1) 453 endif 454 else 455 ngr=0 ! afq -- I do not allow subgrid ocean melt 429 456 endif 430 457 bmelt(i,j)=(ngr/4.)*bmloc+(1-ngr/4.)*bmelt(i,j)
Note: See TracChangeset
for help on using the changeset viewer.
