Changeset 392


Ignore:
Timestamp:
03/23/23 14:56:53 (13 months ago)
Author:
dumas
Message:

use only in subroutine diffusiv

File:
1 edited

Legend:

Unmodified
Added
Removed

  • branches/GRISLIv3/SOURCES/diffusiv-polyn-0.6.f90

    r76 r392  
    1212!! 
    1313!< 
    14       subroutine diffusiv() 
    15  
    16 !     =============================================================== 
    17 !     **     DIFFUSIVITIES                    18 mars 2006   *** 
    18 !     glissement avec reserve d'eau 29 Avril 97 
    19 !     choix sur le type de discretisation (locale iloc=1 )     
    20 ! 
    21 !     Modification Vince --> 2 fev 95 
    22 !     pour couplage avec Shelf 
    23 ! 
    24 !     Modification C. Dumas Dec 99 
    25 !     DDX, SA,S2A,... ont une dimension en plus (loi de deformation) 
    26 ! 
    27 !     Modifications  de Cat du 11 juin 2000 
    28 !     Modif Christophe 24 janv 2002 : passage au fortran 90 
    29 !     pour optimisation de la routine 
    30 !     Modif Vincent dec 2003 : on utilise plus frontmx/y > commentaires 
    31 ! 
    32 !     Modif Cat mars 2003 
    33 !                  modif mars 2007 : moyennes de S2a avec S2a_mx et S2a_my 
    34  
    35 !     =============================================================== 
    36 !$ USE OMP_LIB 
    37 USE module3D_phy 
    38  
    39 use module_choix 
    40  
    41 implicit none 
    42  
    43 REAL :: ulgliss,ultot,uldef,glenexp 
    44 REAL :: INV_4DX ! inverse de dx pour eviter les divisions =1/(4*dx) 
    45 REAL :: INV_4DY ! inverse de dy pour eviter les divisions =1/(4*dy) 
    46  
    47 if (itracebug.eq.1)  call tracebug(' Entree dans routine diffusiv') 
    48  
    49 !     la valeur de ff est donnee dans initial 
    50 !     limitation de la vitesse de glissement de déformation et totale 
    51 ulgliss = V_limit 
    52 ultot   = V_limit 
    53 uldef   = 2000. 
    54  
    55  
    56 INV_4DX=1/(4*dx) 
    57 INV_4DY=1/(4*dy) 
    58  
    59 ! initialisation de la diffusion 
    60 !$OMP PARALLEL 
    61 !$OMP WORKSHARE 
    62 diffmx(:,:)=0. 
    63 diffmy(:,:)=0. 
    64 !$OMP END WORKSHARE 
    65 !$OMP END PARALLEL 
    66  
    67 ! calcul de SDX, SDY et de la pente au carre en mx et en my : 
    68  
    69 call slope_surf 
    70  
    71 call sliding     ! au sens vitesse de glissement 
    72  
    73 !  le glissement est maintenant dans un module a part choisi dans le module choix 
    74 !  pour l'instant seules les lois Heino (loigliss=4) et Bindshadler(loigliss=2) 
    75 !  sont programmees. 
    76  
    77 !      ddbx et ddby termes dus au glissement  
    78 !      relation avec la vitesse de glissement ubx et uby  
    79 !      ubx=-ddbx*sdx et uby=-ddby*sdy 
    80 !      ------------------------------------------------- 
    81  
    82 ! Calcul de la vitesse de glissement, qu'elle soit calculee par diagno ou par sliding 
    83 !------------------------------------------------------------------------------------ 
    84 !$OMP PARALLEL 
    85 !$OMP WORKSHARE 
    86 where (flgzmx(:,:)) 
    87     ubx(:,:) = uxflgz(:,:) 
    88 elsewhere 
    89     ubx(:,:) = ddbx(:,:)* (-sdx(:,:))     ! on pourrait mettre une limitation 
    90 endwhere 
    91  
    92 where (flgzmy(:,:)) 
    93     uby(:,:) = uyflgz(:,:) 
    94 elsewhere 
    95     uby(:,:) = ddby(:,:)* (-sdy(:,:)) 
    96 endwhere 
    97 !$OMP END WORKSHARE 
    98 !$OMP END PARALLEL 
    99 if (itracebug.eq.1) call tracebug(' diffusiv  apres calcul glissement') 
    100  
    101  
    102  
    103 ! Deformation SIA : Calcul de ddy  et ddx pour chaque valeur de iglen 
    104 !------------------------------------------------------------------ 
    105 do iglen=n1poly,n2poly 
    106    glenexp=max(0.,(glen(iglen)-1.)/2.) 
    107   !$OMP PARALLEL  
     14subroutine diffusiv 
     15 
     16  !     =============================================================== 
     17  !     **     DIFFUSIVITIES                    18 mars 2006   *** 
     18  !     glissement avec reserve d'eau 29 Avril 97 
     19  !     choix sur le type de discretisation (locale iloc=1 )     
     20  ! 
     21  !     Modification Vince --> 2 fev 95 
     22  !     pour couplage avec Shelf 
     23  ! 
     24  !     Modification C. Dumas Dec 99 
     25  !     DDX, SA,S2A,... ont une dimension en plus (loi de deformation) 
     26  ! 
     27  !     Modifications  de Cat du 11 juin 2000 
     28  !     Modif Christophe 24 janv 2002 : passage au fortran 90 
     29  !     pour optimisation de la routine 
     30  !     Modif Vincent dec 2003 : on utilise plus frontmx/y > commentaires 
     31  ! 
     32  !     Modif Cat mars 2003 
     33  !                  modif mars 2007 : moyennes de S2a avec S2a_mx et S2a_my 
     34 
     35  !     =============================================================== 
     36  !$ USE OMP_LIB 
     37  use module3D_phy, only:nx,ny,dx,dy,V_limit,iglen,diffmx,diffmy,flgzmx,flgzmy,ubx,uby,& 
     38       uxflgz,uyflgz,ddbx,ddby,sdx,sdy,flotmx,flotmy,slope2mx,slope2my,hmx,hmy,uxdef,uydef,& 
     39       uxbar,uybar,H,flot 
     40  use runparam, only: itracebug 
     41  use param_phy_mod, only:rog 
     42  use deformation_mod_2lois, only:n1poly,n2poly,glen,ddx,ddy,s2a_mx,s2a_my 
     43  use module_choix, only:sliding 
     44 
     45  implicit none 
     46 
     47  integer :: i,j 
     48  real :: ulgliss,ultot,uldef,glenexp 
     49  real :: INV_4DX ! inverse de dx pour eviter les divisions =1/(4*dx) 
     50  real :: INV_4DY ! inverse de dy pour eviter les divisions =1/(4*dy) 
     51 
     52  if (itracebug.eq.1)  call tracebug(' Entree dans routine diffusiv') 
     53 
     54  !     la valeur de ff est donnee dans initial 
     55  !     limitation de la vitesse de glissement de déformation et totale 
     56  ulgliss = V_limit 
     57  ultot   = V_limit 
     58  uldef   = 2000. 
     59 
     60 
     61  INV_4DX=1/(4*dx) 
     62  INV_4DY=1/(4*dy) 
     63 
     64  ! initialisation de la diffusion 
     65  !$OMP PARALLEL 
     66  !$OMP WORKSHARE 
     67  diffmx(:,:)=0. 
     68  diffmy(:,:)=0. 
     69  !$OMP END WORKSHARE 
     70  !$OMP END PARALLEL 
     71 
     72  ! calcul de SDX, SDY et de la pente au carre en mx et en my : 
     73 
     74  call slope_surf 
     75 
     76  call sliding     ! au sens vitesse de glissement 
     77 
     78  !  le glissement est maintenant dans un module a part choisi dans le module choix 
     79  !  pour l'instant seules les lois Heino (loigliss=4) et Bindshadler(loigliss=2) 
     80  !  sont programmees. 
     81 
     82  !      ddbx et ddby termes dus au glissement  
     83  !      relation avec la vitesse de glissement ubx et uby  
     84  !      ubx=-ddbx*sdx et uby=-ddby*sdy 
     85  !      ------------------------------------------------- 
     86 
     87  ! Calcul de la vitesse de glissement, qu'elle soit calculee par diagno ou par sliding 
     88  !------------------------------------------------------------------------------------ 
     89  !$OMP PARALLEL 
     90  !$OMP WORKSHARE 
     91  where (flgzmx(:,:)) 
     92     ubx(:,:) = uxflgz(:,:) 
     93  elsewhere 
     94     ubx(:,:) = ddbx(:,:)* (-sdx(:,:))     ! on pourrait mettre une limitation 
     95  endwhere 
     96 
     97  where (flgzmy(:,:)) 
     98     uby(:,:) = uyflgz(:,:) 
     99  elsewhere 
     100     uby(:,:) = ddby(:,:)* (-sdy(:,:)) 
     101  endwhere 
     102  !$OMP END WORKSHARE 
     103  !$OMP END PARALLEL 
     104  if (itracebug.eq.1) call tracebug(' diffusiv  apres calcul glissement') 
     105 
     106 
     107 
     108  ! Deformation SIA : Calcul de ddy  et ddx pour chaque valeur de iglen 
     109  !------------------------------------------------------------------ 
     110  do iglen=n1poly,n2poly 
     111     glenexp=max(0.,(glen(iglen)-1.)/2.) 
     112     !$OMP PARALLEL  
     113     !$OMP DO 
     114     do j=1,ny 
     115        do i=1,nx 
     116           if (.not.flotmy(i,j)) then !  On calcule quand la deformation même pour les ice streams 
     117              ddy(i,j,iglen)=((slope2my(i,j)**  &  ! slope2my calcule en debut de routine 
     118                   glenexp)*(rog)**glen(iglen)) & 
     119                   *hmy(i,j)**(glen(iglen)+1) 
     120           endif 
     121        end do 
     122     end do 
     123     !$OMP END DO 
     124     !$OMP END PARALLEL 
     125  end do 
     126 
     127 
     128 
     129  do iglen=n1poly,n2poly 
     130     glenexp=max(0.,(glen(iglen)-1.)/2.) 
     131     !$OMP PARALLEL 
     132     !$OMP DO 
     133     do j=1,ny   
     134        do i=1,nx 
     135           if (.not.flotmx(i,j)) then      ! bug y->x corrige le 5 avril 2008 ( 
     136              ddx(i,j,iglen)=((slope2mx(i,j)**  & ! slope2mx calcule en debut de routine 
     137                   glenexp)*(rog)**glen(iglen)) & 
     138                   *hmx(i,j)**(glen(iglen)+1) 
     139           endif 
     140        end do 
     141     end do 
     142     !$OMP END DO 
     143     !$OMP END PARALLEL 
     144  end do 
     145 
     146 
     147  ! vitesse due a la déformation----------------------------------------------- 
     148  !$OMP PARALLEL 
    108149  !$OMP DO 
    109    do j=1,ny 
    110       do i=1,nx 
    111          if (.not.flotmy(i,j)) then !  On calcule quand la deformation même pour les ice streams 
    112             ddy(i,j,iglen)=((slope2my(i,j)**  &  ! slope2my calcule en debut de routine 
    113                  glenexp)*(rog)**glen(iglen)) & 
    114                  *hmy(i,j)**(glen(iglen)+1) 
    115          endif 
    116       end do 
    117    end do 
    118    !$OMP END DO 
    119    !$OMP END PARALLEL 
    120 end do 
    121  
    122  
    123  
    124 do iglen=n1poly,n2poly 
    125    glenexp=max(0.,(glen(iglen)-1.)/2.) 
    126    !$OMP PARALLEL 
    127    !$OMP DO 
    128    do j=1,ny   
    129       do i=1,nx 
    130          if (.not.flotmx(i,j)) then      ! bug y->x corrige le 5 avril 2008 ( 
    131             ddx(i,j,iglen)=((slope2mx(i,j)**  & ! slope2mx calcule en debut de routine 
    132                  glenexp)*(rog)**glen(iglen)) & 
    133                  *hmx(i,j)**(glen(iglen)+1) 
    134          endif 
    135       end do 
    136    end do 
    137    !$OMP END DO 
    138    !$OMP END PARALLEL 
    139 end do 
    140  
    141  
    142 ! vitesse due a la déformation----------------------------------------------- 
    143 !$OMP PARALLEL 
    144 !$OMP DO 
    145 ydef: do j=2,ny 
    146        do i=1,nx  
     150  ydef: do j=2,ny 
     151     do i=1,nx  
    147152 
    148153        if (.not.flotmy(i,j)) then !  On calcule la deformation même pour les ice streams 
    149154 
    150          uydef(i,j)=0. 
    151          diffmy(i,j)=0. 
    152  
    153          do iglen=n1poly,n2poly 
    154            diffmy(i,j)=diffmy(i,j)+ddy(i,j,iglen)*(s2a_my(i,j,1,iglen)) 
    155          end do 
    156  
    157          uydef(i,j)=diffmy(i,j)*(-sdy(i,j))    ! partie deformation 
    158          uydef(i,j)=min(uydef(i,j),uldef) 
    159          uydef(i,j)=max(uydef(i,j),-uldef) 
    160  
    161          diffmy(i,j)=diffmy(i,j)+ddby(i,j)     ! pour utilisation comme diffusion, a multiplier par H 
    162  
    163          uybar(i,j)=uby(i,j)+uydef(i,j) 
    164  
    165  
    166          else    ! points flottants 
    167             uydef(i,j)=0.   ! normalement uxbed est attribue           
    168          endif 
    169     end do 
    170    end do ydef 
    171 !$OMP END DO 
    172  
    173 !$OMP DO 
    174 xdef: do j=1,ny 
    175        do i=2,nx  
     155           uydef(i,j)=0. 
     156           diffmy(i,j)=0. 
     157 
     158           do iglen=n1poly,n2poly 
     159              diffmy(i,j)=diffmy(i,j)+ddy(i,j,iglen)*(s2a_my(i,j,1,iglen)) 
     160           end do 
     161 
     162           uydef(i,j)=diffmy(i,j)*(-sdy(i,j))    ! partie deformation 
     163           uydef(i,j)=min(uydef(i,j),uldef) 
     164           uydef(i,j)=max(uydef(i,j),-uldef) 
     165 
     166           diffmy(i,j)=diffmy(i,j)+ddby(i,j)     ! pour utilisation comme diffusion, a multiplier par H 
     167 
     168           uybar(i,j)=uby(i,j)+uydef(i,j) 
     169 
     170 
     171        else    ! points flottants 
     172           uydef(i,j)=0.   ! normalement uxbed est attribue           
     173        endif 
     174     end do 
     175  end do ydef 
     176  !$OMP END DO 
     177 
     178  !$OMP DO 
     179  xdef: do j=1,ny 
     180     do i=2,nx  
    176181 
    177182        if (.not.flotmx(i,j)) then !  On calcule la deformation même pour les ice streams 
    178183 
    179          uxdef(i,j)=0. 
    180          diffmx(i,j)=0. 
    181  
    182          do iglen=n1poly,n2poly 
    183            diffmx(i,j)=diffmx(i,j)+ddx(i,j,iglen)*(s2a_mx(i,j,1,iglen)) 
    184          end do 
    185  
    186          uxdef(i,j)=diffmx(i,j)*(-sdx(i,j))    ! partie deformation 
    187          uxdef(i,j)=min(uxdef(i,j),uldef) 
    188          uxdef(i,j)=max(uxdef(i,j),-uldef) 
    189  
    190  
    191          diffmx(i,j)=diffmx(i,j)+ddbx(i,j)     ! pour utilisation comme diffusion, 
    192                                                ! a multiplier par H 
    193          uxbar(i,j)=ubx(i,j)+uxdef(i,j) 
    194  
    195  
    196          else    ! points flottants 
    197             uxdef(i,j)=0.   ! normalement uxbed est attribue           
    198          endif 
    199     end do 
    200    end do xdef 
    201 !$OMP END DO 
    202  
    203 if (itracebug.eq.1) call tracebug(' diffusiv  avant limit') 
    204  
    205  
    206 ! limitation par ultot---------------------------------------------------------- 
    207  
    208 ! la limitation selon x 
    209 !$OMP WORKSHARE 
    210 where(.not.flgzmx(:,:)) 
    211    uxbar(:,:)=min(uxbar(:,:),ultot) 
    212    uxbar(:,:)=max(uxbar(:,:),-ultot) 
    213 end where 
    214  
    215 ! la limitation selon y 
    216 where(.not.flgzmy(:,:)) 
    217    uybar(:,:)=min(uybar(:,:),ultot) 
    218    uybar(:,:)=max(uybar(:,:),-ultot) 
    219 end where 
    220 !$OMP END WORKSHARE 
    221  
    222 !    cas particulier des sommets ou il ne reste plus de neige. 
    223 !$OMP DO 
    224 do j=2,ny-1 
    225    do i=2,nx-1 
    226       if ((H(i,j).le.1.).and.(.not.flot(i,j))) then 
    227          uxbar(i+1,j)=min(uxbar(i+1,j),ultot)  ! n'agit que si ux ->  
    228          uxbar(i,j)=max(uxbar(i,j),-ultot)     ! n'agit que si ux <- 
    229          uybar(i,j+1)=min(uybar(i,j+1),ultot) 
    230          uybar(i,j)=max(uybar(i,j),-ultot) 
    231       endif 
    232    end do 
    233 end do 
    234 !$OMP END DO 
    235 !$OMP END PARALLEL 
    236 if (itracebug.eq.1)  call tracebug(' fin diffusiv ') 
    237  
    238 return 
     184           uxdef(i,j)=0. 
     185           diffmx(i,j)=0. 
     186 
     187           do iglen=n1poly,n2poly 
     188              diffmx(i,j)=diffmx(i,j)+ddx(i,j,iglen)*(s2a_mx(i,j,1,iglen)) 
     189           end do 
     190 
     191           uxdef(i,j)=diffmx(i,j)*(-sdx(i,j))    ! partie deformation 
     192           uxdef(i,j)=min(uxdef(i,j),uldef) 
     193           uxdef(i,j)=max(uxdef(i,j),-uldef) 
     194 
     195 
     196           diffmx(i,j)=diffmx(i,j)+ddbx(i,j)     ! pour utilisation comme diffusion, 
     197           ! a multiplier par H 
     198           uxbar(i,j)=ubx(i,j)+uxdef(i,j) 
     199 
     200 
     201        else    ! points flottants 
     202           uxdef(i,j)=0.   ! normalement uxbed est attribue           
     203        endif 
     204     end do 
     205  end do xdef 
     206  !$OMP END DO 
     207 
     208  if (itracebug.eq.1) call tracebug(' diffusiv  avant limit') 
     209 
     210 
     211  ! limitation par ultot---------------------------------------------------------- 
     212 
     213  ! la limitation selon x 
     214  !$OMP WORKSHARE 
     215  where(.not.flgzmx(:,:)) 
     216     uxbar(:,:)=min(uxbar(:,:),ultot) 
     217     uxbar(:,:)=max(uxbar(:,:),-ultot) 
     218  end where 
     219 
     220  ! la limitation selon y 
     221  where(.not.flgzmy(:,:)) 
     222     uybar(:,:)=min(uybar(:,:),ultot) 
     223     uybar(:,:)=max(uybar(:,:),-ultot) 
     224  end where 
     225  !$OMP END WORKSHARE 
     226 
     227  !    cas particulier des sommets ou il ne reste plus de neige. 
     228  !$OMP DO 
     229  do j=2,ny-1 
     230     do i=2,nx-1 
     231        if ((H(i,j).le.1.).and.(.not.flot(i,j))) then 
     232           uxbar(i+1,j)=min(uxbar(i+1,j),ultot)  ! n'agit que si ux ->  
     233           uxbar(i,j)=max(uxbar(i,j),-ultot)     ! n'agit que si ux <- 
     234           uybar(i,j+1)=min(uybar(i,j+1),ultot) 
     235           uybar(i,j)=max(uybar(i,j),-ultot) 
     236        endif 
     237     end do 
     238  end do 
     239  !$OMP END DO 
     240  !$OMP END PARALLEL 
     241  if (itracebug.eq.1)  call tracebug(' fin diffusiv ') 
     242 
     243  return 
    239244end subroutine diffusiv 
    240245 
Note: See TracChangeset for help on using the changeset viewer.