Changeset 264 for trunk

Timestamp:

06/18/19 11:50:42 (5 years ago)

Author:

dumas

Message:

Option to use ocean thermal forcing snapshots (bmelt_time=2) for ISMIP simulations

File:

• trunk/SOURCES/bmelt-ismip6-param_mod.f90 (modified) (10 diffs)

trunk/SOURCES/bmelt-ismip6-param_mod.f90

@@ -17 +17 @@
 module  bmelt_ismip6_param_mod
 
-  use module3D_phy,only: nx,ny,dx,dy,ro,rofresh,row,cl,S,H,sealevel_2d,flot,bmelt,dirnameinp,num_param,num_rep_42
+  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
   ! note: the geom. (nx,ny,dx,dy) come from module_geoplace
   ! note: the densities come from param_phy_mod
@@ -40 +40 @@
   integer :: Nbasin                                 !< number of oceanic basin
   
+  ! For initMIP, we need to read an anomaly of bmelt...:
+  integer                            :: bmelt_time        !< pour selectionner le type de calcul de bmelt  
+  double precision,dimension(nx,ny)  :: bmelt_anom        !< anomalie de bmelt pour exp initMIP abmb
+  double precision,dimension(nx,ny,nzoc)  :: TF_0              !> thermal_forcing initial
+  double precision,dimension(:,:,:,:),allocatable  :: TF_anom  !< anomalie de thermal_forcing pour exp ISMIP
+  real,dimension(:),allocatable      :: time_snap         !> date des snapshots  indice : nb de nb_snap
+  real                               :: time_depart_snaps !> temps du debut premier snapshot
+  integer                            :: nb_snap           !> nombre de snapshots
 contains
 
@@ -49 +57 @@
 
 
-    real*8, dimension(:),     pointer :: tab1d => null()   !< 2d array real pointer, needed for netcdf readings
-    real*8, dimension(:,:),   pointer :: tab2d => 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 :: tab1d => null()   !< 2d array real pointer, needed for netcdf readings
+    real*8, dimension(:,:),     pointer :: tab2d => 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 :: tab4d => null()   !< 4d array real pointer, needed for netcdf readings
     
     character(len=100) :: file_inputs               !< files read
@@ -57 +66 @@
     character(len=100) :: file_basinNumbers         !< files read
     character(len=100) :: file_coef                 !< files read
-
-    integer i,j !juste pour les verifs
+    character(len=100) :: file_bmelt_anom           !< files read
+
+    integer            :: i,j                       !juste pour les verifs
+    integer            :: err                       ! recuperation d'erreur
 
     namelist/bmelt_ismip6_param/file_TF,file_basinNumbers,file_coef,gamma0
+    namelist/bmelt_anom_initMIP/file_bmelt_anom,nb_snap,time_depart_snaps,bmelt_time
     
     rewind(num_param)        ! loop back at the beginning of the param_list.dat file
@@ -100 +112 @@
        thermal_forcing(:,:,:) = 3.5
     endwhere
-    
+    TF_0(:,:,:) = thermal_forcing(:,:,:)
+    
+    ! ______ Anomalies...
+  
+    rewind(num_param)        ! pour revenir au debut du fichier param_list.dat
+    read(num_param,bmelt_anom_initMIP)
+
+    write(num_rep_42,'(A)')'!  module bmelt_ismip6_param                                            '
+    write(num_rep_42,bmelt_anom_initMIP)
+    write(num_rep_42,'(A)')'! file_bmelt_anom   = fichier anomalie bmelt                            '
+    write(num_rep_42,'(A)')'! nb_snap           = nombre de snapshots                               '
+    write(num_rep_42,'(A)')'! time_depart_snaps = debut du forçage                                  '
+    write(num_rep_42,'(A)')'! bmelt_time        = 0:fixe, 1:anomalies, 2:interp snapshots           '
+    write(num_rep_42,'(A)')'!_______________________________________________________________________' 
+
+    file_bmelt_anom  = trim(dirnameinp)//trim(file_bmelt_anom)
+
+    if ( bmelt_time == 1 ) then
+      call Read_Ncdf_var('abmb',file_bmelt_anom,tab2d)
+      bmelt_anom (:,:) = tab2d(:,:)
+    elseif (bmelt_time == 2) then ! fichier 3D de TF_anom : lecture des snapshots
+! allocation dynamique de time_snap
+      if (allocated(time_snap)) then 
+        deallocate(time_snap,stat=err)
+        if (err/=0) then
+          print *,"Erreur à la desallocation de time_snap",err
+          stop 
+        end if
+      end if
+
+      allocate(time_snap(nb_snap),stat=err)
+      if (err/=0) then
+        print *,"erreur a l'allocation du tableau time_snap ",err
+        print *,"nb_snap = ",nb_snap
+        stop 
+      end if
+! allocation dynamique de TF_anom
+      if (allocated(TF_anom)) then 
+        deallocate(TF_anom,stat=err)
+        if (err/=0) then
+          print *,"Erreur à la desallocation de TF_anom",err
+          stop 
+        end if
+      end if
+
+      allocate(TF_anom(nx,ny,nzoc,nb_snap),stat=err)
+      if (err/=0) then
+        print *,"erreur a l'allocation du tableau TF_anom ",err
+        print *,"nx,ny,nb_snap = ",nx,',',ny,',',nb_snap
+        stop 
+      end if
+    ! lecture de tann_snap
+      call Read_Ncdf_var('thermal_forcing',file_bmelt_anom,tab4D)
+      TF_anom(:,:,:,:) = tab4D(:,:,:,:)
+
+! lecture de time_snap
+      call Read_Ncdf_var('time',file_bmelt_anom,tab1D)
+      time_snap(:) = tab1D(:)
+      time_snap(:) = time_snap(:) + time_depart_snaps
+    endif
   end subroutine init_bmelt
 
@@ -115 +186 @@
     real*8  :: bmloc
     
+    real    :: coefanomtime  ! pour initMIP abmb
+    
     allocate( mean_TF(Nbasin), IS_area(Nbasin)  )
+    
+    if (bmelt_time == 0) then
+      coefanomtime = 0.
+    else if (bmelt_time == 1) then
+      coefanomtime = min ( real(time/40.) , 1. )
+    else if (bmelt_time == 2) then
+      call TF_ISMIP_RCM
+    end if
 
     mean_TF(:) = 0.d0
@@ -176 +257 @@
           if ( TF_draft(i,j) .gt. -5.d0 ) then !floating points
              if ((mean_TF(basinNumber(i,j)+1).gt.-9999.d0).and.(H(i,j).gt.200.)) then !should use Hcoup
-                bmelt(i,j) = gamma0 * cste * ( TF_draft(i,j) + deltaT_basin(i,j) )* abs( mean_TF(basinNumber(i,j)+1) + deltaT_basin(i,j) )  
+                bmelt(i,j) = gamma0 * cste * ( TF_draft(i,j) + deltaT_basin(i,j) )* abs( mean_TF(basinNumber(i,j)+1) + deltaT_basin(i,j) ) 
                 !if (bmelt(i,j).lt.-0.2) then
                 !   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)
@@ -183 +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 )
              endif
+             if (bmelt_time.eq.1) bmelt(i,j) = bmelt(i,j) + (coefanomtime * bmelt_anom(i,j))
           endif
        enddo
@@ -209 +291 @@
              endif
              bmelt(i,j)=(ngr/4.)*bmloc+(1-ngr/4.)*bmelt(i,j)
+             if (bmelt_time.eq.1) bmelt(i,j) = bmelt(i,j) + ((ngr/4.) * coefanomtime * bmelt_anom(i,j))
           endif
        enddo
@@ -216 +299 @@
   end subroutine bmeltshelf
   
+subroutine TF_ISMIP_RCM                ! calcule le thermal forcing a partir des snapshots d'anomalies
+
+  implicit none
+  integer             :: k, k_snap               ! pour calculer les indices de temps
+  real :: time_prec
+  real,dimension(nx,ny,nzoc) :: thermal_forcing_time  ! pour calcul Thermal forcing
+
+! calcul TF a partir fichier snapshots TF_ISMIP_RCM
+! Calcule le mass balance d'apres un fichier de snapshots
+
+! calcule thermal_forcing_time par interpolation entre deux snapshots
+! avant prend la valeur de reference
+! apres prend la derniere valeur
+! en general les snapshots vont de 1995 a 2100 (time_snap de 0 Ã  105)
+  if(time.lt.time_snap(1)) then              ! time avant le forcage
+     thermal_forcing_time(:,:,:) = TF_0(:,:,:)
+     k_snap       = 1
+!     S_ref(:,:)   = S(:,:)                   ! du coup sera la surface de reference avant le forcage
+     time_prec    = time
+  else if (time.ge.time_snap(nb_snap)) then  ! time apres le forcage
+     thermal_forcing_time(:,:,:) =  TF_0(:,:,:) + TF_anom(:,:,:,nb_snap)
+     if (abs(time-time_prec-1.).lt.dt) then   ! 
+        time_prec = time_prec + 1
+     endif
+  else                                       ! cas general
+     do k = 1 , nb_snap-1
+        if((time.ge.time_snap(k)).and.(time.lt.time_snap(k+1))) then ! entre k et k+1 
+           thermal_forcing_time(:,:,:) = TF_0(:,:,:) + TF_anom(:,:,:,k) + (TF_anom(:,:,:,k+1)-TF_anom(:,:,:,k)) *   &
+                (time-time_snap(k))/(time_snap(k+1)-time_snap(k))
+! exactement sur le snapshot et avec un ecart 1 an par rapport au precedent stockage
+!           write(6,*) 'time,tests',k,time,time-time_snap(k),time-time_prec-1.
+           if ((abs(time-time_snap(k)).le.dt).and.(abs(time-time_prec-1.).lt.dt)) then    
+              k_snap    = k
+              time_prec = time_snap(k)     ! time_prec est le temps du precedent
+           endif
+           exit
+        endif
+     end do
+  endif
+
+  thermal_forcing(:,:,:) = thermal_forcing_time(:,:,:)
+!  print*,'TF_ISMIP_RCM time_snap TF', k, time_snap(k), time_snap(k+1), thermal_forcing_time(163,89,1), thermal_forcing_time(194,131,1),TF_0(163,89,1),TF_0(194,131,1)
+
+end subroutine TF_ISMIP_RCM
+