Changeset 167


Ignore:
Timestamp:
12/13/17 18:09:05 (7 years ago)
Author:
dumas
Message:

module beta_iter_vitbil_mod is now used by dragging_beta_iter_vitbil_mod

Location:
trunk/SOURCES/Draggings_modules
Files:
1 added
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/SOURCES/Draggings_modules/beta_iter_vitbil_mod.f90

    r70 r167  
    2020  implicit none 
    2121 
    22  
    23   real*8, dimension(:,:),   pointer  :: tab => null() !< tableau 2d real ecrit dans le fichier pour lect ncdf 
    24  
    2522  real, dimension(nx,ny) :: Vcol_x           !< uniquement pour compatibilite 
    2623  real, dimension(nx,ny) :: Vcol_y           !< 
    … …  
    5047  character(len=100)     :: Umag_bil_file           !< fichier qui contient les donnees 
    5148  real                   :: time_iter               !< temps de debut des iterations 
     49  real                   :: time_iter_end           !< temps de fin des iterations 
     50  real                   :: time_reiter             !< nbr annees entre 2 iterations calcul beta 
    5251  integer                :: nb_iter_vitbil          !< nombre d'iteratiosn avant de changer de pas de temps 
    5352  real                   :: coef_iter_vitbil        !< coefficient <= 1 (rapport vitesses) 
    … …  
    6059contains 
    6160 
    62   !----------------------------------------------------------------------------------------- 
    63   !> SUBROUTINE: init_dragging  
    64   !! Cette routine fait l'initialisation du dragging. 
    65   !< 
    66   subroutine init_dragging      ! Cette routine fait l'initialisation du dragging. 
    67     ! nouvelle version : lit les fichiers nc 
    68     implicit none 
    69     character(len=100) :: beta_c_file               !  beta on centered grid 
    70     character(len=100) :: betamx_file               !  beta on staggered grid mx 
    71     character(len=100) :: betamy_file               !  beta on staggered grid mx 
    72  
    73     namelist/beta_prescr/beta_c_file,beta_limgz,beta_min,beta_mult 
    74  
    75     if (itracebug.eq.1)  call tracebug(' Prescr_beta       subroutine init_dragging') 
    76  
    77     iter_beta=0 
    78  
    79     ! lecture des parametres du run                     
    80     !-------------------------------------------------------------------- 
    81  
    82     rewind(num_param)        ! pour revenir au debut du fichier param_list.dat 
    83 428 format(A) 
    84     read(num_param,beta_prescr) 
    85  
    86     write(num_rep_42,428) '!___________________________________________________________'  
    87     write(num_rep_42,428) '!  read beta on centered grid' 
    88     write(num_rep_42,beta_prescr)                     
    89     write(num_rep_42,428) '!  beta_file : nom des fichier qui contiennent les beta_c' 
    90     write(num_rep_42,428) '!  above beta_limgz, gzmx is false' 
    91     write(num_rep_42,428) '!  if grounded, beta_min minimum value ' 
    92     write(num_rep_42,428) '!  beta_mult : coefficient just after reading'  
    93     write(num_rep_42,*) 
    94     write(num_rep_42,428) '!___________________________________________________________'  
    95  
    96  
    97     beta_c_file = trim(dirnameinp)//trim(beta_c_file) 
    98     betamx_file = trim(dirnameinp)//trim(betamx_file) 
    99     betamy_file = trim(dirnameinp)//trim(betamy_file) 
    100     betamax = beta_limgz 
    101     betamax_2d(:,:) = betamax 
    102  
    103 !    call lect_input(1,'beta_c',1,drag_centre,beta_c_file,'') 
    104     write(6,*)  beta_c_file 
    105  
    106     call Read_Ncdf_var('z',beta_c_file,tab) 
    107     drag_centre(:,:) = tab(:,:) 
    108  
    109  
    110     ! multiplication par le coefficient beta_mult 
    111  
    112     drag_centre(:,:) = drag_centre(:,:)*beta_mult 
    113  
    114     where (mk_init(:,:).eq.-2)                                             ! iles a probleme 
    115        drag_centre(:,:) = 2.* abs(beta_limgz) 
    116     end where 
    117  
    118     where(flot(:,:)) 
    119        drag_centre(:,:) = 0. 
    120     end where 
    121  
    122  
    123     call  beta_c2stag(nx,ny,drag_mx,drag_my,drag_centre)         ! redistribute on staggered grid 
    124  
    125     if  (beta_limgz.gt.0.) then 
    126        beta_centre(:,:) =  drag_centre(:,:) 
    127        betamx(:,:)      =  drag_mx(:,:) 
    128        betamy(:,:)      =  drag_my(:,:) 
    129  
    130     else if (beta_limgz.lt.0.) then                              ! attention sans doute pour plastic 
    131  
    132        beta_centre(:,:) = - beta_limgz 
    133        betamx(:,:)      = - beta_limgz 
    134        betamy(:,:)      = - beta_limgz 
    135        drag_centre(:,:) = - beta_limgz  
    136        drag_mx(:,:)     = - beta_limgz 
    137        drag_my(:,:)     = - beta_limgz  
    138        beta_limgz = 1. 
    139     end if 
    140  
    141  
    142  
    143     call beta_min_value(beta_min)      !  valeur mini que peut avoir beta (en Pa an m-1) 
    144     ! on peut envisager une valeur ~ 10  
    145     ! rappel : beta doit être positif. 
    146  
    147  
    148  
    149     !    call beta_c2stag(nx,ny,drag_mx,drag_my,drag_centre)                  ! centered distributed 
    150     !                                                                         ! on staggered 
    151     ! 
    152  
    153  
    154     beta_centre(:,:) =  drag_centre(:,:)                                 ! surtout pour sorties 
    155  
    156  
    157  
    158     where (drag_mx(:,:).ge.beta_limgz) 
    159        mstream_mx(:,:) = 0 
    160        betamx(:,:)     = beta_limgz 
    161        drag_mx(:,:)    = beta_limgz 
    162     elsewhere 
    163        mstream_mx(:,:) = 1 
    164        betamx(:,:)     = drag_mx(:,:) 
    165     endwhere 
    166  
    167     where (drag_my(:,:).ge.beta_limgz) 
    168        mstream_my(:,:) = 0 
    169        betamy(:,:)     = beta_limgz 
    170        drag_my(:,:)    = beta_limgz 
    171     elsewhere 
    172        mstream_my(:,:) = 1 
    173        betamy(:,:)     = drag_my(:,:) 
    174     endwhere 
    175  
    176     if (itracebug.eq.1)  call tracebug(' Prescr_beta  apres lecture') 
    177     mstream_mx(:,:) = 0 
    178     mstream_my(:,:) = 0 
    179  
    180  
    181     call gzm_beta_prescr 
    182  
    183     call init_beta_iter_vitbil 
    184     return 
    185  
    186   end subroutine init_dragging 
    187   !________________________________________________________________________________ 
    188  
    18961!----------------------------------------------------------------------------------------- 
    19062!< subroutine : init_beta_iter_vitbil 
    … …  
    19365subroutine init_beta_iter_vitbil 
    19466 
    195  
    196   namelist/beta_iter_vitbil/time_iter,nb_iter_vitbil,coef_iter_vitbil,Umag_bil_file 
     67  implicit none 
     68   
     69  real*8, dimension(:,:),   pointer  :: tab => null() !< tableau 2d real ecrit dans le fichier pour lect ncdf 
     70 
     71    namelist/beta_iter_vitbil/time_iter,nb_iter_vitbil,coef_iter_vitbil,Umag_bil_file,time_reiter 
    19772 
    19873 
    … …  
    20984    write(num_rep_42,428) '!  nb_iter_vitbil: nombre diteratation a faire avant de changer de pas de temps '  
    21085    write(num_rep_42,*)   '!  coef_iter_vitbil : coefficient pour le rapport vitesses  ' 
     86    write(num_rep_42,*)   '!  time_reiter : nombre d annees entre 2 iterations' 
    21187    write(num_rep_42,428) '!___________________________________________________________'  
    21288 
    21389    time_iter = time_iter + tbegin       ! time_iter est le temps de relaxation (de l'ordre de 5 ans) 
    21490                                         ! ajouter tbegin pour ne pas dependre du temps de depart. 
    215  
     91    time_iter_end = time_iter + 1.       ! 20. ans en version std Cat et Seb 
     92    print*,'debug init_beta_iter_vitbil time_iter=', time_iter, time_iter_end, time 
    21693 
    21794    Umag_bil_file = trim(dirnameinp)//trim(Umag_bil_file) 
    … …  
    227104 
    228105 
    229   !________________________________________________________________________________ 
    230  
    231   !> SUBROUTINE: dragging  
    232   !! Defini la localisation des streams et le frottement basal 
    233   !< 
    234   !------------------------------------------------------------------------- 
    235   subroutine dragging   ! defini la localisation des streams et le frottement basal 
    236  
    237     implicit none 
    238  
    239  
    240  
    241     if (itracebug.eq.1)  call tracebug(' beta_prescr         subroutine dragging') 
    242  
    243  
    244     betamx(:,:)=drag_mx(:,:) 
    245     betamy(:,:)=drag_my(:,:) 
    246  
    247  
    248     call  gzm_beta_prescr        ! determine gz et flgz et met a 0 le beta des points flottants 
    249  
    250  
    251     return 
    252   end subroutine dragging 
    253  
    254   !---------------------------------------------------------------------------------- 
    255   !>SUBROUTINE: gzm_beta_prescr 
    256   !! Calcul de gzmx 
    257   !< 
    258  
    259   subroutine gzm_beta_prescr     ! attribue flgzmx et gzmx (et my) 
    260      
    261     logical ::  test_Tx          ! test if there is a basal melting point in the surrounding 
    262     logical ::  test_Ty          ! test if there is a basal melting point in the surrounding 
    263     logical ::  test_beta_x      ! test if there is a low drag point in the surrounding 
    264     logical ::  test_beta_y      ! test if there is a low drag point in the surrounding 
    265  
    266     real    ::  beta_forc_fleuve ! below this value -> ice stream 
    267  
    268 !    beta_forc_fleuve = 5.e3 
    269     beta_forc_fleuve = beta_limgz 
    270  
    271     ! calcul de gzmx  
    272  
    273  
    274     gzmx(:,:)=.true. 
    275     gzmy(:,:)=.true. 
    276     flgzmx(:,:)=.false. 
    277     flgzmy(:,:)=.false. 
    278  
    279  
    280     ! points flottants : flgzmx mais pas gzmx 
    281     do j=2,ny 
    282        do i=2,nx 
    283           if (flot(i,j).and.(flot(i-1,j))) then 
    284              flotmx(i,j)=.true. 
    285              flgzmx(i,j)=.true. 
    286              gzmx(i,j)=.false. 
    287              betamx(i,j)=0. 
    288           end if 
    289           if (flot(i,j).and.(flot(i,j-1))) then 
    290              flotmy(i,j)=.true. 
    291              flgzmy(i,j)=.true. 
    292              gzmy(i,j)=.false. 
    293              betamy(i,j)=0. 
    294           end if 
    295        end do 
    296     end do 
    297  
    298     if (itracebug.eq.1)  call tracebug(' apres criteres flot') 
    299  
    300     if (itracebug.eq.1) write(num_tracebug,*) 'gzmx', gzmx(127,330)  
    301  
    302     !--------- autres criteres 
    303  
    304     ! points poses 
    305     gzmx(:,:)=gzmx(:,:).and.(betamx(:,:).lt.beta_limgz)   !  Pas de calcul pour les points 
    306     gzmy(:,:)=gzmy(:,:).and.(betamy(:,:).lt.beta_limgz)   !  au fort beta 
    307  
    308  
    309  
    310     ! critere (lache) sur la temperature 
    311     do j = 2, ny-1 
    312        do i =2, nx-1  
    313  
    314 !   test s'il y a un point tempere dans les environs 
    315           test_Tx = any( ibase( i-1:i   , j-1:j+1 ).gt.1) 
    316           test_Ty = any( ibase( i-1:i+1 , j-1:j   ).gt.1) 
    317  
    318 !   test s'il y a un point low drag dans les environs 
    319           test_beta_x = any( drag_centre( i-1:i   , j-1:j+1 ) .lt. beta_forc_fleuve ) 
    320           test_beta_y = any( drag_centre( i-1:i+1 , j-1:j   ) .lt. beta_forc_fleuve ) 
    321  
    322 !   critere pour gzmx 
    323  
    324 ! Attention : les deux lignes suivants sont en commentaire pour voir l'effet 
    325  
    326 !          gzmx(i,j) =  gzmx(i,j) .and. (test_Tx.or.test_beta_x) 
    327 !          gzmy(i,j) =  gzmy(i,j) .and. (test_Ty.or.test_beta_y) 
    328  
    329        end do 
    330     end do 
    331  
    332  
    333     if (itracebug.eq.1)  call tracebug(' apres autres criteres ') 
    334     if (itracebug.eq.1) write(num_tracebug,*) 'gzmx', gzmx(127,330)  
    335  
    336     flgzmx(:,:) = flgzmx(:,:) .or. gzmx(:,:) 
    337     flgzmy(:,:) = flgzmy(:,:) .or. gzmy(:,:) 
    338  
    339     where ((.not.flotmx(:,:)).and.(.not.gzmx(:,:))) 
    340        betamx(:,:) = beta_limgz 
    341     end where 
    342  
    343     where ((.not.flotmy(:,:)).and.(.not.gzmy(:,:))) 
    344        betamy(:,:) = beta_limgz 
    345     end where 
    346  
    347         fleuvemx(:,:)=gzmx(:,:) 
    348         fleuvemy(:,:)=gzmy(:,:) 
    349     
    350   end subroutine gzm_beta_prescr 
     106 
    351107 
    352108!______________________________________________________________________________________ 
    … …  
    394150!< fait des iterations pour affiner la valeur de beta_centre pour s'ajuster sur vitbil 
    395151!----------------------------------------------------------------------------------------- 
    396 subroutine beta_iter_vitbil 
    397  
     152subroutine beta_iter_vitbil(m_iter) 
     153 
     154implicit none 
     155integer :: m_iter ! indice bloucle iteration 
     156 
     157 
     158! calcul des vitesses cibles : 
     159!if ((time.eq.time_iter).and.(time.gt.time_reiter).and.(m_iter.eq.1)) then 
     160if ((abs(time-time_iter).lt.dtmin).and.(time.gt.time_reiter).and.(m_iter.eq.1)) then 
     161  Umag_direct(:,:)    = (((uxbar(:,:)+eoshift(uxbar(:,:),shift=1,boundary=0.,dim=1))**2+ &        ! moyenne 
     162                          (uybar(:,:)+eoshift(uybar(:,:),shift=1,boundary=0.,dim=2))**2)**0.5)*0.5 
     163  Umag_vitbil(:,:)= min(Umag_direct(:,:) * min(1.e3,H(:,:)/H0(:,:)), 5000.) 
     164  print*,'*****debug beta_iter_vitbil calcul de Umag_vitbil',time,m_iter 
     165  do j=1,ny 
     166    do i=1,nx 
     167      write(266,*) Umag_vitbil(i,j) 
     168    enddo 
     169  enddo 
     170endif 
    398171 
    399172! calcule la vitesse centree venant du calcul direct 
    … …  
    485258drag_my(:,:)     = betamy(:,:) 
    486259 
    487 i=319 
    488 j=113 
    489 write(266,'(2(i0,1x),3x,9(es12.4,1x))') nb_umag,nb_udef,time,J_Umag,J_Udef,Umag_direct(i,j),Umag_vitbil(i,j),Uslmag_direct(i,j),Uslid_vitbil(i,j),beta_centre(i,j),betamx(i,j),drag_centre(i,j) 
    490260 
    491261end subroutine beta_iter_vitbil 
Note: See TracChangeset for help on using the changeset viewer.