!> \file lect_GrIce2sa_gen_nc.f90 !!Module pour la lecture de la topographie !< !> \namespace lect_topo_green_gen !! Module pour la lecture de la topographie !! \author Cat !! \date 31 oct 2011 !! @note Used module !! @note - use module3D_phy !! @note - use param_phy_mod !< module lect_topo_green_gen use module3D_phy use param_phy_mod use interface_input use io_netcdf_grisli character(len=100) :: topo_surf !< surface file name character(len=100) :: correc_surf !< ice-real correction file name ! character(len=100) :: topo_thick !< thickness file name character(len=100) :: topo_bed !< bedrock file name character(len=100) :: mask_grounded !< mask file name character(len=100) :: longitude !< longitude file name character(len=100) :: latitude !< latitude file name character(len=100) :: heatflux !< geothermal heat flux logical :: corr_firn !< eventuellement correction firn (pour l'instant surface dejà corrigee) real,dimension(nx,ny) :: work_tab character(len=100) :: file_ncdf !< fichier netcdf issue des fichiers .dat contains subroutine input_topo namelist/topo_groen_gen/topo_surf,correc_surf,topo_bed,mask_grounded,longitude,latitude,heatflux if (itracebug.eq.1) call tracebug(' Entree dans routine input_topo') 428 format(A) rewind(num_param) ! pour revenir au debut du fichier param_list.dat read(num_param,topo_groen_gen) write(num_rep_42,428)'!___________________________________________________________' write(num_rep_42,428)'! module lect_topo_ant_gen ' write(num_rep_42,topo_groen_gen) write(num_rep_42,428)'!___________________________________________________________' if ((trim(correc_surf).eq.'no').or.(trim(correc_surf).eq.'NO')) then corr_firn=.false. else corr_firn=.true. end if topo_surf = trim(dirnameinp)//trim(topo_surf) correc_surf = trim(dirnameinp)//trim(correc_surf) ! topo_thick = trim(dirnameinp)//trim(topo_thick) topo_bed = trim(dirnameinp)//trim(topo_bed) mask_grounded = trim(dirnameinp)//trim(mask_grounded) longitude = trim(dirnameinp)//trim(longitude) latitude = trim(dirnameinp)//trim(latitude) heatflux = trim(dirnameinp)//trim(heatflux) ! Transformation des fichiers .dat en .nc file_ncdf= trim(dirnameinp)//trim(runname)//'.nc' ! lecture adaptee aux fichiers ZBL.dat ou netcdf ou grd call lect_input(1,'Bsoc',1,Bsoc,topo_bed,file_ncdf) ! socle Bsoc0(:,:) = Bsoc(:,:) call lect_input(1,'S',1,S,topo_surf,file_ncdf) ! surface ! call lect_input(1,'H',1,H,topo_thick,file_ncdf) ! epaisseur calculee ensuite if (itracebug.eq.1) call tracebug(' avant corr firn') if (corr_firn) then ! density correction call lect_input(1,'work_tab1',1,work_tab,correc_surf,file_ncdf) where (H(:,:).gt.0.) work_tab(:,:)=min(H(:,:),work_tab(:,:)) S(:,:) = S(:,:) - work_tab(:,:) ! corrected surface H(:,:) = H(:,:) - work_tab(:,:) ! corrected thickness end where end if if (itracebug.eq.1) call tracebug(' apres corr firn') H(:,:)=max(0.,H(:,:)) H0(:,:)= H(:,:) S0(:,:)= S(:,:) mk0(:,:)=0. mk_init(:,:)=0. call lect_input(1,'work_tab2',2,work_tab,mask_grounded,file_ncdf) ! masque ! definition du Mk_init dans ice2sea !__________________________________________ ! outcrops -> eventuellement lu dans un autre fichier ! ocean libre -> 0 ! toundra -> 1 ! pose avec glace -> 2 ! epaisseur fixe -> 3 masque pour ne pas traiter ces points ! ice shelves -> 4 mk_init(:,:) = nint(work_tab(:,:)) if (itracebug.eq.1) call tracebug(' avant worktab') ! make consistency of input data files (flotation, thickness ...) ! Rappel : dans param_phy_mod ! coef_Sflot = (Row-Ro)/Row S = coef_Sflot * H + sealevel pour les shelves ! coef_Bflot = -Ro/Row B = coef_Bflot * H + sealevel pour les shelves do j=1,ny do i=1,nx if (mk_init(i,j).eq.0) then ! ocean H(i,j) = 1. ! valeur bidon pour l'ocean S(i,j) = coef_Sflot * H(i,j) B(i,j) = coef_Bflot * H(i,j) Bsoc(i,j) = -100. ! on surcreuse le socle flot(i,j) = .true. else if (mk_init(i,j).eq.1) then ! toundra S(i,j) = max(Bsoc(i,j),0.) H(i,j) = S(i,j) - Bsoc(i,j) B(i,j) = Bsoc(i,j) flot(i,j) = .false. else if (mk_init(i,j).eq.2) then ! glace H(i,j) = S(i,j) - Bsoc(i,j) B(i,j) = Bsoc(i,j) flot(i,j) = .false. else if (mk_init(i,j).eq.3) then ! zones a ne pas toucher ! le traitement "fixe" se fait dans prescribe_mod H(i,j) = S(i,j) - Bsoc(i,j) B(i,j) = Bsoc(i,j) flot(i,j) = .false. else if (mk_init(i,j).eq.4) then ! ice shelf H(i,j) = S(i,j) / coef_Sflot B(i,j) = S(i,j) - H(i,j) Bsoc(i,j) = min (B(i,j)-20.,Bsoc(i,j)) ! on surcreuse le socle si besoin flot(i,j)=.true. end if end do end do ! garder les cartes debut de run Bsoc0(:,:) = Bsoc(:,:) H0(:,:)= H(:,:) S(:,:) = B(:,:) +H(:,:) S0(:,:)= S(:,:) if (itracebug.eq.1) call tracebug(' apres consistency') ! pour l'Antarctique masque mko vrai partout (version 2006) MK0(:,:)=1 ! mis apres la lecture cptr call lect_input(1,'Xlong',1,Xlong,longitude,file_ncdf) !call lect_datfile(nx,ny,Xlong,1,longitude) ! lecture des coordonnées geographiques call lect_input(1,'Ylat',1,Ylat,latitude,file_ncdf) !call lect_datfile(nx,ny,Ylat,1,latitude) ! if (itracebug.eq.1) call tracebug(' apres lect lon/lat') ! pour recuperer le xmin xmax ymin ymax, on s'inspire de ce que Cat a fait pour l'antarc : ! write(6,*) 'topo_surf',topo_surf ! call Read_Ncdf_var('x',trim(topo_surf),Tab1D) ! lit le tableau x ! xmin = minval(Tab1D(:)) ! xmax = maxval(Tab1d(:)) ! call Read_Ncdf_var('y',trim(topo_surf),Tab1D) ! lit le tableau y ! ymin = minval(Tab1D(:)) ! ymax = maxval(Tab1d(:)) open(22,file=trim(dirnameinp)//'long_ZBL_CISM_5km.dat') ! juste pour avoir l'en tete read(22,*) i,i,xmin,xmax,ymin,ymax xmin=xmin/1000. ymin=ymin/1000. xmax=xmax/1000. ymax=ymax/1000. close(22) ! lecture du fichier de reference pour le calcul du niveau des mers : etat actuel ! pas de version correcte pour l'instant -> valeurs initiales S_sealev(:,:) = S0(:,:) H_sealev(:,:) = H0(:,:) B_sealev(:,:) = Bsoc(:,:) M_sealev(:,:) = Mk0(:,:) if (itracebug.eq.1) call tracebug(' apres lect entete') call lect_input(1,'ghf',1,ghf,heatflux,file_ncdf) ! flux geothermique ! for the purposes of SeaRISE, GTHF should be capped at a value of -0.07 w/m^2, ! and the field provided in the netCDF file does not include these updated values ! ghf(:,:)=min(ghf(:,:),70.) ! pour sea rise ! pour passer les flux des W/m2 au J/m2/an ghf(:,:)=-SECYEAR*ghf(:,:) ! les flux sont donnes en W/m2, correction pour passer a l'annee ghf0(:,:)=ghf(:,:) if (itracebug.eq.1) call tracebug(' a la fin d input_topo') end subroutine input_topo end module lect_topo_green_gen