Changeset 486


Ignore:
Timestamp:
02/13/24 18:13:56 (5 months ago)
Author:
aquiquet
Message:

Cleaning branch: some cleaning and checking of passive tracers

Location:
branches/GRISLIv3
Files:
1 added
5 edited

Legend:

Unmodified
Added
Removed

  • branches/GRISLIv3/SOURCES/Makefile.grisli.inc

    r478 r486  
    594594        $(mod_communs) \ 
    595595        $(mod_clim_tof) \ 
    596         $(mod_no_tracers) \ 
     596        $(mod_tracers) \ 
    597597        $(mod_ell) $(Liste_Ant40) \ 
    598598        $(diagnoshelf) \ 
    … …  
    606606        $(mod_communs) \ 
    607607        $(mod_clim_tof) \ 
    608         $(mod_no_tracers) \ 
     608        $(mod_tracers) \ 
    609609        $(mod_ell) $(Liste_Ant40) \ 
    610610        $(diagnoshelf) \ 

  • branches/GRISLIv3/SOURCES/celltest_tracer.f90

    r481 r486  
    2020 
    2121      integer :: kk   ! indice vertical pour definition de E, mais on veut conserver la valeur de k 
    22 !      real,dimension(nz) :: E   ! vertical coordinate in ice, scaled to H zeta 
    23  
    24       E(1)=0. 
    25       E(NZ)=1. 
    26       do KK=1,NZ 
    27          if ((KK.ne.1).and.(KK.ne.NZ)) E(KK)=(KK-1.)/(NZ-1.) 
    28       end do 
    29  
    3022 
    3123      if(v_x>0.0) then 

  • branches/GRISLIv3/SOURCES/interpolate_tracer.f90

    r481 r486  
    44  use module3d_phy, only:e,H 
    55  use tracer_vars, only:tdepk,tdep,time_max_accu,accucumul,xdepk,ydepk,xdep,ydep 
    6 !cdc  use climat_perturb_mois_mod ! on en a besoin, notamment pour NFT 
    76 
    87  implicit none 
    … …  
    1615  
    1716  integer :: kk   ! indice vertical pour definition de E, mais on veut conserver la valeur de k 
    18 !  real,dimension(nz) :: E   ! vertical coordinate in ice, scaled to H zeta  
    19    
     17  
    2018  !! Note that input (im,jm,km) refer to points on (ij,jj) grid 
    2119  !! whereas CellTest search was performed on (i,j) grid 
    … …  
    7371        al  = abs(tdepk(kc,ic,jc) - a_tmp) 
    7472        ac  = abs(tdepk(kc,ic,jc) - (a_xy + fzz*fprime + fzz**2*fsec/2.)) 
    75         if (ac<(1.3*al)) then    
     73        if (ac<(1.3*al)) then 
    7674           tdep(kc,ic,jc) = a_xy + fzz*fprime + fzz**2*fsec/2. 
    7775        else 
    … …  
    8381  else 
    8482      
    85      index1  = floor(-a_xy/100.)           ! gives the youngest: ind1+1 < tdepk < ind1 
     83     index1  = floor(-a_xy/100.)           ! gives the youngest: ind1+1 < tdepk < ind1 
    8684     fra1 = real((-a_xy/100.) - index1) 
    8785     accucumul1 = real(  -fra1*(accucumul(index1)-accucumul(index1+1))  +  accucumul(index1) ) 
    … …  
    127125        L4: do 
    128126           index0 = index0 + 1 
    129            if ((index0>=NFT-1).or.(accucumul(index0)  < f_k )) exit L4   
     127           if ((index0>=NFT-1).or.(accucumul(index0)  < f_k )) exit L4 
    130128           ! index < ind(t) < index-1,  t(i+1) < t < t(i) 
    131129        end do L4 
    … …  
    133131         
    134132        !========================================= 
    135      else                               !  fz> 0.5, probably downwards 
     133     else                               !  fz> 0.5, probably downwards 
    136134         
    137135        if (km<nz-1) then 
    … …  
    164162      
    165163     if ((index0==1).or.(index0==NFT-1)) then 
    166         fzz=2. 
     164        fzz=2. 
    167165     else   
    168166        fzz = real( (accucumul(index0) - f_k)/ (accucumul(index0)-accucumul(index0+1)) ) 

  • branches/GRISLIv3/SOURCES/tracer_mod.f90

    r465 r486  
    2626 
    2727module tracer_mod 
    28  
    29   use module3d_phy, only: e,time,dtt,num_param,num_rep_42,acc,H,num_forc,& 
    30        ux,uy,uzr,bmelt,S,flot,bm 
    31   use geography, only: nx,ny,nz,dx,dy,dirforcage 
    32   use runparam, only: itracebug,num_tracebug,xmin,ymin,tbegin,tend 
    33   ! module de declaration de variables pour les traceurs 
    34   use tracer_vars, only:file_tr_dat,file_tr_out,file_tr_dep,file_tr_dat,coeft_tra,rappact_tra,nij,& 
    35        xgrid,ygrid,xdepk,ydepk,tdepk,nft_tra,tdate_tra,tpert_tra,accucumul,type_accum,time_max_accu,& 
    36        uxsave,uysave,uzrsave,xdep,ydep,tdep,freezeon,xdep_out,ydep_out,tdep_out 
    37 !cdc  use climat_perturb_mois_mod ! on en a besoin, notamment pour NFT 
    38 ! afq, new version of tracer, does not require climat_perturb 
    3928   
    4029  implicit none 
    41  
    42 ! variables de travail 
    43   integer :: i_ij, i_jj, ip, jp, im, jm, km      ! indice dans les matrices des traceurs 
    44   real    :: v_xij, v_yij, v_zij                 ! vitesses 'cumulees' 
    4530 
    4631! pour chaque point de grille on recherche ou etait la particule au pas de temps precedent 
    … …  
    6752subroutine init_tracer 
    6853 
     54  use module3d_phy, only: e,time,dtt,num_param,num_rep_42,acc,H,num_forc 
     55  use geography, only: nx,ny,nz,dx,dy,dirforcage 
     56  use runparam, only: itracebug,num_tracebug,xmin,ymin,tbegin,tend 
     57  use tracer_vars, only:file_tr_dat,coeft_tra,rappact_tra,nij,& 
     58       xgrid,ygrid,xdepk,ydepk,tdepk,nft_tra,tdate_tra,tpert_tra,accucumul,type_accum,time_max_accu,& 
     59       uxsave,uysave,uzrsave,xdep,ydep,tdep,freezeon 
     60 
     61implicit none 
     62 
    6963! formats pour les ecritures dans 42 
    7064428 format(A) 
     65 
     66! variables de travail 
     67  integer :: i_ij, i_jj 
    7168 
    7269  real :: lambdab  ! fusion basale 'artificielle' pour age des couches  
    … …  
    8481 
    8582  ! on lit les parametres par namelist dorenavant 
    86   namelist/tracer/file_tr_dat,file_tr_out,file_tr_dep,lambdab,agemax,coefT_tra,rappact_tra,filforc,pert_type  
     83  namelist/tracer/file_tr_dat,lambdab,agemax,coefT_tra,rappact_tra,filforc,pert_type  
    8784 
    8885 
    … …  
    109106  print* 
    110107  print*, 'tracer .dat file: ', file_tr_dat 
    111   print*, 'tracer .out file: ', file_tr_out 
    112   print*, 'tracer dep. file: ', file_tr_dep 
    113108 
    114109 
    … …  
    116111  write(num_rep_42,428) '&tracer                                    ! module tracer_mod' 
    117112  write(num_rep_42,'(A,A)')   'tracer .dat file : ', file_tr_dat 
    118   write(num_rep_42,'(A,A)')   'tracer .out file : ', file_tr_out 
    119   write(num_rep_42,'(A,A)')   'tracer dep. file : ', file_tr_dep 
    120113  write(num_rep_42,'(A,A)')   'lambdab : ', lambdab 
    121114  write(num_rep_42,'(A,A)')   'agemax : ', agemax 
    … …  
    302295subroutine tracer 
    303296 
     297  use geography, only: nx,ny,nz 
     298  use module3D_phy, only: E,time,dtt,ux,uy,uzr,h,s,flot,bmelt,bm 
     299  use runparam, only: itracebug,num_tracebug,tend 
     300  use tracer_vars, only: xgrid,ygrid,xdepk,ydepk,tdepk,xdep,ydep,tdep,xdep_out,ydep_out,tdep_out, & 
     301                         uxsave,uysave,uzrsave,freezeon,nft_tra 
     302 
     303 
     304implicit none  
     305 
     306! variables de travail 
     307  integer :: i_ij, i_jj , ip, jp , im, jm, km      ! indice dans les matrices des traceurs 
     308  real    :: v_xij, v_yij, v_zij                   ! vitesses 'cumulees' 
     309 
     310  real,parameter   :: petit = 0.1 ! afq fev2024 - a smaller number breaks the test which is weird... 
     311 
    304312  integer :: i,j,k 
    305   real,dimension(nz) :: E   ! vertical coordinate in ice, scaled to H zeta 
    306  
    307   E(1)=0. 
    308   E(NZ)=1. 
    309   do K=1,NZ 
    310      if ((K.ne.1).and.(K.ne.NZ)) E(K)=(K-1.)/(NZ-1.) 
    311   end do 
    312313 
    313314  If (Itracebug.Eq.1) Write(num_tracebug,*)' Entree Dans Routine tracer at time=',Time 
    … …  
    353354!        tdep(:,1,i_jj) = time_ 
    354355!        tdep(:,nx-2,i_jj) = time_ 
    355         tdep(:,1,i_jj) = time     
    356         tdep(:,nx-2,i_jj) = time  
     356        tdep(:,1,i_jj) = real(time) 
     357        tdep(:,nx-2,i_jj) = real(time) 
    357358     end do 
    358359      !=========== north and south boundaries 
    … …  
    365366!       tdep(:,i_ij,1) =  time_ 
    366367!       tdep(:,i_ij,ny-2) =  time_ 
    367         tdep(:,i_ij,1) =  time    
    368         tdep(:,i_ij,ny-2) =  time 
     368        tdep(:,i_ij,1) =  real(time) 
     369        tdep(:,i_ij,ny-2) =  real(time) 
    369370     end do 
    370371 
    … …  
    380381            xdep(:,i_ij,i_jj) = xgrid(i) 
    381382            ydep(:,i_ij,i_jj) = ygrid(j) 
    382             tdep(:,i_ij,i_jj) =  time   
     383            tdep(:,i_ij,i_jj) =  real(time) 
    383384            cycle 
    384385          end if 
    … …  
    411412                ! *1.001        proven to be dumb 
    412413!               if((tdep(k,i_ij,i_jj)<=time_  -  1500000.) .or. (tdep(k,i_ij,i_jj)>=tend+1.0)) then 
    413                if((tdep(k,i_ij,i_jj)<=time  -  1500000.) .or. (tdep(k,i_ij,i_jj)>=tend+1.0)) then  
     414               if((tdep(k,i_ij,i_jj)<=real(time)  -  1500000.) .or. (tdep(k,i_ij,i_jj)>=tend+1.0)) then  
    414415                  print *,"tdep(",k,",",i_ij,",",i_jj,") range error at time=",time_tra 
    415416                   print *,"tdep(k,i_ij,i_jj) =",tdep(k,i_ij,i_jj), v_xij, v_yij, v_zij 
    … …  
    420421                end if 
    421422 
    422              else if (k==1 .and. bm(i,j)>0.0) then        !! surface mass balance 
     423             else if (k==1 .and. bm(i,j)>0.0) then  !! surface mass balance 
    423424              xdep(k,i_ij,i_jj) = xgrid(i) 
    424425              ydep(k,i_ij,i_jj) = ygrid(j) 
    425               tdep(k,i_ij,i_jj) =  time   
    426             else if ( e_tra < e(1) ) then               !! rapid thickening by 
    427               xdep(k,i_ij,i_jj) = xgrid(i)                              !! surface accumulation 
     426              tdep(k,i_ij,i_jj) =  real(time) 
     427            else if ( e_tra < e(1) ) then       !! rapid thickening by 
     428              xdep(k,i_ij,i_jj) = xgrid(i)      !! surface accumulation 
    428429              ydep(k,i_ij,i_jj) = ygrid(j) 
    429               tdep(k,i_ij,i_jj) =  time    
     430              tdep(k,i_ij,i_jj) =  real(time) 
    430431           else 
    431432                call celltest(time_tra,i,j,k,v_xij, v_yij, v_zij, x_tra, y_tra, e_tra, & 
    432                                 im,jm,km,fx,fy,fz) 
    433                 !! convert seach results from (i,j) grid to (i_ij,i_jj) tracer grid 
     433                                im,jm,km,fx,fy,fz) 
     434                !! convert seach results from (i,j) grid to (i_ij,i_jj) tracer grid 
    434435                !print*,'i,j,k avt interp',i,j,k 
    435436 
    436437                call interpolate(i_ij, i_jj, k, im-1, jm-1, km, fx,fy,fz, e_tra, nft_tra)  
    437438 
    438                 if((tdep(k,i_ij,i_jj)<=time  -  1500000.) .or. (tdep(k,i_ij,i_jj)> time   )) then 
     439                if((tdep(k,i_ij,i_jj)<=real(time)  -  1500000.) .or. (tdep(k,i_ij,i_jj)> real(time)+petit)) then 
    439440                   print *,"tdep(",k,",",i_ij,",",i_jj,") range error at time=",time_tra 
    440441                   print *,"tdep(k,i_ij,i_jj) =",tdep(k,i_ij,i_jj), v_xij, v_yij, v_zij 
    441                    print *, bmelt(i,j), s(i,j), h(i,j), e_tra, tdepk(k,i_ij,i_jj), tdepk(k,i,j) 
    442                    print *, flot(im,jm), flot(im+1,jm), flot(im+1,jm+1), flot(im,jm+1)  
    443                    if (k>1) tdep(k,i_ij,i_jj) = tdepk(k-1,i_ij,i_jj) - real(h(i,j)*(k-1)**2/30.) 
    444                    print *, tdep(k,i_ij,i_jj)  
     442                   print *, bmelt(i,j), s(i,j), h(i,j), e_tra, tdepk(k,i_ij,i_jj), tdepk(k,i,j) 
     443                   print *, flot(im,jm), flot(im+1,jm), flot(im+1,jm+1), flot(im,jm+1)  
     444                   if (k>1) tdep(k,i_ij,i_jj) = tdepk(k-1,i_ij,i_jj) - real(h(i,j)*(k-1)**2/30.) 
     445                   print *, tdep(k,i_ij,i_jj)  
    445446                end if 
    446447 
    … …  
    460461   end if !========== if deltracer 
    461462 
    462    if (abs(time-tend)<19) then    
     463   if (abs(real(time)-tend)<19) then    
    463464      open(505,file="freezeonmap.out",status="replace") 
    464465      print*,"save freezeon file" 
    … …  
    516517 
    517518 
    518  
    519  
    520  
    521  
    522  
    523  
    524  
    525  
    526  
    527  
    528  
    529  
    530 !======================================================================== 
    531 !======================================================================== 
    532 subroutine get_date_time(yr, mon, iday, hr, minu, sec) 
    533  
    534   implicit none 
    535   integer, intent(out) :: yr, mon, iday, hr, minu, sec 
    536  
    537   integer, dimension(8) :: values 
    538   character :: date*8, chtime*10, zone*5 
    539       
    540   call date_and_time(date, chtime, zone, values) 
    541       
    542   yr   = values(1) 
    543   mon  = values(2) 
    544   iday = values(3) 
    545   hr   = values(5) 
    546   minu  = values(6) 
    547   sec  = values(7) 
    548  
    549 end subroutine get_date_time 
    550  
    551519!======================================================================== 
    552520subroutine read_tr_dat 

  • branches/GRISLIv3/SOURCES/tracer_vars_mod.f90

    r465 r486  
    1717   
    1818  integer, parameter :: unit_tr_dat=572 
    19   integer, parameter :: unit_tr_out=573 
    20   integer, parameter :: unit_tr_dep=574 
    2119   
    2220  integer, parameter :: maxspectimes=31 
    … …  
    6866  real , dimension(nz) :: lnzeta 
    6967 
    70   character (LEN=60) :: file_tr_dat, file_tr_out, file_tr_dep 
     68  character (LEN=60) :: file_tr_dat 
    7169   
    7270  real , dimension(nx,ny,nz) :: uxsave 
Note: See TracChangeset for help on using the changeset viewer.