Changeset 366


Ignore:
Timestamp:
02/22/23 16:28:12 (14 months ago)
Author:
dumas
Message:

ablation_mod.f90 : SMB can be computed using ITM with pdd_type=3

Location:
trunk/SOURCES
Files:
5 edited

Legend:

Unmodified
Added
Removed

  • trunk/SOURCES/3D-physique-gen_mod.f90

    r333 r366  
    301301  real,dimension(nx,ny) :: PRECIP       !< precipitation 
    302302  real,dimension(nx,ny) :: PRECIP0      !< initial precipitation  (used in 'heminord') 
     303  real,dimension(nx,ny,12) :: prmois      !< monthly total precipitations 
     304  real,dimension(nx,ny,12) :: sfmois      !< monthly snowfall 
     305  real,dimension(nx,ny,12) :: rfmois      !< monthly rainfall 
    303306  real,dimension(nx,ny) :: PHID         !< flux de chaleur lie a la deformation et glissement basal 
    304307  real,dimension(nx,ny) :: chalgliss_maj !< chaleur de glissement seulement 
    … …  
    320323  real,dimension(nx,ny) :: SLOPE2mx     !< = Sdx**2 + Sdymx**2 '>' 
    321324  real,dimension(nx,ny) :: SLOPE2my     !< = Sdy**2 + Sdxmy**2 '^' 
     325  real,dimension(nx,ny,12) :: swmois    !< monthly SW radiations at TOA 
    322326  double precision,dimension(nx,ny) :: S0           !< altitude actuelle de la surface 
    323327!afq -- replaced by sealevel_2D:  real,dimension(nx,ny) :: slv          !< niveau de flottaison (sealevel et lakes) 

  • trunk/SOURCES/Grismip6_files/module_choix-grismip6.f90

    r262 r366  
    4444!use climat_forcage_mois_mod ! forcage mensuel GCM 1 Snapshot Fev 2015 
    4545!use climat_perturb_mod   ! climat perturbe a reverifier Dec 2015 
    46 use climat_Grice2sea_years_mod  ! climat force par fichier SMB directement (grice2sea) 
     46!use climat_Grice2sea_years_mod  ! climat force par fichier SMB directement (grice2sea) 
    4747!use climat_Grice2sea_years_perturb_mod ! climat force par fichier SMB directement (grice2sea) + index temperature carotte de glace 
    4848!use climat_InitMIP_years_perturb_mod ! climat pour experiences initMIP 
    4949 
    50 !use climat_forcage_mod 
     50use climat_forcage_mod 
    5151!use climat_synthes_mod 
    5252!use climat_profil_mod 
    5353!use climat_regions_delta 
    5454 
    55 !use ablation_mod ! calcul de l'ablation (PDD ou autre methode) 
    56 use no_ablation  ! pas de calcul de l'ablation => lecture fichier SMB (necessaire avec climat_Grice2sea_years_mod) 
     55use ablation_mod ! calcul de l'ablation (PDD ou autre methode) 
     56!use no_ablation  ! pas de calcul de l'ablation => lecture fichier SMB (necessaire avec climat_Grice2sea_years_mod) 
    5757 
    5858! pas de differences locales de niveau marin 
    … …  
    6565!--------------Choix isostasie---------------------- 
    6666! use isostasie_mod  ! module permettant de calculer la deflexion isostasique 
    67 use noisostasie_mod ! module pour ne pas avoir d'isostasie 
     67use isostasie_mod ! module pour ne pas avoir d'isostasie 
    6868 
    6969 

  • trunk/SOURCES/ablation_mod.f90

    r330 r366  
    66! pdd_type=1  ! pdd fausto 
    77! pdd_type=2  ! pdd Tarasov 
     8! pdd_type=3  ! ITM 
    89module ablation_mod 
    910 
    … …  
    3334  REAL, dimension(nx,ny) :: PDDCT2               ! ct for PDD calculation 
    3435  integer :: methode_abl=0                       ! selection methode pdd (0) ou van den Berg 2008 (1) 
     36  ! ITM variables 
     37  REAL :: c ! Short-wave radiation and sensible heat flux constant 
     38  REAL :: lambda ! Long-wave radiation coefficient 
     39 
     40  REAL, DIMENSION(nx,ny) :: c_2D             
     41  REAL, DIMENSION(nx,ny) :: lambda_2D   
    3542 
    3643contains 
    … …  
    5966  write(num_rep_42,'(A,F12.6)')    'sigma_ice = ',sigma_ice 
    6067  write(num_rep_42,*)'/'                       
    61   write(num_rep_42,428) '! pdd_type : 0 reeh, 1 Fausto, 2 Tarasov' 
     68  write(num_rep_42,428) '! pdd_type : 0 reeh, 1 Fausto, 2 Tarasov, 3 ITM' 
    6269  write(num_rep_42,428) '! annual : T = annuel, F = mensuel' 
    6370  write(num_rep_42,428)'! Cice and Csnow, melting factors for ice and snow ' 
    … …  
    7481  PDDCT=DTP/sigma_ice/SQRT(2.*PI_L)/NYEAR*365. 
    7582  PDDCT2=DTP/sigma_ice/SQRT(2.*PI_L)/NYEAR2*365. 
     83   
     84  ! namelist specifique ITM : 
     85  namelist/ablation_ITM/c,lambda,csi 
     86 
     87  if (pdd_type.EQ.3) then ! ITM => lecture de la namelist ablation_ITM 
     88     rewind(num_param)        ! pour revenir au debut du fichier param_list.dat 
     89     read(num_param,ablation_ITM) 
     90     write(num_rep_42,428)'!___________________________________________________________________'  
     91     write(num_rep_42,428) '&ablation_ITM                             ! module ablation_ITM_mod' 
     92     write(num_rep_42,'(A,F12.6)')    'c = ',c 
     93     write(num_rep_42,'(A,F12.6)')    'lambda = ',lambda 
     94     write(num_rep_42,*)'/' 
     95     write(num_rep_42,428) '! c is the short-wave radiation and sensible heat flux constant.' 
     96     write(num_rep_42,428)'! lambda is the long-wave radiation coefficient.' 
     97     write(num_rep_42,428)'! csi proportion of melted water that can refreeze ' 
     98     write(num_rep_42,*) 
     99 
     100     ! initialisation des variables 2D : 
     101     c_2D(:,:)=c 
     102     lambda_2D(:,:)=lambda 
     103  endif 
     104      
     105   
    76106end subroutine init_ablation 
    77107 
    … …  
    84114 
    85115 
    86   USE module3d_phy,only:Tjuly,Tann,Tmois,acc,pdd,TS,Tshelf,precip,BM,Abl,S,dice,cl 
     116  USE module3d_phy,only:Tjuly,Tann,Tmois,acc,pdd,TS,Tshelf,precip,BM,Abl,S,dice,cl, & 
     117       swmois,ROFRESH,H,prmois 
    87118 
    88119  IMPLICIT NONE 
    … …  
    90121  real, dimension(nx,ny) :: pds, simax, pdsi, sif 
    91122  integer :: i,j,k,mo,nday 
    92   real :: summ 
     123  real, dimension(nx,ny) :: summ 
    93124  real :: temp 
    94125  real,dimension(365) :: TT             !< air temperature yearly cycle, for PDD 
    … …  
    96127  REAL, DIMENSION(nx,ny) :: pr_ice_eq, totmelt, snowmelt, cpsurf 
    97128  REAL, DIMENSION(nx,ny) :: refr2, refreezed_ice 
    98  
    99  
    100 ! pour pdd fausto calcul Cice_2D, csnow_2D, csi_2D,sigma_ice_2D 
     129  ! ITM variables : 
     130  REAL :: alphas=0.85 
     131  REAL :: coef_reduc=0.3 
     132  REAL :: epsilon=10.**-6. 
     133  REAL, DIMENSION(nx,ny,12) :: alb 
     134  REAL, DIMENSION(nx,ny,12) :: ft, snow, rain 
     135  REAL, DIMENSION(nx,ny) :: alphag 
     136 
     137   
     138  ! pour pdd fausto calcul Cice_2D, csnow_2D, csi_2D,sigma_ice_2D 
    101139  IF (pdd_type.eq.1) THEN 
    102140      WHERE (Tjuly(:,:).GE.10.) 
    … …  
    110148          Csnow_2D(:,:)=3. 
    111149      ENDWHERE 
    112 ! pour etre en m/an : 
     150      ! pour etre en m/an : 
    113151      Cice_2D(:,:)=Cice_2D(:,:)/1000. 
    114152      Csnow_2D(:,:)=Csnow_2D(:,:)/1000. 
    115 !cdc version std   
     153      ! version std   
    116154      sigma_ice_2D(:,:)=1.574+0.0012224*S(:,:) 
    117 ! calcul de S22, PDDCT et PDDCT2 
     155      ! calcul de S22, PDDCT et PDDCT2 
    118156      S22(:,:)=0.5/SIGMA_ICE_2D(:,:)/SIGMA_ICE_2D(:,:) 
    119157      PDDCT(:,:)=DTP/SIGMA_ICE_2D(:,:)/SQRT(2.*PI_L)/NYEAR*365. 
    120158      PDDCT2(:,:)=DTP/SIGMA_ICE_2D(:,:)/SQRT(2.*PI_L)/NYEAR2*365. 
    121159 
    122 ! calcul du tx de regel 
     160      ! calcul du tx de regel 
    123161      WHERE (S(:,:).LE.800) 
    124162          CSI_2D(:,:)=0. 
    … …  
    145183 
    146184      sigma_ice_2D(:,:)=sigma_ice 
    147   ENDIF 
    148  
    149  
    150  
    151   IF (annual) THEN 
     185       
     186   ELSEIF (pdd_type.EQ.3) THEN ! ITM 
     187      ! calcul fraction de neige/fration de pluie 
     188      WHERE (Tmois(:,:,:).LT.-10.) 
     189         ft(:,:,:) = 1. 
     190      ELSEWHERE (Tmois(:,:,:).GT.7.) 
     191         ft(:,:,:) = 0. 
     192      ELSEWHERE 
     193         ft(:,:,:) = (7.-Tmois(:,:,:))/17. 
     194      ENDWHERE 
     195      snow(:,:,:) = ft(:,:,:)*prmois(:,:,:) 
     196      rain(:,:,:) = (1.-ft(:,:,:))*prmois(:,:,:) 
     197      summ(:,:) = 0. 
     198 
     199      WHERE (H(:,:).GT.0.) 
     200         alphag(:,:) = 0.4 
     201      ELSEWHERE  
     202         alphag(:,:) = 0.2 
     203      ENDWHERE 
     204 
     205      ! calcul de l'albedo 
     206      DO mo=1,12 
     207         WHERE (snow(:,:,mo).LT.epsilon) 
     208            alb(:,:,mo) = alphag(:,:) 
     209         ELSEWHERE 
     210            alb(:,:,mo) = amax1(alphas-coef_reduc*rain(:,:,mo)/snow(:,:,mo)*(alphas-alphag(:,:)),alphag(:,:)) 
     211         ENDWHERE 
     212      ENDDO 
     213      ! calcul de la fonte totale a partir de la formule de l'ITM (van der Berg et al. 2008) en m water equivalent.  
     214      DO mo=1,12 
     215         summ(:,:) = summ(:,:) + amax1(((1.-alb(:,:,mo))*swmois(:,:,mo)+c_2D(:,:)+lambda_2D(:,:)*Tmois(:,:,mo))*(3600.*24.*30.)/(ROFRESH*CL),0.) 
     216      END DO 
     217   ENDIF 
     218 
     219 
     220   IF (pdd_type.LE.2) THEN ! PDD : reconstruction du cycle diurne et mensuel 
     221   IF (annual) THEN 
    152222!     (*** positive degree days Tann et Tjuly***) 
    153223      DO i=1,nx 
    154         DO j=1,ny 
    155           summ=0.0 
    156           month_reconstr: DO nday=1,nyear   ! reconstitution du cycle annuel 
    157             tt(nday)=tann(i,j)+(tjuly(i,j)-tann(i,j))*COS(pyg*nday)  
    158             temp=0.0 
    159             k=1 
    160             average_day_reconstr: DO WHILE ((temp.LE.tt(nday)+2.5*sigma_ice_2D(i,j)).AND.k.LE.50) ! pdd d'un jour par mois     
    161               summ=summ+temp*EXP(-(temp-tt(nday))*(temp-tt(nday))*s22(i,j))  
    162               temp=temp+dtp 
    163               k=k+1 
    164             END DO average_day_reconstr 
    165           END DO month_reconstr 
    166           pdd(i,j)=summ*pddct(i,j)  
    167         END DO 
     224         DO j=1,ny 
     225            summ(i,j)=0.0 
     226            month_reconstr: DO nday=1,nyear   ! reconstitution du cycle annuel 
     227               tt(nday)=tann(i,j)+(tjuly(i,j)-tann(i,j))*COS(pyg*nday)  
     228               temp=0.0 
     229               k=1 
     230               average_day_reconstr: DO WHILE ((temp.LE.tt(nday)+2.5*sigma_ice_2D(i,j)).AND.k.LE.50) ! pdd d'un jour par mois     
     231                  summ(i,j)=summ(i,j)+temp*EXP(-(temp-tt(nday))*(temp-tt(nday))*s22(i,j))  
     232                  temp=temp+dtp 
     233                  k=k+1 
     234               END DO average_day_reconstr 
     235            END DO month_reconstr 
     236            pdd(i,j)=summ(i,j)*pddct(i,j)  
     237         END DO 
    168238      END DO 
    169   ELSE            ! pdd mensuel 
     239   ELSE            ! pdd mensuel 
    170240      DO i=1,nx 
    171         DO j=1,ny 
    172           summ=0.0 
     241         DO j=1,ny 
     242            summ(i,j)=0.0 
    173243! On a deja calcule Tmois(i,j,m) : temperature moyenne de chaque mois 
    174           month: DO mo=1,12 ! boucle sur les mois 
    175             temp=0.0          ! variable d'integration 
    176             k=1 
    177             average_day: DO WHILE ((temp.LE.Tmois(i,j,mo)+2.5*sigma_ice_2D(i,j)).AND.k.LE.50) ! pdd d'un jour par mois     
    178               summ=summ+temp*EXP(-((temp-Tmois(i,j,mo))*(temp-Tmois(i,j,mo)))*s22(i,j)) 
    179               temp=temp+dtp     ! pdtp pas d'integration 
    180               k= k+1 
    181             END DO average_day 
    182           END DO month          !  boucle sur le mois 
    183           pdd(i,j)=summ*pddct(i,j)   ! pdd pour toute l'annee 
    184         END DO 
     244            month: DO mo=1,12 ! boucle sur les mois 
     245               temp=0.0          ! variable d'integration 
     246               k=1 
     247               average_day: DO WHILE ((temp.LE.Tmois(i,j,mo)+2.5*sigma_ice_2D(i,j)).AND.k.LE.50) ! pdd d'un jour par mois     
     248                  summ(i,j)=summ(i,j)+temp*EXP(-((temp-Tmois(i,j,mo))*(temp-Tmois(i,j,mo)))*s22(i,j)) 
     249                  temp=temp+dtp     ! pdtp pas d'integration 
     250                  k= k+1 
     251               END DO average_day 
     252            END DO month          !  boucle sur le mois 
     253            pdd(i,j)=summ(i,j)*pddct(i,j)   ! pdd pour toute l'annee 
     254         END DO 
    185255      END DO 
    186   ENDIF 
     256   ENDIF 
     257   ENDIF 
    187258 
    188259 
    189260!     calcul du Bilan de masse 
    190   IF (methode_abl.EQ.0) THEN 
    191      IF (pdd_type.EQ.1) THEN 
    192         PDS(:,:)=ACC(:,:)/Csnow_2D(:,:) 
    193         SIMAX(:,:)=ACC(:,:)*CSI_2D(:,:) 
    194         PDSI(:,:)=SIMAX(:,:)/Cice_2D(:,:) 
     261   IF (methode_abl.EQ.0) THEN 
     262      IF (pdd_type.EQ.1) THEN 
     263         PDS(:,:)=ACC(:,:)/Csnow_2D(:,:) 
     264         SIMAX(:,:)=ACC(:,:)*CSI_2D(:,:) 
     265         PDSI(:,:)=SIMAX(:,:)/Cice_2D(:,:) 
    195266! avec regel de 60% puis fonte (2 premiers where) : 
    196267!      WHERE (PDD(:,:).LE.CSI_2D*PDS(:,:)) 
    … …  
    202273!          SIF(:,:)=SIMAX(:,:) 
    203274!      endwhere 
    204         WHERE (PDD(:,:).LE.PDS(:,:))  ! test avec regel de 60% progressif 
    205            BM(:,:)=ACC(:,:)-PDD(:,:)*Csnow_2D*(1-CSI_2D(:,:)) 
    206            SIF(:,:)=PDD(:,:)*Csnow_2D(:,:)*(1-CSI_2D(:,:)) 
    207         endwhere 
    208         WHERE ((PDS(:,:).LT.PDD(:,:)).AND.(PDD(:,:).LE.PDS(:,:)+PDSI(:,:))) 
    209            BM(:,:)=SIMAX(:,:)-(PDD(:,:)-PDS(:,:))*Cice_2D 
    210            SIF(:,:)=SIMAX(:,:) 
    211         endwhere 
    212         WHERE (PDS(:,:)+PDSI(:,:).LE.PDD(:,:)) 
    213            BM(:,:)=(PDS(:,:)+PDSI(:,:)-PDD(:,:))*Cice_2D 
    214            SIF(:,:)=SIMAX(:,:) 
    215         endwhere 
    216      ELSEIF (PDD_type.EQ.2) THEN  ! PDD Tarasov 
    217         PDS(:,:)=ACC(:,:)/Csnow_2D(:,:) 
    218         pr_ice_eq(:,:) = amax1(0.,((PRECIP(:,:)/DICE)-ACC(:,:)))  ! precipe liquide (ice equivalent) 
    219  
    220         WHERE (PDD(:,:).LE.PDS(:,:)) 
    221            totmelt(:,:) = Csnow_2D(:,:)*PDD(:,:) 
    222         ELSEWHERE 
    223            totmelt(:,:) = Csnow_2D(:,:)*PDS(:,:) + Cice_2D(:,:)*(PDD(:,:)-PDS(:,:)) 
    224         ENDWHERE 
    225  
    226         snowmelt(:,:) = amin1(totmelt(:,:),ACC(:,:)) 
     275         WHERE (PDD(:,:).LE.PDS(:,:))  ! test avec regel de 60% progressif 
     276            BM(:,:)=ACC(:,:)-PDD(:,:)*Csnow_2D*(1-CSI_2D(:,:)) 
     277            SIF(:,:)=PDD(:,:)*Csnow_2D(:,:)*(1-CSI_2D(:,:)) 
     278         endwhere 
     279         WHERE ((PDS(:,:).LT.PDD(:,:)).AND.(PDD(:,:).LE.PDS(:,:)+PDSI(:,:))) 
     280            BM(:,:)=SIMAX(:,:)-(PDD(:,:)-PDS(:,:))*Cice_2D 
     281            SIF(:,:)=SIMAX(:,:) 
     282         endwhere 
     283         WHERE (PDS(:,:)+PDSI(:,:).LE.PDD(:,:)) 
     284            BM(:,:)=(PDS(:,:)+PDSI(:,:)-PDD(:,:))*Cice_2D 
     285            SIF(:,:)=SIMAX(:,:) 
     286         endwhere 
     287      ELSEIF (PDD_type.EQ.2) THEN  ! PDD Tarasov 
     288         PDS(:,:)=ACC(:,:)/Csnow_2D(:,:) 
     289         pr_ice_eq(:,:) = amax1(0.,((PRECIP(:,:)/DICE)-ACC(:,:)))  ! precipe liquide (ice equivalent) 
     290 
     291         WHERE (PDD(:,:).LE.PDS(:,:)) 
     292            totmelt(:,:) = Csnow_2D(:,:)*PDD(:,:) 
     293         ELSEWHERE 
     294            totmelt(:,:) = Csnow_2D(:,:)*PDS(:,:) + Cice_2D(:,:)*(PDD(:,:)-PDS(:,:)) 
     295         ENDWHERE 
     296 
     297         snowmelt(:,:) = amin1(totmelt(:,:),ACC(:,:)) 
    227298 
    228299! Deux formules possibles pour la capacité calorifique (en J/kg.K): 
    229300! Formule 1 correspond à celle figurant dans prop-thermiques_mod.f90 
    230301! Formule 2 : celle de Tarasov 
    231         cpsurf(:,:) = 2115.3+7.79293*TANN(:,:) ! Formule Catherine 
    232 !      cpsurf(:,:) = 152.5 + 7.122*TANN(:,:) ! Formule Tarasov 
     302         cpsurf(:,:) = 2115.3+7.79293*TANN(:,:) ! Formule Catherine 
     303!        cpsurf(:,:) = 152.5 + 7.122*TANN(:,:) ! Formule Tarasov 
    233304 
    234305!        where(snowmelt(:,:).lt.ACC(:,:)) 
    235         refr2(:,:) = 2.2*(ACC(:,:)-snowmelt(:,:))-(cpsurf(:,:)/CL)*amin1(TANN(:,:),0.) 
    236         refreezed_ice(:,:) = amin1(pr_ice_eq(:,:)+snowmelt(:,:),refr2(:,:)) 
     306         refr2(:,:) = 2.2*(ACC(:,:)-snowmelt(:,:))-(cpsurf(:,:)/CL)*amin1(TANN(:,:),0.) 
     307         refreezed_ice(:,:) = amin1(pr_ice_eq(:,:)+snowmelt(:,:),refr2(:,:)) 
    237308!        elsewhere 
    238309!           refr2(:,:) = -(cpsurf(:,:)/CL)*amin1(TANN(:,:),0.) 
    … …  
    240311!        endwhere 
    241312 
    242         bm(:,:) = ACC(:,:)+refreezed_ice(:,:)-totmelt(:,:) 
    243         SIF(:,:)=refreezed_ice(:,:) 
    244  
    245      ELSE ! pdd standard reeh 
     313         bm(:,:) = ACC(:,:)+refreezed_ice(:,:)-totmelt(:,:) 
     314         SIF(:,:)=refreezed_ice(:,:) 
     315 
     316      ELSEIF (PDD_type.EQ.0) THEN  ! PDD standard reeh 
    246317!       (* Positive degrees required to melt the snow layer *) 
    247         PDS(:,:)=ACC(:,:)/Csnow_2D(:,:) 
     318         PDS(:,:)=ACC(:,:)/Csnow_2D(:,:) 
    248319!       (* Maximum amount of super. ice that can be formed *) 
    249         SIMAX(:,:)=ACC(:,:)*CSI_2D(:,:) 
     320         SIMAX(:,:)=ACC(:,:)*CSI_2D(:,:) 
    250321!       (* Pos. degrees required to melt the superimposed ice *) 
    251         PDSI(:,:)=SIMAX(:,:)/Cice_2D(:,:) 
     322         PDSI(:,:)=SIMAX(:,:)/Cice_2D(:,:) 
    252323! avec regel de 60% puis fonte (2 premiers where) : 
    253         WHERE (PDD(:,:).LE.CSI_2D*PDS(:,:)) 
    254            BM(:,:)=ACC(:,:) 
    255            SIF(:,:)=PDD(:,:)*Csnow_2D(:,:) 
    256         endwhere 
    257         WHERE ((CSI_2D*PDS(:,:).LT.PDD(:,:)).AND.(PDD(:,:).LE.PDS(:,:))) 
    258            BM(:,:)=ACC(:,:)+SIMAX(:,:)-PDD(:,:)*Csnow_2D 
    259            SIF(:,:)=SIMAX(:,:) 
    260         endwhere 
    261 !       WHERE (PDD(:,:).LE.PDS(:,:))                  ! test avec regel de 60% progressif| 
     324         WHERE (PDD(:,:).LE.CSI_2D*PDS(:,:)) 
     325            BM(:,:)=ACC(:,:) 
     326            SIF(:,:)=PDD(:,:)*Csnow_2D(:,:) 
     327         endwhere 
     328         WHERE ((CSI_2D*PDS(:,:).LT.PDD(:,:)).AND.(PDD(:,:).LE.PDS(:,:))) 
     329            BM(:,:)=ACC(:,:)+SIMAX(:,:)-PDD(:,:)*Csnow_2D 
     330            SIF(:,:)=SIMAX(:,:) 
     331         endwhere 
     332!        WHERE (PDD(:,:).LE.PDS(:,:))                  ! test avec regel de 60% progressif| 
    262333!           BM(:,:)=ACC(:,:)-PDD(:,:)*Csnow_2D*(1-CSI_2D)   ! remplace les 2 premiers where    | 
    263334!           SIF(:,:)=PDD(:,:)*Csnow_2D(:,:)*(1-CSI_2D(:,:)) !----------------------------------| 
    264 !      endwhere 
    265         WHERE ((PDS(:,:).LT.PDD(:,:)).AND.(PDD(:,:).LE.PDS(:,:)+PDSI(:,:))) 
    266            BM(:,:)=SIMAX(:,:)-(PDD(:,:)-PDS(:,:))*Cice_2D 
    267            SIF(:,:)=SIMAX(:,:) 
    268         endwhere 
    269         WHERE (PDS(:,:)+PDSI(:,:).LE.PDD(:,:)) 
    270            BM(:,:)=(PDS(:,:)+PDSI(:,:)-PDD(:,:))*Cice_2D 
    271            SIF(:,:)=SIMAX(:,:) 
    272         endwhere 
    273      ENDIF 
    274   ELSEIF ( methode_abl.EQ.1 ) THEN ! Insolation Temperature Melt equation, van den Berg, 2008 
    275      SIMAX(:,:)=ACC(:,:)*CSI 
    276      SIF(:,:)=SIMAX(:,:) 
    277   ELSE 
    278      print*,'ablation.f90 : pb methode calcul BM' 
    279      print*,'ATTENTION methode_abl > 1 NON COMPATIBLE AVEC iLOVECLIM'  
    280      stop 
    281   ENDIF 
    282  
    283 ! calcul de la temperature de surface (utilisee dans icetemp) : 
    284   TS(:,:)=(TANN(:,:)+26.6*SIF(:,:)) 
    285   TS(:,:)=min(0.0,TS(:,:)) 
    286   tshelf(:,:)=TS(:,:) 
    287   Abl(:,:)=BM(:,:)-Acc(:,:) 
    288  
    289 END SUBROUTINE ABLATION 
     335!        ENDWHERE 
     336         WHERE ((PDS(:,:).LT.PDD(:,:)).AND.(PDD(:,:).LE.PDS(:,:)+PDSI(:,:))) 
     337            BM(:,:)=SIMAX(:,:)-(PDD(:,:)-PDS(:,:))*Cice_2D 
     338            SIF(:,:)=SIMAX(:,:) 
     339         ENDWHERE 
     340         WHERE (PDS(:,:)+PDSI(:,:).LE.PDD(:,:)) 
     341            BM(:,:)=(PDS(:,:)+PDSI(:,:)-PDD(:,:))*Cice_2D 
     342            SIF(:,:)=SIMAX(:,:) 
     343         ENDWHERE 
     344      ELSEIF (PDD_type.EQ.3) THEN  ! ITM (van der Berg et al. 2008) en m water equivalent. 
     345         totmelt(:,:) = summ(:,:)/DICE ! annual melt (m ice equivalent) 
     346         ! calcul du regel 
     347         ! taux de regel (Tarasov and Peltier, 2002) 
     348         pr_ice_eq(:,:) = amax1(0.,((PRECIP(:,:)/DICE)-ACC(:,:)))  ! precipe liquide (ice equivalent) 
     349         snowmelt(:,:) = amin1(totmelt(:,:),ACC(:,:)) 
     350 
     351         ! Deux formules possibles pour la capacité calorifique (en J/kg.K): 
     352         cpsurf(:,:) = 2115.3+7.79293*TANN(:,:) ! Capacité calorifique (en J/kg.K) - Formule Catherine 
     353         !  cpsurf(:,:) = 152.5 + 7.122*TANN(:,:) ! Capacité calorifique (en J/kg.K) - Formule Tarasov 
     354 
     355         refr2(:,:) = 2.2*(ACC(:,:)-snowmelt(:,:))-(cpsurf(:,:)/CL)*amin1(TANN(:,:),0.) 
     356         refreezed_ice(:,:) = amin1(pr_ice_eq(:,:)+snowmelt(:,:),refr2(:,:)) 
     357 
     358         SIMAX(:,:)=refreezed_ice(:,:) 
     359         SIF(:,:)=SIMAX(:,:) 
     360 
     361         ! calcul du bilan de masse 
     362         BM(:,:) = ACC(:,:)+SIMAX(:,:)-totmelt(:,:) 
     363      ELSE 
     364         print*,'ablation_mod.f90 : pb methode calcul BM' 
     365         print*,'ATTENTION PDD_type > 3 NON IMPLEMENTE'  
     366         stop 
     367      ENDIF 
     368   ELSEIF ( methode_abl.EQ.1 ) THEN ! Insolation Temperature Melt equation, van den Berg, 2008 
     369      SIMAX(:,:)=ACC(:,:)*CSI 
     370      SIF(:,:)=SIMAX(:,:) 
     371   ELSE 
     372      print*,'ablation.f90 : pb methode calcul BM' 
     373      print*,'ATTENTION methode_abl > 1 NON COMPATIBLE AVEC iLOVECLIM'  
     374      stop 
     375   ENDIF 
     376 
     377   ! calcul de la temperature de surface (utilisee dans icetemp) : 
     378   TS(:,:)=(TANN(:,:)+26.6*SIF(:,:)) 
     379   TS(:,:)=min(0.0,TS(:,:)) 
     380   tshelf(:,:)=TS(:,:) 
     381   Abl(:,:)=BM(:,:)-Acc(:,:) 
     382 
     383 END SUBROUTINE ablation 
    290384 
    291385end module ablation_mod 

  • trunk/SOURCES/climat_forcage_mod.f90

    r347 r366  
    66! possibilite d'utiliser autant de snapshots que l'on veut (ntr) 
    77! fichiers de forcage, Snapshot et valeurs de lapse rate definis dans fichier param 
    8   use module3d_phy,only:nx,ny,sealevel,sealevel_2d,S,Tmois,Tann,Tjuly,acc,num_forc,num_param,num_rep_42,tafor,coefbmshelf,PI,ro,time,dirforcage,dirnameinp 
     8  use module3d_phy,only:nx,ny,sealevel,sealevel_2d,S,Tmois,Tann,Tjuly,acc,swmois,prmois,num_forc,num_param,num_rep_42,tafor,coefbmshelf,PI,ro,time,dirforcage,dirnameinp 
    99  use netcdf 
    1010  use io_netcdf_grisli 
    … …  
    2323 
    2424  character(len=100) :: clim_ref_file ! climat de reference 
    25   character(len=70), dimension(5) :: forcage_file ! fichier de forcage (climat+topo) => dimension >= ntr 
     25  character(len=100), dimension(12) :: forcage_file ! fichier de forcage (climat+topo) => dimension >= ntr 
    2626!  character(len=100) :: forcage_file2 ! fichier de forcage 2 (climat+topo) 
    2727 
    … …  
    3232  real    :: r_atmvar              !< a ratio to account fast atmos variability between the 2 snapshots (in %: 0=>no fast var,1=>only fast var) 
    3333 
    34  
    35  
    3634!  character(len=200) :: filtr(ntr)   ! fichiers de forcage 
    3735 
    … …  
    4341                     ! 2 for steady state using snap 2 ... 
    4442                     ! 0 for transient 
     43 
     44  integer :: smbcomp ! 0 for PDD 
     45                     ! 1 for ITM 
    4546   
    4647! variables pour climato mensuelle :   
    4748  real,dimension(nx,ny,12) :: t2m0   !< initial monthly air temperature (climato) 
     49  real,dimension(nx,ny,12) :: sw0   !< initial monthly downward SW radiations at the surface (climato) 
    4850  real,dimension(nx,ny,12) :: pr0    !< initial monthly precip (climato) 
    4951  real,dimension(nx,ny) :: orog0     !< topo de la climato 
    … …  
    5153! climat des snapshots : 
    5254  real,dimension(:,:,:,:), allocatable :: t2m_ss  !< t2m monthly for every GCM snapshot  
     55  real,dimension(:,:,:,:), allocatable :: sw_ss  !< sw radiations monthly for every GCM snapshot 
    5356  real,dimension(:,:,:,:), allocatable :: pr_ss   !< precip monthly for every GCM snapshot 
    5457  real,dimension(:,:,:), allocatable :: orog_ss    !< topo for every GCM snapshot  
    5558   
    5659  real,dimension(:,:,:,:), allocatable :: t2m_sstopo0 !< t2m GCM snapshot on orog0 
     60  real,dimension(:,:,:,:), allocatable :: sw_sstopo0  !< sw GCM snapshot on orog0 
    5761  real,dimension(:,:,:,:), allocatable :: pr_sstopo0  !< pr GCM snapshot on orog0 
    58  
    5962 
    6063!~   real,dimension(nx,ny,12,ntr) :: t2m_ss  !< t2m monthly for every GCM snapshot  
    … …  
    141144  enddo 
    142145 
     146  if (smbcomp .EQ. 1) then 
     147     ! lecture du climat de reference : climato mensuelle 
     148     call Read_Ncdf_var('rsds',trim(DIRNAMEINP)//TRIM(clim_ref_file),tab3d) 
     149     sw0(:,:,:)=tab3d(:,:,:) 
     150 
     151     ! lecture des fichiers snapshots pour n'importe quel geoplace 
     152     do k=1,ntr 
     153        call Read_Ncdf_var('rsds',trim(DIRNAMEINP)//'Snapshots-GCM/'//trim(forcage_file(k)),tab3d) 
     154        sw_ss(:,:,:,k)=tab3d(:,:,:) 
     155     end do 
     156 
     157     ! calcul de sw sur la topo des obs : 
     158     do m=1,12 
     159        do k=1,ntr 
     160           !sw_sstopo0(:,:,m,k)=sw_ss(:,:,m,k) 
     161           sw_sstopo0(:,:,m,k)=sw_ss(:,:,m,k)+sw_ss(:,:,m,k)*0.00006*(orog0(:,:)-orog_ss(:,:,k)) ! Adapted from Robinson et al., 2011 
     162        enddo 
     163     enddo 
     164 
     165  end if 
    143166 
    144167! lecture des fichiers d'evolution pour tout geoplace 
    … …  
    214237subroutine init_forclim 
    215238 
    216   real,dimension(5) :: ttr_temp         !< date des tranches (snapshots)  
     239  real,dimension(12) :: ttr_temp         !< date des tranches (snapshots)  
    217240  integer :: err                  !< recuperation erreur 
    218241  integer :: k 
    219242   
    220   namelist/clim_forcage/clim_ref_file,ntr,ttr_temp,forcage_file,typerun,lapserate,rappact,filforc,pertbmb,coeft,coefbmb,r_atmvar 
     243  namelist/clim_forcage/clim_ref_file,ntr,ttr_temp,forcage_file,typerun,lapserate,rappact,filforc,pertbmb,coeft,coefbmb,r_atmvar,smbcomp 
    221244 
    222245  rewind(num_param)        ! pour revenir au debut du fichier param_list.dat 
    … …  
    255278   end if 
    256279  end if 
     280 
     281  if (.not.allocated(sw_ss)) THEN 
     282   allocate(sw_ss(nx,ny,12,ntr),stat=err) 
     283   if (err/=0) then 
     284      print *,"Erreur à l'allocation du tableau t2m_ss",err 
     285      stop 4 
     286   end if 
     287  end if 
    257288   
    258289  if (.not.allocated(pr_ss)) THEN 
    … …  
    274305  if (.not.allocated(t2m_sstopo0)) THEN 
    275306   allocate(t2m_sstopo0(nx,ny,12,ntr),stat=err) 
     307   if (err/=0) then 
     308      print *,"Erreur à l'allocation du tableau t2m_sstopo0",err 
     309      stop 4 
     310   end if 
     311  end if 
     312 
     313  if (.not.allocated(sw_sstopo0)) THEN 
     314   allocate(sw_sstopo0(nx,ny,12,ntr),stat=err) 
    276315   if (err/=0) then 
    277316      print *,"Erreur à l'allocation du tableau t2m_sstopo0",err 
    … …  
    333372  real :: TEMP                       !< calcul du nbr de jour < psolid 
    334373  real,dimension(nx,ny) :: deltaZs   ! diff temp par rapport a topo orog0  
    335   real,dimension(nx,ny,12) :: prmois ! precip sur topo Zs 
    336    
    337374  real (kind=kind(0.d0)) ::  timeloc 
    338375   
    … …  
    451488    prmois(:,:,m)= pr0(:,:,m)*((pr_sstopo0(:,:,m,itr)*(1-coeftime) + pr_sstopo0(:,:,m,itr+1)*coeftime)/ pr_sstopo0(:,:,m,ntr)) * exp(rappact*(deltaZs(:,:))) 
    452489    !$OMP END WORKSHARE 
    453  
    454490! Temp et precip fonction de coeftime : 
    455491  else if (itr.LT.ntr) then 
    … …  
    466502enddo  
    467503 
    468  
     504if (smbcomp .EQ. 1) then 
     505   ! correction d'altitude 
     506   do m=1,12  
     507      if ((itr.LT.ntr-1).and.(ntr.ge.3)) then 
     508         ! if the is more than 2 snapshots : time climate is mean between 2 snapshots and anomaly with snapshot ntr (ctrl)   
     509         !$OMP WORKSHARE 
     510         swmois(:,:,m)= (sw_sstopo0(:,:,m,itr)*(1-coeftime) + sw_sstopo0(:,:,m,itr+1)*coeftime) - sw_sstopo0(:,:,m,ntr) + sw0(:,:,m)+sw0(:,:,m)*0.00006*(Zs(:,:)-orog0(:,:)) 
     511         !$OMP END WORKSHARE 
     512      ! Temp et precip fonction de coeftime : 
     513      elseif (itr.LT.ntr) then 
     514         !$OMP WORKSHARE 
     515         swmois(:,:,m)= (sw_sstopo0(:,:,m,itr) - sw_sstopo0(:,:,m,itr+1))*(1-coeftime) + sw0(:,:,m)+sw0(:,:,m)*0.00006*(Zs(:,:)-orog0(:,:)) 
     516         !$OMP END WORKSHARE 
     517      else ! itr=ntr (time>tsnapmax) 
     518         !$OMP WORKSHARE 
     519         !swmois(:,:,m) = sw0(:,:,m) 
     520         swmois(:,:,m) = sw0(:,:,m)+sw0(:,:,m)*0.00006*(Zs(:,:)-orog0(:,:)) 
     521         !$OMP END WORKSHARE 
     522      endif 
     523   enddo 
     524end if 
    469525 
    470526  !$OMP WORKSHARE 
    … …  
    472528  Tjuly=Tmois(:,:,7) 
    473529   
    474  
    475  
    476530 
    477531! calcul de l'accumulation de neige : 

  • trunk/SOURCES/climat_transient_GCM_mod.f90

    r353 r366  
    77! but not the one included in the climate anomalies forcing fields 
    88   
    9   use module3d_phy,only:nx,ny,sealevel_2d,S,S0,Tmois,Tann,Tjuly,acc,num_forc,num_param,num_rep_42,ro,coefbmshelf,dt,time,dirforcage,dirnameinp 
     9  use module3d_phy,only:nx,ny,sealevel_2d,S,S0,Tmois,Tann,Tjuly,acc,prmois,num_forc,num_param,num_rep_42,ro,coefbmshelf,dt,time,dirforcage,dirnameinp 
    1010  use netcdf 
    1111  use io_netcdf_grisli 
    … …  
    1717  real,dimension(:),allocatable       :: time_snap    !< snapshots dates   indexes: nb_snap 
    1818 
    19   real,dimension(nx,ny,12)    :: t2m0         !< reference period temperature 
    20   real,dimension(nx,ny,12)    :: pr0          !< reference period temperature 
    21   real,dimension(nx,ny,12) :: prmois ! precip sur topo Zs 
    22   real,dimension(nx,ny) :: orog0     !< topo de la climato 
    23    
    24   real :: coef_prec_unit     !< conversion factor for the precip 
     19  real,dimension(nx,ny,12)    :: t2m0            !< reference period temperature 
     20  real,dimension(nx,ny,12)    :: pr0             !< reference period temperature 
     21  real,dimension(nx,ny)       :: orog0           !< topo de la climato 
     22   
     23  real :: coef_prec_unit                         !< conversion factor for the precip 
    2524  real :: psolid=2.                              !< temp limit between liquid and solid precip 
    2625  real :: lapserate                              !< vertical temperature gradient 
Note: See TracChangeset for help on using the changeset viewer.