Changeset 461


Ignore:
Timestamp:
01/12/24 16:26:19 (4 months ago)
Author:
aquiquet
Message:

Cleaning branch: update of bmelt-beckmann to follow the trunk

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/GRISLIv3/SOURCES/bmelt-beckmann-gcm_mod.f90

    r449 r461  
    3535  real*8 :: K_t 
    3636  real*8 :: n_tour 
    37   real*8, dimension (:,:,:), allocatable :: temp_ocean , mask_oce2, temp_ocean_2 !< thermal forcing            , input 
     37  real*8, dimension (:,:,:), allocatable :: temp_ocean     !< thermal forcing            , input 
    3838  real*8, dimension (:,:,:), allocatable :: salinity_ocean !< thermal forcing            , input 
    3939   
     
    4141       
    4242!  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 
    4346  integer, dimension(nx,ny) :: bassin 
    4447  real*8, dimension (:), allocatable :: zoc             !< depth of oceanic levels    , input 
     
    4750  character(len=200) :: TO_file, SO_file, Bassin_file   ! fichiers de forcage 
    4851  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) 
    4954   
    5055contains 
     
    6368    character(len=100) :: SO_file                   !< files read 
    6469    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) 
    6572    integer            :: i,j,k                
    6673    integer :: n,imin,imax,jmin,jmax 
     
    7481    integer :: nmax ! to avoid inifite loop  
    7582    integer :: err 
    76     namelist/bmelt_beckmann_gcm_mod/K_t,TO_file,SO_file,Bassin_file,nbassins,bmelt_empty 
     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 
    7785     
    7886    rewind(num_param)        ! loop back at the beginning of the param_list.dat file 
     
    112120       end if 
    113121    end if 
    114     if (.not.allocated(mask_oce2)) then 
    115       allocate(mask_oce2(nx,ny,nzoc),stat=err) 
    116       if (err/=0) then 
    117           print *,"erreur a l'allocation du tableau mask_oce2",err 
    118           stop 4 
    119        end if 
    120     end if 
    121     if (.not.allocated(temp_ocean_2)) then 
    122       allocate(temp_ocean_2(nx,ny,nzoc),stat=err) 
    123       if (err/=0) then 
    124           print *,"erreur a l'allocation du tableau temp_ocean_2",err 
    125           stop 4 
    126        end if 
    127     end if 
    128122    if (.not.allocated(salinity_ocean)) then 
    129123      allocate(salinity_ocean(nx,ny,nzoc),stat=err) 
     
    151145      if (err/=0) then 
    152146          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 
    153161          stop 4 
    154162       end if 
     
    163171    call Read_Ncdf_var('so',file_inputs,tab3d) 
    164172    salinity_ocean(:,:,:)  = tab3d(:,:,:) 
    165      
     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 
    166182    file_inputs=TRIM(DIRNAMEINP)//'Snapshots-GCM/'//TRIM(Bassin_file) 
    167183    call Read_Ncdf_var('mask',file_inputs,tab2d) 
     
    226242   where (temp_ocean(:,:,:).gt.1.e10) 
    227243      mask_gcm_init(:,:,:) = 0 
    228         endwhere 
     244   endwhere 
    229245    
    230246! mask_gcm = 1 si valeur dans gcm, 0 si pas de valeur  => tableau update lorsqu'on rempli la grille    
     
    232248   where (temp_ocean(:,:,:).gt.1.e10) 
    233249      mask_gcm(:,:,:) = 0 
    234         endwhere 
     250   endwhere 
    235251   !$OMP END WORKSHARE   
    236252    
     
    352368    !$OMP WORKSHARE 
    353369    ice_draft(:,:) = S(:,:)-H(:,:)-sealevel_2d(:,:) 
     370    mean_TF(:) = 0.d0 
     371    IS_area(:) = 0.d0 
    354372    !$OMP END WORKSHARE 
    355373 
     
    370388                kinf = ksup - 1 
    371389                if ( ice_draft(i,j) .gt. zoc(1) ) then  
    372                   TO_draft(i,j) = temp_ocean(i,j,1) 
    373                   SO_draft(i,j) = salinity_ocean(i,j,1) 
    374                   T_freez(i,j) = 0.0939 - 0.057 * salinity_ocean(i,j,1) + 0.000764 * zoc(1)  
     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)  
    375393                elseif ( ice_draft(i,j) .lt. zoc(nzoc) ) then 
    376                   TO_draft(i,j) = temp_ocean(i,j,nzoc) 
    377                   SO_draft(i,j) = salinity_ocean(i,j,nzoc) 
    378                   T_freez(i,j) = 0.0939 - 0.057 * salinity_ocean(i,j,nzoc) + 0.000764 * zoc(nzoc)  
     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)  
    379397                else 
    380398                  !TO_draft(i,j) = (   (zoc(ksup)-ice_draft(i,j)) * temp_ocean(i,j,kinf)   & 
     
    382400                  !SO_draft(i,j) = (   (zoc(ksup)-ice_draft(i,j)) * salinity_ocean(i,j,kinf)   & 
    383401                  !    &                   + (ice_draft(i,j)-zoc(kinf)) * salinity_ocean(i,j,ksup) ) / (zoc(ksup)-zoc(kinf)) 
    384                   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)))  
    385                   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))) 
    386                   T_freez(i,j) = 0.0939 - 0.057 * SO_draft(i,j) + 0.000764 * ice_draft(i,j) 
     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) 
    387410                endif 
    388411                 
     
    401424    enddo 
    402425    !$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 
    403432     
    404433    ! 3  - Calculation of melting rate : 
     
    409438    do j=1,ny 
    410439       do i=1,nx 
    411           if ( TO_draft(i,j) .gt. -9000.d0 ) then ! floating points 
     440          !if ( TO_draft(i,j) .gt. -9000.d0 ) then ! floating points 
     441          if ( flot(i,j) ) then ! floating points 
    412442            !if (mask_bassin_vide(bassin(i,j)) .or. (bassin(i,j).eq.1)) then ! bassins without any GCM data or bassin 1 
    413443            !if (mask_bassin_vide(bassin(i,j)) .or. (bassin(i,j).eq.1)) then !bassins without any GCM data 
     
    415445               bmelt(i,j) = bmelt_empty 
    416446            else  ! Bassins with GCM data 
    417                bmelt(i,j) = coef_OM * ( TO_draft(i,j) - T_freez(i,j) )* abs( TO_draft(i,j) - T_freez(i,j) ) 
     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 
    418452            endif 
    419453          else 
     
    427461    do j=2,ny-1 
    428462       do i=2,nx-1 
    429           if ( TO_draft(i,j) .le. -9000.d0 ) then !grounded points 
    430              ngr=0 
    431              bmloc=0.d0 
    432              !if (TO_draft(i+1,j) .le. -9000.d0 ) then 
    433              if (flot(i+1,j)) then 
    434                 ngr=ngr+1 
    435                 bmloc=bmloc+bmelt(i+1,j) 
     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 
    436490             endif 
    437              !if (TO_draft(i-1,j) .le. -9000.d0 ) then 
    438              if (flot(i-1,j)) then 
    439                 ngr=ngr+1 
    440                 bmloc=bmloc+bmelt(i-1,j) 
    441              endif 
    442              !if (TO_draft(i,j+1) .le. -9000.d0 ) then 
    443              if (flot(i,j+1)) then 
    444                 ngr=ngr+1 
    445                 bmloc=bmloc+bmelt(i,j+1) 
    446              endif 
    447              !if (TO_draft(i,j-1) .le. -9000.d0 ) then 
    448              if (flot(i,j-1)) then 
    449                 ngr=ngr+1 
    450                 bmloc=bmloc+bmelt(i,j-1) 
    451              endif 
    452              ngr=0   ! afq -- I do not allow subgrid ocean melt 
    453491             bmelt(i,j)=(ngr/4.)*bmloc+(1-ngr/4.)*bmelt(i,j)  
    454492          endif 
Note: See TracChangeset for help on using the changeset viewer.