Changeset 482 for branches


Ignore:
Timestamp:
02/09/24 17:08:19 (4 months ago)
Author:
aquiquet
Message:

Cleaning branch: some routines with explicit arguments in New-remplimat

Location:
branches/GRISLIv3/SOURCES/New-remplimat
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • branches/GRISLIv3/SOURCES/New-remplimat/diagno-L2_mod.f90

    r481 r482  
    1616real, dimension(nx,ny) :: uyb1ramollo 
    1717real, dimension(nx,ny) :: pvi_keep 
    18  
    19 !cdc transfere dans module3d pour compatibilite avec furst_schoof_mod 
    20 !cdc  integer, dimension(nx,ny) :: imx_diag 
    21 !cdc  integer, dimension(nx,ny) :: imy_diag 
    2218 
    2319integer :: nneigh=max(int(400000./dx),1) !nb. points around the grline kept for back force comp. in case we reduce the extent 
     
    151147!$OMP END WORKSHARE 
    152148!$OMP END PARALLEL 
    153   ! avec couplage thermomecanique 
    154 !  write(166,*) ' apres call calc_pvi',itour_pvi 
    155149 
    156150else 
     
    158152end if 
    159153 
    160   call imx_imy_nx_ny         ! pour rempli_L2 : calcule les masques imx et imy qui  
     154  call imx_imy_nx_ny(imx_diag,imy_diag,flgzmx,flgzmy)         ! pour rempli_L2 : calcule les masques imx et imy qui  
    161155 
    162156  diagno_grline=sum(gr_line(2:(nx-1),2:(ny-1))) 
    163157  if ((Schoof.ge.1).and.(diagno_grline.gt.0)) then 
    164        call imx_imy_nx_ny_reduce(1) !afq: provides imx_reduce and imy_reduce 
     158       call imx_imy_nx_ny_reduce(1,flot,imx_diag,imy_diag,gr_line) !afq: provides imx_reduce and imy_reduce 
    165159  endif 
    166  
    167 !cdc debug Schoof !!!!!!!!!!!!   
    168 !~   do j=1,ny 
    169 !~              do i=1,nx 
    170 !~                      write(578,*) uxbar(i,j) 
    171 !~                      write(579,*) uybar(i,j) 
    172 !~              enddo 
    173 !~      enddo    
    174    
    175   !if (Schoof.eq.1.and.nt.GT.15000) then ! flux grounding line Schoof 
    176 ! afq -- below:  if (Schoof.eq.1) then ! flux grounding line Schoof 
    177 ! afq -- below:     call interpol_glflux ! calcul flux GL + interpolation sur voisins 
    178 ! afq -- below:  endif 
    179  
    180 !~       do j=1,ny 
    181 !~              do i=1,nx 
    182 !~                      write(588,*) uxbar(i,j) 
    183 !~                      write(589,*) uybar(i,j) 
    184 !~              enddo 
    185 !~      enddo    
    186 !~      print*,'ecriteure termineee !!!!!!' 
    187 !~      read(*,*) 
    188160 
    189161  ! donnent les cas de conditions aux limites 
     
    202174  nyd1=1 
    203175  nyd2=ny 
    204  
    205   !call rempli_L2(1,nx,1,ny,uxbar,uybar,uxb1,uyb1,imx_diag,imy_diag,ifail_diagno)     
    206   !nxd1=15 
    207   !nxd2=19 
    208   !nyd1=30 
    209   !nyd2=34 
    210  
    211   !nxd1=35 
    212   !nxd2=60 
    213   !nyd1=35 
    214   !nyd2=60 
    215  
    216176        
    217   !if (Schoof.eq.1.and.nt.GT.15000) then ! flux grounding line Schoof avec calcul de la back force par shelf ramollo 
    218177  if ((Schoof.ge.1).and.(diagno_grline.gt.0)) then ! flux grounding line Schoof avec calcul de la back force par shelf ramollo 
    219178 
     
    224183    pvi_keep(:,:)=pvi(:,:)        
    225184    where (flot(:,:).and.H(:,:).GT.2.) 
    226 !       pvi(:,:)=1.e5     
    227185       pvi(:,:)=pvimin 
    228186    endwhere 
     
    250208 
    251209    if (ifail_diagno_ramollo.gt.0) then 
    252 !       write(6,*) ' Probleme resolution systeme L2. ramollo ifail=',ifail_diagno_ramollo 
    253 !       STOP 
    254210       write(*,*) ' Probleme resolution systeme L2. ramollo ifail=',ifail_diagno_ramollo 
    255211       write(*,*) '          ... we go on anyway!' 
    256212    endif 
    257 !~   do j=1,ny 
    258 !~              do i=1,nx 
    259 !~                if (sqrt(uxb1(i,j)**2+ uyb1(i,j)**2).gt.0..and..not.flot(i,j)) then 
    260 !~                              write(1034,*) sqrt(uxb1(i,j)**2+ uyb1(i,j)**2) / sqrt(uxb1ramollo(i,j)**2 + uyb1ramollo(i,j)**2) 
    261 !~                else 
    262 !~                              write(1034,*) 1. 
    263 !~                      endif 
    264 !~              enddo            
    265 !~      enddo 
    266  
    267 !~   print*,'apres calcul rempli_L2' 
    268 !~   read(*,*) 
    269213     
    270214    call interpol_glflux ! calcul flux GL + interpolation sur voisins 
     
    503447debug_3D(:,:,27)=pvi(:,:) 
    504448!$OMP END WORKSHARE 
    505 ! attention run 35 
    506 !--------------------  
    507 !!$if (time.gt.10.) then 
    508 !!$   where (flot(:,:)) 
    509 !!$      pvi(:,:)=pvimin 
    510 !!$   end where 
    511 !!$end if 
    512449 
    513450!  calcul de la viscosite integree au milieu des mailles (pvm) 
     
    528465!------------------------------------------------------------------ 
    529466 
    530 subroutine imx_imy_nx_ny 
    531  
    532 use module3d_phy, only: imx_diag,imy_diag,flgzmx,flgzmy  
    533           
     467subroutine imx_imy_nx_ny(imx_diagl,imy_diagl,flgzmxl,flgzmyl) 
     468 
    534469implicit none 
     470 
     471logical,dimension(:,:), intent(in)    :: flgzmxl,flgzmyl 
     472integer,dimension(:,:), intent(inout) :: imx_diagl,imy_diagl 
    535473 
    536474! definition des masques 
     
    555493!$OMP PARALLEL 
    556494!$OMP WORKSHARE 
    557 imx_diag(:,:)=0 
    558 imy_diag(:,:)=0 
     495imx_diagl(:,:)=0 
     496imy_diagl(:,:)=0 
    559497 
    560498! a l'interieur du domaine 
    561499!------------------------- 
    562500 
    563 where (flgzmx(:,:)) 
    564    imx_diag(:,:)=2             ! shelf ou stream         
     501where (flgzmxl(:,:)) 
     502   imx_diagl(:,:)=2             ! shelf ou stream         
    565503elsewhere 
    566    imx_diag(:,:)=1             ! vitesse imposee 
     504   imx_diagl(:,:)=1             ! vitesse imposee 
    567505end where 
    568506 
    569 where (flgzmy(:,:))       
    570    imy_diag(:,:)=2             ! shelf ou stream 
     507where (flgzmyl(:,:))       
     508   imy_diagl(:,:)=2             ! shelf ou stream 
    571509elsewhere 
    572    imy_diag(:,:)=1             ! vitesse imposee 
     510   imy_diagl(:,:)=1             ! vitesse imposee 
    573511end where 
    574512 
    575513! bord sud  
    576 imx_diag(:,1)=-1 
    577 imy_diag(:,2)=-1 
     514imx_diagl(:,1)=-1 
     515imy_diagl(:,2)=-1 
    578516 
    579517! bord nord 
    580 imx_diag(:,ny)=-3 
    581 imy_diag(:,ny)=-3 
     518imx_diagl(:,ny)=-3 
     519imy_diagl(:,ny)=-3 
    582520 
    583521! bord Est 
    584 imx_diag(1,:)=0    ! hors domaine a cause des mailles alternees 
    585 imx_diag(2,:)=-4 
    586 imy_diag(1,:)=-4 
     522imx_diagl(1,:)=0    ! hors domaine a cause des mailles alternees 
     523imx_diagl(2,:)=-4 
     524imy_diagl(1,:)=-4 
    587525 
    588526! bord West 
    589 imx_diag(nx,:)=-2 
    590 imy_diag(nx,:)=-2 
     527imx_diagl(nx,:)=-2 
     528imy_diagl(nx,:)=-2 
    591529 
    592530! Coins 
    593 imx_diag(2,1)=-41       ! SW 
    594 imy_diag(1,2)=-41 
    595  
    596 imx_diag(nx,1)=-12      ! SE 
    597 imy_diag(nx,2)=-12 
    598  
    599 imx_diag(nx,ny)=-23     ! NE 
    600 imy_diag(nx,ny)=-23 
    601  
    602 imx_diag(2,ny)=-34      ! NW 
    603 imy_diag(1,ny)=-34 
     531imx_diagl(2,1)=-41       ! SW 
     532imy_diagl(1,2)=-41 
     533 
     534imx_diagl(nx,1)=-12      ! SE 
     535imy_diagl(nx,2)=-12 
     536 
     537imx_diagl(nx,ny)=-23     ! NE 
     538imy_diagl(nx,ny)=-23 
     539 
     540imx_diagl(2,ny)=-34      ! NW 
     541imy_diagl(1,ny)=-34 
    604542 
    605543! hors domaine 
    606 imx_diag(1,:)=0     ! hors domaine a cause des mailles alternees 
    607 imy_diag(:,1)=0     ! hors domaine a cause des mailles alternees 
     544imx_diagl(1,:)=0     ! hors domaine a cause des mailles alternees 
     545imy_diagl(:,1)=0     ! hors domaine a cause des mailles alternees 
    608546!$OMP END WORKSHARE 
    609547!$OMP END PARALLEL 
     
    612550 
    613551 
    614 subroutine imx_imy_nx_ny_reduce(choix) 
    615  
    616 use module3d_phy, only: flot,imx_diag,imy_diag,gr_line  
    617           
     552subroutine imx_imy_nx_ny_reduce(choix,flot,imx_diag,imy_diag,gr_line) 
     553 
    618554implicit none 
    619555 
     
    621557!       We simply compute the velocities around (nneigh) the grounding line 
    622558 
    623   integer, intent(in) :: choix 
     559  integer,                 intent(in)    :: choix 
     560  logical, dimension(:,:), intent(in)    :: flot 
     561  integer, dimension(:,:), intent(in)    :: imx_diag,imy_diag,gr_line 
     562 
    624563 
    625564  integer i,j,nvx,nvy 
     
    717656end subroutine mismip_boundary_cond 
    718657 
    719 !end module diagno_L2_mod 
    720658end module diagno_mod 
  • branches/GRISLIv3/SOURCES/New-remplimat/eq_ellipt_sgbsv_mod-0.2.f90

    r478 r482  
    5757 
    5858!---------------------------- 
    59 subroutine  redim_matrice 
     59subroutine  redim_matrice(mmatl) 
    6060 
    6161! allocation dynamique de mmat 
     
    6363implicit none 
    6464 
     65real,dimension(:,:), allocatable, intent(inout) :: mmatl 
     66 
    6567integer ::  err 
    6668 
    67 if (allocated(mmat)) then  
    68    deallocate(mmat,stat=err) 
     69if (allocated(mmatl)) then  
     70   deallocate(mmatl,stat=err) 
    6971     if (err/=0) then 
    7072        print *,"Erreur à la desallocation de mmat",err 
     
    7476 
    7577 
    76 allocate(mmat(ldb,nblig),stat=err) 
     78allocate(mmatl(ldb,nblig),stat=err) 
    7779if (err/=0) then 
    7880   print *,"erreur a l'allocation du tableau mmat",err 
     
    8082   stop  
    8183end if 
    82  
    83  
    84 !write(6,*) "allocation mmat ok, ldb,nblig",ldb,nblig 
    8584 
    8685end subroutine redim_matrice 
     
    364363! si la matrice mmat  change de taille on la re allocate. 
    365364 
    366 if ((ldb.ne.ldb_old).and.(nblig.ne.nblig_old)) call redim_matrice 
     365if ((ldb.ne.ldb_old).and.(nblig.ne.nblig_old)) call redim_matrice(mmat) 
    367366 
    368367 
  • branches/GRISLIv3/SOURCES/New-remplimat/remplimat-shelves-tabTu.f90

    r481 r482  
    131131!$OMP END PARALLEL 
    132132 
    133 ! calcul de Hmx_oppos et Hmy_oppos dans le cas de fronts 
    134 ! ce calcul pourrait être dans diagno et faire passer hmx_oppos par module 
    135  
    136 ! 4 avril 2012 : Finalement il me semble que c'est faux, il suffit de garder le Hmx habituel. 
    137 ! cat 
    138  
    139 !!$do j=ny1,ny2 
    140 !!$   do i=nx1,nx2 
    141 !!$ 
    142 !!$      il=max(1,i-1) 
    143 !!$      jl=max(1,j-1) 
    144 !!$ 
    145 !!$ 
    146 !!$         if ((H(i,j).gt.1.).and.(H(il,j).le.1.)) then   ! Bord West : on garde H(i,j) 
    147 !!$            Hmx_oppos(i,j)=H(i,j) 
    148 !!$ 
    149 !!$         else if ((H(i,j).le.1.).and.(H(il,j).gt.1.)) then   ! Bord Est : on garde H(il,j) 
    150 !!$            Hmx_oppos(i,j)=H(il,j) 
    151 !!$ 
    152 !!$         end if 
    153 !!$ 
    154 !!$         if ((H(i,j).gt.1.).and.(H(i,jl).le.1.)) then   ! Bord Sud : on garde H(i,j) 
    155 !!$            Hmy_oppos(i,j)=H(i,j) 
    156 !!$ 
    157 !!$         else if ((H(i,j).le.1.).and.(H(i,jl).gt.1.)) then   ! Bord Nord on garde H(i,jl) 
    158 !!$            Hmy_oppos(i,j)=H(i,jl) 
    159 !!$ 
    160 !!$         end if 
    161 !!$ 
    162 !!$ 
    163 !!$      end do 
    164 !!$   end do 
    165  
    166  
    167 ! mailles mineures en dehors du domaine 
    168 !imy(:,ny1) = 0 
    169 !imx(nx1,:) = 0            
    170  
    171133 
    172134! limite la valeur de frotmx a betamax 
     
    175137call limit_frotm 
    176138 
    177  
    178139! remplissage des Tuij, ....  
    179140!---------------------------- 
    180141call rempli_Tuij 
    181  
    182 ! debug : impression des Tuij en un point 
    183 !write(166,*) 
    184 !write(166,*) 'impression des Tuij' 
    185  
    186 !j = 12 
    187  
    188 !do i = 244,245 
    189  
    190 !!$do i = 244,246 
    191 !!$   write(166,801) i,imx(i,j),frotmx(i,j),hmx_oppos(i,j),sdx(i,j) 
    192 !!$   do k= 1,-1,-1 
    193 !!$      write(166,800)  k,Tu(i,j,-1:1,k),opposx(i,j) 
    194 !!$   end do 
    195 !!$   write(166,*) 
    196 !!$   do k= 1,-1,-1 
    197 !!$      write(166,800)  k,Tv(i,j,-1:1,k) 
    198 !!$   end do 
    199 !!$   write(166,*) 
    200 !!$   write(166,*) 'impression des Svij' 
    201 !!$   write(166,801) i,imy(i,j),frotmy(i,j),hmy_oppos(i,j),sdy(i,j) 
    202 !!$   do k= 1,-1,-1 
    203 !!$      write(166,800)  k,Su(i,j,-1:1,k),opposy(i,j) 
    204 !!$   end do 
    205 !!$   write(166,*) 
    206 !!$   do k= 1,-1,-1 
    207 !!$      write(166,800)  k,Sv(i,j,-1:1,k) 
    208 !!$   end do 
    209 !!$    
    210 !!$ 
    211 !!$ 
    212 !!$800 format(' k=',(i0,1x),': ',4(es12.4,1x) ) 
    213 !!$801 format(' i,im,frotm,hm,sd   ',2(i0.1,x),3(es12.4,1x)) 
    214 !!$   write(166,*) 
    215 !!$end do 
    216 !!$write(166,*) 
    217  
    218  
    219 ! debug : differents termes de l equation MacAyeal 
    220 !--------------------------------------------------- 
    221  
    222  
    223 !!$do j=ny1+1,ny2-1  
    224 !!$   do i=nx1+1,nx2-1 
    225 !!$ 
    226 !!$      if (gzmx(i,j)) then   
    227 !!$         do jl=-1,1 
    228 !!$            do il=-1,1 
    229 !!$               debug_3D(i,j,82) =debug_3D(i,j,82)+(Tu(i,j,il,jl)*Vcol_x(i+il,j+jl) &  
    230 !!$                    +Tv(i,j,il,jl)*Vcol_y(i+il,j+jl))              
    231 !!$ 
    232 !!$            end do 
    233 !!$         end do 
    234 !!$!         debug_3D(i,j,83)=frotmx(i,j)*dx2*Vcol_x(i,j) 
    235 !!$      end if 
    236 !!$   end do 
    237 !!$end do 
    238  
    239  
    240 !!$do j=ny1,ny2 
    241 !!$   do i=nx1,nx2 
    242 !!$      debug_3D(i,j,84) = imx(i,j) 
    243 !!$   end do 
    244 !!$end do 
    245  
    246 !   call rempli_Tuij                    ! pour reutiliser le nouveau beta 
    247  
    248142 
    249143! Pour determiner les colonnes Mu, Mv, Nu, Nv 
     
    268162         if(abs(scal).lt.1.e-30) then 
    269163            write(num_tracebug,*)'                      i,j,imx,pvi',i,j,imx(i,j),pvi(i,j) 
    270 !            write(num_tracebug,*) Tu(i,j,:,:) 
    271164         end if 
    272165 
     
    275168         Tv(i,j,:,:)=Tv(i,j,:,:)/scal                ! Tu(:,:,0,0) est la diagonale 
    276169         opposx(i,j)=opposx(i,j)/scal 
    277  
    278 !        if ((i.eq.200).and.(j.eq.251)) then  
    279 !            write(num_tracebug,*)'      i,j,imx',i,j,imx(i,j) 
    280 !            write(num_tracebug,*)'     pvi,scal',pvi(i,j),scal 
    281 !            write(num_tracebug,*) Tu(i,j,:,-1) 
    282 !            write(num_tracebug,*) Tu(i,j,:,0) 
    283 !            write(num_tracebug,*) Tu(i,j,:,1) 
    284 !            write(num_tracebug,*)'       opposx',opposx(i,j) 
    285 !         end if 
    286  
    287  
    288170 
    289171      end if 
     
    293175         if(abs(scal).lt.1.e-30) then 
    294176            write(num_tracebug,*)'                      i,j,Svij,imy,pvi',i,j,imy(i,j),pvi(i,j) 
    295 !            write(num_tracebug,*) Sv(i,j,:,:) 
    296177         end if 
    297178 
     
    339220!$OMP END PARALLEL 
    340221 
    341 ! pour appeler la routine de resolution du systeme lineaire 
    342  
    343 !!$i = 244 
    344 !!$do j=11,11,-1 
    345 !!$ 
    346 !!$write(166,*) 'j=',j 
    347 !!$   write(166,803) i,uxprec(i,j),i+1,uxprec(i+1,j),i+2,uxprec(i+2,j) 
    348 !!$   write(166,802) i,uxnew(i,j),i+1,uxnew(i+1,j),i+2,uxnew(i+2,j) 
    349 !!$802 format(3('uxnew ','   i=',(i0,1x,es12.4,1x) )) 
    350 !!$803 format(3('uxprec','   i=',(i0,1x,es12.4,1x) )) 
    351 !!$write(166,*) 
    352  
    353 !!$end do 
    354  
    355 ! limitation des nouvelles vitesses 
    356  
    357 !!$uxnew(:,:)=max(-3900.,uxnew(:,:)) 
    358 !!$uxnew(:,:)=min( 3900.,uxnew(:,:)) 
    359 !!$uynew(:,:)=max(-3900.,uynew(:,:)) 
    360 !!$uynew(:,:)=min( 3900.,uynew(:,:)) 
    361 !!$ 
    362  
    363  
    364  
    365  
    366222888 return 
    367  
    368223 
    369224 
     
    695550real, intent(out) :: opposx 
    696551 
    697  
    698 !  Tu(i,j,0,0)    =  4.*pvi(i-1,j) 
    699 !  Tu(i,j,1,0)    = -4.*pvi(i-1,j) 
    700  
    701552  Tu(0,0)    =  4.*pvi              ! quand le bord est epais, il faut prendre  
    702553  Tu(1,0)    = -4.*pvi              ! le noeud du shelf pas celui en H=1 
    703                                              ! s'il est fin, pas de difference 
    704  
    705 !!$  Tv(i-1,j,-1,0) =  2.*pvi(i-1,j)         ! ???? 
    706 !!$  Tv(i-1,j,-1,1) = -2.*pvi(i-1,j) 
     554                                    ! s'il est fin, pas de difference 
    707555 
    708556  Tv(-1,0) =  2.*pvi_imoinsun 
     
    12271075!$OMP END WORKSHARE 
    12281076!$OMP END PARALLEL 
    1229  
    1230 ! sortie Netcdf pour verifier ok_umat 
    1231 !~ where (ok_umat(:,:)) 
    1232 !~    debug_3D(:,:,39)=1 
    1233 !~ elsewhere 
    1234 !~    debug_3D(:,:,39)=0 
    1235 !~ end where 
    1236 !~  
    1237 !~ where (ok_vmat(:,:)) 
    1238 !~    debug_3D(:,:,40)=1 
    1239 !~ elsewhere 
    1240 !~    debug_3D(:,:,40)=0 
    1241 !~ end where 
    1242 !~  
    1243 !~ ! sortie Netcdf pour verifier ghost 
    1244 !~ where (ghost_x(:,:)) 
    1245 !~    debug_3D(:,:,41)=1 
    1246 !~ elsewhere 
    1247 !~    debug_3D(:,:,41)=0 
    1248 !~  
    1249 !~ end where 
    1250 !~  
    1251 !~ where (ghost_y(:,:)) 
    1252 !~    debug_3D(:,:,42)=1 
    1253 !~ elsewhere 
    1254 !~    debug_3D(:,:,42)=0 
    1255 !~ end where 
    12561077 
    12571078if (itracebug.eq.1)  call tracebug(' Sortie  Okmat') 
     
    13741195         do il=-1,1 
    13751196            betamx(i,j) = betamx(i,j)+(Tu(i,j,il,jl)*uxgiven(i+il,j+jl) &  
    1376                                       +Tv(i,j,il,jl)*uygiven(i+il,j+jl))  
    1377  
    1378 !            debug_3D(i,j,82)=debug_3D(i,j,82) + (Tu(i,j,il,jl)*uxgiven(i+il,j+jl) &  
    1379 !                                      +Tv(i,j,il,jl)*uygiven(i+il,j+jl))  
    1380  
     1197                                      +Tv(i,j,il,jl)*uygiven(i+il,j+jl)) 
    13811198         end do 
    13821199      end do 
    1383  
    1384 !      debug_3D(i,j,82)=debug_3D(i,j,82)/dx/dx   ! longitudinal stress computed from velocities 
    1385  
    13861200 
    13871201!  fonctionne même quand frotmx n'est pas nul 
     
    14081222   
    14091223      betamy(i,j)=-opposy(i,j) 
    1410 !      debug_3D(i,j,36) = opposy(i,j)/dx/dx      ! driving stress ou pression hydrostatique (Pa) 
    14111224      do jl=-1,1 
    14121225         do il=-1,1 
     
    14151228                                      +Sv(i,j,il,jl)*uygiven(i+il,j+jl)) 
    14161229 
    1417 !            debug_3D(i,j,83)=debug_3D(i,j,83) + (Su(i,j,il,jl)*uxgiven(i+il,j+jl) &  
    1418 !                                      +Sv(i,j,il,jl)*uygiven(i+il,j+jl)) 
    1419  
    1420           
    14211230         end do 
    14221231      end do 
    1423  
    1424 !      debug_3D(i,j,83)=debug_3D(i,j,83) /dx/dx  ! longitudinal stress computed from velocities 
    14251232 
    14261233!  fonctionne même quand frotmx n'est pas nul 
     
    14791286!$OMP END PARALLEL 
    14801287 
    1481  
    1482 ! les noeuds negatifs veulent dire qu'il faudrait ajouter un moteur pour aller aussi vite 
    1483 ! que les vitesses de bilan.  
    1484 ! il faudrait repartir plutot que mettre a 0 
    1485  
    1486 ! a mettre ailleurs 
    1487 ! betamx(:,:)=max(betamx(:,:),0.) 
    1488 ! betamy(:,:)=max(betamy(:,:),0.) 
    1489 ! beta_centre(:,:)=max(beta_centre(:,:),0.) 
    1490  
    1491 !  iter_beta=2 
    1492  
    1493  
    14941288if (itracebug.eq.1)  call tracebug(' Fin calc_beta') 
    14951289end subroutine calc_beta 
Note: See TracChangeset for help on using the changeset viewer.