Changeset 485 for branches


Ignore:
Timestamp:
02/11/24 16:02:32 (4 months ago)
Author:
aquiquet
Message:

Cleaning branch: some minor cleaning + remove local isostasy computation (flag nlith)

Location:
branches/GRISLIv3/SOURCES
Files:
9 edited

Legend:

Unmodified
Added
Removed

  • branches/GRISLIv3/SOURCES/flottab2-0.7.f90

    r481 r485  
    281281    !  ------------------------------------------------------------- 
    282282 
    283     if (nt.ge.2) then   ! pour ne pas faire ce calcul lors du premier passage 
     283    if (nt.ge.2) then   ! pour ne pas faire ce calcul lors du premier passage -- afq fev2024, pb restart?? 
    284284       !$OMP WORKSHARE 
    285285       uxbar(:,:)=uxs1(:,:) 

  • branches/GRISLIv3/SOURCES/initial-phy-2.f90

    r465 r485  
    2222       sealevel_2d,secyear,isynchro 
    2323  use geography, only: geoplace 
    24   use runparam, only :runname,itracebug,num_tracebug,tbegin,tend,dirsource,dirnameout,dttest,& 
    25        nt 
     24  use runparam, only :runname,itracebug,num_tracebug,tbegin,tend,dirsource,dirnameout,nt 
    2625 
    2726 
    … …  
    151150  dt=dtmin           ! sera réattribue a chaque pas de temps 
    152151  ntmax=90000000     ! nombre de tours maxi dans la boucle temps 
    153   dttest=dtmin       ! sert dans plein de tests  
    154152  time=tbegin  
    155153  nt=-1 

  • branches/GRISLIv3/SOURCES/iso_declar_mod-0.3.f90

    r181 r485  
    2222 
    2323implicit none 
    24  
    25 INTEGER ::  NLITH    !< defini dans isostasie  permet choisir modele lithosphere 
    26                      !< 0 enfoncement local, 1 enfoncement regional 
    2724 
    2825INTEGER :: NBED      !< permet de choisir quel type d'isostasie est utilise 

  • branches/GRISLIv3/SOURCES/isostasie_mod-0.3.f90

    r481 r485  
    3232  subroutine init_iso ! routine qui permet d'initialiser les variables  
    3333 
    34     use iso_declar,only: nlith,dt_iso,tausoc,dl,rl,lbloc,we,charge 
     34    use iso_declar,only: dt_iso,tausoc,dl,rl,lbloc,we,charge 
    3535    use module3D_phy, only: err,h0,bsoc0,sealevel_2d 
    3636    use geography, only: geoplace,nx,ny,dx,dy 
    … …  
    4545 
    4646    NBED=1 
    47     NLITH=1 
    4847 
    4948    !      temps de reaction en annees 
    50     if (NBED.eq.1) tausoc=3000 
     49    tausoc=3000 
    5150 
    52     if (NLITH.eq.1) then 
    53        !        DL lithosphere flexural rigidity (N.m) 
    54        DL=9.87E24 
    55        !        radius of relative stiffness  (metre) 
    56        RL=131910. 
    57        !        LBLOC, 400 km de part et d'autre 
    58        LBLOC=int((400000.-0.1)/DX)+1 
    59        !         LBLOC=int((480000.-0.1)/DX)+1 
    60        !         LBLOC=int(400000./DX)+1 
     51    !        DL lithosphere flexural rigidity (N.m) 
     52    DL=9.87E24 
     53    !        radius of relative stiffness  (metre) 
     54    RL=131910. 
     55    !        LBLOC, 400 km de part et d'autre 
     56    LBLOC=int((400000.-0.1)/DX)+1 
    6157 
    6258 
    63        !******** allocation des tableaux en fonction de la valeur de LBLOC ***** 
     59    !******** allocation des tableaux en fonction de la valeur de LBLOC ***** 
    6460 
    65        if (.not.allocated(WE)) then 
    66           allocate(WE(-LBLOC:LBLOC,-LBLOC:LBLOC),stat=err) 
    67           if (err/=0) then 
    68              print *,"Erreur a l'allocation du tableau WE",err 
    69              stop 4 
    70           end if 
     61    if (.not.allocated(WE)) then 
     62       allocate(WE(-LBLOC:LBLOC,-LBLOC:LBLOC),stat=err) 
     63       if (err/=0) then 
     64          print *,"Erreur a l'allocation du tableau WE",err 
     65          stop 4 
    7166       end if 
    72        !    PRINT*,'shape(we)',SHAPE(we),'lbloc',lbloc 
     67    end if 
    7368 
    74        if (.not.allocated(CHARGE)) then 
    75           allocate(CHARGE(1-LBLOC:NX+LBLOC,1-LBLOC:NY+LBLOC),stat=err) 
    76           if (err/=0) then 
    77              print *,"Erreur a l'allocation du tableau CHARGE",err 
    78              stop 4 
    79           end if 
     69    if (.not.allocated(CHARGE)) then 
     70       allocate(CHARGE(1-LBLOC:NX+LBLOC,1-LBLOC:NY+LBLOC),stat=err) 
     71       if (err/=0) then 
     72          print *,"Erreur a l'allocation du tableau CHARGE",err 
     73          stop 4 
    8074       end if 
    81        !********* FIN DE L'ALLOCATION DES TABLEAUX *********** 
    82        call tab_litho 
    8375    end if 
     76    !********* FIN DE L'ALLOCATION DES TABLEAUX *********** 
     77     
     78    call tab_litho 
    8479 
    8580    !******** initialisation de CHARGE *********** 
    … …  
    9186       do j=1,ny 
    9287          if (ro*h0(i,j).ge.-row*(BSOC0(i,j)-sealevel_2d(i,j))) then 
    93  
    9488             !        glace ou terre 
    9589             charge(i,j)=rog*h0(i,j) 
    96  
    9790          else 
    9891             !        ocean 
    … …  
    125118 
    126119    !     deflection de la lithosphere 
    127     if (nlith.eq.1) then 
     120    call litho 
    128121 
    129        call litho 
    130  
    131        do j=1,ny 
    132           do i=1,nx 
    133              w0(i,j)=w1(i,j) 
    134           end do 
    135        end do 
    136     else 
    137        do j=1,ny 
    138           do i=1,nx 
    139              w0(i,j)=charge(i,j)/romg 
    140           end do 
    141        end do 
    142     end if 
     122    w0(:,:)=w1(:,:) 
    143123 
    144124    dt_iso = 50. 

  • branches/GRISLIv3/SOURCES/main3D-0.4-40km.f90

    r481 r485  
    134134                          ice,bm,bmelt,ablbord,ablbord_dtt,dt,    & 
    135135                          s,h,b,bsoc,flot,mk,mk0,uxbar,uybar,hwater,time,timemax,ndebug,ndebug_max 
    136   use runparam, only: nt,tbegin,dtprofile,dtcpt,dirnameout,itracebug 
    137   use geography, only: nx,ny,geoplace 
     136  use runparam, only: nt,tbegin,dirnameout,itracebug 
    138137  use deformation_mod_2lois, only:n1poly,n2poly 
    139138  use bilan_eau_mod, only: init_bilan_eau 
    140139  use module_choix, only: forclim,ablation,bmeltshelf,calving,flow_general,flowlaw 
    141   !  module_choix donne acces a tous les modules 
    142   !  de declaration des packages 
    143140  use flottab_mod, only: flottab 
    144   use sorties_ncdf_grisli, only: iglob_ncdf,testsort_time_ncdf,init_sortie_ncdf,testsort_time_ncdf, & 
    145                                  sortie_ncdf_cat 
    146   use util_recovery, only: dtout 
     141  use sorties_ncdf_grisli, only: iglob_ncdf,testsort_time_ncdf,init_sortie_ncdf,sortie_ncdf_cat 
    147142   
    148 !  use track_debug  
    149  
    150143  implicit none 
    151144 
    152   integer :: i,j 
    153  
    154145  if (itracebug.eq.1)  call tracebug(' Entree dans routine grisli_init') 
    155   !      switch pour passer ou non par T lliboutry calcule => 0, ne passe pas, 
    156   !      1 ou 2 passe (se met a 0 tout seul si on prend un fichier .cptr) 
    157  
    158   nt=-1   ! utilisee dans initialisation flottab 
    159   !     sortie profile tous les dtprofile 
    160   DTPROFILE=50000. 
    161   !     ----------------------------------fin des modifs run les plus usuelles 
     146 
    162147  DIRNAMEOUT='../RESULTATS/' 
    163   !DIRNAMEOUT='./' 
    164148 
    165149  call initial  ! routine qui appel toutes les routines d'initialisation 
    166  
    167  
    168   !      call init_sortie_ncdf 
    169   !      call sortie_ncdf_cat 
    170   !      STOP 
    171  
    172   !     compteur tous les DTCPT 
    173   DTCPT=dtout 
    174150 
    175151  !------------------------------ INITIALISATION ---------------------------- 
    176152  ! 
    177 ! ecriture netcdf apres initialisation 
    178  
    179  
    180  
     153  ! ecriture netcdf apres initialisation 
    181154  call testsort_time_ncdf(dble(tbegin)) 
    182155  if (iglob_ncdf .EQ. 1) call sortie_ncdf_cat 
    183156 
    184   ! test vincent car certains H(i,j)=0 dans fichier de reprise 
    185   do j=1,ny 
    186      do i=1,nx 
    187         H(i,j)=max(0.,H(i,j)) 
    188      enddo 
    189   enddo 
    190  
    191  
    192157  call forclim                   !  initialisation BM et TS          
    193158  call ablation 
    194159 
    195160 
    196  
    197161  !  -----------                  CALCULATION OF INITIAL TEMPERATURES  
    198162 
    … …  
    200164 
    201165           call masque(flot,mk,mk0,itracebug) 
    202  
    203            call Neffect() 
    204  
     166           call Neffect() 
    205167           call flottab() 
    206  
    207            call Neffect() 
    208  
    209  
    210            !          call vitbilan_lect   ! routine de lecture des vitesses de bilan 
    211            !       ======================================================== 
    212  
     168           call Neffect() 
    213169           call lineartemp() 
    214  
    215170           call bmelt_grounded  
    216            call  bmeltshelf 
    217  
     171           call bmeltshelf 
    218172 
    219173           call flow_general 
    220  
    221174           do iglen=n1poly,n2poly 
    222175              call flowlaw(iglen) 
    … …  
    232185           call SIA_velocities() 
    233186 
    234   else  ! tcpt     on reprend un fichier compteur (ICOMPTEUR.eq.1) 
     187  else  !     on reprend un fichier compteur (ICOMPTEUR.eq.1) 
    235188 
    236189     time=tbegin       ! prend le temps du compteur 
    237  
    238190 
    239191     call masque(flot,mk,mk0,itracebug) 
    … …  
    242194     call flottab() 
    243195     call masque(flot,mk,mk0,itracebug) 
    244  
    245      do i=1,nx 
    246         do j=1,ny 
    247            if (S(i,j).lt.0) then 
    248               print*,i,j,S(i,j) 
    249               goto 11115 
    250            endif 
    251         enddo 
    252      enddo 
    253 11115 continue 
    254  
    255      call  bmeltshelf ! afq --  
    256  
    257      !       ======================================================== 
     196     call bmeltshelf 
     197 
    258198     call flow_general 
    259  
    260199     do iglen=n1poly,n2poly 
    261200        call flowlaw(iglen) 
    262201     end do 
    263  
    264202 
    265203     call Neffect() 
    … …  
    272210  !     fin du test sur icompteur 
    273211 
    274   !      call init_sortie_ncdf 
    275   !      call sortie_ncdf_cat 
    276  
    277212  call flottab() 
    278213  call Neffect() 
    … …  
    280215 
    281216  if (icompteur.eq.0) then 
    282      do i=1,nx 
    283         do j=1,ny 
    284            if (.not.flot(i,j)) then 
    285               B(i,j) = Bsoc(i,j) 
    286               Uxbar(i,j) = 0. 
    287               Uybar(i,j) = 0. 
    288            end if 
    289         end do 
    290      end do 
     217     where(.not.flot(:,:)) 
     218        B(:,:) = Bsoc(:,:) 
     219        Uxbar(:,:) = 0. 
     220        Uybar(:,:) = 0. 
     221     endwhere 
    291222  endif 
    292223 
    293   do i=2,nx-1 
    294      do j=2,ny-1 
    295         hwater(i,j)=max(hwater(i,j),0.) 
    296      enddo 
    297   enddo 
    298224  timemax=time 
    299225  isynchro=1 

  • branches/GRISLIv3/SOURCES/runparam_mod.f90

    r468 r485  
    1818  real ::  TBEGIN                       !< temps debut 
    1919  real ::  TEND                         !< temps fin 
    20   real ::  dttest                       !< plus petit pas de temps=dtmin, sert dans plein de tests 
    2120  character(len=40)  :: DIRNAMEOUT      !< declaration du dir de sortie, pour les resultats 
    2221  character(len=8)   :: RUNNAME         !< 
    … …  
    2524  integer :: itracebug                  !< si 1 fait une impression au début des routines 
    2625  integer :: num_tracebug               !< numero de l'unite itracebug 
    27   real dtsortie 
    28   real dtcpt 
    29   real dtprofile 
    3026 
    3127end module runparam 

  • branches/GRISLIv3/SOURCES/steps_time_loop.f90

    r481 r485  
    2626       bilan_eau 
    2727  use bilan_flux_mod, only: bilan_flux_output 
    28   !  module_choix donne acces aux modules interchangeables 
    2928  use module_choix, only: shortoutput 
    3029 
    … …  
    3938  !============================= 
    4039 
     40  ! afq fev2024, a verifier : - cas nt < 0 - pourquoi ? a garder ? 
     41  !                           - variantes ispinup - toutes utiles ? 
    4142 
    4243  spinup_icethick: if (ispinup.eq.0.or.ispinup.eq.2.or.nt.lt.nt_init) then  
    … …  
    113114  call testsort_time_ncdf(time) 
    114115  if (iglob_ncdf .EQ. 1) call sortie_ncdf_cat 
    115   !  if ((isynchro.eq.1).or.(nt.eq.1).and.(mod(abs(TIME),1.).lt.dtmin)) then 
    116116  if ((nt.eq.1).or.((isynchro.eq.1).and.(mod(abs(TIME),1.).lt.dtmin))) then 
    117117     call shortoutput 
    … …  
    125125  endif 
    126126 
    127   !   sortie compteur tous les dtcpt ans  
     127  !   sortie compteur tous les xx ans  
    128128  !------------------------------------------------------------------ 
    129129  !iout == 1 sortie cptr 
    … …  
    177177 
    178178  integer :: m 
    179   logical :: shelf_vitbil = .true.   ! si vrai on prend les vitesses de bilan pour le shelf 
    180179  integer :: nt_init_tm = 0          ! number of loops for initialisation of thermo mechanical 
    181180  integer :: iter_visco              ! number of iterations for ssa viscosity 
    182181  real ::  test_iter_diag            ! test sur les vitesses pour iterations diagnostiques 
    183182 
    184   if (ispinup.le.1) shelf_vitbil = .false.         ! general case, ice shelves velocities are computed by diagnoshelf 
    185  
    186183  timemax=time 
     184 
    187185  if (itracebug.eq.1) write(num_tracebug,*) 'entree dans step_thermomeca', & 
    188186       '   nt=',nt 
    … …  
    196194  !   debut d'un bloc appels tous les dtt (pas de temps long- asynchrone) ! 
    197195  !-------------   -------------------------------------------------------! 
    198  
    199196 
    200197 
    … …  
    208205     ! because eustatic sea level is updated in forclim, we need to call the RSL here 
    209206     call rsl         ! provide the map of the local sea level (sealevel_2d) 
    210  
    211207 
    212208 
    … …  
    248244     ! update values in the structures Geom_g, Temperature_g, ... 
    249245 
    250      calc_temp:     if (geoplace(1:5).ne.'mism3') then 
    251         if (itracebug.eq.1)  call tracebug('avant appel icetemp') 
    252         call icetemp 
    253         ! subroutines pour le calcul de la fusion basale 
    254         call bmelt_grounded 
    255         call bmeltshelf 
    256  
    257      endif calc_temp 
     246     if (itracebug.eq.1)  call tracebug('avant appel icetemp') 
     247     call icetemp 
     248     ! subroutines pour le calcul de la fusion basale 
     249     call bmelt_grounded 
     250     call bmeltshelf 
    258251 
    259252 
    … …  
    309302 
    310303 
    311  
    312304     call Neffect() 
    313305 
    … …  
    327319 
    328320     isync_2:    if (isynchro.eq.1) then                 
    329  
    330         !les 3 lignes ci-dessous pour avoir les champs avant diagno (en cas de pb colineaire) 
    331         !call init_sortie_ncdf 
    332         !call sortie_ncdf_cat 
    333321 
    334322        if (icompteur.eq.1) then    ! reprise de tout le vecteur d'etat     

  • branches/GRISLIv3/SOURCES/taubed-0.3.f90

    r467 r485  
    4848  USE param_phy_mod, only: ro,row,rog,rowg,rom 
    4949  USE isostasie_mod, only:w0,w1 
    50   USE iso_declar,only: nlith,lbloc,charge,dt_iso,tausoc ! module de declaration des variables de l'isostasie 
     50  USE iso_declar,only: lbloc,charge,dt_iso,tausoc ! module de declaration des variables de l'isostasie 
    5151 
    5252  implicit none 
    … …  
    5555 
    5656  !    ********* calcul de W1 l'enfoncement d'equilibre au temps t 
    57   ! NLITH est defini dans isostasie et permet le choix du modele d'isostasie 
    5857 
    59    
    60   if (NLITH.eq.1) then 
    61      !       avec rigidite de la lithosphere 
    62      !$OMP PARALLEL 
    63      !$OMP DO 
    64      do J=1,NY  
    65         do I=1,NX 
    66            if (RO*H(I,J).ge.-ROW*(BSOC(I,J)-SEALEVEL_2D(i,j))) then 
    67               !           glace ou terre  
    68               CHARGE(I,J)=ROG*H(I,J) 
    69            else 
    70               !           ocean 
    71               CHARGE(I,J)=-ROWG*(BSOC(I,J)-SEALEVEL_2D(i,j)) 
    72            endif 
    73         end do 
     58  !       avec rigidite de la lithosphere 
     59  !$OMP PARALLEL 
     60  !$OMP DO 
     61  do J=1,NY  
     62     do I=1,NX 
     63        if (RO*H(I,J).ge.-ROW*(BSOC(I,J)-SEALEVEL_2D(i,j))) then 
     64           !           glace ou terre  
     65           CHARGE(I,J)=ROG*H(I,J) 
     66        else 
     67           !           ocean 
     68           CHARGE(I,J)=-ROWG*(BSOC(I,J)-SEALEVEL_2D(i,j)) 
     69        endif 
    7470     end do 
    75      !$OMP END DO 
     71  end do 
     72  !$OMP END DO 
    7673 
    77      ! il faut remplir CHARGE dans les parties a l'exterieur de la grille : 
    78      ! a l'exterieur de la grille CHARGE est egale a la valeur sur le bord 
    79      !$OMP DO 
    80      do J=1,NY 
    81         CHARGE(1-LBLOC:0,J)=CHARGE(1,J)      ! bord W 
    82         CHARGE(NX+1:NX+LBLOC,J)=CHARGE(NX,J) ! bord E 
    83      end do 
    84      !$OMP END DO 
    85      !$OMP DO 
    86      do I=1,NX 
    87         CHARGE(I,1-LBLOC:0)=CHARGE(I,1)      ! bord S 
    88         CHARGE(I,NY+1:NY+LBLOC)=CHARGE(I,NY) ! bord N 
    89      end do 
    90      !$OMP END DO 
     74  ! il faut remplir CHARGE dans les parties a l'exterieur de la grille : 
     75  ! a l'exterieur de la grille CHARGE est egale a la valeur sur le bord 
     76  !$OMP DO 
     77  do J=1,NY 
     78     CHARGE(1-LBLOC:0,J)=CHARGE(1,J)      ! bord W 
     79     CHARGE(NX+1:NX+LBLOC,J)=CHARGE(NX,J) ! bord E 
     80  end do 
     81  !$OMP END DO 
     82  !$OMP DO 
     83  do I=1,NX 
     84     CHARGE(I,1-LBLOC:0)=CHARGE(I,1)      ! bord S 
     85     CHARGE(I,NY+1:NY+LBLOC)=CHARGE(I,NY) ! bord N 
     86  end do 
     87  !$OMP END DO 
    9188 
    92      ! valeur dans les quatres coins exterieurs au domaine        
    93           !$OMP WORKSHARE 
    94      CHARGE(1-LBLOC:0,1-LBLOC:0)=CHARGE(1,1)           ! coin SW 
    95      CHARGE(1-LBLOC:0,NY+1:NY+LBLOC)=CHARGE(1,NY)      ! coin NW 
    96      CHARGE(NX+1:NX+LBLOC,1-LBLOC:0)=CHARGE(NX,1)      ! coin SE 
    97      CHARGE(NX+1:NX+LBLOC,NY+1:NY+LBLOC)=CHARGE(NX,NY) ! coin NE 
    98      !$OMP END WORKSHARE 
    99      !$OMP END PARALLEL  
    100      call litho 
    101  
    102   else 
    103      !     enfoncement local 
    104      !$OMP PARALLEL 
    105      !$OMP DO 
    106      do J=1,NY 
    107         do I=1,NX 
    108            if (RO*H(I,J).ge.-ROW*(BSOC(I,J)-SEALEVEL_2D(i,j))) then 
    109               !           glace ou terre  
    110               W1(I,J)=RO/ROM*H(I,J) 
    111            else 
    112               !           ocean 
    113               W1(I,J)=-ROW/ROM*(BSOC(I,J)-SEALEVEL_2D(i,j)) 
    114            endif 
    115         end do 
    116      end do 
    117      !$OMP END DO 
    118      !$OMP END PARALLEL 
    119   endif 
     89  ! valeur dans les quatres coins exterieurs au domaine        
     90  !$OMP WORKSHARE 
     91  CHARGE(1-LBLOC:0,1-LBLOC:0)=CHARGE(1,1)           ! coin SW 
     92  CHARGE(1-LBLOC:0,NY+1:NY+LBLOC)=CHARGE(1,NY)      ! coin NW 
     93  CHARGE(NX+1:NX+LBLOC,1-LBLOC:0)=CHARGE(NX,1)      ! coin SE 
     94  CHARGE(NX+1:NX+LBLOC,NY+1:NY+LBLOC)=CHARGE(NX,NY) ! coin NE 
     95  !$OMP END WORKSHARE 
     96  !$OMP END PARALLEL  
     97  call litho 
    12098 
    12199  !     decroissance exponentielle de l'enfoncement 

  • branches/GRISLIv3/SOURCES/util_recovery.f90

    r481 r485  
    4747    end do comment1 
    4848 
    49     ! lecture de dtcpt 
     49    ! lecture de frequence des restarts 
    5050    read(num_file,*) dtout 
    5151    read(num_file,*)       !saut de la ligne "-------" 
Note: See TracChangeset for help on using the changeset viewer.