Changeset 263 for trunk


Ignore:
Timestamp:
06/18/19 11:48:41 (5 years ago)
Author:
dumas
Message:

Option to use climate snapshots f(massb_time=2) for ISMIP simulations

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/SOURCES/Ant16_files/climat_InitMIP_years_perturb_mod.f90

    r237 r263  
    2424real :: T_lapse_rate              !< pour la temperature 
    2525 
    26 ! anomalies: smb and bmelt 
    27 real,dimension(nx,ny)    :: smb_anom            !> bilan de masse des anomalies indices : nx,ny 
    28  
     26! anomalies: smb and Tann 
     27real,dimension(:,:,:),allocatable   :: smb_anom    !> bilan de masse des anomalies indices : nx,ny,tsnap 
     28real,dimension(:,:,:),allocatable   :: tann_anom   !> Tann anomalies indices : nx,ny,tsnap 
     29real,dimension(:),allocatable       :: time_snap   !> date des snapshots  indice : nb de nb_snap 
     30real,dimension(:,:,:),allocatable   :: grad_bm     !> gradient vertical de smb indices : nx,ny,tsnap 
     31real,dimension(nx,ny)               :: lapse_rates !> gradient vertical de temperature 
    2932 
    3033! declaration pour le bilan de masse et le bmelt 
    31 real,dimension(nx,ny)                :: bm_0              !> bilan de masse initial 
    32  
    33 real                                 :: coef_smb_unit     ! pour corriger l'unite 
    34  
    35 real,dimension(nx,ny)                :: TA0               !> Temp annuelle sur S0 
     34real,dimension(nx,ny)    :: bm_0              !> bilan de masse initial 
     35 
     36real                     :: coef_smb_unit     ! pour corriger l'unite 
     37real                     :: coef_smb_anom_unit ! facteur correction unité anomalies de smb 
     38real,dimension(nx,ny)    :: TA0               !> Temp annuelle sur S0 
    3639 
    3740 
    3841! aurel, pour climat type perturb: 
    39 integer nft             !> nombre de lignes a lire dans le fichier forcage climatique 
    40 real,dimension(:),allocatable :: tdate          !> time for climate forcing 
    41 real,dimension(:),allocatable :: tpert          !> temperature for climate forcing 
    42 real,dimension(:),allocatable :: spert          !> sea surface perturbation 
    43 real :: coefT                    !> pour modifier l'amplitude de la perturb. T 
    44 character(len=120) :: filforc    !> nom du fichier forcage 
    45 integer :: pertsmb               !> boolean: do we modify the smb? 
    46 real    :: rapsmb                !> if we modify the smb this is the equivalent of rappact 
     42!integer nft             !> nombre de lignes a lire dans le fichier forcage climatique 
     43!real,dimension(:),allocatable :: tdate          !> time for climate forcing 
     44!real,dimension(:),allocatable :: tpert          !> temperature for climate forcing 
     45!real,dimension(:),allocatable :: spert          !> sea surface perturbation 
     46!real :: coefT                    !> pour modifier l'amplitude de la perturb. T 
     47!character(len=120) :: filforc    !> nom du fichier forcage 
     48!integer :: pertsmb               !> boolean: do we modify the smb? 
     49!real    :: rapsmb                !> if we modify the smb this is the equivalent of rappact 
    4750 
    4851 
    4952! pour les lectures ncdf 
    50 real*8, dimension(:,:),   pointer    :: tab               !< tableau 2d real pointer 
     53real*8, dimension(:), pointer        :: tab1D             !< tableau 1d real pointer 
     54real*8, dimension(:,:),   pointer    :: tab2D             !< tableau 2d real pointer 
     55real*8, dimension(:,:,:), pointer    :: tab3D             !< tableau 3d real pointer 
    5156Character(len=10), dimension(3)      :: dimtestname       !> pour la sortie test netcdf 
    5257integer                              :: ncidloc           !> pour la sortie test netcdf 
    … …  
    5459 
    5560integer                              :: massb_time        !< pour selectionner le type de calcul de smb 
    56 ! smb fixe ou en appliquant les anomalies 
     61 
     62integer                              :: nb_snap           !> nombre de snapshots 
     63 
    5764 
    5865contains 
    … …  
    7481 
    7582  !aurel for Tafor: 
    76   character(len=8) :: control      !label to check clim. forc. file (filin) is usable 
    77   character(len=80):: filin 
     83  character(len=8)   :: control      !label to check clim. forc. file (filin) is usable 
     84  character(len=80)  :: filin 
     85   
     86  real               :: time_depart_snaps  !> temps du debut premier snapshot 
    7887 
    7988  namelist/clim_smb_T_gen/smb_file,coef_smb_unit,temp_annual_file 
    80   namelist/smb_anom_initMIP/file_smb_anom,massb_time 
     89  namelist/smb_anom_initMIP/file_smb_anom,coef_smb_anom_unit,nb_snap,time_depart_snaps,massb_time 
     90   
     91!!!!!  namelist/clim_snap/nb_snap,file_smb_snap,massb_time_snap 
    8192 
    8293428 format(A) 
    … …  
    8596 
    8697  write(num_rep_42,428)'!________________________________________________________________'  
    87   write(num_rep_42,428)'!  module climat_InitMIP_years_mod lecture climat ref          ' 
     98  write(num_rep_42,428)'!  module climat_InitMIP_years_perturb_mod lecture climat ref          ' 
    8899  write(num_rep_42,clim_smb_T_gen) 
    89100  write(num_rep_42,428)'! smb_file          = fichier SMB (kg/m2/an)                     ' 
    … …  
    96107  smb_file  = trim(dirnameinp)//trim(smb_file) 
    97108 
    98   call Read_Ncdf_var('smb',smb_file,tab) 
    99  
    100  
    101   bm(:,:)  = tab(:,:) * coef_smb_unit 
     109  call Read_Ncdf_var('smb',smb_file,tab2D) 
     110 
     111 
     112  bm(:,:)  = tab2D(:,:) * coef_smb_unit 
    102113 
    103114  acc(:,:) = 0. 
    … …  
    115126  temp_annual_file = trim(dirnameinp)//trim(temp_annual_file) 
    116127 
    117   call Read_Ncdf_var('Tann',temp_annual_file,tab) 
    118   Tann(:,:)  = tab(:,:) 
     128  call Read_Ncdf_var('Tann',temp_annual_file,tab2D) 
     129  Tann(:,:)  = tab2D(:,:) 
    119130 
    120131  ta0(:,:)   = Tann(:,:) 
    121132  Tjuly(:,:) = Tann(:,:) 
    122  
    123  
    124   ! ______ Anomalies... 
    125    
     133  bm_0(:,:) = bm(:,:) 
     134   
     135  sealevel=0. 
     136  tafor=0. 
     137  sealevel_2d(:,:) = sealevel 
     138 
     139  ! ______ Anomalies de SMB (simule asmb initMIP)  
    126140  rewind(num_param)        ! pour revenir au debut du fichier param_list.dat 
    127141  read(num_param,smb_anom_initMIP) 
    128142 
    129   write(num_rep_42,428)'!_______________________________________________________________________'  
    130   write(num_rep_42,428)'!  module climat_InitMIP_years_mod                                      ' 
     143  write(num_rep_42,428)'!____________________________________________________________________________'  
     144  write(num_rep_42,428)'!  module climat_InitMIP_years_perturb_mod                                   ' 
    131145  write(num_rep_42,smb_anom_initMIP) 
    132   write(num_rep_42,428)'! file_smb_anom     = fichier anomalie SMB de GCM                       ' 
    133   write(num_rep_42,428)'! massb_time        = 0:fixe, 1:anomalies                               ' 
    134   write(num_rep_42,428)'!_______________________________________________________________________'  
    135  
    136  
     146  write(num_rep_42,428)'! file_smb_anom = fichier anomalie SMB de GCM                                ' 
     147  write(num_rep_42,428)'! coef_smb_anom_unit = coef passage m glace/an  (1/910 ou 1/918)             ' 
     148  write(num_rep_42,428)'! nb_snap       = nombre de snapshots                                        ' 
     149  write(num_rep_42,428)'! time_depart_snaps = debut du forçage                                       ' 
     150  write(num_rep_42,428)'! massb_time = 0:fixe, 1:anomalies, 2:interp snapshots, 3:snapsh+interp vert ' 
     151  write(num_rep_42,428)'!____________________________________________________________________________' 
     152   
    137153  file_smb_anom  = trim(dirnameinp)//trim(file_smb_anom) 
    138   call Read_Ncdf_var('asmb',file_smb_anom,tab) 
    139   smb_anom (:,:) = Tab(:,:) !* coef_smb_unit already in m/yr 
    140    
    141   bm_0(:,:) = bm(:,:) 
    142    
    143   filin=trim(dirforcage)//trim(filforc) 
    144  
    145   open(num_forc,file=filin,status='old') 
    146    
    147   read(num_forc,*) control,nft 
    148    
    149     ! Determination of file size (line nb), allocation of perturbation array 
    150    
    151   if (control.ne.'nb_lines') then 
    152      write(6,*) filin,'indiquer le nb de ligne en debut de fichier:' 
    153      write(6,*) 'le nb de lignes et le label de control nb_lines' 
    154      stop 
    155   endif 
     154   
     155  if (massb_time == 1) then ! fichier 2D de smb (exp initMIP asmb) 
     156! allocation dynamique de smb_snap 
     157      if (allocated(smb_anom)) then  
     158        deallocate(smb_anom,stat=err) 
     159        if (err/=0) then 
     160          print *,"Erreur à la desallocation de smb_anom",err 
     161          stop  
     162        end if 
     163      end if 
     164 
     165      allocate(smb_anom(nx,ny,1),stat=err) 
     166      if (err/=0) then 
     167        print *,"erreur a l'allocation du tableau smb_anom ",err 
     168        print *,"nx,ny,1 = ",nx,',',ny,',',1 
     169        stop  
     170      end if 
     171      call Read_Ncdf_var('asmb',file_smb_anom,tab2D) 
     172      smb_anom(:,:,1) = Tab2D(:,:) * coef_smb_anom_unit 
     173  elseif (massb_time >= 2) then ! fichier 3D de smb : lecture des snapshots 
     174! allocation dynamique de time_snap 
     175      if (allocated(time_snap)) then  
     176        deallocate(time_snap,stat=err) 
     177        if (err/=0) then 
     178          print *,"Erreur à la desallocation de time_snap",err 
     179          stop  
     180        end if 
     181      end if 
     182 
     183      allocate(time_snap(nb_snap),stat=err) 
     184      if (err/=0) then 
     185        print *,"erreur a l'allocation du tableau time_snap ",err 
     186        print *,"nb_snap = ",nb_snap 
     187        stop  
     188      end if 
     189! allocation dynamique de smb_snap 
     190      if (allocated(smb_anom)) then  
     191        deallocate(smb_anom,stat=err) 
     192        if (err/=0) then 
     193          print *,"Erreur à la desallocation de smb_anom",err 
     194          stop  
     195        end if 
     196      end if 
     197 
     198      allocate(smb_anom(nx,ny,nb_snap),stat=err) 
     199      if (err/=0) then 
     200        print *,"erreur a l'allocation du tableau smb_anom ",err 
     201        print *,"nx,ny,nb_snap = ",nx,',',ny,',',nb_snap 
     202        stop  
     203      end if 
     204 
     205! allocation dynamique de tann_snap 
     206      if (allocated(tann_anom)) then  
     207        deallocate(tann_anom,stat=err) 
     208        if (err/=0) then 
     209          print *,"Erreur à la desallocation de tann_anom",err 
     210          stop  
     211        end if 
     212      end if 
     213 
     214      allocate(tann_anom(nx,ny,nb_snap),stat=err) 
     215      if (err/=0) then 
     216        print *,"erreur a l'allocation du tableau tann_anom ",err 
     217        print *,"nx,ny,nb_snap = ",nx,',',ny,',',nb_snap 
     218        stop  
     219      end if 
     220   
     221! lecture de smb_snap 
     222      call Read_Ncdf_var('smb_anomaly',file_smb_anom,tab3D) 
     223      smb_anom(:,:,:) = tab3D(:,:,:) * coef_smb_anom_unit 
     224   
     225! lecture de tann_snap 
     226      call Read_Ncdf_var('ts_anomaly',file_smb_anom,tab3D) 
     227      tann_anom(:,:,:) = tab3D(:,:,:) 
     228 
     229! lecture de time_snap 
     230      call Read_Ncdf_var('time',file_smb_anom,tab1D) 
     231      time_snap(:) = tab1D(:) 
     232      time_snap(:) = time_snap(:) + time_depart_snaps 
    156233     
    157     ! Dimensionnement des tableaux tdate, .... 
    158  
    159   if (.not.allocated(tdate)) then 
    160      allocate(tdate(nft),stat=err) 
    161      if (err/=0) then 
    162         print *,"erreur a l'allocation du tableau Tdate",err 
    163         stop 4 
    164      end if 
    165   end if 
    166  
    167   if (.not.allocated(spert)) then 
    168      allocate(spert(nft),stat=err) 
    169      if (err/=0) then 
    170         print *,"erreur a l'allocation du tableau Spert",err 
    171         stop 4 
    172      end if 
    173   end if 
    174    
    175   if (.not.allocated(tpert)) then 
    176      allocate(tpert(nft),stat=err) 
    177      if (err/=0) then 
    178         print *,"erreur a l'allocation du tableau Tpert",err 
    179         stop 4 
    180      end if 
    181   end if 
    182    
    183   do i=1,nft 
    184      read(num_forc,*) tdate(i),spert(i),tpert(i) 
    185   end do 
    186   close(num_forc) 
    187    
    188   tpert(:)=tpert(:)*coefT 
    189    
     234      if (massb_time == 2) then ! lecture lapse rate 
     235      ! A coder....... 
     236      endif 
     237       
     238  endif      
    190239 
    191240    
    … …  
    200249  implicit none 
    201250  namelist/lapse_rates/T_lapse_rate 
    202   namelist/clim_pert_massb/coefT,filforc,pertsmb,rapsmb 
     251!  namelist/clim_pert_massb/coefT,filforc,pertsmb,rapsmb 
    203252   
    204253  rewind(num_param)        ! pour revenir au debut du fichier param_list.dat 
    … …  
    217266  write(num_rep_42,428)'!________________________________________________________________'  
    218267 
    219    
    220   rewind(num_param)   
    221   read(num_param,clim_pert_massb) 
    222    
    223   write(num_rep_42,428)'!___________________________________________________________' 
    224   write(num_rep_42,428) '&clim_pert                      ! module climat_perturb_mod ' 
    225   write(num_rep_42,*) 
    226   write(num_rep_42,*) 'coefT        = ', coefT 
    227   write(num_rep_42,'(A,A)') ' filforc      = ', filforc 
    228   write(num_rep_42,*) 'pertsmb        = ', pertsmb 
    229   write(num_rep_42,*) 'rapsmb         = ', rapsmb 
    230   write(num_rep_42,*)'/' 
    231   write(num_rep_42,*) 
    232  
    233  
    234 ! appelle la routine de lecture des smb annuels 
    235   call input_clim 
    236  
     268  rewind(num_param) 
    237269  return 
     270   
    238271end subroutine init_forclim 
    239272 
    … …  
    242275!!   
    243276!!  Routine qui permet le calcul climatique au cours du temps 
    244 !!  @note Au temps considere (time) attribue les scalaires 
    245 !!  @note   - tafor : forcage en temperature 
    246 !!  @note   - sealevel : forcage niveau des mers 
    247 !!  @note   - coefbmelt : forcage fusion basale ice shelves  
    248277!> 
     278 
    249279subroutine forclim               !  au temps considere (time)  
    250280 
    … …  
    253283  integer i,ift 
    254284  real             :: coefanomtime 
    255  
    256   coefanomtime = min ( real(time/40.) , 1. )  
    257   if (massb_time.eq.0) then 
    258      bm(:,:) = bm_0(:,:) 
    259   else if (massb_time.eq.1) then 
    260      bm(:,:) = bm_0(:,:) + coefanomtime * smb_anom(:,:) 
     285   
     286   
     287  if (massb_time == 0) then 
     288    bm(:,:) = bm_0(:,:) 
     289  elseif (massb_time == 1) then 
     290    coefanomtime = min ( real(time/40.) , 1. )  
     291    bm(:,:) = bm_0(:,:) + coefanomtime * smb_anom(:,:,1) 
     292  elseif (massb_time >= 2) then 
     293    call massb_ISMIP_RCM 
    261294  end if 
    262295   
    263   if(time.lt.tdate(1)) then 
    264      tafor=tpert(1) 
    265      sealevel=spert(1) 
    266      ift=1 
    267  
    268   else if (time.ge.tdate(nft)) then 
    269      tafor=tpert(nft) 
    270      sealevel=spert(nft) 
    271      ift=nft 
    272             
    273   else 
    274      do i=1,nft-1 
    275         if((time.ge.tdate(i)).and.(time.lt.tdate(i+1))) then ! entre i et i+1 : cas general 
    276            tafor=tpert(i)+(tpert(i+1)-tpert(i))*       & 
    277                 (time-tdate(i))/(tdate(i+1)-tdate(i)) 
    278            sealevel=spert(i)+(spert(i+1)-spert(i))*    & 
    279                 (time-tdate(i))/(tdate(i+1)-tdate(i)) 
    280            ift=i 
    281            goto 100 
     296  Ts(:,:)    = Tann(:,:) 
     297 
     298        
     299! coefmshelf est un coefficient qui fait varier bmgrz et bmshelf en fonction de tafor 
     300  coefbmshelf=1. 
     301       
     302end subroutine forclim 
     303 
     304subroutine massb_ISMIP_RCM                ! calcule le mass balance  
     305 
     306  implicit none 
     307  integer             :: k, k_snap               ! pour calculer les indices de temps 
     308  real :: time_prec 
     309  real,dimension(nx,ny) :: bm_time, tann_time ! pour calcul Bm et Tann 
     310 
     311! calcul smb a partir fichier snapshots massb_ISMIP_RCM 
     312! Calcule le mass balance d'apres un fichier de snapshots 
     313 
     314! calcule bm_time par interpolation entre deux snapshots 
     315! avant prend la valeur de reference 
     316! apres prend la derniere valeur 
     317! en general les snapshots vont de 1995 a 2100 (time_snap de 0 à 105) 
     318  if(time.lt.time_snap(1)) then              ! time avant le forcage 
     319     bm_time(:,:) = bm_0(:,:) 
     320     tann_time(:,:) = ta0(:,:) 
     321     k_snap       = 1 
     322!     S_ref(:,:)   = S(:,:)                   ! du coup sera la surface de reference avant le forcage 
     323     time_prec    = time 
     324  else if (time.ge.time_snap(nb_snap)) then  ! time apres le forcage 
     325     bm_time(:,:) =  bm_0(:,:) + smb_anom (:,:,nb_snap) 
     326     tann_time(:,:) = ta0(:,:) + tann_anom(:,:,nb_snap) 
     327     if (abs(time-time_prec-1.).lt.dt) then   !  
     328        time_prec = time_prec + 1 
     329     endif 
     330  else                                       ! cas general 
     331     do k = 1 , nb_snap-1 
     332        if((time.ge.time_snap(k)).and.(time.lt.time_snap(k+1))) then ! entre k et k+1  
     333           bm_time(:,:) = bm_0(:,:) + smb_anom(:,:,k) + (smb_anom(:,:,k+1)-smb_anom(:,:,k)) *   & 
     334                (time-time_snap(k))/(time_snap(k+1)-time_snap(k)) 
     335           tann_time(:,:) = ta0(:,:) + tann_anom(:,:,k) + (tann_anom(:,:,k+1)-tann_anom(:,:,k)) *   & 
     336                (time-time_snap(k))/(time_snap(k+1)-time_snap(k)) 
     337! exactement sur le snapshot et avec un ecart 1 an par rapport au precedent stockage 
     338!           write(6,*) 'time,tests',k,time,time-time_snap(k),time-time_prec-1. 
     339           if ((abs(time-time_snap(k)).le.dt).and.(abs(time-time_prec-1.).lt.dt)) then     
     340              k_snap    = k 
     341              time_prec = time_snap(k)     ! time_prec est le temps du precedent 
     342           endif 
     343           exit 
    282344        endif 
    283345     end do 
    284346  endif 
    285 100 continue 
    286  
    287   sealevel_2d(:,:) = sealevel 
    288  
    289   Tann (:,:) = Ta0 (:,:) + T_lapse_rate * (S(:,:)-S0(:,:)) +Tafor 
    290   Ts(:,:)    = Tann(:,:) 
    291  
    292      ! aurel marion dufresne: we might want to decrease the SMB during glacials..? 
    293   if (pertsmb.eq.1) then 
    294      bm(:,:) = bm_0(:,:) * exp( rapsmb *(Tann(:,:)-Ta0(:,:))) 
    295      if (Tafor.lt.0.) then 
    296         where(bm(:,:).lt.0.) bm(:,:)=min(bm(:,:)-Tafor*0.05,1.) !10 degrees less give 0.5 meter more ? 
    297      end if 
    298   end if 
    299          
    300      ! coefmshelf est un coefficient qui fait varier bmgrz et bmshelf en fonction de tafor 
    301   coefbmshelf=(1.+tafor/10.)       ! coefbmshelf=0 pour tafor=-10deg 
    302   coefbmshelf=max(coefbmshelf,0.)   
    303   coefbmshelf=min(coefbmshelf,2.) 
    304        
    305 end subroutine forclim 
     347 
     348  if (massb_time == 2) then ! pas d'interpolation verticale 
     349     bm(:,:) = bm_time(:,:) 
     350     Tann(:,:) = tann_time(:,:) 
     351  else if (massb_time == 3) then ! interpolation verticale 
     352     ! ajuste bm en fonction du temps et du gradient 
     353!     bm(:,:) = bm_time(:,:) + grad_bm(:,:,????) * (S(:,:) - S0(:,:)) 
     354     tann(:,:)= tann_time(:,:) + lapse_rates(:,:) * (S(:,:) - S0(:,:)) 
     355! finir de coder le cas lapse rates 
     356!     Tann(:,:) = tann_time(:,:) + lapse_rates(:,:)*(S(:,:) - S0(:,:)) 
     357!     write(6,897) time, time_prec, icum, i_moy 
     358!897  format('test temps smb   ',2(f0.3,1x),2(i0,1x)) 
     359  endif 
     360 
     361!~ ! garde les 10 dernieres annees et calcule la moyenne 
     362!~   if (icum.eq.1) then                           ! stockage dans le tableau bm_ref_10 
     363!~      do k = 9,1,-1 
     364!~         bm_ref_10(:,:,k+1) = bm_ref_10(:,:,k)   ! on decale tous les elements 
     365!~      end do 
     366!~      bm_ref_10(:,:,k+1) = bm(:,:)               ! le plus recent est en position 1 
     367!~      i_moy              = i_moy +1              ! compte combien il y en a pour la moyenne 
     368!~      i_moy              = min(i_moy,10) 
     369!~      bm_ref_t(:,:)   = 0. 
     370!~      do k = 1,i_moy 
     371!~         bm_ref_t(:,:) = bm_ref_t(:,:) + bm_ref_10(:,:,k) 
     372!~      end do 
     373!~      bm_ref_t(:,:)    = bm_ref_t(:,:)/i_moy 
     374!~      write(6,898) time, time_prec, icum, i_moy 
     375!~ 898  format('cumul pour gradient  ',2(f0.3,1x),2(i0,1x)) 
     376!~   end if 
     377end subroutine massb_ISMIP_RCM 
     378 
     379 
    306380 
    307381 
Note: See TracChangeset for help on using the changeset viewer.