Changeset 389


Ignore:
Timestamp:
03/23/23 09:18:11 (13 months ago)
Author:
dumas
Message:

use only in module equat_adv_diff_2D_vect

File:
1 edited

Legend:

Unmodified
Added
Removed

  • branches/GRISLIv3/SOURCES/conserv-mass-adv-diff_sept2009_mod.f90

    r285 r389  
    1818!< 
    1919module equat_adv_diff_2D_vect                          ! Cat nouvelle mouture juin 2009 
    20 use module3D_phy 
    21 use reso_adv_diff_2D_vect 
    22 use prescribe_H 
    23  
    24 implicit none 
    25  
    26 real, dimension(nx,ny) :: advmx           !< partie advective en x 
    27 real, dimension(nx,ny) :: advmy           !< partie advective en y 
    28 real, dimension(nx,ny) :: dmx             !< partie advective en x 
    29 real, dimension(nx,ny) :: dmy             !< partie advective en y 
    30 double precision, dimension(nx,ny) :: VieuxH          !< ancienne valeur de H 
    31 real, dimension(nx,ny) :: bilmass         !< bilan de masse pour la colonne 
    32 logical, dimension(nx,ny) :: zonemx       !< pour separer advection-diffusion 
    33 logical,  dimension(nx,ny) :: zonemy      !< pour separer advection-diffusion 
    34 real :: adv_frac                          !< fraction du flux traitee en advection 
    35 !real, parameter :: Hmin=1.001             !< Hmin pour être considere comme point ice  
    36 integer :: itesti 
    37 integer :: itour 
     20  use module3D_phy, only:nx,ny,V_limit,num_param,num_rep_42,dx1,dx,mk0,i_Hp,Hp,H,uxbar,uybar,testdiag,& 
     21       dtmax,dt,dtmin,time,dtt,isynchro,diffmx,diffmy,sdx,sdy,hmx,hmy,flgzmx,flgzmy,flot,tabtest,timemax,& 
     22       marine,dtdx2,dtdx,geoplace,bm,bmelt,igrdline,ibmelt_inv,ice,ablbord,hdot 
     23  use runparam, only:itracebug,num_tracebug,tgrounded 
     24  use reso_adv_diff_2D_vect, only:init_reso_adv_diff_2D,resol_adv_diff_2D_vect 
     25  use prescribe_H, only:prescribe_fixed_points,prescribe_paleo_gl_shelf,prescribe_present_H_gl,& 
     26       prescribe_present_H_gl_bmelt,init_prescribe_H 
     27  use toy_retreat_mod, only:time_step_recul 
     28 
     29  implicit none 
     30 
     31  real, dimension(nx,ny) :: advmx           !< partie advective en x 
     32  real, dimension(nx,ny) :: advmy           !< partie advective en y 
     33  real, dimension(nx,ny) :: dmx             !< partie advective en x 
     34  real, dimension(nx,ny) :: dmy             !< partie advective en y 
     35  double precision, dimension(nx,ny) :: VieuxH          !< ancienne valeur de H 
     36  real, dimension(nx,ny) :: bilmass         !< bilan de masse pour la colonne 
     37  logical, dimension(nx,ny) :: zonemx       !< pour separer advection-diffusion 
     38  logical,  dimension(nx,ny) :: zonemy      !< pour separer advection-diffusion 
     39  real :: adv_frac                          !< fraction du flux traitee en advection 
     40  !real, parameter :: Hmin=1.001             !< Hmin pour être considere comme point ice  
     41  integer :: itesti 
     42  integer :: itour 
    3843 
    3944contains 
    4045 
    41 !> SUBROUTINE: init_icethick 
    42 !! Definition et initialisation des parametres qui gerent la repartition advection diffusion 
    43 !! @return adv_frac 
    44  
    45 subroutine init_icethick 
    46  
    47 implicit none 
    48  
    49 namelist/mass_conserv/adv_frac,V_limit   
    50  
    51  
    52 rewind(num_param)        ! pour revenir au debut du fichier param_list.dat 
    53 read(num_param,mass_conserv) 
    54  
    55 write(num_rep_42,mass_conserv) 
    56 write(num_rep_42,*)'! Conservation de la masse avec equation advection-diffusion ' 
    57 write(num_rep_42,*)'! la repartition depend de adv_frac' 
    58 write(num_rep_42,*)'! >1  -> advection seule' 
    59 write(num_rep_42,*)'! 0   -> diffusion seule' 
    60 write(num_rep_42,*)'! 0<*<1   -> fraction de l advection' 
    61 write(num_rep_42,*)'! -1 -> zones diffusion + zones advection' 
    62 write(num_rep_42,*)'! V_limit depend de la calotte : ' 
    63 write(num_rep_42,*)'! typiquement 3000 m/an en Antarctique, 10000 m/an au Groenland' 
    64 write(num_rep_42,*)'!___________________________________________________________'  
    65  
    66  
    67 call init_reso_adv_diff_2D      
    68 call init_prescribe_H        ! initialize grounding line mask 
    69  
    70 Dx1=1./Dx 
    71 itour=0 
    72  
    73 end subroutine init_icethick 
    74  
    75 !------------------------------------------------------------------ 
    76 !> SUBROUTINE: icethick3 
    77 !! Routine pour le calcul de l'epaisseur de glace  
    78 !< 
    79 subroutine icethick3 
    80  
    81 !$ USE OMP_LIB 
    82  
    83 implicit none 
    84 real,dimension(nx,ny) :: Dminx,Dminy 
    85 real,dimension(nx,ny) :: Uxdiff,Uydiff           ! vitesse due a la diffusion 
    86 real aux                                         ! pour le calcul du critere 
    87 real maxdia                                      ! sur le pas de temps 
    88 integer imaxdia,jmaxdia 
    89  
    90  
    91 if (itracebug.eq.1) call tracebug(" Entree dans routine icethick") 
    92  
    93 ! Copie de H sur vieuxH         
    94 ! -------------------- 
    95 !$OMP PARALLEL 
    96 !$OMP WORKSHARE 
    97 where (mk0(:,:).eq.0) 
    98    vieuxH(:,:)=0. 
    99 elsewhere (i_Hp(:,:).eq.1)                 ! previous prescribed thickness 
    100    vieuxH(:,:)=Hp(:,:)                     ! H may have been changed by calving 
    101 elsewhere 
    102    vieuxH(:,:)=H(:,:) 
    103 end where 
    104 !$OMP END WORKSHARE 
    105  
    106 ! mk0 est le masque à ne jamais dépasser 
    107 ! mk0=0  ->  pas de glace  
    108 ! mk0=-1 ->  on garde la même epaisseur 
    109 ! pour l'Antarctique masque mko=1 partout 
    110  
    111  
    112 Maxdia = -1.0 
    113 imaxdia=0 
    114 jmaxdia=0 
    115  
    116 ! attention limitation ! 
    117 !$OMP WORKSHARE  
    118 uxbar(:,:)=max(uxbar(:,:),-V_limit) 
    119 uxbar(:,:)=min(uxbar(:,:),V_limit) 
    120 uybar(:,:)=max(uybar(:,:),-V_limit) 
    121 uybar(:,:)=min(uybar(:,:),V_limit) 
    122 !$OMP END WORKSHARE 
    123  
    124 ! le pas de temps est choisi pour rendre la matrice advective diagonale dominante 
    125 !$OMP END PARALLEL 
    126 !$OMP PARALLEL 
    127 !$OMP DO PRIVATE(aux) 
    128 !do i=2,nx-1 
    129 !   do j=2,ny-1 
    130 do j=2,ny-1 
    131   do i=2,nx-1 
    132       aux = (abs(uxbar(i,j)) + abs(uxbar(i+1,j))) & 
    133            + (abs(uybar(i,j)) + abs(uybar(i,j+1))) 
    134       aux=aux*dx1*0.5 
    135       !$OMP CRITICAL 
    136       if (aux.gt.maxdia) then 
    137          maxdia=aux 
    138 !         imaxdia=i 
    139 !         jmaxdia=j 
    140       endif 
    141       !$OMP END CRITICAL 
    142    end do 
    143 end do 
    144 !$OMP END DO 
    145 !$OMP END PARALLEL 
    146  
    147 if (maxdia.ge.(testdiag/dtmax)) then 
    148    dt = testdiag/Maxdia 
    149    dt=max(dt,dtmin) 
    150 else 
    151    dt = dtmax 
    152 end if 
    153 !~ print*,'testdiag/Maxdia, ',time, dt   
    154 ! next_time ajuste le pas de temps pour faciliter la synchronisation  
    155 ! avec le pas de temps long (dtt) 
    156  
    157 call next_time(time,dt,dtt,dtmax,dtmin,isynchro,itracebug,num_tracebug) 
    158  
    159  
    160 ! calcul diffusivite - vitesse  
    161 !-----------------------------------------------------------------  
    162 !$OMP PARALLEL 
    163 !$OMP WORKSHARE 
    164 Uxdiff(:,:)=diffmx(:,:)*(-sdx(:,:))   ! vitesse qui peut s'exprimer en terme diffusif 
    165 Uydiff(:,:)=diffmy(:,:)*(-sdy(:,:))   ! dans les where qui suivent elle est limitee a uxbar 
    166                                       ! pour eviter des problemes dans le cas 0< adv_frac <1  
    167 dmx(:,:)=diffmx(:,:) 
    168 dmy(:,:)=diffmy(:,:) 
    169  
    170 ! schema amont pour la diffusion en x 
    171 where (Uxbar(:,:).ge.0.) 
    172   dmx(:,:)=diffmx(:,:)*eoshift(H(:,:),shift=-1,boundary=0d0,dim=1) 
    173   uxdiff(:,:)=min(uxdiff(:,:),uxbar(:,:))       ! pour tenir compte limitation utot 
    174 elsewhere 
    175   dmx(:,:)=diffmx(:,:)*H(:,:) 
    176   uxdiff(:,:)=max(uxdiff(:,:),uxbar(:,:))       ! pour tenir compte limitation utot 
    177 end where 
    178  
    179  
    180  
    181 ! schema amont pour la diffusion en y 
    182 where (uybar(:,:).ge.0.) 
    183   dmy(:,:)=diffmy(:,:)*eoshift(H(:,:),shift=-1,boundary=0d0,dim=2) 
    184   uydiff(:,:)=min(uydiff(:,:),uybar(:,:))       ! pour tenir compte limitation utot 
    185 elsewhere 
    186   dmy(:,:)=dmy(:,:)*H(:,:) 
    187   uydiff(:,:)=max(uydiff(:,:),uybar(:,:))       ! pour tenir compte limitation utot 
    188 end where 
    189 !$OMP END WORKSHARE 
    190  
    191  
    192 ! tests pour la répartition advection - diffusion 
    193 !------------------------------------------------ 
    194  
    195 if (adv_frac.gt.1.) then ! advection seulement 
    196    !$OMP WORKSHARE 
    197    advmx(:,:)=Uxbar(:,:) 
    198    advmy(:,:)=Uybar(:,:) 
    199    Dmx(:,:)=0. 
    200    Dmy(:,:)=0. 
    201    !$OMP END WORKSHARE 
    202 else if (abs(adv_frac).lt.1.e-8) then   ! diffusion seulement 
    203    !$OMP WORKSHARE 
    204    advmx(:,:)=0. 
    205    advmy(:,:)=0. 
    206    !$OMP END WORKSHARE 
    207    ! D reste la valeur au dessus) 
    208  
    209 else if ((adv_frac.ge.1.e-8).and.(adv_frac.le.1.)) then        ! advection-diffusion) 
    210  
    211 ! selon x -------------------------------- 
    212  
    213 ! advection = adv_frac* tout ce qui n'est pas diffusion 
    214 ! diffusion = ce qui peut être exprime en diffusion  
    215 !             + une partie (1-adv_frac) de l'advection 
    216    !$OMP WORKSHARE 
    217    where ((abs(sdx(:,:)).gt.1.e-8).and.(Hmx(:,:).gt.2.))      ! cas general 
    218        advmx(:,:)=(Uxbar(:,:)-Uxdiff(:,:))                   ! tout ce qui n'est pas diffusion 
    219        advmx(:,:)=advmx(:,:)*adv_frac                         ! la fraction adv_frac du precedent 
    220        Dminx(:,:)=-(Uxbar(:,:)-advmx(:,:))/sdx(:,:)          ! ce qui reste exprime en diffusivite 
    221                                                             ! a multiplier par H 
    222  
    223 ! schema amont pour la diffusivite 
    224        where (uxbar(:,:).ge.0.) 
    225           dmx(:,:)=dminx(:,:)*eoshift(H(:,:),shift=-1,boundary=0d0,dim=1) 
     46  !> SUBROUTINE: init_icethick 
     47  !! Definition et initialisation des parametres qui gerent la repartition advection diffusion 
     48  !! @return adv_frac 
     49 
     50  subroutine init_icethick 
     51 
     52    implicit none 
     53 
     54    namelist/mass_conserv/adv_frac,V_limit   
     55 
     56 
     57    rewind(num_param)        ! pour revenir au debut du fichier param_list.dat 
     58    read(num_param,mass_conserv) 
     59 
     60    write(num_rep_42,mass_conserv) 
     61    write(num_rep_42,*)'! Conservation de la masse avec equation advection-diffusion ' 
     62    write(num_rep_42,*)'! la repartition depend de adv_frac' 
     63    write(num_rep_42,*)'! >1  -> advection seule' 
     64    write(num_rep_42,*)'! 0   -> diffusion seule' 
     65    write(num_rep_42,*)'! 0<*<1   -> fraction de l advection' 
     66    write(num_rep_42,*)'! -1 -> zones diffusion + zones advection' 
     67    write(num_rep_42,*)'! V_limit depend de la calotte : ' 
     68    write(num_rep_42,*)'! typiquement 3000 m/an en Antarctique, 10000 m/an au Groenland' 
     69    write(num_rep_42,*)'!___________________________________________________________'  
     70 
     71 
     72    call init_reso_adv_diff_2D      
     73    call init_prescribe_H        ! initialize grounding line mask 
     74 
     75    Dx1=1./Dx 
     76    itour=0 
     77 
     78  end subroutine init_icethick 
     79 
     80  !------------------------------------------------------------------ 
     81  !> SUBROUTINE: icethick3 
     82  !! Routine pour le calcul de l'epaisseur de glace  
     83  !< 
     84  subroutine icethick3 
     85 
     86    !$ USE OMP_LIB 
     87 
     88    implicit none 
     89    integer :: i,j 
     90    real,dimension(nx,ny) :: Dminx,Dminy 
     91    real,dimension(nx,ny) :: Uxdiff,Uydiff           ! vitesse due a la diffusion 
     92    real aux                                         ! pour le calcul du critere 
     93    real maxdia                                      ! sur le pas de temps 
     94    integer imaxdia,jmaxdia 
     95 
     96 
     97    if (itracebug.eq.1) call tracebug(" Entree dans routine icethick") 
     98 
     99    ! Copie de H sur vieuxH         
     100    ! -------------------- 
     101    !$OMP PARALLEL 
     102    !$OMP WORKSHARE 
     103    where (mk0(:,:).eq.0) 
     104       vieuxH(:,:)=0. 
     105    elsewhere (i_Hp(:,:).eq.1)                 ! previous prescribed thickness 
     106       vieuxH(:,:)=Hp(:,:)                     ! H may have been changed by calving 
     107    elsewhere 
     108       vieuxH(:,:)=H(:,:) 
     109    end where 
     110    !$OMP END WORKSHARE 
     111 
     112    ! mk0 est le masque à ne jamais dépasser 
     113    ! mk0=0  ->  pas de glace  
     114    ! mk0=-1 ->  on garde la même epaisseur 
     115    ! pour l'Antarctique masque mko=1 partout 
     116 
     117 
     118    Maxdia = -1.0 
     119    imaxdia=0 
     120    jmaxdia=0 
     121 
     122    ! attention limitation ! 
     123    !$OMP WORKSHARE  
     124    uxbar(:,:)=max(uxbar(:,:),-V_limit) 
     125    uxbar(:,:)=min(uxbar(:,:),V_limit) 
     126    uybar(:,:)=max(uybar(:,:),-V_limit) 
     127    uybar(:,:)=min(uybar(:,:),V_limit) 
     128    !$OMP END WORKSHARE 
     129 
     130    ! le pas de temps est choisi pour rendre la matrice advective diagonale dominante 
     131    !$OMP END PARALLEL 
     132    !$OMP PARALLEL 
     133    !$OMP DO PRIVATE(aux) 
     134    !do i=2,nx-1 
     135    !   do j=2,ny-1 
     136    do j=2,ny-1 
     137       do i=2,nx-1 
     138          aux = (abs(uxbar(i,j)) + abs(uxbar(i+1,j))) & 
     139               + (abs(uybar(i,j)) + abs(uybar(i,j+1))) 
     140          aux=aux*dx1*0.5 
     141          !$OMP CRITICAL 
     142          if (aux.gt.maxdia) then 
     143             maxdia=aux 
     144             !         imaxdia=i 
     145             !         jmaxdia=j 
     146          endif 
     147          !$OMP END CRITICAL 
     148       end do 
     149    end do 
     150    !$OMP END DO 
     151    !$OMP END PARALLEL 
     152 
     153    if (maxdia.ge.(testdiag/dtmax)) then 
     154       dt = testdiag/Maxdia 
     155       dt=max(dt,dtmin) 
     156    else 
     157       dt = dtmax 
     158    end if 
     159    !~ print*,'testdiag/Maxdia, ',time, dt   
     160    ! next_time ajuste le pas de temps pour faciliter la synchronisation  
     161    ! avec le pas de temps long (dtt) 
     162 
     163    call next_time(time,dt,dtt,dtmax,dtmin,isynchro,itracebug,num_tracebug) 
     164 
     165 
     166    ! calcul diffusivite - vitesse  
     167    !-----------------------------------------------------------------  
     168    !$OMP PARALLEL 
     169    !$OMP WORKSHARE 
     170    Uxdiff(:,:)=diffmx(:,:)*(-sdx(:,:))   ! vitesse qui peut s'exprimer en terme diffusif 
     171    Uydiff(:,:)=diffmy(:,:)*(-sdy(:,:))   ! dans les where qui suivent elle est limitee a uxbar 
     172    ! pour eviter des problemes dans le cas 0< adv_frac <1  
     173    dmx(:,:)=diffmx(:,:) 
     174    dmy(:,:)=diffmy(:,:) 
     175 
     176    ! schema amont pour la diffusion en x 
     177    where (Uxbar(:,:).ge.0.) 
     178       dmx(:,:)=diffmx(:,:)*eoshift(H(:,:),shift=-1,boundary=0d0,dim=1) 
     179       uxdiff(:,:)=min(uxdiff(:,:),uxbar(:,:))       ! pour tenir compte limitation utot 
     180    elsewhere 
     181       dmx(:,:)=diffmx(:,:)*H(:,:) 
     182       uxdiff(:,:)=max(uxdiff(:,:),uxbar(:,:))       ! pour tenir compte limitation utot 
     183    end where 
     184 
     185 
     186 
     187    ! schema amont pour la diffusion en y 
     188    where (uybar(:,:).ge.0.) 
     189       dmy(:,:)=diffmy(:,:)*eoshift(H(:,:),shift=-1,boundary=0d0,dim=2) 
     190       uydiff(:,:)=min(uydiff(:,:),uybar(:,:))       ! pour tenir compte limitation utot 
     191    elsewhere 
     192       dmy(:,:)=dmy(:,:)*H(:,:) 
     193       uydiff(:,:)=max(uydiff(:,:),uybar(:,:))       ! pour tenir compte limitation utot 
     194    end where 
     195    !$OMP END WORKSHARE 
     196 
     197 
     198    ! tests pour la répartition advection - diffusion 
     199    !------------------------------------------------ 
     200 
     201    if (adv_frac.gt.1.) then ! advection seulement 
     202       !$OMP WORKSHARE 
     203       advmx(:,:)=Uxbar(:,:) 
     204       advmy(:,:)=Uybar(:,:) 
     205       Dmx(:,:)=0. 
     206       Dmy(:,:)=0. 
     207       !$OMP END WORKSHARE 
     208    else if (abs(adv_frac).lt.1.e-8) then   ! diffusion seulement 
     209       !$OMP WORKSHARE 
     210       advmx(:,:)=0. 
     211       advmy(:,:)=0. 
     212       !$OMP END WORKSHARE 
     213       ! D reste la valeur au dessus) 
     214 
     215    else if ((adv_frac.ge.1.e-8).and.(adv_frac.le.1.)) then        ! advection-diffusion) 
     216 
     217       ! selon x -------------------------------- 
     218 
     219       ! advection = adv_frac* tout ce qui n'est pas diffusion 
     220       ! diffusion = ce qui peut être exprime en diffusion  
     221       !             + une partie (1-adv_frac) de l'advection 
     222       !$OMP WORKSHARE 
     223       where ((abs(sdx(:,:)).gt.1.e-8).and.(Hmx(:,:).gt.2.))      ! cas general 
     224          advmx(:,:)=(Uxbar(:,:)-Uxdiff(:,:))                   ! tout ce qui n'est pas diffusion 
     225          advmx(:,:)=advmx(:,:)*adv_frac                         ! la fraction adv_frac du precedent 
     226          Dminx(:,:)=-(Uxbar(:,:)-advmx(:,:))/sdx(:,:)          ! ce qui reste exprime en diffusivite 
     227          ! a multiplier par H 
     228 
     229          ! schema amont pour la diffusivite 
     230          where (uxbar(:,:).ge.0.) 
     231             dmx(:,:)=dminx(:,:)*eoshift(H(:,:),shift=-1,boundary=0d0,dim=1) 
     232          elsewhere 
     233             dmx(:,:)=dminx(:,:)*H(:,:) 
     234          end where 
     235 
     236 
     237       elsewhere   ! zones trop plates ou trop fines 
     238          Dmx(:,:)=0. 
     239          advmx(:,:)=Uxbar(:,:) 
     240       end where 
     241 
     242 
     243       ! selon y -------------------------------- 
     244 
     245       ! advection = adv_frac* tout ce qui n'est pas diffusion 
     246       ! diffusion = ce qui peut être exprime en diffusion  
     247       !             + une partie (1-adv_frac) de l'advection 
     248 
     249       where ((abs(sdy(:,:)).gt.1.e-8).and.(Hmy(:,:).gt.2.))      ! cas general 
     250          advmy(:,:)=(Uybar(:,:)-Uydiff(:,:))                    ! tout ce qui n'est pas diffusion 
     251          advmy(:,:)=advmy(:,:)*adv_frac                         ! la fraction adv_frac du precedent 
     252          Dminy(:,:)=-(Uybar(:,:)-advmy(:,:))/sdy(:,:)          ! ce qui reste exprime en diffusivite 
     253          ! a multiplier par H 
     254 
     255          ! schema amont pour la diffusivite 
     256          where (uybar(:,:).ge.0.) 
     257             dmy(:,:)=dminy(:,:)*eoshift(H(:,:),shift=-1,boundary=0d0,dim=2) 
     258          elsewhere 
     259             dmy(:,:)=dminy(:,:)*H(:,:) 
     260          end where 
     261       elsewhere   ! zones trop plates ou trop fines 
     262          Dmy(:,:)=0. 
     263          advmy(:,:)=Uybar(:,:) 
     264       end where 
     265       !$OMP END WORKSHARE 
     266 
     267    else if (adv_frac.lt.0) then                  ! diffusif dans zones SIA, advectif ailleurs 
     268       !$OMP WORKSHARE 
     269       zonemx(:,:)=flgzmx(:,:) 
     270       zonemy(:,:)=flgzmy(:,:) 
     271 
     272       ! selon x -------------- 
     273       where (zonemx(:,:)) 
     274          advmx(:,:)=Uxbar(:,:) 
     275          Dmx(:,:)=0. 
    226276       elsewhere 
    227           dmx(:,:)=dminx(:,:)*H(:,:) 
     277          advmx(:,:)=0. 
    228278       end where 
    229279 
    230  
    231    elsewhere   ! zones trop plates ou trop fines 
    232       Dmx(:,:)=0. 
    233       advmx(:,:)=Uxbar(:,:) 
    234    end where 
    235  
    236  
    237 ! selon y -------------------------------- 
    238  
    239 ! advection = adv_frac* tout ce qui n'est pas diffusion 
    240 ! diffusion = ce qui peut être exprime en diffusion  
    241 !             + une partie (1-adv_frac) de l'advection 
    242  
    243    where ((abs(sdy(:,:)).gt.1.e-8).and.(Hmy(:,:).gt.2.))      ! cas general 
    244        advmy(:,:)=(Uybar(:,:)-Uydiff(:,:))                    ! tout ce qui n'est pas diffusion 
    245        advmy(:,:)=advmy(:,:)*adv_frac                         ! la fraction adv_frac du precedent 
    246        Dminy(:,:)=-(Uybar(:,:)-advmy(:,:))/sdy(:,:)          ! ce qui reste exprime en diffusivite 
    247                                                             ! a multiplier par H 
    248  
    249 ! schema amont pour la diffusivite 
    250        where (uybar(:,:).ge.0.) 
    251           dmy(:,:)=dminy(:,:)*eoshift(H(:,:),shift=-1,boundary=0d0,dim=2) 
     280       ! selon y -------------- 
     281       where (zonemy(:,:)) 
     282          advmy(:,:)=Uybar(:,:) 
     283          Dmy(:,:)=0. 
    252284       elsewhere 
    253           dmy(:,:)=dminy(:,:)*H(:,:) 
     285          advmy(:,:)=0. 
    254286       end where 
    255    elsewhere   ! zones trop plates ou trop fines 
    256       Dmy(:,:)=0. 
    257       advmy(:,:)=Uybar(:,:) 
    258    end where 
    259    !$OMP END WORKSHARE 
    260  
    261 else if (adv_frac.lt.0) then                  ! diffusif dans zones SIA, advectif ailleurs 
    262    !$OMP WORKSHARE 
    263    zonemx(:,:)=flgzmx(:,:) 
    264    zonemy(:,:)=flgzmy(:,:) 
    265  
    266 ! selon x -------------- 
    267    where (zonemx(:,:)) 
    268     advmx(:,:)=Uxbar(:,:) 
    269     Dmx(:,:)=0. 
    270    elsewhere 
    271     advmx(:,:)=0. 
    272    end where 
    273  
    274 ! selon y -------------- 
    275    where (zonemy(:,:)) 
    276     advmy(:,:)=Uybar(:,:) 
    277     Dmy(:,:)=0. 
    278    elsewhere 
    279     advmy(:,:)=0. 
    280    end where 
    281    !$OMP END WORKSHARE 
    282 end if 
    283  
    284 !$OMP WORKSHARE 
    285 where (flot(:,:) ) 
    286    tabtest(:,:)=1. 
    287 elsewhere 
    288    tabtest(:,:)=0. 
    289 end where 
    290 !$OMP END WORKSHARE 
    291 !$OMP END PARALLEL 
    292 !------------------------------------------------------------------detect 
     287       !$OMP END WORKSHARE 
     288    end if 
     289 
     290    !$OMP WORKSHARE 
     291    where (flot(:,:) ) 
     292       tabtest(:,:)=1. 
     293    elsewhere 
     294       tabtest(:,:)=0. 
     295    end where 
     296    !$OMP END WORKSHARE 
     297    !$OMP END PARALLEL 
     298    !------------------------------------------------------------------detect 
    293299!!$ tabtest(:,:)=dmx(:,:) 
    294300!!$ 
    … …  
    301307!!$   write(6,*) ' pas d asymetrie sur dmx avant resol pour time=',time,itesti 
    302308!!$end if 
    303 !----------------------------------------------------------- fin detect   
    304  
    305  
    306 !------------------------------------------------------------------detect 
     309    !----------------------------------------------------------- fin detect   
     310 
     311 
     312    !------------------------------------------------------------------detect 
    307313!!$ tabtest(:,:)=dmy(:,:) 
    308314!!$ 
    … …  
    315321!!$   write(6,*) ' pas d asymetrie sur dmy avant resol pour time=',time,itesti 
    316322!!$end if 
    317 !----------------------------------------------------------- fin detect   
    318  
    319 !------------------------------------------------------------------detect 
     323    !----------------------------------------------------------- fin detect   
     324 
     325    !------------------------------------------------------------------detect 
    320326!!$ tabtest(:,:)=advmx(:,:) 
    321327!!$ 
    … …  
    328334!!$   write(6,*) ' pas d asymetrie sur advdmx avant resol pour time=',time,itesti 
    329335!!$end if 
    330 !----------------------------------------------------------- fin detect   
    331  
    332 !------------------------------------------------------------------detect 
     336    !----------------------------------------------------------- fin detect   
     337 
     338    !------------------------------------------------------------------detect 
    333339!!$ tabtest(:,:)=advmy(:,:) 
    334340!!$ 
    … …  
    341347!!$   write(6,*) ' pas d asymetrie sur advmy avant resol pour time=',time,itesti 
    342348!!$end if 
    343 !----------------------------------------------------------- fin detect   
    344  
    345 ! nouveau calcul de dt 
    346 aux=maxval( (abs(advmx(:,:))+abs(advmy(:,:)))/dx+(abs(dmx(:,:))+abs(dmy(:,:)))/dx/dx) 
    347  
    348 ! write(166,*) 'critere pour dt',time-dt,testdiag/aux,dt 
    349  
    350  
    351 if (aux.gt.1.e-20) then 
    352    if (testdiag/aux.lt.dt) then 
    353       time=time-dt 
    354       dt=min(testdiag/aux,dtmax)   !cdc *4. 
    355 !~       print*,'aux*4, ',time, dt, aux, testdiag/aux 
    356                         if (itracebug.eq.1) call tracebug("aux tres petit,deuxieme appel a next_time dt*4") 
    357       call next_time(time,dt,dtt,dtmax,dtmin,isynchro,itracebug,num_tracebug) 
    358    else 
    359 !~       print*,'testdiag/Maxdia, ',time, dt   
    360    end if 
    361 end if 
    362  
    363  
    364 timemax=time 
    365  
    366 ! etait avant dans step 
    367 if (time.lt.TGROUNDED) then 
    368    MARINE=.false. 
    369 else 
    370    MARINE=.true. 
    371 endif 
    372 ! fin vient de step 
    373  
    374 ! les variables dtdx et dtdx2 sont globales 
    375  Dtdx2=Dt/(Dx**2) 
    376  dtdx=dt/dx 
    377  
    378  
    379 if (geoplace(1:5).eq.'mism3') then       ! pas de flux sur les limites laterales 
    380    dmy(:,2)       = 0. 
    381    dmy(:,3)       = 0. 
    382    dmy(:,ny-1)    = 0. 
    383    dmy(:,ny)      = 0. 
    384    advmy(:,2)     = 0. 
    385    advmy(:,3)     = 0. 
    386    advmy(:,ny-1)  = 0. 
    387    advmy(:,ny)    = 0. 
    388    uybar(:,2)     = 0 
    389    uybar(:,3)     = 0 
    390    uybar(:,ny-1)  = 0 
    391    uybar(:,ny)    = 0    
    392 end if 
    393  
    394  
    395 !debug_3D(:,:,45)=dmx(:,:) 
    396 !debug_3D(:,:,46)=dmy(:,:) 
    397 !debug_3D(:,:,47)=advmx(:,:) 
    398 !debug_3D(:,:,48)=advmy(:,:) 
    399  
    400 !$OMP PARALLEL 
    401 !$OMP WORKSHARE 
    402 bilmass(:,:)=bm(:,:)-bmelt(:,:)                       ! surface and bottom mass balance 
    403 !$OMP END WORKSHARE 
    404 !$OMP END PARALLEL 
    405  
    406 ! diverses precription de l'epaisseur en fonction de l'objectif 
    407 !------------------------------------------------------------------- 
    408  
    409  
    410 call prescribe_fixed_points            ! ceux qui sont fixes pour tout le run (bords de grille, regions) 
    411  
    412  
    413  if ((igrdline.eq.1)) then               ! present grounding line 
    414  
    415 !    if ((time.lt.time_gl(1)).or.(nt.lt.2)) then      ATTENTION test seulement si version prescribe-H_mod.f90 
    416    if (ibmelt_inv.eq.0) then 
    417       call prescribe_present_H_gl       ! prescribe present thickness at the grounding line 
    418    else if (ibmelt_inv.eq.1) then ! inversion bmelt 
    419       call prescribe_present_H_gl_bmelt ! prescribe only grounding line points 
    420    else 
    421       print*,'ibmelt_inv value must be 0 or 1 !!!' 
    422       stop 
    423    endif 
    424 !    else 
    425 !       call prescribe_retreat 
    426 !    endif 
    427  end if 
    428  
    429  if ((igrdline.eq.2)) then               ! paleo grounding line 
    430     call prescribe_paleo_gl_shelf 
    431  end if 
    432  
    433  if ((igrdline.eq.3)) then   
    434 !    where (flot(:,:)) 
    435 !       mk_gr(:,:) = 0 
    436 !    elsewhere 
    437 !       mk_gr(:,:) = 1 
    438 !    end where 
    439     if (itracebug.eq.1) call tracebug(" avant time_step_recul") 
    440     call time_step_recul                 ! version prescribe-H-i2s_mod.f90 avec use proto_recul_mod  
    441     if (itracebug.eq.1) call tracebug(" apres time_step_recul") 
    442 !    call prescr_ice2sea_retreat         ! version prescribe-H_mod.f90 
    443  endif 
    444  
    445  
    446 !if (time.ge.t_break) then                ! action sur les ice shelves 
    447 !   call melt_ice_shelves                 ! ATTENTION version prescribe-H_mod.f90 
    448 !   call break_all_ice_shelves 
    449 !end if 
     349    !----------------------------------------------------------- fin detect   
     350 
     351    ! nouveau calcul de dt 
     352    aux=maxval( (abs(advmx(:,:))+abs(advmy(:,:)))/dx+(abs(dmx(:,:))+abs(dmy(:,:)))/dx/dx) 
     353 
     354    ! write(166,*) 'critere pour dt',time-dt,testdiag/aux,dt 
     355 
     356 
     357    if (aux.gt.1.e-20) then 
     358       if (testdiag/aux.lt.dt) then 
     359          time=time-dt 
     360          dt=min(testdiag/aux,dtmax)   !cdc *4. 
     361          !~       print*,'aux*4, ',time, dt, aux, testdiag/aux 
     362          if (itracebug.eq.1) call tracebug("aux tres petit,deuxieme appel a next_time dt*4") 
     363          call next_time(time,dt,dtt,dtmax,dtmin,isynchro,itracebug,num_tracebug) 
     364       else 
     365          !~       print*,'testdiag/Maxdia, ',time, dt   
     366       end if 
     367    end if 
     368 
     369 
     370    timemax=time 
     371 
     372    ! etait avant dans step 
     373    if (time.lt.TGROUNDED) then 
     374       MARINE=.false. 
     375    else 
     376       MARINE=.true. 
     377    endif 
     378    ! fin vient de step 
     379 
     380    ! les variables dtdx et dtdx2 sont globales 
     381    Dtdx2=Dt/(Dx**2) 
     382    dtdx=dt/dx 
     383 
     384 
     385    if (geoplace(1:5).eq.'mism3') then       ! pas de flux sur les limites laterales 
     386       dmy(:,2)       = 0. 
     387       dmy(:,3)       = 0. 
     388       dmy(:,ny-1)    = 0. 
     389       dmy(:,ny)      = 0. 
     390       advmy(:,2)     = 0. 
     391       advmy(:,3)     = 0. 
     392       advmy(:,ny-1)  = 0. 
     393       advmy(:,ny)    = 0. 
     394       uybar(:,2)     = 0 
     395       uybar(:,3)     = 0 
     396       uybar(:,ny-1)  = 0 
     397       uybar(:,ny)    = 0    
     398    end if 
     399 
     400 
     401    !debug_3D(:,:,45)=dmx(:,:) 
     402    !debug_3D(:,:,46)=dmy(:,:) 
     403    !debug_3D(:,:,47)=advmx(:,:) 
     404    !debug_3D(:,:,48)=advmy(:,:) 
     405 
     406    !$OMP PARALLEL 
     407    !$OMP WORKSHARE 
     408    bilmass(:,:)=bm(:,:)-bmelt(:,:)                       ! surface and bottom mass balance 
     409    !$OMP END WORKSHARE 
     410    !$OMP END PARALLEL 
     411 
     412    ! diverses precription de l'epaisseur en fonction de l'objectif 
     413    !------------------------------------------------------------------- 
     414 
     415 
     416    call prescribe_fixed_points            ! ceux qui sont fixes pour tout le run (bords de grille, regions) 
     417 
     418 
     419    if ((igrdline.eq.1)) then               ! present grounding line 
     420 
     421       !    if ((time.lt.time_gl(1)).or.(nt.lt.2)) then      ATTENTION test seulement si version prescribe-H_mod.f90 
     422       if (ibmelt_inv.eq.0) then 
     423          call prescribe_present_H_gl       ! prescribe present thickness at the grounding line 
     424       else if (ibmelt_inv.eq.1) then ! inversion bmelt 
     425          call prescribe_present_H_gl_bmelt ! prescribe only grounding line points 
     426       else 
     427          print*,'ibmelt_inv value must be 0 or 1 !!!' 
     428          stop 
     429       endif 
     430       !    else 
     431       !       call prescribe_retreat 
     432       !    endif 
     433    end if 
     434 
     435    if ((igrdline.eq.2)) then               ! paleo grounding line 
     436       call prescribe_paleo_gl_shelf 
     437    end if 
     438 
     439    if ((igrdline.eq.3)) then   
     440       !    where (flot(:,:)) 
     441       !       mk_gr(:,:) = 0 
     442       !    elsewhere 
     443       !       mk_gr(:,:) = 1 
     444       !    end where 
     445       if (itracebug.eq.1) call tracebug(" avant time_step_recul") 
     446       call time_step_recul                 ! version prescribe-H-i2s_mod.f90 avec use proto_recul_mod  
     447       if (itracebug.eq.1) call tracebug(" apres time_step_recul") 
     448       !    call prescr_ice2sea_retreat         ! version prescribe-H_mod.f90 
     449    endif 
     450 
     451 
     452    !if (time.ge.t_break) then                ! action sur les ice shelves 
     453    !   call melt_ice_shelves                 ! ATTENTION version prescribe-H_mod.f90 
     454    !   call break_all_ice_shelves 
     455    !end if 
    450456 
    451457 
    … …  
    461467!!$end if 
    462468 
    463 !debug_3D(:,:,87)=S(:,:)-debug_3D(:,:,88) 
    464  
    465 !debug_3D(:,:,86)=Hp(:,:) 
    466 !debug_3D(:,:,85)=i_Hp(:,:) 
     469    !debug_3D(:,:,87)=S(:,:)-debug_3D(:,:,88) 
     470 
     471    !debug_3D(:,:,86)=Hp(:,:) 
     472    !debug_3D(:,:,85)=i_Hp(:,:) 
    467473 
    468474!!$where (i_hp(:,:).eq.1) 
    … …  
    471477 
    472478    !$OMP WORKSHARE 
    473 !cdc    where (flot(:,:)) 
    474 !cdc 1m       where (H(:,:).gt.max(Hmin,Hmin+BM(:,:)-Bmelt(:,:))) 
    475 !cdc      where (H(:,:).gt.0.) 
    476 !cdc          ice(:,:)=1 
    477 !cdc       elsewhere 
    478 !cdc          ice(:,:)=0 
    479 !cdc       end where 
    480 !cdc    elsewhere 
    481        where (H(:,:).gt.0.) 
    482           ice(:,:)=1 
    483        elsewhere 
    484           ice(:,:)=0 
    485        end where 
    486 !cdc    end where 
    487     !$OMP END WORKSHARE 
    488      
    489 ! Appel a la routine de resolution resol_adv_diff_2D_vect 
    490 !------------------------------------------------------------ 
    491 ! On ecrit les vitesses sous la forme 
    492 ! Ux = -Dmx * sdx + advmx 
    493  
    494 ! Dmx, Dmy sont les termes diffusifs de la vitesse 
    495 ! advmx, advmy sont les termes advectifs. 
    496 ! la repartition entre les deux depend de adv_frac. 
    497  
    498 ! i_Hp est un tableau des points ou l'epaisseur est imposee a la valeur Hp 
    499 ! bilmass est le bilan de masse (surface et fond) 
    500 ! vieuxH est la precedente valeur de H 
    501 ! H est le retour de la nouvelle valeur. 
    502  
    503  
    504 call resol_adv_diff_2D_vect (dble(Dmx),dble(Dmy),dble(advmx),dble(advmy),dble(Hp),i_Hp,dble(bilmass),dble(vieuxH),H) 
    505  
    506 !$OMP PARALLEL 
    507 !$OMP DO 
    508 ! on met à 0 l'épaisseur sur les bords si nécessaire 
    509 do i=1,nx 
    510 !cdc 1m if (H(i,1).gt.max(Hmin,Hmin+BM(i,1)-Bmelt(i,1))) then 
    511         if (H(i,1).gt.0.) then 
    512                 ice(i,1)=1 
    513   else 
    514     ice(i,1)=0 
    515     H(i,1)=0. 
    516   endif   
    517 !cdc 1m if (H(i,ny).gt.max(Hmin,Hmin+BM(i,ny)-Bmelt(i,ny))) then 
    518         if (H(i,ny).gt.0.) then 
    519                 ice(i,ny)=1 
    520   else 
    521     ice(i,ny)=0 
    522     H(i,ny)=0. 
    523   endif 
    524 enddo 
    525 !$OMP END DO 
    526 !$OMP DO 
    527 do j=1,ny 
    528 !cdc 1m if (H(1,j).gt.max(Hmin,Hmin+BM(1,j)-Bmelt(1,j))) then 
    529         if (H(1,j).gt.0.) then 
    530                 ice(1,j)=1 
    531   else 
    532     ice(1,j)=0 
    533     H(1,j)=0. 
    534   endif   
    535 !cdc 1m if (H(nx,j).gt.max(Hmin,Hmin+BM(nx,j)-Bmelt(nx,j))) then 
    536         if (H(nx,j).gt.0.) then 
    537                 ice(nx,j)=1 
    538   else 
    539     ice(nx,j)=0 
    540     H(nx,j)=0. 
    541   endif 
    542 enddo 
    543 !$OMP END DO 
    544  
    545 ! remise a 0 des epaisseurs negatives, on garde la difference dans ablbord 
    546  
    547 !$OMP WORKSHARE 
    548 where (H(:,:).lt.0.) 
    549 !where ((H(:,:).lt.0.).OR.(H(:,:).GT.0..AND.H(:,:).LT.min(1.,vieuxH(:,:)).AND.flot(:,:))) 
    550    ablbord(:,:)=H(:,:)  ! ici ablbord = Hcons_mass = Hadv + Bm - Bmelt 
    551 elsewhere 
    552    ablbord(:,:)=0. 
    553 endwhere 
    554  
    555 H(:,:)=max(H(:,:),0.)       ! pas d'epaisseur negative 
    556  
    557 ! eventuellement retirer apres spinup ou avoir un cas serac 
    558 Hdot(:,:)=(H(:,:)-vieuxH(:,:))/dt 
    559 !$OMP END WORKSHARE 
    560  
    561 if (ibmelt_inv.eq.1) then ! inversion du bmelt avec igrdline=1 
    562    !$OMP WORKSHARE 
    563    where (i_Hp(:,:).eq.1) 
    564       H(:,:)=Hp(:,:) 
    565    endwhere 
    566    !$OMP END WORKSHARE 
    567 endif 
    568  
    569 !$OMP WORKSHARE 
    570 ! si mk0=-1, on garde l'epaisseur precedente 
    571 ! en attendant un masque plus propre 
    572  
    573 !~ where(mk0(:,:).eq.-1) 
    574 !~    H(:,:)=vieuxH(:,:) 
    575 !~    Hdot(:,:)=0. 
    576 !~ end where 
    577 !$OMP END WORKSHARE 
    578 !$OMP END PARALLEL 
    579  
    580 if (itracebug.eq.1)  call tracebug(' Fin routine icethick') !, maxval(H) ',maxval(H(:,:)) 
    581  
    582 end subroutine icethick3 
     479    !cdc    where (flot(:,:)) 
     480    !cdc 1m       where (H(:,:).gt.max(Hmin,Hmin+BM(:,:)-Bmelt(:,:))) 
     481    !cdc      where (H(:,:).gt.0.) 
     482    !cdc          ice(:,:)=1 
     483    !cdc       elsewhere 
     484    !cdc          ice(:,:)=0 
     485    !cdc       end where 
     486    !cdc    elsewhere 
     487    where (H(:,:).gt.0.) 
     488       ice(:,:)=1 
     489    elsewhere 
     490       ice(:,:)=0 
     491    end where 
     492    !cdc    end where 
     493    !$OMP END WORKSHARE 
     494 
     495    ! Appel a la routine de resolution resol_adv_diff_2D_vect 
     496    !------------------------------------------------------------ 
     497    ! On ecrit les vitesses sous la forme 
     498    ! Ux = -Dmx * sdx + advmx 
     499 
     500    ! Dmx, Dmy sont les termes diffusifs de la vitesse 
     501    ! advmx, advmy sont les termes advectifs. 
     502    ! la repartition entre les deux depend de adv_frac. 
     503 
     504    ! i_Hp est un tableau des points ou l'epaisseur est imposee a la valeur Hp 
     505    ! bilmass est le bilan de masse (surface et fond) 
     506    ! vieuxH est la precedente valeur de H 
     507    ! H est le retour de la nouvelle valeur. 
     508 
     509 
     510    call resol_adv_diff_2D_vect (dble(Dmx),dble(Dmy),dble(advmx),dble(advmy),dble(Hp),i_Hp,dble(bilmass),dble(vieuxH),H) 
     511 
     512    !$OMP PARALLEL 
     513    !$OMP DO 
     514    ! on met à 0 l'épaisseur sur les bords si nécessaire 
     515    do i=1,nx 
     516       !cdc 1m  if (H(i,1).gt.max(Hmin,Hmin+BM(i,1)-Bmelt(i,1))) then 
     517       if (H(i,1).gt.0.) then 
     518          ice(i,1)=1 
     519       else 
     520          ice(i,1)=0 
     521          H(i,1)=0. 
     522       endif 
     523       !cdc 1m  if (H(i,ny).gt.max(Hmin,Hmin+BM(i,ny)-Bmelt(i,ny))) then 
     524       if (H(i,ny).gt.0.) then 
     525          ice(i,ny)=1 
     526       else 
     527          ice(i,ny)=0 
     528          H(i,ny)=0. 
     529       endif 
     530    enddo 
     531    !$OMP END DO 
     532    !$OMP DO 
     533    do j=1,ny 
     534       !cdc 1m  if (H(1,j).gt.max(Hmin,Hmin+BM(1,j)-Bmelt(1,j))) then 
     535       if (H(1,j).gt.0.) then 
     536          ice(1,j)=1 
     537       else 
     538          ice(1,j)=0 
     539          H(1,j)=0. 
     540       endif 
     541       !cdc 1m  if (H(nx,j).gt.max(Hmin,Hmin+BM(nx,j)-Bmelt(nx,j))) then 
     542       if (H(nx,j).gt.0.) then 
     543          ice(nx,j)=1 
     544       else 
     545          ice(nx,j)=0 
     546          H(nx,j)=0. 
     547       endif 
     548    enddo 
     549    !$OMP END DO 
     550 
     551    ! remise a 0 des epaisseurs negatives, on garde la difference dans ablbord 
     552 
     553    !$OMP WORKSHARE 
     554    where (H(:,:).lt.0.) 
     555       !where ((H(:,:).lt.0.).OR.(H(:,:).GT.0..AND.H(:,:).LT.min(1.,vieuxH(:,:)).AND.flot(:,:))) 
     556       ablbord(:,:)=H(:,:)  ! ici ablbord = Hcons_mass = Hadv + Bm - Bmelt 
     557    elsewhere 
     558       ablbord(:,:)=0. 
     559    endwhere 
     560 
     561    H(:,:)=max(H(:,:),0.)       ! pas d'epaisseur negative 
     562 
     563    ! eventuellement retirer apres spinup ou avoir un cas serac 
     564    Hdot(:,:)=(H(:,:)-vieuxH(:,:))/dt 
     565    !$OMP END WORKSHARE 
     566 
     567    if (ibmelt_inv.eq.1) then ! inversion du bmelt avec igrdline=1 
     568       !$OMP WORKSHARE 
     569       where (i_Hp(:,:).eq.1) 
     570          H(:,:)=Hp(:,:) 
     571       endwhere 
     572       !$OMP END WORKSHARE 
     573    endif 
     574 
     575    !$OMP WORKSHARE 
     576    ! si mk0=-1, on garde l'epaisseur precedente 
     577    ! en attendant un masque plus propre 
     578 
     579    !~ where(mk0(:,:).eq.-1) 
     580    !~    H(:,:)=vieuxH(:,:) 
     581    !~    Hdot(:,:)=0. 
     582    !~ end where 
     583    !$OMP END WORKSHARE 
     584    !$OMP END PARALLEL 
     585 
     586    if (itracebug.eq.1)  call tracebug(' Fin routine icethick') !, maxval(H) ',maxval(H(:,:)) 
     587 
     588  end subroutine icethick3 
    583589 
    584590end module equat_adv_diff_2D_vect 
Note: See TracChangeset for help on using the changeset viewer.