!> \file climat_GrIce2sea_mod.f90 ! forcage avec BM !! Module ou les variations temporelles des variables climatiques !! sont directement imposees !< !> \namespace climat_Grice2sea_mod !! Module ou les variations temporelles des variables climatiques !! sont directement imposees !! \author Cat !! \date 31 oct !! @note Used modules: !! @note - use module3D_phy !< module climat_Grice2sea_years_perturb_mod use module3d_phy,only: nx,ny,S,S0,H,Bsoc,acc,abl,BM,Tann,Tjuly,Ts,time,dt,num_param,num_rep_42,num_forc,dirforcage,dirnameinp,tafor,time,sealevel !use lect_climref_Ice2sea use netcdf use io_netcdf_grisli implicit none real :: T_lapse_rate !< pour la temperature ! declarations specifiques year by year real,dimension(:),allocatable :: time_snap !> date des snapshots indice : nb de nb_snap) real,dimension(:,:,:),allocatable :: smb_snap !> bilan de masse des snapshots indices : nx,ny,nb_snap real :: ecart_snap !> ecart temporel entre les snapshots real :: time_depart_snaps !> temps du debut premier snapshot integer :: nb_snap !> nombre de snapshots ! declaration pour le bilan de masse real,dimension(nx,ny) :: bm_0 !> bilan de masse initial equ. S_0ref de Tamsin real,dimension(nx,ny) :: bm_time !> bilan de masse interpole entre 2 snapshots (~ S_t^RCM) ! calcul du gradient real,dimension(nx,ny) :: bm_ref_t !> bilan de masse de reference au temps t real,dimension(nx,ny,10) :: bm_ref_10 !> tableau des bm_ref pour moyenne 10 ans real :: time_prec !> temps du precedent snapshot integer :: icum !> pour le test stockage dans bm_ref_10 integer :: i_moy !> nombre de snapshots stockes real,dimension(nx,ny) :: grad_bm !> gradient du bilan de masse real,dimension(nx,ny) :: S_ref !> surface de reference real :: grad_N_smb_neg !> SMB < 0 North real :: grad_N_smb_pos !> SMB >= 0 North real :: grad_S_smb_neg !> SMB < 0 South real :: grad_S_smb_pos !> SMB >= 0 South real :: coef_smb_unit ! pour corriger l'unite real,dimension(nx,ny) :: TA0 !> Temp annuelle sur S0 ! aurel, pour climat type perturb: integer nft !< nombre de lignes a lire dans le fichier forcage climatique real,dimension(:),allocatable :: tdate !< time for climate forcing real,dimension(:),allocatable :: tpert !< temperature for climate forcing real,dimension(:),allocatable :: spert !< sea surface perturbation real :: coefT !< pour modifier l'amplitude de la perturb. T character(len=120) :: filforc !< nom du fichier forcage ! pour les lectures ncdf real*8, dimension(:,:,:), pointer :: tab3D !< tableau 3d real pointer real*8, dimension(:,:), pointer :: tab !< tableau 2d real pointer Character(len=10), dimension(3) :: dimtestname !> pour la sortie test netcdf integer :: ncidloc !> pour la sortie test netcdf integer :: status !> pour la sortie test netcdf integer :: massb_time !< pour selectionner le type de calcul de smb ! On a deux routines : celle avec un seul fichier (donnees observees) : massb_Ice2sea_fixe ! Celle avec le bilan de masse sur plusieurs snapshots (annuels par ex.) entre lesquels on interpole. contains !-------------------------------------------------------------------------------- !> SUBROUTINE: input_clim !! Routine qui permet d'initialiser les variations temporelles des variables climatiques !> !-------------------------------------------------------------------------------- subroutine input_clim character(len=100) :: smb_file ! fichier smb character(len=100) :: temp_annual_file ! temperature annuelles character(len=100) :: file_smb_snap !> nom du fichier dans lequel sont les snapshots integer :: err ! recuperation d'erreur integer :: i,j,k !aurel for Tafor: character(len=8) :: control !label to check clim. forc. file (filin) is usable character(len=80):: filin namelist/clim_smb_T_gen/smb_file,coef_smb_unit,temp_annual_file namelist/clim_snap/nb_snap,time_depart_snaps,ecart_snap,file_smb_snap,massb_time 428 format(A) rewind(num_param) ! pour revenir au debut du fichier param_list.dat read(num_param,clim_smb_T_gen) write(num_rep_42,428)'!________________________________________________________________' write(num_rep_42,428)'! module climat_Grice2sea_years_mod lecture climat ref ' write(num_rep_42,clim_smb_T_gen) write(num_rep_42,428)'! smb_file = fichier SMB (kg/m2/an) ' write(num_rep_42,428)'! coef_smb_unit = coef passage m glace/an (1/910 ou 1/918) ' write(num_rep_42,428)'! temp_annual_file = Temp moy annuelle (°C) ' write(num_rep_42,428)'!________________________________________________________________' ! smb : surface mass balance smb_file = trim(dirnameinp)//trim(smb_file) call Read_Ncdf_var('smb',smb_file,tab) ! call lect_input(3,'smb',1,bm,smb_file,trim(dirnameinp)//trim(runname)//'.nc') bm(:,:) = tab(:,:) * coef_smb_unit ! where ((H(:,:).lt.1.).and.(Bsoc(:,:).gt.0.)) ! bm(:,:) = bm(:,:) - 5. ! pour faire un masque a l'exterieur du Groenland actuel ! end where !cdc test debug Hemin15 et Greeneem15 ! where (bm(:,:).lt.-1000) bm(:,:)=0. acc(:,:) = 0. abl(:,:) = 0. where (bm(:,:).gt.0.) acc(:,:) = bm(:,:) ! accumulation quand positif elsewhere abl(:,:) = - bm(:,:) ! ablation quand negatif end where ! surface temperature Tann temp_annual_file = trim(dirnameinp)//trim(temp_annual_file) call Read_Ncdf_var('Tann',temp_annual_file,tab) Tann(:,:) = tab(:,:) ! call lect_input(3,'Tann',1,Tann,temp_annual_file,trim(dirnameinp)//trim(runname)//'.nc') !cdc test debug Hemin15 ! where (Tann(:,:).lt.-100.) Tann(:,:)=10. ta0(:,:) = Tann(:,:) Tjuly(:,:) = Tann(:,:) rewind(num_param) ! pour revenir au debut du fichier param_list.dat read(num_param,clim_snap) ! formats pour les ecritures dans 42 rewind(num_param) ! pour revenir au debut du fichier param_list.dat read(num_param,clim_snap) write(num_rep_42,428)'!_______________________________________________________________________' write(num_rep_42,428)'! module climat_Grice2sea_years_mod ' write(num_rep_42,clim_snap) write(num_rep_42,428)'! nb_snap = nombre de snapshots ' write(num_rep_42,428)'! time_depart_snaps = debut du forçage ' write(num_rep_42,428)'! ecart_snap = ecart entre les snapshots ' write(num_rep_42,428)'! file_smb_snap = fichier serie temp anomalie SMB de GCM ' write(num_rep_42,428)'! massb_time = 0:fixe, 1:interp snapshots, 2:snapsh+interp vert ' write(num_rep_42,428)'!_______________________________________________________________________' if (massb_time == 1) then ! 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 ! remplissage de time_snap !write(6,*) 'time_snap' do i=1,nb_snap time_snap(i) = time_depart_snaps + ecart_snap * (i-1) ! write(6,*) i,time_snap(i) end do ! allocation dynamique de smb_snap if (allocated(smb_snap)) then deallocate(smb_snap,stat=err) if (err/=0) then print *,"Erreur à la desallocation de smb_snap",err stop end if end if allocate(smb_snap(nx,ny,nb_snap),stat=err) if (err/=0) then print *,"erreur a l'allocation du tableau smb_snap ",err print *,"nx,ny,nb_snap = ",nx,',',ny,',',nb_snap stop end if ! lecture de smb_snap file_smb_snap = trim(dirnameinp)//trim(file_smb_snap) call Read_Ncdf_var('z',file_smb_snap,tab3D) smb_snap (:,:,:) = Tab3D(:,:,:) * coef_smb_unit ! ce sont des anomalies : ajoute les valeurs de reference do k = 1,nb_snap do j = 1,ny do i = 1,nx smb_snap (i,j,k) = smb_snap(i,j,k) + bm(i,j) end do end do end do ! copie la valeur de reference dans bm_0 bm_0(:,:) = bm(:,:) endif filin=trim(dirforcage)//trim(filforc) open(num_forc,file=filin,status='old') read(num_forc,*) control,nft ! Determination of file size (line nb), allocation of perturbation array if (control.ne.'nb_lines') then write(6,*) filin,'indiquer le nb de ligne en debut de fichier:' write(6,*) 'le nb de lignes et le label de control nb_lines' stop endif ! Dimensionnement des tableaux tdate, .... if (.not.allocated(tdate)) then allocate(tdate(nft),stat=err) if (err/=0) then print *,"erreur a l'allocation du tableau Tdate",err stop 4 end if end if if (.not.allocated(spert)) then allocate(spert(nft),stat=err) if (err/=0) then print *,"erreur a l'allocation du tableau Spert",err stop 4 end if end if if (.not.allocated(tpert)) then allocate(tpert(nft),stat=err) if (err/=0) then print *,"erreur a l'allocation du tableau Tpert",err stop 4 end if end if do i=1,nft read(num_forc,*) tdate(i),spert(i),tpert(i) end do close(num_forc) tpert(:)=tpert(:)*coefT ! ecriture de verification !file_smb_snap = 'test_sortie_smb.nc' !file_smb_snap = trim(dirnameout)//trim(file_smb_snap) !ncidloc = 501 !call lect_netcdf_type !write(6,*) 'ncdf_type' !if (ncdf_type.eq.32) then ! status = nf90_create(TRIM(file_smb_snap),NF90_WRITE,ncidloc) ! ouverture du fichier !else if (ncdf_type.eq.64) then ! status = nf90_create(trim(file_smb_snap),and(nf90_write,nf90_64bit_offset),ncidloc) ! r2d2 !end if ! status = nf90_close(ncidloc) ! fermeture ! call write_ncdf_dim('x',trim(file_smb_snap),nx) ! dimensions des tableaux ! call write_ncdf_dim('y',trim(file_smb_snap),ny) ! call write_ncdf_dim('time',trim(file_smb_snap),nb_snap) !Tab3D(:,:,:) = smb_snap (:,:,:) !dimtestname(1) = 'x' !dimtestname(2) = 'y' !dimtestname(3) = 'time' !call write_ncdf_var(trim('smb'),dimtestname,trim(file_smb_snap),tab3D,'double') end subroutine input_clim !-------------------------------------------------------------------------------- !> SUBROUTINE: init_forclim !! Routine qui permet d'initialiser les variables climatiques au cours du temps !> subroutine init_forclim implicit none namelist/lapse_rates/T_lapse_rate namelist/clim_pert_massb/coefT,filforc rewind(num_param) ! pour revenir au debut du fichier param_list.dat read(num_param,lapse_rates) ! formats pour les ecritures dans 42 428 format(A) rewind(num_param) ! pour revenir au debut du fichier param_list.dat read(num_param,lapse_rates) write(num_rep_42,428)'!________________________________________________________________' write(num_rep_42,428)'! module climat_Grice2sea_years_mod ' write(num_rep_42,lapse_rates) write(num_rep_42,428)'!T_lapse_rate = lapse rate temp annuelle ' write(num_rep_42,428)'!________________________________________________________________' rewind(num_param) read(num_param,clim_pert_massb) write(num_rep_42,428)'!___________________________________________________________' write(num_rep_42,428) '&clim_pert ! module climat_perturb_mod ' write(num_rep_42,*) write(num_rep_42,*) 'coefT = ', coefT write(num_rep_42,'(A,A)') ' filforc = ', filforc write(num_rep_42,*)'/' write(num_rep_42,*) ! appelle la routine de lecture des smb annuels call input_clim if (massb_time == 1) then ! lecture gradients smb call init_grad_smb endif return end subroutine init_forclim !-------------------------------------------------------------------------------- !> SUBROUTINE: forclim !! !! Routine qui permet le calcul climatique au cours du temps !! @note Au temps considere (time) attribue les scalaires !! @note - tafor : forcage en temperature !! @note - sealevel : forcage niveau des mers !! @note - coefbmelt : forcage fusion basale ice shelves !> subroutine forclim ! au temps considere (time) implicit none integer i,j,ift select case (massb_time) case(0) ! smb fixe et Tann avec lapse rate + Tafor if(time.lt.tdate(1)) then tafor=tpert(1) sealevel=spert(1) ift=1 else if (time.ge.tdate(nft)) then tafor=tpert(nft) sealevel=spert(nft) ift=nft else do i=1,nft-1 if((time.ge.tdate(i)).and.(time.lt.tdate(i+1))) then ! entre i et i+1 : cas general tafor=tpert(i)+(tpert(i+1)-tpert(i))* & (time-tdate(i))/(tdate(i+1)-tdate(i)) sealevel=spert(i)+(spert(i+1)-spert(i))* & (time-tdate(i))/(tdate(i+1)-tdate(i)) ift=i goto 100 endif end do endif 100 continue Tann (:,:) = Ta0 (:,:) + T_lapse_rate * (S(:,:)-S0(:,:)) +Tafor Ts(:,:) = Tann(:,:) case(1) call massb_Ice2sea_RCM case default print *, 'Numero de massb invalide dans climat_Grice2sea_years_mod' stop end select end subroutine forclim subroutine massb_Ice2sea_RCM ! calcule le mass balance implicit none integer :: k_snap ! pour calculer les indices de temps real :: time_bis ! pour repliquer les dernieres annees integer :: k ! calcul smb a partir fichier snapshots massb_Ice2sea_RCM ! Calcule le mass balance d'apres un fichier de snapshots ! avec la temperature parametree ! surface temperature et accumulation Tann (:,:) = Ta0 (:,:) + T_lapse_rate * (S(:,:)-S0(:,:)) Ts(:,:) = Tann(:,:) ! calcule bm_time par interpolation entre deux snapshots ! avant prend la valeur de reference ! apres prend la derniere valeur ! en general les snapshots vont de 2000 a 2200 if(time.lt.time_snap(1)) then ! time avant le forcage bm_time(:,:) = bm_0(:,:) k_snap = 1 S_ref(:,:) = S(:,:) ! du coup sera la surface de reference avant le forcage icum = 0 i_moy = 0 bm_ref_t(:,:)= bm_0(:,:) ! bilan de masse de reference time_prec = time else if (time.ge.time_snap(nb_snap)) then ! time apres le forcage k_snap = nb_snap bm_time(:,:) = smb_snap (:,:,k_snap) if (abs(time-time_prec-1.).lt.dt) then ! time_prec = time_prec + 1 icum = 1 else icum = 0 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 bm_time(:,:) = smb_snap(:,:,k)+(smb_snap(:,:,k+1)-smb_snap(:,:,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 icum = 1 time_prec = time_snap(k) ! time_prec est le temps du precedent ! snapshot garde else icum = 0 endif exit endif end do endif call grad_smb !-----------------------------> A faire if (massb_time == 1) then ! pas d'interpolation verticale bm(:,:) = bm_time(:,:) else if (massb_time == 2) then ! interpolation verticale ! ajuste bm en fonction du temps et du gradient bm(:,:) = bm_time(:,:) + grad_bm(:,:) *(S(:,:) - S_ref(:,:)) write(6,897) time, time_prec, icum, i_moy 897 format('test temps smb ',2(f0.3,1x),2(i0,1x)) endif ! garde les 10 dernieres annees et calcule la moyenne if (icum.eq.1) then ! stockage dans le tableau bm_ref_10 do k = 9,1,-1 bm_ref_10(:,:,k+1) = bm_ref_10(:,:,k) ! on decale tous les elements end do bm_ref_10(:,:,k+1) = bm(:,:) ! le plus recent est en position 1 i_moy = i_moy +1 ! compte combien il y en a pour la moyenne i_moy = min(i_moy,10) bm_ref_t(:,:) = 0. do k = 1,i_moy bm_ref_t(:,:) = bm_ref_t(:,:) + bm_ref_10(:,:,k) end do bm_ref_t(:,:) = bm_ref_t(:,:)/i_moy write(6,898) time, time_prec, icum, i_moy 898 format('cumul pour gradient ',2(f0.3,1x),2(i0,1x)) end if end subroutine massb_Ice2sea_RCM !------------------------------------------------------------------------------ !> initialise le calcul du gradient vertical de smb subroutine init_grad_smb use module3D_phy implicit none character(len=120) :: file_grad_smb ! nom du fichier gradients de Tamsin character(len=40) :: fin_ligne ! fin de la ligne real :: grad namelist/grad_smb/file_grad_smb ! Dans lequel sont : ! grad_N_smb_pos,grad_N_smb_neg,grad_S_smb_pos,grad_S_smb_neg,lim_lat rewind(num_param) ! pour revenir au debut du fichier param_list.dat read(num_param,grad_smb) ! formats pour les ecritures dans 42 428 format(A) rewind(num_param) ! pour revenir au debut du fichier param_list.dat read(num_param,grad_smb) write(num_rep_42,428)'!________________________________________________________________' write(num_rep_42,428)'! gradient smb climat_Grice2sea_years_mod ' write(num_rep_42,grad_smb) write(num_rep_42,428)'!grad_smb = fichier gradient SMB ' write(num_rep_42,428)'!________________________________________________________________' file_grad_smb=trim(dirnameinp)//'SMB-H-Feedback/gradients_18_07_2012/'//trim(file_grad_smb) open(622,file=file_grad_smb) do i=1,4 read(622,'(f9.4,A)') grad,fin_ligne write(6,*) grad,fin_ligne if (index(fin_ligne, "North").ne.0) then ! North est dans la ligne if (index(fin_ligne, "<").ne.0) then ! smb negatif grad_N_smb_neg = grad else if (index(fin_ligne, ">=").ne.0) then ! smb positif grad_N_smb_pos = grad else write(6,*) 'probleme lecture North fichier smb', file_grad_smb STOP end if else if (index(fin_ligne, "South").ne.0) then !South est dans la ligne if (index(fin_ligne, "<").ne.0) then ! smb negatif grad_S_smb_neg = grad else if (index(fin_ligne, ">=").ne.0) then ! smb positif grad_S_smb_pos = grad else write(6,*) 'probleme lecture South fichier smb ', file_grad_smb STOP end if else write(6,*) 'probleme lecture ni North ni South fichier smb ', file_grad_smb write(6,*) 'fin_ligne',fin_ligne,' index North', index(fin_ligne, "North"),' index South', index(fin_ligne, "South") STOP end if end do write(6,*) 'coefficients lus ' write(6,*) 'grad_N_smb_pos', grad_N_smb_pos write(6,*) 'grad_N_smb_neg', grad_N_smb_neg write(6,*) 'grad_S_smb_pos', grad_S_smb_pos write(6,*) 'grad_S_smb_neg', grad_S_smb_neg grad_N_smb_pos = coef_smb_unit * grad_N_smb_pos grad_N_smb_neg = coef_smb_unit * grad_N_smb_neg grad_S_smb_pos = coef_smb_unit * grad_S_smb_pos grad_S_smb_neg = coef_smb_unit * grad_S_smb_neg write(6,*) 'coefficients *coef_smb_unit ' write(6,*) 'grad_N_smb_pos', grad_N_smb_pos write(6,*) 'grad_N_smb_neg', grad_N_smb_neg write(6,*) 'grad_S_smb_pos', grad_S_smb_pos write(6,*) 'grad_S_smb_neg', grad_S_smb_neg end subroutine init_grad_smb !------------------------------------------------------------------------------ !> Calcule le gradient vertical de smb subroutine grad_smb use module3D_phy implicit none do j = 1,ny do i = 1,nx if (Ylat(i,j).gt.77.) then ! region nord if (bm_ref_t(i,j).lt. 0.) then ! smb negatif grad_bm (i,j) = grad_N_smb_neg else if (bm_ref_t(i,j).ge. 0.) then ! smb positif grad_bm (i,j) = grad_N_smb_pos end if else if (Ylat(i,j).le.77.) then ! region sud if (bm_ref_t(i,j).lt. 0.) then ! smb negatif grad_bm (i,j) = grad_S_smb_neg else if (bm_ref_t(i,j).ge. 0.) then ! smb positif grad_bm (i,j) = grad_S_smb_pos end if end if end do end do end subroutine grad_smb end module climat_Grice2sea_years_perturb_mod