Changeset 73


Ignore:
Timestamp:
06/22/16 15:43:40 (8 years ago)
Author:
dumas
Message:

OpenMP parallelization in Temperature-routines, deformation_mod_2lois and resol_adv_diff_2D-sept2009

Location:
trunk/SOURCES
Files:
7 edited

Legend:

Unmodified
Added
Removed

  • trunk/SOURCES/Temperature-routines/advec_icetemp.f90

    r24 r73  
    1414!> 
    1515 
    16 Subroutine  Advec_icetemp(Ii,Jj,Prodq_m,Cro_m,Advecx_m,Advecy_m,Advec_m,Ct_m) 
     16Subroutine  Advec_icetemp(Nz,Nzm,Nn,Prodq_m,Cro_m,Advecx_m,Advecy_m,Advec_m,Ct_m,Iadvec_w,Iadvec_e,Iadvec_s,Iadvec_n, & 
     17        Uxij,Uxi_plus_un,Uyij,Uyj_plus_un,Tij,Ti_moins_un,Ti_plus_un,Tj_moins_un,Tj_plus_un,Hij,Tsij,Aa,Bb,Cc,Rr, & 
     18        Dx11,Dou,DTT,Dee,Uzrij,Dttdx) 
    1719 
    18   use icetemp_declar  
     20!  use icetemp_declar  
    1921 
    2022  Implicit None 
    2123    !< Arguments 
    22   Integer,                     Intent(In)  ::  Ii,Jj          !< Indice Selon X Et Y 
    23   Real, Dimension(Nz),         Intent(In)   :: Prodq_m        !< Tableau 1d Vert. heat production 
    24   Real, Dimension(Nz),         Intent(In)   :: Cro_m          !< Tableau 1d Vert. heat capcity 
    25   Real, Dimension(Nz),         Intent(Inout):: Advecx_m       !< Advection selon x 
    26   Real, Dimension(Nz),         Intent(Inout):: Advecy_m       !< Advection selon y 
    27   Real, Dimension(Nz),         Intent(Inout):: Advec_m        !< Advection total 
    28   Real, Dimension(Nz),         Intent(In)   :: Ct_m           !< Tableau 1d Vert.  thermal cond. 
     24  Integer,                     intent(in)  ::  Nz,Nzm,Nn      !< Taille des tableaux 
     25  Real, Dimension(Nz),         intent(in)   :: Prodq_m        !< Tableau 1d Vert. heat production 
     26  Real, Dimension(Nz),         intent(in)   :: Cro_m          !< Tableau 1d Vert. heat capcity 
     27  Real, Dimension(Nz),         intent(inout):: Advecx_m       !< Advection selon x 
     28  Real, Dimension(Nz),         intent(inout):: Advecy_m       !< Advection selon y 
     29  Real, Dimension(Nz),         intent(inout):: Advec_m        !< Advection total 
     30  Real, Dimension(Nz),         intent(in)   :: Ct_m           !< Tableau 1d Vert.  thermal cond. 
     31   
     32  Integer, intent(in) :: Iadvec_w,Iadvec_e,Iadvec_s,Iadvec_n 
     33  real, dimension(Nz), Intent(in) :: Uxij, Uxi_plus_un 
     34  real, dimension(Nz), Intent(in) :: Uyij, Uyj_plus_un 
     35  real, dimension(Nz+Nzm), Intent(inout) :: Tij 
     36  real, dimension(Nz+Nzm), Intent(in) :: Ti_moins_un,Ti_plus_un,Tj_moins_un,Tj_plus_un 
     37  real, intent(in) :: Hij 
     38  real,intent(in) :: Tsij 
     39  Real,Dimension(Nn),intent(inout) :: Aa    !< Work Arrays For Tridiag  !Dim Nn 
     40  Real,Dimension(Nn),intent(inout) :: Bb    !< Work Arrays For Tridiag  !Dim Nn 
     41  Real,Dimension(Nn),intent(inout) :: Cc    !< Work Arrays For Tridiag  !Dim Nn 
     42  Real,Dimension(Nn),intent(inout) :: Rr    !< Work Arrays For Tridiag  !Dim Nn 
     43  Real, intent(in) :: Dx11 
     44  Real, intent(inout) :: Dou 
     45  real, intent(in) :: DTT 
     46  real, intent(in) :: Dee 
     47  real, dimension(Nz) :: Uzrij 
     48  real :: Dttdx 
     49   
     50 
     51   
     52  ! variables locales 
     53  Integer :: K 
     54  Real :: Dzz,Dah,Ct_haut,Ct_bas 
    2955 
    3056  
    … …  
    3561     !         Avection Selon X  (advecx_m en general > 0) 
    3662     !         ---------------- 
    37      Advecx_m(K) = Iadvec_w(Ii,Jj) * Ux(Ii,Jj,K) *          & ! ux west if upwind 
    38           (T(Ii,Jj,K) - T(Ii-1,Jj,K))  & ! west T gradient 
     63     Advecx_m(K) = Iadvec_w * Uxij(K) *          & ! ux west if upwind 
     64          (Tij(K) - Ti_moins_un(K))  & ! west T gradient 
    3965             
    40           + Iadvec_e(Ii+1,Jj)* Ux(Ii+1,Jj,K)*      & ! ux east if upwind 
    41           (T(Ii+1,Jj,K)-T(Ii,Jj,K))     ! east T gradient 
     66          + Iadvec_e * Uxi_plus_un(K)*      & ! ux east if upwind 
     67          (Ti_plus_un(K)-Tij(K))     ! east T gradient 
    4268 
    4369     Advecx_m(K) =  Advecx_m(K) * Dx11                                        !Dx11=1/Dx 
    … …  
    4672     !         Avection Selon Y 
    4773     !         ---------------- 
    48      Advecy_m(K) =   Iadvec_s(Ii,Jj) * Uy(Ii,Jj,K) *        & ! uy sud if upwind 
    49           (T(Ii,Jj,K) - T(Ii,Jj-1,K))& ! south T gradient 
     74     Advecy_m(K) =   Iadvec_s * Uyij(K) *        & ! uy sud if upwind 
     75          (Tij(K) - Tj_moins_un(K))& ! south T gradient 
    5076           
    51           + Iadvec_n(Ii,Jj+1) * Uy(Ii,Jj+1,K) *   & ! uy nord is upwind 
    52           (T(Ii,Jj+1,K)-T(Ii,Jj,K))    ! north T gradient 
     77          + Iadvec_n * Uyj_plus_un(K) *   & ! uy nord is upwind 
     78          (Tj_plus_un(K)-Tij(K))    ! north T gradient 
    5379 
    5480     Advecy_m(K) =  Advecy_m(K) * Dx11 
    … …  
    6288 
    6389  !      -----------------------------Cas General (H>10m) 
    64   thick_ice:  If (H(Ii,Jj).Gt.10.) Then 
     90  thick_ice:  If (Hij.Gt.10.) Then 
    6591 
    6692     !        Variables De Calcul Dans La Glace 
    6793     ! dou = dtt/dz^2 
    68      Dou=Dtt/Dee/Dee/ H(Ii,Jj) / H(Ii,Jj)   
     94     Dou=Dtt/Dee/Dee/ Hij / Hij   
    6995 
    7096     ! Dah = dtt/dz 
    71      Dah=Dtt /Dee /H(Ii,Jj) 
     97     Dah=Dtt /Dee /Hij 
    7298 
    7399     ! thermal conductivity at mid point just below the surface                    
    … …  
    75101 
    76102     ! surface temperature : cannot go above 0 celsius 
    77      T(Ii,Jj,1)=Min(0.,Ts(Ii,Jj)) 
     103     Tij(1)=Min(0.,Tsij) 
    78104 
    79105     ! surface boundary condition 
    … …  
    81107     Bb(1)=1. 
    82108     Cc(1)=0. 
    83      Rr(1)=T(Ii,Jj,1) 
     109     Rr(1)=Tij(1) 
    84110 
    85111     Do K=2,Nz-1 
    … …  
    91117 
    92118        ! Advection Verticale Centree 
    93         Aa(K) = -Dzz * Ct_haut - Uzr(Ii,Jj,K) * Dah / 2.  ! lower diag 
    94         Cc(K) = -Dzz * Ct_bas  + Uzr(Ii,Jj,K) * Dah / 2.  ! upper diag 
     119        Aa(K) = -Dzz * Ct_haut - Uzrij(K) * Dah / 2.  ! lower diag 
     120        Cc(K) = -Dzz * Ct_bas  + Uzrij(K) * Dah / 2.  ! upper diag 
    95121        Bb(K) =  1.+ Dzz * (Ct_haut+Ct_bas)                                ! diag 
    96122         
    97         Rr(K) =  T(Ii,Jj,K) + Prodq_m(K) * Dtt   ! vector 
     123        Rr(K) =  Tij(K) + Prodq_m(K) * Dtt   ! vector 
    98124        Rr(K) =  Rr(K)- Dttdx * (Advec_m(K)) 
    99125 

  • trunk/SOURCES/Temperature-routines/icetemp_declar_mod.f90

    r71 r73  
    2222  Integer :: Nfracq   !< Exposant Fracq 
    2323  Integer :: Iq       !< Choix Du Type De Routine Chaleur 
    24   Real :: Sx,Sy,Sx2,Sy2,Deh22,Tss 
    25   Real ::  Dou,Dah,Duu,Dzz,Dzi,Chalbed,Ct_bas,Ct_haut 
     24  Real :: Sx,Sy,Sx2,Sy2,Deh22 ! ,Tss 
     25  Real ::  Dou,Duu,Chalbed ! ,Dah,Dzz,Dzi,Ct_bas,Ct_haut 
    2626 
    2727  Real,Parameter :: Acof1=-0.0575  !< Pour La Temperature De L'Eau De Mer 
    … …  
    4040  ! _______________ 
    4141 
    42   Real,allocatable,Dimension(:) :: Aa    !< Work Arrays For Tridiag  !Dim Nn 
    43   Real,allocatable,Dimension(:) :: Bb    !< Work Arrays For Tridiag  !Dim Nn 
    44   Real,allocatable,Dimension(:) :: Cc    !< Work Arrays For Tridiag  !Dim Nn 
    45   Real,allocatable,Dimension(:) :: Rr    !< Work Arrays For Tridiag  !Dim Nn 
    46   Real,allocatable,Dimension(:) :: Hh    !< Work Arrays For Tridiag  !Dim Nn 
     42  Real,Dimension(Nn) :: Aa    !< Work Arrays For Tridiag  !Dim Nn 
     43  Real,Dimension(Nn) :: Bb    !< Work Arrays For Tridiag  !Dim Nn 
     44  Real,Dimension(Nn) :: Cc    !< Work Arrays For Tridiag  !Dim Nn 
     45  Real,Dimension(Nn) :: Rr    !< Work Arrays For Tridiag  !Dim Nn 
     46  Real,Dimension(Nn) :: Hh    !< Work Arrays For Tridiag  !Dim Nn 
    4747 
    48   Real,allocatable,Dimension(:) :: Tdot   !< Temperature Time Derivative*Dtt   
    49   Real,allocatable,Dimension(:) :: Abis,Bbis,Cbis,Rbis,Hbis                    
    50   Real,allocatable,Dimension(:) :: Ee        !< Vertical Coordinate In Ice, Scaled To H Zeta   
     48  Real,Dimension(Nz+Nzm) :: Tdot   !< Temperature Time Derivative*Dtt   
     49!  Real,allocatable,Dimension(:) :: Abis,Bbis,Cbis,Rbis,Hbis                    
     50  Real,Dimension(Nz) :: Ee        !< Vertical Coordinate In Ice, Scaled To H Zeta   
    5151 
    5252  ! Tableaux De Travail 2d  
    5353  ! ___________________________ 
    54   Real,allocatable,Dimension(:,:):: Tbmer    !< Temperature De La Mer A La Base De L'Ice Shelf   
    55   Real,allocatable,Dimension(:,:) :: Chalglissx,Chalglissy !< Chaleur De Glissement              
    56   Integer,allocatable,Dimension(:,:) :: Iadvec_w,Iadvec_e,Iadvec_s,Iadvec_n    
     54  Real,Dimension(Nx,Ny):: Tbmer    !< Temperature De La Mer A La Base De L'Ice Shelf   
     55  Real,Dimension(Nx,Ny) :: Chalglissx,Chalglissy !< Chaleur De Glissement              
     56  Integer,Dimension(Nx,Ny) :: Iadvec_w,Iadvec_e,Iadvec_s,Iadvec_n    
    5757  ! Pour L'Advection         
    5858 
    5959  ! Tableaux De Travail 3d 
    6060  ! __________________________ 
    61   Real,allocatable,Dimension(:,:,:,:) :: Chalx, Chaly         !Dim Nx,Ny,Nz,N1poly:N2poly 
     61  Real,Dimension(Nx,Ny,Nz,size(Btt,4)) :: Chalx, Chaly         !Dim Nx,Ny,Nz,N1poly:N2poly 
    6262  ! Utilise Pour Calculer La Chaleur De Deformation Selon Xx Yy Zz Et Xy 
    63   Real,allocatable,Dimension(:,:,:) :: Ffx,Ffy                                        
    64   Real,allocatable,Dimension(:,:,:) :: T3d_new                       
    65   Real,allocatable,Dimension(:,:,:) :: Chal2_x, Chal2_y, Chal2_z, Chal2_xy 
    66   Real,allocatable,Dimension(:,:,:) :: Chaldef_maj   !< Chaleur De Deformation                    
    67   Real,allocatable,Dimension(:,:,:) :: Advecx,Advecy,Advec 
     63  Real,Dimension(Nx,Ny,size(Btt,4)) :: Ffx,Ffy                                        
     64  Real,Dimension(Nx,Ny,Nz+Nzm) :: T3d_new                       
     65  Real,Dimension(Nx,Ny,Nz) :: Chal2_x, Chal2_y, Chal2_z, Chal2_xy 
     66  Real,Dimension(Nx,Ny,Nz) :: Chaldef_maj   !< Chaleur De Deformation                    
     67  Real,Dimension(Nx,Ny,Nz) :: Advecx,Advecy,Advec 
    6868  Real,Dimension(Nx,Ny,Nz) :: Cp            !< Specific Heat Capacity (J/(M-3)/K)=Ro Cp   
    6969  Real,Dimension(Nx,Ny,Nz) :: Ct            !< Thermal Conductivity (J/M/K/A)             

  • trunk/SOURCES/Temperature-routines/icetemp_mod.f90

    r29 r73  
    6666 
    6767  subroutine icetemp 
    68  
     68    !$ use OMP_LIB  
    6969    use Icetemp_declar 
    7070 
    7171 
    7272    Implicit None  
    73  
    74     ! Allocation de la memoire pour les tableaux dans icetemp_declar 
    75  
    76     If (.Not. Allocated(Aa)) then  
    77        ! tab Nn 
    78        Allocate (Aa(Nn),Bb(Nn),Cc(Nn),Rr(Nn),Hh(Nn))   
    79        ! tab NzNzm_m 
    80        Allocate(Tdot(Nz+Nzm))    
    81        ! tab Nzm+1 
    82        Allocate(Abis(Nzm+1),Bbis(Nzm+1),Cbis(Nzm+1),Rbis(Nzm+1),Hbis(Nzm+1))                   
    83        ! tab nz_m 
    84        Allocate(Ee(Nz))         
    85        !Tab Nx,Ny 
    86        Allocate(Tbmer(Nx,Ny),&    
    87             Chalglissx(Nx,Ny),& 
    88             Chalglissy(Nx,Ny),& 
    89             Iadvec_w(Nx,Ny),& 
    90             Iadvec_e(Nx,Ny),& 
    91             Iadvec_s(Nx,Ny),& 
    92             Iadvec_n(Nx,Ny) )       
    93        ! Tab Nx,Ny,Nz 
    94        Allocate(Chal2_x(Nx,Ny,Nz),& 
    95             Chal2_y(Nx,Ny,Nz),& 
    96             Chal2_z(Nx,Ny,Nz),& 
    97             Chal2_xy(Nx,Ny,Nz),& 
    98             Chaldef_maj(Nx,Ny,Nz),& 
    99             Advecx(Nx,Ny,Nz),& 
    100             Advecy(Nx,Ny,Nz),& 
    101             Advec(Nx,Ny,Nz)) 
    102        ! Tab Nx,Ny,Nz,N1poly:N2poly 
    103        Allocate(Chalx(Nx,Ny,Nz,size(Btt,4)),& 
    104             Chaly(Nx,Ny,Nz,size(Btt,4)))          
    105        ! Tab Nx,Ny,N1poly:N2poly 
    106        Allocate(Ffx( Nx,Ny,size(Btt,4)),& 
    107             Ffy( Nx,Ny,size(Btt,4)))                                         
    108        ! Tab Nx,Ny,NzNzm 
    109        Allocate( T3d_new (Nx,Ny,Nz+Nzm))                   
    110     end If 
    111  
     73     
     74!~   integer :: t1,t2,ir 
     75!~   real                           :: temps, t_cpu_0, t_cpu_1, t_cpu, norme   
     76!~    
     77!~    ! Temps CPU de calcul initial. 
     78!~    call cpu_time(t_cpu_0) 
     79!~    ! Temps elapsed de reference. 
     80!~    call system_clock(count=t1, count_rate=ir) 
     81     
     82     
     83         !$OMP PARALLEL 
     84         !$OMP WORKSHARE 
    11285    ! init 
    11386    Aa=0.     
    … …  
    12497    Chalglissx=0. 
    12598    Chalglissy=0. 
     99     
    126100    !Instruction  
    127101    Dee=1./(Nz-1)  
    … …  
    130104    Ee(1)=0. 
    131105    Ee(Nz)=1. 
     106    !$OMP END WORKSHARE 
     107     
     108    !$OMP DO 
    132109    Do K=1,Nz 
    133110       If ((K.Ne.1).And.(K.Ne.Nz)) Ee(K)=(K-1.)/(Nz-1.) 
    134111    End Do 
     112    !$OMP END DO 
     113    !$OMP END PARALLEL 
    135114 
    136115    !   Variables Dependant Du Pas De Temps Dtt 
    … …  
    156135    Call Qprod_ice(Iq) 
    157136 
    158  
     137   
    159138    ! Dans Le Socle : Les elements de La Matrice Tridiag Sont Invariants Dans L'Espace 
     139    !$OMP PARALLEL 
     140    !$OMP DO 
    160141    Do K=Nz+1,nz+nzm-1 
    161142       Aa(K)=-Ctm 
    … …  
    163144       Cc(K)=-Ctm 
    164145    End Do 
     146         !$OMP END DO 
    165147 
    166148    ! Conditions Aux Limites Pour Le Socle 
    … …  
    174156    !         Advection Selon X 
    175157    !         ----------------- 
    176  
     158    !$OMP WORKSHARE 
    177159    Iadvec_w(:,:)=Nint(Max(Sign(1.0,Uxbar(:,:)),0.0)) 
    178160    Iadvec_e(:,:)=Nint(Abs(Min(Sign(1.0,Uxbar(:,:)),0.0))) 
    … …  
    183165    Iadvec_s(:,:)=Nint(Max(Sign(1.0,Uybar(:,:)),0.0))  
    184166    Iadvec_n(:,:)=Nint(Abs(Min(Sign(1.0,Uybar(:,:)),0.0))) 
    185  
     167    !$OMP END WORKSHARE 
    186168    ! Boucle De Resolution Des Temperatures Colonne Par Colonne 
    187  
     169         !$OMP SINGLE 
    188170    If (Itracebug.Eq.1) Write(num_tracebug,*)' avant appel  boucle Advec_icetemp' 
    189  
     171    !$OMP END SINGLE 
     172 
     173         !$OMP DO PRIVATE(Deh22,Duu,Aa,Bb,Cc,Rr,Hh,Dou) 
    190174    Do J=2,Ny-1 
    191175       Do I=2,Nx-1 
    192176 
    193177 !  If (Itracebug.Eq.1) Write(num_tracebug,*)' avant appel  Advect_icetemp  i,j', i,j 
    194           call Advec_icetemp(I,J,Chaldef_maj(I,J,:),Cp(I,J,:),&  
    195                Advecx(I,J,:),Advecy(I,J,:),Advec(I,J,:),Ct(I,J,:)) 
     178          call Advec_icetemp(Nz,Nzm,Nn,Chaldef_maj(I,J,:),Cp(I,J,:),&  
     179               Advecx(I,J,:),Advecy(I,J,:),Advec(I,J,:),Ct(I,J,:),Iadvec_w(I,J),Iadvec_e(I+1,J),Iadvec_s(I,J),Iadvec_n(I,J+1), & 
     180                                        Ux(i,j,:),Ux(i+1,j,:),Uy(i,j,:),Uy(i,j+1,:),T(i,j,:),T(i-1,j,:),T(i+1,j,:),T(i,j-1,:),T(i,j+1,:), & 
     181                                        H(i,j),Ts(i,j),Aa,Bb,Cc,Rr,Dx11,Dou,DTT,Dee,Uzr(i,j,:),Dttdx) 
    196182 
    197183          ! variables de calcul dans la glace 
    … …  
    201187 
    202188!   If (Itracebug.Eq.1) Write(num_tracebug,*)' avant appel  temp_col  i,j', i,j 
    203           Call temp_col(I,J,T3d_new,Ct(I,J,:),& 
    204                flot(I,J),H(I,J),Tbmer(I,J)) 
     189          Call temp_col(Nz,Nzm,Nn,T3d_new(i,j,:),Ct(i,j,:),& 
     190               flot(i,j),H(i,j),Tbmer(i,j),Ibase(i,j),T(i,j,:),Tpmp(i,j,:),Aa,Bb,Cc,Rr,Hh,Ncond,Dee,Cm,& 
     191               Dzm,Phid(i,j),Ghf(i,j),Ifail1,Tbdot(i,j),DTT,Dou,Ts(i,j),Bmelt(i,j)) 
    205192 
    206193          ! Test Permettant D'Ecrire Les Variables Quand Il Y A Une Erreur 
    … …  
    218205       End Do 
    219206    End Do 
    220  
     207         !$OMP END DO 
     208 
     209         !$OMP DO 
    221210    Do J=1,Ny 
    222211       Do I=1,Nx 
    … …  
    226215       End Do 
    227216    End Do 
     217    !$OMP END DO 
    228218 
    229219    ! Attribution Finale De La Temperature 
    230220    ! Limitation a 3 deg De Variation Par Dtt Et Pas Plus Froid Que -70 
    231  
     221         !$OMP DO COLLAPSE(2) PRIVATE(Tdot) 
    232222    Do K=1,nz+nzm 
    233223       Do J=1,Ny 
    … …  
    250240       End Do 
    251241    End Do 
    252  
    253  
     242    !$OMP END DO 
     243 
     244         !$OMP DO COLLAPSE(2) 
    254245    Do K=2,Nz 
    255246       Do J=1,Ny 
    … …  
    262253       End Do 
    263254    End Do 
     255    !$OMP END DO 
    264256 
    265257    !     Temperature Sur Les Bords : 
    266258    !     ------------------------- 
    267259   
     260    !$OMP DO COLLAPSE(2) 
    268261    Do K=1,Nz 
    269262       Do J=1,Ny 
    … …  
    281274       End Do 
    282275    End Do 
    283     
    284  
     276    !$OMP END DO 
     277 
     278    !$OMP DO COLLAPSE(2)  
    285279    Do K=1,Nz 
    286280       Do J=1,Ny 
    … …  
    294288       End Do 
    295289    End Do 
    296  
    297     Deallocate (Aa,Bb,Cc,Rr,Hh,Tdot,Abis,Bbis,Cbis,Rbis,Hbis,Ee,Tbmer,            &  
    298                 Chalglissx,Chalglissy,Iadvec_w,Iadvec_e,                          &  
    299                 Iadvec_s,Iadvec_n,Chal2_x,Chal2_y,Chal2_z,Chal2_xy,Chaldef_maj,   & 
    300                 Advecx,Advecy,Advec,Chalx,Chaly,Ffx,Ffy,T3d_new)         
     290    !$OMP END DO 
     291    !$OMP END PARALLEL 
     292  
     293!~   ! Temps elapsed final 
     294!~   call system_clock(count=t2, count_rate=ir) 
     295!~   temps=real(t2 - t1,kind=4)/real(ir,kind=4) 
     296!~   ! Temps CPU de calcul final 
     297!~   call cpu_time(t_cpu_1) 
     298!~   t_cpu = t_cpu_1 - t_cpu_0 
     299   
     300!~   ! Impression du resultat. 
     301!~   print '(//,3X,"Valeurs de nx et ny : ",I5,I5/,              & 
     302!~            & 3X,"Temps elapsed       : ",1PE10.3," sec.",/, & 
     303!~            & 3X,"Temps CPU           : ",1PE10.3," sec.",/, & 
     304!~            & 3X,"Norme (PB si /= 0)  : ",1PE10.3,//)', & 
     305!~            nx,ny,temps,t_cpu,norme 
     306 
    301307 
    302308    If (Itracebug.Eq.1) Write(num_tracebug,*)' fin routine icetemp' 

  • trunk/SOURCES/Temperature-routines/temp_col.f90

    r24 r73  
    1414!> 
    1515 
    16 Subroutine  Temp_col(Ii,Jj,Newtempcol_m,Ct_m,Flot_m,H_m,Tbmer_m) 
    17  
    18   Use Icetemp_declar 
     16Subroutine  Temp_col(Nz,Nzm,Nn,Newtempcol,Ctij,Flotij,Hij,Tbmerij,Ibase,Tij,Tpmpij,Aa,Bb,Cc,Rr,Hh,Ncond,Dee,Cm, & 
     17                                Dzm,Phidij,Ghfij,Ifail1,Tbdotij,DTT,Dou,Tsij,Bmeltij) 
     18 
     19!  Use Icetemp_declar 
    1920  Use Tridiagmod 
    2021   
    2122   
    2223  Implicit None 
    23   integer :: Ii,Jj 
    24   Real,Dimension(Nx,Ny,Nz+Nzm),Intent(Out)  :: Newtempcol_m        !< Tableau 1d Vert. Des Temperatures En Sortie 
    25   Real,Dimension(Nz),Intent(In)             :: Ct_m                !< Tableau 1d Vert. Conductivite Thermique 
    26   Logical,Intent(In)                        :: Flot_m              !< Vrai Si Flot_mtant (Test D'Archimede) 'O' 
    27   Real,Intent(In)                           :: H_m                 !< Ice Thickness  'O' 
    28   Real,Intent(In)                           :: Tbmer_m             !< Temperature De La Mer A La Base De L'Ice Shelf 
    29  
    30  
    31   If (H_m.Gt.10.) Then 
     24  integer,intent(in) :: Nz,Nzm,Nn 
     25  Real,Dimension(Nz+Nzm),Intent(Out)  :: Newtempcol        !< Tableau 1d Vert. Des Temperatures En Sortie 
     26  Real,Dimension(Nz),Intent(In)       :: Ctij                !< Tableau 1d Vert. Conductivite Thermique 
     27  Logical,Intent(In)                  :: Flotij              !< Vrai Si Flotijtant (Test D'Archimede) 'O' 
     28  Real,Intent(In)                     :: Hij                 !< Ice Thickness  'O' 
     29  Real,Intent(In)                     :: Tbmerij             !< Temperature De La Mer A La Base De L'Ice Shelf 
     30  integer,intent(inout)               :: Ibase 
     31  real,dimension(Nz+Nzm),intent(inout) :: Tij             !< Temperature sur la colonne 
     32  real,dimension(Nz),intent(in)       :: Tpmpij               !< pressure melting point temperature in ice sheet 'o' 
     33  Real,Dimension(Nn),intent(inout) :: Aa    !< Work Arrays For Tridiag  !Dim Nn 
     34  Real,Dimension(Nn),intent(inout) :: Bb    !< Work Arrays For Tridiag  !Dim Nn 
     35  Real,Dimension(Nn),intent(inout) :: Cc    !< Work Arrays For Tridiag  !Dim Nn 
     36  Real,Dimension(Nn),intent(inout) :: Rr    !< Work Arrays For Tridiag  !Dim Nn 
     37  Real,Dimension(Nn),intent(inout) :: Hh    !< Work Arrays For Tridiag  !Dim Nn 
     38  Integer,intent(in)               :: Ncond !< Switch : Ncond=1 Thermal Conductivity In Mantle 
     39  real, intent(in) :: Dee 
     40  real,intent(in) :: Cm 
     41  real,intent(in) :: Dzm 
     42  real,intent(inout) :: Phidij     !< flux de chaleur lie a la deformation et glissement basal 
     43  real,intent(in) :: Ghfij      !< flux geothermique 
     44  integer,intent(out) :: Ifail1 !< Diagnostique erreur dans Tridiag 
     45  real, intent(out) :: Tbdotij  !< variation in time of basal temperature 
     46  real, intent(in) :: DTT 
     47  Real, intent(inout) :: Dou 
     48  real, intent(in) :: Tsij        !< surface ice temperature  'o' 
     49  real, intent(out) :: Bmeltij    !< fonte basale 
     50   
     51   
     52   
     53  ! variables locales 
     54  integer :: K 
     55  Real :: Dzi,Tss 
     56  Real,Dimension(Nzm+1) :: Abis,Bbis,Cbis,Rbis,Hbis    
     57 
     58  If (Hij.Gt.10.) Then 
    3259 
    3360     ! Conditions Aux Limites a La Base De La Glace 
    … …  
    3562 
    3663     ! Base Froide  
    37      If ( ((Ibase(Ii,Jj).Eq.1).Or.(Ibase(Ii,Jj).Eq.4) & 
    38           .Or. ((Ibase(Ii,Jj).Eq.5).And.(T(Ii,Jj,Nz).Lt.Tpmp(Ii,Jj,Nz)))) & 
    39           .And..Not.Flot_m) Then 
    40  
    41         Ibase(Ii,Jj)=1 
     64     If ( ((Ibase.Eq.1).Or.(Ibase.Eq.4) & 
     65          .Or. ((Ibase.Eq.5).And.(Tij(Nz).Lt.Tpmpij(Nz)))) & 
     66          .And..Not.Flotij) Then 
     67 
     68        Ibase=1 
    4269        Bb(Nz)=1. 
    4370 
    4471        If (Ncond.Eq.1) Then ! Avec Socle 
    45            Dzi=H_m*Dee*Cm/Ct_m(Nz) 
     72           Dzi=Hij*Dee*Cm/Ctij(Nz) 
    4673           Aa(Nz)=-Dzm/(Dzm+Dzi) 
    4774           Cc(Nz)=-Dzi/(Dzm+Dzi) 
    48            Rr(Nz)=H_m*Dee*Phid(Ii,Jj)*Dzm/Ct_m(Nz)/(Dzm+Dzi) 
     75           Rr(Nz)=Hij*Dee*Phidij*Dzm/Ctij(Nz)/(Dzm+Dzi) 
    4976        Else                 ! Sans Socle 
    5077           Aa(Nz)=-1. 
    5178           Cc(Nz)=0. 
    52            Rr(Nz)=-(Ghf(Ii,Jj)-Phid(Ii,Jj))/Ct_m(Nz)*H_m*Dee 
     79           Rr(Nz)=-(Ghfij-Phidij)/Ctij(Nz)*Hij*Dee 
    5380        Endif 
    5481 
    … …  
    5784        ! Base Temperee Ou Shelf 
    5885        ! ----------------------                
    59         If (Ibase(Ii,Jj).Eq.5)  Ibase(Ii,Jj)=2 
    60         Ibase(Ii,Jj)=Max(Ibase(Ii,Jj),2) 
     86        If (Ibase.Eq.5)  Ibase=2 
     87        Ibase=Max(Ibase,2) 
    6188        Aa(Nz)=0. 
    6289        Bb(Nz)=1. 
    6390        Cc(Nz)=0. 
    6491 
    65         If (.Not.Flot_m) Then 
    66            Rr(Nz)=Tpmp(Ii,Jj,Nz) 
    67         Else 
    68            Rr(Nz)=Tbmer_m 
     92        If (.Not.Flotij) Then 
     93           Rr(Nz)=Tpmpij(Nz) 
     94        Else 
     95           Rr(Nz)=Tbmerij 
    6996        End If 
    7097 
    … …  
    74101     If (Ncond.Eq.1) Then ! Avec Socle 
    75102        Do K=Nz+1,nz+nzm-1 
    76            Rr(K)=T(Ii,Jj,K) 
     103           Rr(K)=Tij(K) 
    77104        End Do 
    78105        Call Tridiag (Aa,Bb,Cc,Rr,Hh,nz+nzm,Ifail1) 
    … …  
    87114     If (Ncond.Eq.0) Then 
    88115        Do K=Nz+1,nz+nzm 
    89            Hh(K)=Hh(Nz)-Dzm*(K-Nz)*Ghf(Ii,Jj)/Cm 
     116           Hh(K)=Hh(Nz)-Dzm*(K-Nz)*Ghfij/Cm 
    90117        End Do 
    91118     Endif 
    92119 
    93120     Do K=1,nz+nzm 
    94         Newtempcol_m(Ii,Jj,K)=Hh(K) 
     121        Newtempcol(K)=Hh(K) 
    95122     End Do 
    96      Tbdot(Ii,Jj)=(Newtempcol_m(Ii,Jj,Nz)-T(Ii,Jj,Nz))/Dtt 
     123     Tbdotij=(Newtempcol(Nz)-Tij(Nz))/Dtt 
    97124 
    98125     !      ------------------------------------------------- Glace Tres Fine (H<10m Ou H=0) 
    99   Else If (H_m.Le.10.) Then 
     126  Else If (Hij.Le.10.) Then 
    100127 
    101128     !       Pour Eviter Des Problemes Lorsque H Est Inferieur A 10 M 
    … …  
    104131 
    105132     ! ........................................ Avec Calcul Dans Le Socle 
    106      If (Ncond.Eq.1.And..Not.Flot_m) Then    
    107         If (H_m.Gt.0.) Then 
     133     If (Ncond.Eq.1.And..Not.Flotij) Then    
     134        If (Hij.Gt.0.) Then 
    108135           ! Gradient Dans Le Socle 
    109            Dou=(T(Ii,Jj,Nz+1)-T(Ii,Jj,Nz))/Dzm*Cm      
    110            Dou=Dou/Ct_m(Nz)*Dee*H_m 
    111  
    112            Tss=Min(0.,Ts(Ii,Jj)) 
    113            Do K=1,Nz 
    114               Newtempcol_m(Ii,Jj,K)=Tss+Dou*(K-1.) 
     136           Dou=(Tij(Nz+1)-Tij(Nz))/Dzm*Cm      
     137           Dou=Dou/Ctij(Nz)*Dee*Hij 
     138 
     139           Tss=Min(0.,Tsij) 
     140           Do K=1,Nz 
     141              Newtempcol(K)=Tss+Dou*(K-1.) 
    115142           End Do 
    116143 
    117144        Else 
    118145           ! Pas De Glace 
    119            Tss=Ts(Ii,Jj) 
    120            Do K=1,Nz 
    121               Newtempcol_m(Ii,Jj,K)=Tss 
     146           Tss=Tsij 
     147           Do K=1,Nz 
     148              Newtempcol(K)=Tss 
    122149           End Do 
    123150        End If 
    124151 
    125152        ! Calcul Dans Le Socle Meme S'Il N'Y A Pas De Glace 
    126         Rr(Nz)=Newtempcol_m(Ii,Jj,Nz) 
     153        Rr(Nz)=Newtempcol(Nz) 
    127154        Aa(Nz)=0 
    128155        Bb(Nz)=1 
    … …  
    130157 
    131158        Do K=Nz+1,nz+nzm-1 
    132            Rr(K)=T(Ii,Jj,K) 
     159           Rr(K)=Tij(K) 
    133160        End Do 
    134161 
    … …  
    150177     Else 
    151178        ! Calotte Posee                     
    152         If ((H_m.Gt.0.).And..Not.Flot_m) Then 
    153            Dou=-Ghf(Ii,Jj)/Ct_m(Nz)*Dee*H_m 
    154            Tss=Min(0.,Ts(Ii,Jj)) 
    155            Do K=1,Nz 
    156               Newtempcol_m(Ii,Jj,K)=Tss+Dou*(K-1.) 
     179        If ((Hij.Gt.0.).And..Not.Flotij) Then 
     180           Dou=-Ghfij/Ctij(Nz)*Dee*Hij 
     181           Tss=Min(0.,Tsij) 
     182           Do K=1,Nz 
     183              Newtempcol(K)=Tss+Dou*(K-1.) 
    157184           End Do 
    158185 
    159186           !  Shelf 
    160         Else If ((H_m.Gt.0.).And.Flot_m) Then 
    161            Tss=Min(0.,Ts(Ii,Jj)) 
    162            Dou=(Tbmer_m-Tss)*Dee 
    163            Do K=1,Nz 
    164               Newtempcol_m(Ii,Jj,K)=Tss+Dou*(K-1.) 
     187        Else If ((Hij.Gt.0.).And.Flotij) Then 
     188           Tss=Min(0.,Tsij) 
     189           Dou=(Tbmerij-Tss)*Dee 
     190           Do K=1,Nz 
     191              Newtempcol(K)=Tss+Dou*(K-1.) 
    165192           End Do 
    166193 
    167194           ! Pas De Glace 
    168195        Else 
    169            Tss=Ts(Ii,Jj)       
    170            Do K=1,Nz 
    171               Newtempcol_m(Ii,Jj,K)=Tss 
     196           Tss=Tsij       
     197           Do K=1,Nz 
     198              Newtempcol(K)=Tss 
    172199           End Do 
    173200        End If 
    174201 
    175202        ! Temperature Dans Le Socle, Lineaire Avec Le Gradient Ghf 
    176         If (.Not.Flot_m) Then 
     203        If (.Not.Flotij) Then 
    177204           Do K=Nz+1,nz+nzm 
    178               Hh(K)=Newtempcol_m(Ii,Jj,Nz)-Dzm*(K-Nz)*Ghf(Ii,Jj)/Cm 
     205              Hh(K)=Newtempcol(Nz)-Dzm*(K-Nz)*Ghfij/Cm 
    179206           End Do 
    180207        Else 
    181208           Do K=Nz+1,nz+nzm 
    182               Hh(K)=Tbmer_m-Dzm*(K-Nz)*Ghf(Ii,Jj)/Cm 
     209              Hh(K)=Tbmerij-Dzm*(K-Nz)*Ghfij/Cm 
    183210           End Do 
    184211        End If 
    … …  
    189216 
    190217     Do K=Nz+1,nz+nzm 
    191         Newtempcol_m(Ii,Jj,K)=Hh(K) 
     218        Newtempcol(K)=Hh(K) 
    192219     End Do 
    193220 
    194      Tbdot(Ii,Jj)=(Newtempcol_m(Ii,Jj,Nz)-T(Ii,Jj,Nz))/Dtt 
    195      Bmelt(Ii,Jj)=0. 
    196      Ibase(Ii,Jj)=5 
    197      Phid(Ii,Jj)=0. 
     221     Tbdotij=(Newtempcol(Nz)-Tij(Nz))/Dtt 
     222     Bmeltij=0. 
     223     Ibase=5 
     224     Phidij=0. 
    198225 
    199226 

  • trunk/SOURCES/deformation_mod_2lois.f90

    r71 r73  
    468468!debug_3d(:,:,90+iiglen-1) = s2a(:,:,1,iiglen) 
    469469 
    470  
    471 !~   ! Temps elapsed final 
    472 !~   call system_clock(count=t2, count_rate=ir) 
    473 !~   temps=real(t2 - t1,kind=4)/real(ir,kind=4) 
    474 !~   ! Temps CPU de calcul final 
    475 !~   call cpu_time(t_cpu_1) 
    476 !~   t_cpu = t_cpu_1 - t_cpu_0 
    477   
    478 !~   ! Impression du resultat. 
    479 !~   print '(//,3X,"Valeurs de nx et ny : ",I5,I5/,              & 
    480 !~            & 3X,"Temps elapsed       : ",1PE10.3," sec.",/, & 
    481 !~            & 3X,"Temps CPU           : ",1PE10.3," sec.",/, & 
    482 !~            & 3X,"Norme (PB si /= 0)  : ",1PE10.3,//)', & 
    483 !~            nx,ny,temps,t_cpu,norme 
    484  
    485470end subroutine flowlaw 
    486471 

  • trunk/SOURCES/flottab2-0.7.f90

    r72 r73  
    105105 
    106106 
    107 !~   integer :: t1,t2,ir 
    108 !~   real                           :: temps, t_cpu_0, t_cpu_1, t_cpu, norme   
    109    
    110 !~    ! Temps CPU de calcul initial. 
    111 !~    call cpu_time(t_cpu_0) 
    112 !~    ! Temps elapsed de reference. 
    113 !~    call system_clock(count=t1, count_rate=ir) 
    114  
    115107 
    116108if (itracebug.eq.1)  call tracebug(' Entree dans routine flottab') 
    … …  
    725717!fin de routine flottab2 
    726718!print*, 'front',front(50,30),ice(50,30),flotmx(i,j),uxbar(i,j) 
    727  
    728 !~   ! Temps elapsed final 
    729 !~   call system_clock(count=t2, count_rate=ir) 
    730 !~   temps=real(t2 - t1,kind=4)/real(ir,kind=4) 
    731 !~   ! Temps CPU de calcul final 
    732 !~   call cpu_time(t_cpu_1) 
    733 !~   t_cpu = t_cpu_1 - t_cpu_0 
    734   
    735 !~   ! Impression du resultat. 
    736 !~   print '(//,3X,"Valeurs de nx et ny : ",I5,I5/,              & 
    737 !~            & 3X,"Temps elapsed       : ",1PE10.3," sec.",/, & 
    738 !~            & 3X,"Temps CPU           : ",1PE10.3," sec.",/, & 
    739 !~            & 3X,"Norme (PB si /= 0)  : ",1PE10.3,//)', & 
    740 !~            nx,ny,temps,t_cpu,norme 
    741719 
    742720end subroutine flottab 

  • trunk/SOURCES/resol_adv_diff_2D-sept2009.f90

    r65 r73  
    101101subroutine resol_adv_diff_2D_vect(Dfx,Dfy,advx,advy,H_presc,i_Hpresc,bil,vieuxH,newH) 
    102102 
     103!$ USE OMP_LIB 
    103104implicit none 
    104105 
    … …  
    148149if (itracebug.eq.1)  call tracebug(' Entree dans routine resolution_diffusion') 
    149150 
    150  
     151!$OMP PARALLEL 
     152!$OMP WORKSHARE 
    151153! calcul de bdx et hdx 
    152154hdx(:,:)=dx1*(vieuxH(:,:)-eoshift(vieuxH(:,:),shift=-1,boundary=0.,dim=1)) 
    … …  
    164166Frelax(:,:)=0. 
    165167DeltaH(:,:)=0. 
    166  
     168!$OMP END WORKSHARE 
    167169 
    168170! schema spatial 
    169171 
    170172if (upwind.eq.1.) then                 !schema amont 
    171  
     173!$OMP WORKSHARE 
    172174   where (Advx(:,:).ge.0.) 
    173175      c_west(:,:)=1. 
    … …  
    185187      c_north(:,:)=1. 
    186188   end where 
    187  
     189!$OMP END WORKSHARE 
    188190else if (upwind.lt.1.) then             ! schema centre 
     191!$OMP WORKSHARE 
    189192      c_west(:,:)=0.5 
    190193      c_east(:,:)=0.5 
    191194      c_south(:,:)=0.5 
    192195      c_north(:,:)=0.5 
     196!$OMP END WORKSHARE       
    193197end if 
    194198 
     199 
    195200! attribution des elements des diagonales 
     201!$OMP DO PRIVATE(frdx,frdy,fraxw,fraxe,frays,frayn) 
    196202do j=2,ny-1 
    197203  do i=2,nx-1 
    … …  
    254260     end do 
    255261  end do 
    256  
    257  
     262!$OMP END DO 
     263 
     264!$OMP WORKSHARE 
    258265where (i_hpresc(:,:) .eq.1)          ! thickness prescribed 
    259266   frelax(:,:) = H_presc(:,:) 
    … …  
    264271   erelax(:,:) = 0. 
    265272end where 
     273!$OMP END WORKSHARE 
     274!$OMP END PARALLEL 
    266275 
    267276stopp = .false. 
    268277ntour=0 
    269278 
    270  
    271  
    272279relax_loop: do while(.not.stopp) 
    273280ntour=ntour+1 
    274  
    275  
     281!$OMP PARALLEL 
     282!$OMP DO PRIVATE(reste) 
    276283do j=2,ny-1 
    277284   do i=2,nx-1 
    … …  
    281288           + crelax(i,j)*newH(i,j))- frelax(i,j) 
    282289 
    283 if (ntour.eq.1)  debug_3D(i,j,49)=(((arelax(i,j)*newH(i-1,j) +drelax(i,j)*newH(i,j-1)) & 
    284            + (brelax(i,j)*newH(i+1,j) + erelax(i,j)*newH(i,j+1))) & 
    285            + crelax(i,j)*newH(i,j)) 
     290!if (ntour.eq.1)  debug_3D(i,j,49)=(((arelax(i,j)*newH(i-1,j) +drelax(i,j)*newH(i,j-1)) & 
     291!           + (brelax(i,j)*newH(i+1,j) + erelax(i,j)*newH(i,j+1))) & 
     292!           + crelax(i,j)*newH(i,j)) 
    286293 
    287294 
    … …  
    290297   end do 
    291298end do 
    292  
    293 debug_3D(:,:,50)=arelax(:,:) 
    294 debug_3D(:,:,51)=brelax(:,:) 
    295 debug_3D(:,:,52)=crelax(:,:) 
    296 debug_3D(:,:,53)=drelax(:,:) 
    297 debug_3D(:,:,54)=erelax(:,:) 
    298 debug_3D(:,:,55)=frelax(:,:) 
    299  
    300  
    301  
     299!$OMP END DO 
     300 
     301!debug_3D(:,:,50)=arelax(:,:) 
     302!debug_3D(:,:,51)=brelax(:,:) 
     303!debug_3D(:,:,52)=crelax(:,:) 
     304!debug_3D(:,:,53)=drelax(:,:) 
     305!debug_3D(:,:,54)=erelax(:,:) 
     306!debug_3D(:,:,55)=frelax(:,:) 
     307 
     308 
     309!$OMP WORKSHARE 
    302310newH(:,:)=newH(:,:)-deltaH(:,:) 
    303  
    304  
     311!$OMP END WORKSHARE 
     312!$OMP END PARALLEL 
    305313 
    306314 
    … …  
    310318delh=0 
    311319 
    312  
     320!$OMP PARALLEL 
     321!$OMP DO REDUCTION(+:delh) 
    313322do j=2,ny-1 
    314323   do i=2,nx-1 
    … …  
    316325   end do 
    317326end do 
     327!$OMP END DO 
     328!$OMP END PARALLEL 
    318329 
    319330if (delh.gt.0.) then 
    … …  
    330341 
    331342! thickness at the upwind node 
    332 do j = 1, ny 
    333    do i = 2, nx 
    334       debug_3D(i,j,92) = c_west(i,j) * newH(i-1,j) + c_east(i,j) * newH(i,j) 
    335    end do 
    336 end do 
    337 do j = 2, ny 
    338    do i = 1, nx 
    339       debug_3D(i,j,93) = c_south(i,j) * newH(i,j-1) + c_north(i,j) * newH(i,j) 
    340    end do 
    341 end do 
    342  
    343  if (itracebug.eq.1)  call tracebug(' fin routine resolution_diffusion') 
     343!do j = 1, ny 
     344!   do i = 2, nx 
     345!      debug_3D(i,j,92) = c_west(i,j) * newH(i-1,j) + c_east(i,j) * newH(i,j) 
     346!   end do 
     347!end do 
     348!do j = 2, ny 
     349!   do i = 1, nx 
     350!      debug_3D(i,j,93) = c_south(i,j) * newH(i,j-1) + c_north(i,j) * newH(i,j) 
     351!   end do 
     352!end do 
     353 
     354if (itracebug.eq.1)  call tracebug(' fin routine resolution_diffusion') 
     355 
     356  
    344357return 
    345358 
Note: See TracChangeset for help on using the changeset viewer.