Changeset 393


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

use only in module dragging_param_beta

File:
1 edited

Legend:

Unmodified
Added
Removed

  • branches/GRISLIv3/SOURCES/dragging_param_beta_mod.f90

    r333 r393  
    1818module dragging_param_beta 
    1919 
    20 ! Defini les zones de stream avec : 
    21 ! * un critere sur la pression effective  
    22 ! * un critere de continuite depuis la cote 
    23 ! * un critere sur la courbure du socle (si negatif, vallees)  
    24  
    25 use module3d_phy 
    26 use param_phy_mod 
    27 use interface_input 
    28 use io_netcdf_grisli 
    29 use fake_beta_iter_vitbil_mod 
    30  
    31 implicit none 
    32  
    33 logical,dimension(nx,ny) :: cote 
    34  
    35 real,dimension(nx,ny) :: beta_param ! local variable 
    36  
    37 real :: betamin   ! betamin : (Pa) frottement mini sous les streams 
    38  
    39 real :: beta_slope        ! = A in:     B = A x Neff ** K 
    40 real :: beta_expo         ! = K in:     B = A x Neff ** K 
    41  
    42 real,dimension(nx,ny) :: neff   ! pression effective noeuds majeurs 
    43 real,dimension(nx,ny) :: hwatmx ! hauteur d'eau staggered grille - afq 
    44 real,dimension(nx,ny) :: hwatmy ! hauteur d'eau staggered grille - afq 
    45  
    46 real :: coef_ile      ! coef frottement zones iles   (0.1) 
    47  
    48 real, dimension(nx,ny) :: Vcol_x           !< uniquement pour compatibilite avec spinup cat 
    49 real, dimension(nx,ny) :: Vcol_y           !< uniquement pour compatibilite avec spinup cat 
    50 real, dimension(nx,ny) :: Vsl_x            !< uniquement pour compatibilite avec spinup cat 
    51 real, dimension(nx,ny) :: Vsl_y            !< uniquement pour compatibilite avec spinup cat 
    52 logical :: corr_def = .false.              !< for deformation correction, pour compatibilite beta 
     20  ! Defini les zones de stream avec : 
     21  ! * un critere sur la pression effective  
     22  ! * un critere de continuite depuis la cote 
     23  ! * un critere sur la courbure du socle (si negatif, vallees)  
     24 
     25  use module3d_phy, only:nx,ny,betamax,num_param,num_rep_42,inv_beta,mstream_mx,mstream_my,drag_mx,drag_my,& 
     26       betamax_2d,neffmx,neffmy,niter_nolin,hwater,betamx,betamy,ilemx,ilemy,flot,beta_centre,fleuvemx,fleuvemy,& 
     27       gzmx,gzmy,flgzmx,flgzmy,flotmx,flotmy 
     28  use runparam, only:itracebug 
     29  use param_phy_mod, only: 
     30  use interface_input, only: 
     31  use io_netcdf_grisli 
     32  use fake_beta_iter_vitbil_mod, only:time_iter,time_iter_end,nb_iter_vitbil,beta_iter_vitbil 
     33 
     34  implicit none 
     35 
     36  logical,dimension(nx,ny) :: cote 
     37 
     38  real,dimension(nx,ny) :: beta_param ! local variable 
     39 
     40  real :: betamin   ! betamin : (Pa) frottement mini sous les streams 
     41 
     42  real :: beta_slope        ! = A in:     B = A x Neff ** K 
     43  real :: beta_expo         ! = K in:     B = A x Neff ** K 
     44 
     45  real,dimension(nx,ny) :: neff   ! pression effective noeuds majeurs 
     46  real,dimension(nx,ny) :: hwatmx ! hauteur d'eau staggered grille - afq 
     47  real,dimension(nx,ny) :: hwatmy ! hauteur d'eau staggered grille - afq 
     48 
     49  real :: coef_ile      ! coef frottement zones iles   (0.1) 
     50 
     51  real, dimension(nx,ny) :: Vcol_x           !< uniquement pour compatibilite avec spinup cat 
     52  real, dimension(nx,ny) :: Vcol_y           !< uniquement pour compatibilite avec spinup cat 
     53  real, dimension(nx,ny) :: Vsl_x            !< uniquement pour compatibilite avec spinup cat 
     54  real, dimension(nx,ny) :: Vsl_y            !< uniquement pour compatibilite avec spinup cat 
     55  logical :: corr_def = .false.              !< for deformation correction, pour compatibilite beta 
    5356 
    5457contains 
    55 !------------------------------------------------------------------------------- 
    56  
    57 !> SUBROUTINE: init_dragging 
    58 !! Cette routine fait l'initialisation du dragging. 
    59 !> 
    60 subroutine init_dragging      ! Cette routine fait l'initialisation du dragging. 
    61  
    62 implicit none 
    63  
    64 namelist/drag_param_beta/beta_slope,beta_expo,betamax,betamin,coef_ile 
    65  
    66 if (itracebug.eq.1)  call tracebug(' dragging avec neff et slope') 
    67  
    68  
    69 ! formats pour les ecritures dans 42 
     58  !------------------------------------------------------------------------------- 
     59 
     60  !> SUBROUTINE: init_dragging 
     61  !! Cette routine fait l'initialisation du dragging. 
     62  !> 
     63  subroutine init_dragging      ! Cette routine fait l'initialisation du dragging. 
     64 
     65    implicit none 
     66 
     67    namelist/drag_param_beta/beta_slope,beta_expo,betamax,betamin,coef_ile 
     68 
     69    if (itracebug.eq.1)  call tracebug(' dragging avec neff et slope') 
     70 
     71 
     72    ! formats pour les ecritures dans 42 
    7073428 format(A) 
    7174 
    72 ! lecture des parametres du run                      block drag neff 
    73 !-------------------------------------------------------------------- 
    74  
    75 rewind(num_param)        ! pour revenir au debut du fichier param_list.dat 
    76 read(num_param,drag_param_beta) 
    77  
    78 write(num_rep_42,428)'!___________________________________________________________'  
    79 write(num_rep_42,428) '&drag_param_beta          ! nom du bloc dragging param beta' 
    80 write(num_rep_42,*) 
    81 write(num_rep_42,*) 'beta_slope            = ', beta_slope 
    82 write(num_rep_42,*) 'beta_expo             = ', beta_expo 
    83 write(num_rep_42,*) 'betamax               = ', betamax 
    84 write(num_rep_42,*) 'betamin               = ', betamin 
    85 write(num_rep_42,*)'/'                             
    86 write(num_rep_42,428) '! Slope & expo of beta = - slope x Neff ** expo' 
    87  
    88 inv_beta=0 
    89 !------------------------------------------------------------------- 
    90 ! masque stream 
    91 mstream_mx(:,:)=1 
    92 mstream_my(:,:)=1 
    93  
    94  
    95 ! coefficient permettant de modifier le basal drag. 
    96 drag_mx(:,:)=1. 
    97 drag_my(:,:)=1. 
    98        
    99 betamax_2d(:,:) = betamax 
    100  
    101 !afq -- toblim = 0. !0.7e5 ! afq -- pour les iles, mais est-ce vraiment utile? 
    102  
    103 niter_nolin = 1 
    104  
    105 return 
    106 end subroutine init_dragging 
    107 !________________________________________________________________________________ 
    108  
    109 !> SUBROUTINE: dragging 
    110 !! Defini la localisation des streams et le frottement basal 
    111 !> 
    112 !------------------------------------------------------------------------- 
    113 subroutine dragging   ! defini le frottement basal 
    114 !$ USE OMP_LIB 
    115  
    116 !         les iles n'ont pas de condition neff mais ont la condition toblim 
    117 !         (ce bloc etait avant dans flottab) 
    118  
    119 !$OMP PARALLEL 
    120  
    121 !    -- !$OMP DO 
    122 !    -- do j=2,ny 
    123 !    --    do i=2,nx 
    124 !    --       ilemx(i,j)=ilemx(i,j).and.(abs(rog*Hmx(i,j)*sdx(i,j)).lt.toblim) 
    125 !    --       ilemy(i,j)=ilemy(i,j).and.(abs(rog*Hmy(i,j)*sdy(i,j)).lt.toblim) 
    126 !    --    end do 
    127 !    -- end do 
    128 !    -- !$OMP END DO 
    129  
    130  
    131  
    132 ! calcul de neff (pression effective sur noeuds majeurs) 
    133 if (sum(neffmx(:,:)).le.0.) neffmx(:,:) =1.e8 
    134 if (sum(neffmy(:,:)).le.0.) neffmy(:,:) =1.e8 
    135  
    136 !$OMP DO 
    137 do j=1,ny-1 
    138    do i=1,nx-1 
    139         neff(i,j)=(neffmx(i,j)+neffmx(i+1,j)+neffmy(i,j)+neffmy(i,j+1))/4 
    140    enddo 
    141 enddo 
    142 !$OMP END DO 
    143 !aurel, for the borders: 
    144 !$OMP WORKSHARE 
    145 neff(:,ny)=1.e5 
    146 neff(nx,:)=1.e5 
    147 ! calcul de hwat (staggered) 
    148 hwatmx(:,:)=0. 
    149 hwatmy(:,:)=0. 
    150 !$OMP END WORKSHARE 
    151 !$OMP DO 
    152 do i=2,nx 
    153    hwatmx(i,:)=(hwater(i-1,:)+hwater(i,:))/2. 
    154 enddo 
    155 !$OMP END DO 
    156 !$OMP DO 
    157 do j=2,ny 
    158    hwatmy(:,j)=(hwater(:,j-1)+hwater(:,j))/2. 
    159 enddo 
    160 !$OMP END DO 
    161  
    162 !$OMP WORKSHARE 
    163  
    164 ! new parametrisation of beta on Neff: 
    165 betamx(:,:)= beta_slope*(neffmx(:,:)**beta_expo) 
    166 betamy(:,:)= beta_slope*(neffmy(:,:)**beta_expo) 
    167  
    168 where (ilemx(:,:)) betamx(:,:) = betamx(:,:) * coef_ile 
    169 where (ilemy(:,:)) betamy(:,:) = betamy(:,:) * coef_ile 
    170  
    171  
    172 betamx(:,:)=max(betamx(:,:),betamin) 
    173 betamy(:,:)=max(betamy(:,:),betamin) 
    174 betamx(:,:)=min(betamx(:,:),betamax) 
    175 betamy(:,:)=min(betamy(:,:),betamax) 
    176  
    177 where ( hwatmx(:,:) .lt. 0.5 ) betamx(:,:) = betamax 
    178 where ( hwatmy(:,:) .lt. 0.5 ) betamy(:,:) = betamax 
    179  
    180 !$OMP END WORKSHARE 
    181  
    182 ! points flottants  
    183 !$OMP DO 
    184 do j=2,ny 
    185    do i=2,nx 
    186       if (flot(i,j).and.(flot(i-1,j))) then 
    187          betamx(i,j)=0. 
    188       end if 
    189       if (flot(i,j).and.(flot(i,j-1))) then 
    190          betamy(i,j)=0. 
    191       end if 
    192    end do 
    193 end do 
    194 !$OMP END DO 
    195  
    196  
    197 !$OMP DO 
    198 do j=2,ny-1 
    199    do i=2,nx-1 
    200       beta_centre(i,j) = ((betamx(i,j)+betamx(i+1,j)) & 
    201           + (betamy(i,j)+betamy(i,j+1)))*0.25 
    202    end do 
    203 end do 
    204 !$OMP END DO 
    205 !$OMP END PARALLEL 
    206  
    207 return 
    208 end subroutine dragging 
    209 !________________________________________________________________________________ 
    210  
    211 !> SUBROUTINE: mstream_dragging 
    212 !! Defini la localisation des streams 
    213 !> 
    214 !------------------------------------------------------------------------- 
    215 subroutine mstream_dragging   ! defini la localisation des streams 
    216  
    217 !$OMP PARALLEL 
    218  
    219 !$OMP WORKSHARE 
    220 fleuvemx(:,:)=.false. 
    221 fleuvemy(:,:)=.false. 
    222 cote(:,:)=.false. 
    223 gzmx(:,:)=.true. 
    224 gzmy(:,:)=.true. 
    225 flgzmx(:,:)=.false. 
    226 flgzmy(:,:)=.false. 
    227 !$OMP END WORKSHARE 
    228  
    229 ! detection des cotes 
    230 !$OMP DO 
    231 do  j=2,ny-1   
    232    do i=2,nx-1   
    233       if ((.not.flot(i,j)).and. &  
    234       ((flot(i+1,j)).or.(flot(i,j+1)).or.(flot(i-1,j)).or.(flot(i,j-1)))) then 
    235          cote(i,j)=.true. 
    236       endif 
    237    end do 
    238 end do 
    239 !$OMP END DO 
    240  
    241  
    242 ! calcul de gzmx  
    243  
    244 ! points flottants : flgzmx mais pas gzmx 
    245 !$OMP DO 
    246 do j=2,ny 
    247    do i=2,nx 
    248       if (flot(i,j).and.(flot(i-1,j))) then 
    249          flotmx(i,j)=.true. 
    250          flgzmx(i,j)=.true. 
    251          gzmx(i,j)=.false. 
    252          betamx(i,j)=0. 
    253       end if 
    254       if (flot(i,j).and.(flot(i,j-1))) then 
    255          flotmy(i,j)=.true. 
    256          flgzmy(i,j)=.true. 
    257          gzmy(i,j)=.false. 
    258          betamy(i,j)=0. 
    259       end if 
    260    end do 
    261 end do 
    262 !$OMP END DO 
    263  
    264 !--------- autres criteres 
    265  
    266 ! points poses 
    267 !$OMP WORKSHARE 
    268 gzmx(:,:)=gzmx(:,:).and.(betamx(:,:).lt.betamax)   !  Pas de calcul pour les points 
    269 gzmy(:,:)=gzmy(:,:).and.(betamy(:,:).lt.betamax)   !  au fort beta 
    270  
    271 flgzmx(:,:) = flgzmx(:,:) .or. gzmx(:,:) 
    272 flgzmy(:,:) = flgzmy(:,:) .or. gzmy(:,:) 
    273  
    274 fleuvemx(:,:)=gzmx(:,:) 
    275 fleuvemy(:,:)=gzmy(:,:) 
    276 !$OMP END WORKSHARE 
    277  
    278 !$OMP END PARALLEL 
    279  
    280 return 
    281 end subroutine mstream_dragging 
     75    ! lecture des parametres du run                      block drag neff 
     76    !-------------------------------------------------------------------- 
     77 
     78    rewind(num_param)        ! pour revenir au debut du fichier param_list.dat 
     79    read(num_param,drag_param_beta) 
     80 
     81    write(num_rep_42,428)'!___________________________________________________________'  
     82    write(num_rep_42,428) '&drag_param_beta          ! nom du bloc dragging param beta' 
     83    write(num_rep_42,*) 
     84    write(num_rep_42,*) 'beta_slope            = ', beta_slope 
     85    write(num_rep_42,*) 'beta_expo             = ', beta_expo 
     86    write(num_rep_42,*) 'betamax               = ', betamax 
     87    write(num_rep_42,*) 'betamin               = ', betamin 
     88    write(num_rep_42,*)'/'                             
     89    write(num_rep_42,428) '! Slope & expo of beta = - slope x Neff ** expo' 
     90 
     91    inv_beta=0 
     92    !------------------------------------------------------------------- 
     93    ! masque stream 
     94    mstream_mx(:,:)=1 
     95    mstream_my(:,:)=1 
     96 
     97 
     98    ! coefficient permettant de modifier le basal drag. 
     99    drag_mx(:,:)=1. 
     100    drag_my(:,:)=1. 
     101 
     102    betamax_2d(:,:) = betamax 
     103 
     104    !afq -- toblim = 0. !0.7e5 ! afq -- pour les iles, mais est-ce vraiment utile? 
     105 
     106    niter_nolin = 1 
     107 
     108    return 
     109  end subroutine init_dragging 
     110  !________________________________________________________________________________ 
     111 
     112  !> SUBROUTINE: dragging 
     113  !! Defini la localisation des streams et le frottement basal 
     114  !> 
     115  !------------------------------------------------------------------------- 
     116  subroutine dragging   ! defini le frottement basal 
     117 
     118    !$ USE OMP_LIB 
     119 
     120    !         les iles n'ont pas de condition neff mais ont la condition toblim 
     121    !         (ce bloc etait avant dans flottab) 
     122 
     123     
     124 
     125    !    -- !$OMP DO 
     126    !    -- do j=2,ny 
     127    !    --    do i=2,nx 
     128    !    --       ilemx(i,j)=ilemx(i,j).and.(abs(rog*Hmx(i,j)*sdx(i,j)).lt.toblim) 
     129    !    --       ilemy(i,j)=ilemy(i,j).and.(abs(rog*Hmy(i,j)*sdy(i,j)).lt.toblim) 
     130    !    --    end do 
     131    !    -- end do 
     132    !    -- !$OMP END DO 
     133     
     134    implicit none 
     135    integer :: i,j 
     136 
     137    !$OMP PARALLEL 
     138    ! calcul de neff (pression effective sur noeuds majeurs) 
     139    if (sum(neffmx(:,:)).le.0.) neffmx(:,:) =1.e8 
     140    if (sum(neffmy(:,:)).le.0.) neffmy(:,:) =1.e8 
     141 
     142    !$OMP DO 
     143    do j=1,ny-1 
     144       do i=1,nx-1 
     145          neff(i,j)=(neffmx(i,j)+neffmx(i+1,j)+neffmy(i,j)+neffmy(i,j+1))/4 
     146       enddo 
     147    enddo 
     148    !$OMP END DO 
     149    !aurel, for the borders: 
     150    !$OMP WORKSHARE 
     151    neff(:,ny)=1.e5 
     152    neff(nx,:)=1.e5 
     153    ! calcul de hwat (staggered) 
     154    hwatmx(:,:)=0. 
     155    hwatmy(:,:)=0. 
     156    !$OMP END WORKSHARE 
     157    !$OMP DO 
     158    do i=2,nx 
     159       hwatmx(i,:)=(hwater(i-1,:)+hwater(i,:))/2. 
     160    enddo 
     161    !$OMP END DO 
     162    !$OMP DO 
     163    do j=2,ny 
     164       hwatmy(:,j)=(hwater(:,j-1)+hwater(:,j))/2. 
     165    enddo 
     166    !$OMP END DO 
     167 
     168    !$OMP WORKSHARE 
     169 
     170    ! new parametrisation of beta on Neff: 
     171    betamx(:,:)= beta_slope*(neffmx(:,:)**beta_expo) 
     172    betamy(:,:)= beta_slope*(neffmy(:,:)**beta_expo) 
     173 
     174    where (ilemx(:,:)) betamx(:,:) = betamx(:,:) * coef_ile 
     175    where (ilemy(:,:)) betamy(:,:) = betamy(:,:) * coef_ile 
     176 
     177 
     178    betamx(:,:)=max(betamx(:,:),betamin) 
     179    betamy(:,:)=max(betamy(:,:),betamin) 
     180    betamx(:,:)=min(betamx(:,:),betamax) 
     181    betamy(:,:)=min(betamy(:,:),betamax) 
     182 
     183    where ( hwatmx(:,:) .lt. 0.5 ) betamx(:,:) = betamax 
     184    where ( hwatmy(:,:) .lt. 0.5 ) betamy(:,:) = betamax 
     185 
     186    !$OMP END WORKSHARE 
     187 
     188    ! points flottants  
     189    !$OMP DO 
     190    do j=2,ny 
     191       do i=2,nx 
     192          if (flot(i,j).and.(flot(i-1,j))) then 
     193             betamx(i,j)=0. 
     194          end if 
     195          if (flot(i,j).and.(flot(i,j-1))) then 
     196             betamy(i,j)=0. 
     197          end if 
     198       end do 
     199    end do 
     200    !$OMP END DO 
     201 
     202 
     203    !$OMP DO 
     204    do j=2,ny-1 
     205       do i=2,nx-1 
     206          beta_centre(i,j) = ((betamx(i,j)+betamx(i+1,j)) & 
     207               + (betamy(i,j)+betamy(i,j+1)))*0.25 
     208       end do 
     209    end do 
     210    !$OMP END DO 
     211    !$OMP END PARALLEL 
     212 
     213    return 
     214  end subroutine dragging 
     215  !________________________________________________________________________________ 
     216 
     217  !> SUBROUTINE: mstream_dragging 
     218  !! Defini la localisation des streams 
     219  !> 
     220  !------------------------------------------------------------------------- 
     221  subroutine mstream_dragging   ! defini la localisation des streams 
     222 
     223    implicit none 
     224    integer :: i,j 
     225     
     226    !$OMP PARALLEL 
     227 
     228    !$OMP WORKSHARE 
     229    fleuvemx(:,:)=.false. 
     230    fleuvemy(:,:)=.false. 
     231    cote(:,:)=.false. 
     232    gzmx(:,:)=.true. 
     233    gzmy(:,:)=.true. 
     234    flgzmx(:,:)=.false. 
     235    flgzmy(:,:)=.false. 
     236    !$OMP END WORKSHARE 
     237 
     238    ! detection des cotes 
     239    !$OMP DO 
     240    do  j=2,ny-1   
     241       do i=2,nx-1   
     242          if ((.not.flot(i,j)).and. &  
     243               ((flot(i+1,j)).or.(flot(i,j+1)).or.(flot(i-1,j)).or.(flot(i,j-1)))) then 
     244             cote(i,j)=.true. 
     245          endif 
     246       end do 
     247    end do 
     248    !$OMP END DO 
     249 
     250 
     251    ! calcul de gzmx  
     252 
     253    ! points flottants : flgzmx mais pas gzmx 
     254    !$OMP DO 
     255    do j=2,ny 
     256       do i=2,nx 
     257          if (flot(i,j).and.(flot(i-1,j))) then 
     258             flotmx(i,j)=.true. 
     259             flgzmx(i,j)=.true. 
     260             gzmx(i,j)=.false. 
     261             betamx(i,j)=0. 
     262          end if 
     263          if (flot(i,j).and.(flot(i,j-1))) then 
     264             flotmy(i,j)=.true. 
     265             flgzmy(i,j)=.true. 
     266             gzmy(i,j)=.false. 
     267             betamy(i,j)=0. 
     268          end if 
     269       end do 
     270    end do 
     271    !$OMP END DO 
     272 
     273    !--------- autres criteres 
     274 
     275    ! points poses 
     276    !$OMP WORKSHARE 
     277    gzmx(:,:)=gzmx(:,:).and.(betamx(:,:).lt.betamax)   !  Pas de calcul pour les points 
     278    gzmy(:,:)=gzmy(:,:).and.(betamy(:,:).lt.betamax)   !  au fort beta 
     279 
     280    flgzmx(:,:) = flgzmx(:,:) .or. gzmx(:,:) 
     281    flgzmy(:,:) = flgzmy(:,:) .or. gzmy(:,:) 
     282 
     283    fleuvemx(:,:)=gzmx(:,:) 
     284    fleuvemy(:,:)=gzmy(:,:) 
     285    !$OMP END WORKSHARE 
     286 
     287    !$OMP END PARALLEL 
     288 
     289    return 
     290  end subroutine mstream_dragging 
    282291 
    283292end module dragging_param_beta 
Note: See TracChangeset for help on using the changeset viewer.