Changeset 15


Ignore:
Timestamp:
06/10/15 15:39:11 (9 years ago)
Author:
dumas
Message:

Greeneem15 : nouveau fichier de topo (Bamber 2013 via fichier Tamsin 5km) et climat-forcage-insolation_mod_oneway.f90 pour utilisation PLIOMIP

Location:
trunk
Files:
8 added
8 edited

Legend:

Unmodified
Added
Removed

  • trunk/SOURCES/Fichiers-parametres/LISTE-VAR-NETCDF.dat

    r12 r15  
    288288Nefmx            
    28928968 
    290 1       1       1                                        
     2901       1       1.e-5                                    
    291291----------------------------------------------------------- 
    292292Nefmy            
    29329369 
    294 1       2       1                                        
     2941       2       1.e-5                                    
    295295----------------------------------------------------------- 
    296296posx             

  • trunk/SOURCES/Fichiers-parametres/TEMPS-CPTR-NC.dat

    r4 r15  
    44 
    55~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ fin des commentaires~~~~~~~~~~~ 
    6 100.          Dtcpt 
     61000.          Dtcpt 
    77------------------------------------------------------------------------------------------- 
    8 0           Nbr des temps de sortie cptr ou nc 
    9 20. 
    10 100. 
    11 200 
     815           Nbr des temps de sortie cptr ou nc 
     9-126000. 
    1210-20000.      
     11-10000. 
     1210000. 
    131320000. 
    14 23000. 
    15 22000. 
    161430000. 
    171540000. 
    … …  
    191760000. 
    201870000. 
     1980000. 
     2090000. 
    2121100000. 
    22 -126000.    
    23 140000. 
    2422150000. 
    25 180000. 
    26 -297000.         
    27 -380000.     
    28 -350000.     
    29 -390000.     
    30 -400000.      
    31 -410000.      
    32 -420000. 
    33 -430000. 
    34 -435000. 
    35 -440000. 
    36 -450000. 
    37 -470000. 
    38 -490000. 
     23200000. 
    3924------------------------------------------------------------------------------------------- 

  • trunk/SOURCES/GrIce2sea_files/climat_GrIce2sea_years_mod.f90

    r14 r15  
    112112!  end where 
    113113 
    114 !cdc test debug Hemin15 
    115 !  where (bm(:,:).lt.-20) bm(:,:)=-10. 
     114!cdc test debug Hemin15 et Greeneem15 
     115!  where (bm(:,:).lt.-1000) bm(:,:)=0. 
    116116 
    117117 

  • trunk/SOURCES/Greeneem_files/lect-greeneem_mod.f90

    r13 r15  
    33use module3D_phy 
    44use param_phy_mod   
     5!cdc 
     6use interface_input 
     7use io_netcdf_grisli 
    58 
    69character(len=100) :: topo_dep       ! Topo de départ 
    … …  
    1215!real :: ghf0                             ! flux geothermique constant en mW/m2 
    1316 
     17real, dimension(nx,ny) :: Bsoc_bamber, S_bamber, Bsoc_new 
     18character(len=100) :: topo_surf      !< surface file name 
     19character(len=100) :: topo_bed       !< bedrock file name 
     20character(len=100) :: file_ncdf      !< fichier netcdf issue des fichiers .dat 
     21 
    1422integer :: n1,n2,ndx 
    1523 
    1624contains 
    1725  
    18 subroutine input_topo 
     26  subroutine input_topo 
    1927 
    20 namelist/toponeem/topo_ref,topo_dep,grid_topo,ghf_fich 
     28    namelist/toponeem/topo_ref,topo_dep,grid_topo,ghf_fich 
    2129 
    2230428 format(A) 
    2331 
    24 rewind(num_param)        ! pour revenir au debut du fichier param_list.dat 
     32    rewind(num_param)        ! pour revenir au debut du fichier param_list.dat 
    2533 
    26 read(num_param,toponeem) 
     34    read(num_param,toponeem) 
    2735 
    28 write(num_rep_42,428)'!___________________________________________________________'  
    29 write(num_rep_42,428) '&toponeen                     ! module lect_topo_greeneem' 
    30 write(num_rep_42,'(A,A)')   'topo_ref  =', topo_ref 
    31 write(num_rep_42,'(A,A)')   'topo_dep  =', topo_dep 
    32 write(num_rep_42,'(A,A)')   'grid_topo =', grid_topo 
    33 write(num_rep_42,*)         'ghf_fich  =', ghf_fich 
    34 write(num_rep_42,*)'/'                       
    35 write(num_rep_42,428) '! topo_ref= topo actuelle' 
    36 write(num_rep_42,428) '! ghf0 flux geothermique en mW/m2' 
    37 write(num_rep_42,428) '! grid_topo : fichier i,j,x,y,lon,lat' 
    38 write(num_rep_42,428) '! ghf_fich  : fichier flux geothermique' 
    39 write(num_rep_42,*) 
     36    write(num_rep_42,428)'!___________________________________________________________'  
     37    write(num_rep_42,428) '&toponeen                     ! module lect_topo_greeneem' 
     38    write(num_rep_42,'(A,A)')   'topo_ref  =', topo_ref 
     39    write(num_rep_42,'(A,A)')   'topo_dep  =', topo_dep 
     40    write(num_rep_42,'(A,A)')   'grid_topo =', grid_topo 
     41    write(num_rep_42,*)         'ghf_fich  =', ghf_fich 
     42    write(num_rep_42,*)'/'                       
     43    write(num_rep_42,428) '! topo_ref= topo actuelle' 
     44    write(num_rep_42,428) '! ghf0 flux geothermique en mW/m2' 
     45    write(num_rep_42,428) '! grid_topo : fichier i,j,x,y,lon,lat' 
     46    write(num_rep_42,428) '! ghf_fich  : fichier flux geothermique' 
     47    write(num_rep_42,*) 
    4048 
    4149 
    42 topo_ref=trim(dirnameinp)//trim(topo_ref)  
    43 topo_dep=trim(dirnameinp)//trim(topo_dep)  
    44 grid_topo=trim(dirnameinp)//trim(grid_topo) 
    45 ghf_fich=trim(dirnameinp)//trim(ghf_fich)!trim(ghf0) 
     50    topo_ref=trim(dirnameinp)//trim(topo_ref)  
     51    topo_dep=trim(dirnameinp)//trim(topo_dep)  
     52    grid_topo=trim(dirnameinp)//trim(grid_topo) 
     53    ghf_fich=trim(dirnameinp)//trim(ghf_fich)!trim(ghf0) 
    4654 
    47 sealevel0=0.  ! voir a passer dans le fichier parametre 
     55    sealevel0=0.  ! voir a passer dans le fichier parametre 
    4856 
    49 !     lecture de la topo de référence 
    50 !     ------------------------------- 
    51 ! Cette topo sert a calculer le socle de reference pour l'isostasie 
    52 ! voir init_iso et a avoir une surface de reference pour les temperatures 
     57    !     lecture de la topo de référence 
     58    !     ------------------------------- 
     59    ! Cette topo sert a calculer le socle de reference pour l'isostasie 
     60    ! voir init_iso et a avoir une surface de reference pour les temperatures 
     61    ! lecture adaptee aux fichiers ZBL.dat ou netcdf ou grd 
     62    call lect_input(1,'Bsoc',1,Bsoc0,topo_ref,file_ncdf)    ! socle 
     63    call lect_input(1,'S',1,S0,topo_ref,file_ncdf)          ! surface 
     64    call lect_input(1,'H',1,H0,topo_ref,file_ncdf)          ! epaisseur 
    5365 
    54      open (20,file=topo_ref) 
     66    S0(:,:)=max(S0(:,:),sealevel0)   ! pour etre au niveau des mers : ATTENTION si SEALEV <0 
     67    mk0(:,:) = 1                   ! mk0=0 pour les zones interdites 
    5568 
    56      read(20,*)  
    57      read(20,*)  
    58      read(20,*) n1,n2,ndx 
    59      read(20,*)  
     69    !     lecture de la topo de depart 
     70    !     --------------------------- 
     71    ! lecture adaptee aux fichiers ZBL.dat ou netcdf ou grd 
     72    call lect_input(1,'Bsoc',1,Bsoc,topo_dep,file_ncdf)    ! socle 
     73    call lect_input(1,'S',1,S,topo_dep,file_ncdf)          ! surface 
     74    call lect_input(1,'H',1,H,topo_dep,file_ncdf)          ! epaisseur 
     75    S(:,:)=max(S(:,:),0.)   ! pour etre au niveau des mers : ATTENTION si SEALEV <0 
     76    H(:,:)=max(H(:,:),1.)   ! pour avoir au moins 1 m 
    6077 
    61      if ((n1.eq.nx).and.(n2.eq.ny).and.(ndx.eq.nint(dx/1000.))) then 
    62                                ! le fichier correspond bien a  
    63                                ! la grille definie dans paradim 
    64          do J=1,ny  
    65           do I=1,nx 
    66               read (20,*) S0(I,J),H0(I,J),BSOC0(I,J) 
    67               S0(i,j)=max(S0(i,j),sealevel0)   ! pour etre au niveau des mers : ATTENTION si SEALEV <0 
    68           end do 
    69         end do 
    70    
    71      else 
    72         write(6,*) 'le fichier ne correspond pas a la grille definie' 
    73         STOP 
    74      end if 
    75      close(20) 
     78           
     79    ! calcul de flottaison 
     80    do j=1,ny 
     81       do i=1,nx 
    7682 
    77      mk0(:,:) = 1                   ! mk0=0 pour les zones interdites 
    78       
    79 !     lecture de la topo de depart 
    80 !     --------------------------- 
    81      open (20,file=topo_dep) 
     83          archim = Bsoc(i,j)+H(i,j)*ro/row -sealevel 
    8284 
    83      read(20,*)  
    84      read(20,*)  
    85      read(20,*) n1,n2,ndx 
    86      read(20,*)  
     85          if ((ARCHIM.LT.0.).and.(H(I,J).gt.1.E-3)) then    ! le point flotte 
     86             flot(i,j)=.true. 
     87          else 
     88             flot(i,j)=.false. 
     89          end if 
     90 
     91       end do 
     92    end do 
     93 
     94    ! lecture des coordonnées geographiques 
    8795 
    8896 
    89      if ((n1.eq.nx).and.(n2.eq.ny).and.(ndx.eq.nint(dx/1000.))) then 
    90                                ! le fichier correspond bien a  
    91                                ! la grille definie dans paradim 
    92    
    93          do J=1,ny  
    94           do I=1,nx 
    95               read (20,*) S(I,J),H(I,J),BSOC(I,J) 
    96               S(i,j)=max(S(i,j),0.)   ! pour etre au niveau des mers : ATTENTION si SEALEV <0 
    97               H(i,j)=max(H(i,j),1.)   ! pour avoir au moins 1 m 
     97    ! les fichiers EISMINT etaient ranges i,j,lat,lon 
     98    ! je n'arrive pas a retrouver la projection exacte. 
     99    ! le x est donc calcule en partant de 0, 0 
     100 
     101    open(unit=20,file=grid_topo,iostat=ios) 
     102    read(20,*)  
     103    read(20,*)  
     104    read(20,*) n1,n2,ndx 
     105    read(20,*)  
     106 
     107    if ((n1.eq.nx).and.(n2.eq.ny).and.(ndx.eq.nint(dx/1000.))) then 
     108       ! le fichier correspond bien a  
     109       ! la grille definie dans paradim 
     110 
     111       do k=1,nx*ny 
     112          read(20,*) i,j,xcc(i,j),ycc(i,j),xlong(i,j),ylat(i,j) 
     113       enddo 
     114       ! pour assurer la continuite des temperatures parametrees : 
     115       where(xlong(:,:).lt.100.)     
     116          xlong(:,:)=xlong(:,:)+360. 
     117       end where 
     118 
     119    else 
     120       write(6,*) 'le fichier ne correspond pas a la grille definie' 
     121       STOP 
     122    end if 
     123    close(20) 
     124 
     125    ! xcc est en m 
     126    xmin=xcc(1,1) 
     127    ymin=ycc(1,1) 
     128    do j=1,ny 
     129       do i=1,nx 
     130          xcc(i,j)=xcc(i,j)*1000. 
     131          ycc(i,j)=ycc(i,j)*1000. 
     132       end do 
     133    end do 
     134    xmax=xcc(nx,ny)/1000 
     135    ymax=ycc(nx,ny)/1000 
     136 
     137    ! lecture de la carte de flux geothermique 
     138    open (20,file=ghf_fich) 
     139    read(20,*) n1,n2,ndx 
     140    if ((n1.eq.nx).and.(n2.eq.ny).and.(ndx.eq.nint(dx/1000.))) then 
     141       ! le fichier correspond bien a  
     142       ! la grille definie dans paradim 
     143       do j=1,ny 
     144          do i=1,nx 
     145             read(20,*) ghf(i,j) 
    98146          end do 
    99         end do 
    100    
    101      else 
    102         write(6,*) 'le fichier ne correspond pas a la grille definie' 
    103         STOP 
    104      end if 
    105      close(20) 
    106  
    107 ! calcul de flottaison 
    108 do j=1,ny 
    109    do i=1,nx 
    110  
    111       archim = Bsoc(i,j)+H(i,j)*ro/row -sealevel 
    112  
    113       if ((ARCHIM.LT.0.).and.(H(I,J).gt.1.E-3)) then    ! le point flotte 
    114          flot(i,j)=.true. 
    115       else 
    116          flot(i,j)=.false. 
    117       end if 
    118      
    119    end do 
    120 end do 
    121  
    122 ! lecture des coordonnées geographiques 
     147       end do 
     148       close(20) 
     149    else 
     150       write(6,*) 'le fichier ne correspond pas a la grille definie' 
     151       STOP 
     152    end if 
    123153 
    124154 
    125 ! les fichiers EISMINT etaient ranges i,j,lat,lon 
    126 ! je n'arrive pas a retrouver la projection exacte. 
    127 ! le x est donc calcule en partant de 0, 0 
    128  
    129 open(unit=20,file=grid_topo,iostat=ios) 
    130 read(20,*)  
    131 read(20,*)  
    132 read(20,*) n1,n2,ndx 
    133 read(20,*)  
    134  
    135 if ((n1.eq.nx).and.(n2.eq.ny).and.(ndx.eq.nint(dx/1000.))) then 
    136    ! le fichier correspond bien a  
    137    ! la grille definie dans paradim 
    138     
    139    do k=1,nx*ny 
    140       read(20,*) i,j,xcc(i,j),ycc(i,j),xlong(i,j),ylat(i,j) 
    141    enddo 
    142    ! pour assurer la continuite des temperatures parametrees : 
    143    where(xlong(:,:).lt.100.)     
    144       xlong(:,:)=xlong(:,:)+360. 
    145    end where 
    146     
    147 else 
    148    write(6,*) 'le fichier ne correspond pas a la grille definie' 
    149    STOP 
    150 end if 
    151 close(20) 
    152  
    153 ! xcc est en m 
    154 xmin=xcc(1,1) 
    155 ymin=ycc(1,1) 
    156 do j=1,ny 
    157    do i=1,nx 
    158       xcc(i,j)=xcc(i,j)*1000. 
    159       ycc(i,j)=ycc(i,j)*1000. 
    160    end do 
    161 end do 
    162 xmax=xcc(nx,ny)/1000 
    163 ymax=ycc(nx,ny)/1000 
    164  
    165 ! lecture de la carte de flux geothermique 
    166 open (20,file=ghf_fich) 
    167 read(20,*) n1,n2,ndx 
    168 if ((n1.eq.nx).and.(n2.eq.ny).and.(ndx.eq.nint(dx/1000.))) then 
    169                                ! le fichier correspond bien a  
    170                                ! la grille definie dans paradim 
    171    do j=1,ny 
    172       do i=1,nx 
    173          read(20,*) ghf(i,j) 
    174       end do 
    175    end do 
    176    close(20) 
    177 else 
    178    write(6,*) 'le fichier ne correspond pas a la grille definie' 
    179    STOP 
    180 end if 
    181  
    182  
    183 !  flux geothermique  : pour passer les flux des mW/m2 aux J/m2/an      
    184 ghf(:,:)=-secyear*ghf(:,:) 
    185 !ghf(:,:)=-secyear*50.0/1000.   !ghf0 ! pour flux constant 
     155    !  flux geothermique  : pour passer les flux des mW/m2 aux J/m2/an      
     156    ghf(:,:)=-secyear*ghf(:,:) 
     157    !ghf(:,:)=-secyear*50.0/1000.   !ghf0 ! pour flux constant 
    186158 
    187159 
    188160 
    189 ! calcul des courbures du socle 
    190 call courbure(nx,ny,dx,Bsoc,bidon(:,:,1),bidon(:,:,2),bidon(:,:,3), & 
    191      bidon(:,:,4),socle_cry,bidon(:,:,5)) 
    192 socle_cry(:,:)=socle_cry(:,:)*dx*dx 
    193      
    194 !------------------------------------------------      
    195 end subroutine input_topo 
     161    ! calcul des courbures du socle 
     162    call courbure(nx,ny,dx,Bsoc,bidon(:,:,1),bidon(:,:,2),bidon(:,:,3), & 
     163         bidon(:,:,4),socle_cry,bidon(:,:,5)) 
     164    socle_cry(:,:)=socle_cry(:,:)*dx*dx 
     165 
     166    !------------------------------------------------      
     167  end subroutine input_topo 
    196168 
    197169 

  • trunk/SOURCES/Greeneem_files/module_choix-greeneem.f90

    r14 r15  
    2626!--------------Module climat --------------- 
    2727!use climat_forcage_mois_mod ! forcage mensuel GCM 1 Snapshot Fev 2015 
    28 use climat_Grice2sea_years_mod ! forcage avec SMB type MAR Fev2015 
     28!use climat_Grice2sea_years_mod ! forcage avec SMB type MAR Fev2015 
     29!use climat_forcage_mois_mod ! climat mensuel GCM, possibilite d'utiliser plusieurs snapshots avec fichier d'evolution (Pas de correction topo) 
     30use climat_forcage_insolation_mod ! Methode JB plusieurs snaphots, chacun avec topo + courbe insolation et CO2 
    2931 
    3032! Anciens modules pas encore valides 
    … …  
    3739!!!!!!use climat_regions_delta 
    3840 
    39 !use ablation_mod ! calcul de l'ablation (PDD ou autre methode) 
    40 use no_ablation  ! pas de calcul de l'ablation => lecture fichier SMB (necessaire avec climat_Grice2sea_years_mod) 
     41use ablation_mod ! calcul de l'ablation (PDD ou autre methode) 
     42!use no_ablation  ! pas de calcul de l'ablation => lecture fichier SMB (necessaire avec climat_Grice2sea_years_mod) 
    4143 
    4244! pas de lacs proglaciaires 

  • trunk/SOURCES/Makefile.grisli.inc

    r14 r15  
    4141 
    4242# module de forcage climatique C. Dumas 
    43 mod_clim_tof = climat_forcage_mois_mod.o climat_GrIce2sea_years_mod.o \ 
     43mod_clim_tof = climat_forcage_mois_mod.o climat-forcage-insolation_mod_oneway.o \ 
     44        climat_GrIce2sea_years_mod.o \ 
    4445        ablation_mod.o no_ablation_mod.o  
    4546 
    … …  
    612613        $(FT) climat_forcage_mois_mod.f90 
    613614 
     615climat-forcage-insolation_mod_oneway.o : climat-forcage-insolation_mod_oneway.f90 
     616        $(FT) climat-forcage-insolation_mod_oneway.f90 
     617 
    614618ablation_mod.o : ablation_mod.f90 
    615619        $(FT) ablation_mod.f90 

  • trunk/SOURCES/main3D-0.4-40km.f90

    r10 r15  
    370370  call testsort_time(tbegin) 
    371371  call sortie_hz_multi   ! pour les variables declarees dans 3D 
    372   call hz_output(tbegin) 
     372!  call hz_output(tbegin) 
    373373 
    374374!  call limit_file(1,real(time),dt,tend,dtsortie,dtcpt,testdiag,dtt,runname) 

  • trunk/SOURCES/steps_time_loop.f90

    r10 r15  
    2424 
    2525  use module_choix ! module de choix du type de run 
    26                    !  module_choix donne acces a tous les modules 
    27                    !  de declaration des packages 
     26  !  module_choix donne acces a tous les modules 
     27  !  de declaration des packages 
    2828 
    2929  use sorties_ncdf_grisli 
    … …  
    3131  use interface_icetempmod 
    3232  use diagno_mod  
    33 !  use track_debug  
     33  !  use track_debug  
    3434 
    3535  implicit none 
    … …  
    161161  !   sortie compteur tous les dtcpt ans  
    162162  !------------------------------------------------------------------ 
    163      !iout == 1 sortie cptr 
    164      !iout == 2 sortie nc 
    165     if (itracebug.eq.1)  call tracebug('dans steps_time_loop avant call out_recovery ') 
    166  
    167      call out_recovery(iout) 
     163  !iout == 1 sortie cptr 
     164  !iout == 2 sortie nc 
     165  if (itracebug.eq.1)  call tracebug('dans steps_time_loop avant call out_recovery ') 
     166 
     167  call out_recovery(iout) 
    168168 
    169169  ! end of outputs 
Note: See TracChangeset for help on using the changeset viewer.