Changeset 429 for branches/GRISLIv3

Timestamp:

04/25/23 18:08:16 (15 months ago)

Author:

dumas

Message:

Use only in module spinup_vitbil

File:

• branches/GRISLIv3/SOURCES/spinup_mod.f90 (modified) (27 diffs)

branches/GRISLIv3/SOURCES/spinup_mod.f90

@@ -19 +19 @@
 module spinup_vitbil
 
-  use module3D_phy
-  use param_phy_mod
-  use deformation_mod_2lois
-  use interface_input
-  use io_netcdf_grisli
-  
+  use geography, only: nx,ny,nz
+
   implicit none
 
@@ -36 +32 @@
   real, dimension(nx,ny) :: init_coefmy      !< rescalling coefficient before limitation
 
-! compatibilites avec remplimat
+  ! compatibilites avec remplimat
   logical :: corr_def = .False.               !< for deformation correction (compatibility)
 
   real, dimension(nx,ny) :: Vsl_x             !<
   real, dimension(nx,ny) :: Vsl_y             !<
-  
+
   integer :: type_vitbil ! defini le type de fichier a lire : lect_vitbil ou lect_vitbil_Lebrocq
 
@@ -51 +47 @@
   !<
   subroutine init_spinup
-  namelist/spinup/ispinup,type_vitbil
+
+    use module3D_phy, only: ispinup,num_param,num_rep_42
+    use runparam, only : itracebug
+    
+    namelist/spinup/ispinup,type_vitbil
 
     ! put here which type of spinup
@@ -63 +63 @@
     ! ispinup = 3                          ! fait le calcul des temperatures 
     !                                      ! avec les vitesses de bilan
-  
-
-
-
-! lecture des parametres
-
-rewind(num_param)        ! pour revenir au debut du fichier param_list.dat
-read(num_param,spinup)
-
-write(num_rep_42,*)'!___________________________________________________________'
-write(num_rep_42,*)'!        spinup          module spinup_vitbil               '
-write(num_rep_42,spinup)                   ! pour ecrire les parametres lus
-write(num_rep_42,*)
-write(num_rep_42,*)
-write(num_rep_42,*)'! ispinup = 0     run standard voir no_spinup'
-write(num_rep_42,*)'! ispinup = 1     equilibre temperature avec vitesses grisli voir no_spinup'
-write(num_rep_42,*)
-write(num_rep_42,*)'! ispinup >1      use spinup_vitbil'
-write(num_rep_42,*)'! ispinup = 2     conservation de la masse avec vitesses bilan '
-write(num_rep_42,*)'! ispinup = 3     equilibre temperature avec vitesses bilan'
-write(num_rep_42,*)'! type_vitbil type de lecture des vitesses de bilan 1: lect_vitbil, 2 : lect_vitbil_Lebrocq'
-write(num_rep_42,*)
-
- 
-if (ispinup.eq.0) then
-   write(6,*)' ispinup = 0 et ispinup = 1 doivent etre appeles par no_spinup '
-   write(6,*) 'et il faut rajouter un dragging'
-   stop
-endif
-
-if (ispinup.eq.3) then
-  select case (type_vitbil)
-  case (1)
-!cdc Methode Cat lecture fichier selon x et y sur grille stag
-    call lect_vitbil
-  case(2)
-!cdc Methode Tof lecture fichier Lebrocq vitesse et direction
-    call lect_vitbil_Lebrocq
-  case default
-    print*,'type_vitbil valeur invalide dans spinup_vitbil'
-    stop
-  end select
-endif
+
+
+
+
+    ! lecture des parametres
+
+    rewind(num_param)        ! pour revenir au debut du fichier param_list.dat
+    read(num_param,spinup)
+
+    write(num_rep_42,*)'!___________________________________________________________'
+    write(num_rep_42,*)'!        spinup          module spinup_vitbil               '
+    write(num_rep_42,spinup)                   ! pour ecrire les parametres lus
+    write(num_rep_42,*)
+    write(num_rep_42,*)
+    write(num_rep_42,*)'! ispinup = 0     run standard voir no_spinup'
+    write(num_rep_42,*)'! ispinup = 1     equilibre temperature avec vitesses grisli voir no_spinup'
+    write(num_rep_42,*)
+    write(num_rep_42,*)'! ispinup >1      use spinup_vitbil'
+    write(num_rep_42,*)'! ispinup = 2     conservation de la masse avec vitesses bilan '
+    write(num_rep_42,*)'! ispinup = 3     equilibre temperature avec vitesses bilan'
+    write(num_rep_42,*)'! type_vitbil type de lecture des vitesses de bilan 1: lect_vitbil, 2 : lect_vitbil_Lebrocq'
+    write(num_rep_42,*)
+
+
+    if (ispinup.eq.0) then
+       write(6,*)' ispinup = 0 et ispinup = 1 doivent etre appeles par no_spinup '
+       write(6,*) 'et il faut rajouter un dragging'
+       stop
+    endif
+
+    if (ispinup.eq.3) then
+       select case (type_vitbil)
+       case (1)
+          !cdc Methode Cat lecture fichier selon x et y sur grille stag
+          call lect_vitbil
+       case(2)
+          !cdc Methode Tof lecture fichier Lebrocq vitesse et direction
+          call lect_vitbil_Lebrocq
+       case default
+          print*,'type_vitbil valeur invalide dans spinup_vitbil'
+          stop
+       end select
+    endif
 
     if (itracebug.eq.1)  call tracebug(' fin routine init_spinup de spinup_vitbil')
@@ -111 +111 @@
 
   subroutine lect_vitbil
-
+    
+    use module3D_phy, only: debug_3D,num_rep_42,num_param
+    use io_netcdf_grisli, only: Read_Ncdf_var
+    use geography, only: dirnameinp
+    use runparam, only : itracebug
+    
     character(len=100) :: balance_Ux_file        !  balance velocity file Ux
     character(len=100) :: balance_Uy_file        !  balance velocity file Uy
-          real*8, dimension(:,:),   pointer    :: tab               !< tableau 2d real pointer
-          
+    real*8, dimension(:,:),   pointer    :: tab               !< tableau 2d real pointer
+
     namelist/vitbil_upwind/balance_Ux_file, balance_Uy_file  
 
@@ -135 +140 @@
 
     ! read balance velocities on staggered nodes
-    
+
     balance_Ux_file = trim(dirnameinp)//trim(balance_Ux_file)
     balance_Uy_file = trim(dirnameinp)//trim(balance_Uy_file)
-    
+
     call Read_Ncdf_var('Vcol_x',balance_Ux_file,tab)
     Vcol_x(:,:)=tab(:,:)
@@ -144 +149 @@
     Vcol_y(:,:)=tab(:,:)
 
-!    call lect_input(2,'Vcol',1,Vcol_x,balance_vel_file,trim(dirnameinp)//trim(runname)//'.nc')  !read Vcol_x
-!    call lect_input(2,'Vcol',2,Vcol_y,balance_vel_file,trim(dirnameinp)//trim(runname)//'.nc')  !read Vcol_y
+    !    call lect_input(2,'Vcol',1,Vcol_x,balance_vel_file,trim(dirnameinp)//trim(runname)//'.nc')  !read Vcol_x
+    !    call lect_input(2,'Vcol',2,Vcol_y,balance_vel_file,trim(dirnameinp)//trim(runname)//'.nc')  !read Vcol_y
     !call lect_datfile(nx,ny,Vcol_x,1,balance_vel_file)   ! read Vcol_x   
     !call lect_datfile(nx,ny,Vcol_y,2,balance_vel_file)   ! read Vcol_y 
@@ -153 +158 @@
 
   end subroutine lect_vitbil
-  
+
   subroutine lect_vitbil_Lebrocq
 
-
+    use module3D_phy, only: num_rep_42,mstream_mx,itracebug,num_param,pi
+    use io_netcdf_grisli, only: Read_Ncdf_var
+    use geography, only: dirnameinp
+    use runparam, only : itracebug
+    
     character(len=100) :: vit_balance_file        !  balance velocity file
     character(len=100) :: flowdir_balance_file    !  balance flow direction file
@@ -164 +173 @@
     real, dimension(nx,ny) :: Vx_Lebrocq, Vy_Lebrocq ! vitesses selon x et y
     integer, dimension(nx,ny) :: v_good  !points avec vitesse Lebrocq OK
-    
+    integer :: i,j
+
     namelist/vitbil_Lebrocq/vit_balance_file, flowdir_balance_file  
 
@@ -184 +194 @@
 
     ! read balance velocities on staggered nodes
-    
+
     vit_balance_file = trim(dirnameinp)//trim(vit_balance_file)
     flowdir_balance_file = trim(dirnameinp)//trim(flowdir_balance_file)
-    
+
     call Read_Ncdf_var('z',vit_balance_file,tab)
     V_Lebrocq(:,:)=tab(:,:)
@@ -194 +204 @@
 
     where (V_Lebrocq(:,:).ge.0.)
-      v_good(:,:)=1
+       v_good(:,:)=1
     elsewhere
-      v_good(:,:)=0
+       v_good(:,:)=0
     endwhere
-    
+
     do j=2,ny-1
-      do i=2,nx-1
-        if (V_Lebrocq(i,j).lt.0.) then
-          
-          V_Lebrocq(i,j)= (V_Lebrocq(i-1,j)*v_good(i-1,j) + V_Lebrocq(i+1,j)*v_good(i+1,j) + &
-                         V_Lebrocq(i,j-1)*v_good(i,j-1) + V_Lebrocq(i,j+1)*v_good(i,j+1)) / &
-                         (v_good(i-1,j)+v_good(i+1,j)+v_good(i,j-1)+v_good(i,j-1))
-        endif
-      enddo
+       do i=2,nx-1
+          if (V_Lebrocq(i,j).lt.0.) then
+
+             V_Lebrocq(i,j)= (V_Lebrocq(i-1,j)*v_good(i-1,j) + V_Lebrocq(i+1,j)*v_good(i+1,j) + &
+                  V_Lebrocq(i,j-1)*v_good(i,j-1) + V_Lebrocq(i,j+1)*v_good(i,j+1)) / &
+                  (v_good(i-1,j)+v_good(i+1,j)+v_good(i,j-1)+v_good(i,j-1))
+          endif
+       enddo
     enddo
 
 
     where (V_Lebrocq(:,:).ge.100.)
-      V_Lebrocq(:,:)=100.
+       V_Lebrocq(:,:)=100.
     endwhere
 
-                ! calcul des vitesses selon x et selon y :
-                where (V_Lebrocq(:,:).ge.0.)
-                        Vx_Lebrocq(:,:) = V_Lebrocq(:,:)*sin(flowdir_Lebrocq(:,:))
-                        Vy_Lebrocq(:,:) = V_Lebrocq(:,:)*cos(flowdir_Lebrocq(:,:))
-                elsewhere
-                        Vx_Lebrocq(:,:) = 0.
-                        Vy_Lebrocq(:,:) = 0.
-                endwhere
-
-                
-                ! calcul des vitesses sur les points staggered
-                
-                do j=2,ny
-                        do i=2,nx
-                                Vcol_x(i,j) = (Vx_Lebrocq(i-1,j) + Vx_Lebrocq(i,j))/2.
-                                Vcol_y(i,j) = (Vy_Lebrocq(i,j-1) + Vy_Lebrocq(i,j))/2.
-                        enddo
-                enddo
-
-  end subroutine lect_vitbil_Lebrocq 
+    ! calcul des vitesses selon x et selon y :
+    where (V_Lebrocq(:,:).ge.0.)
+       Vx_Lebrocq(:,:) = V_Lebrocq(:,:)*sin(flowdir_Lebrocq(:,:))
+       Vy_Lebrocq(:,:) = V_Lebrocq(:,:)*cos(flowdir_Lebrocq(:,:))
+    elsewhere
+       Vx_Lebrocq(:,:) = 0.
+       Vy_Lebrocq(:,:) = 0.
+    endwhere
+
+
+    ! calcul des vitesses sur les points staggered
+
+    do j=2,ny
+       do i=2,nx
+          Vcol_x(i,j) = (Vx_Lebrocq(i-1,j) + Vx_Lebrocq(i,j))/2.
+          Vcol_y(i,j) = (Vy_Lebrocq(i,j-1) + Vy_Lebrocq(i,j))/2.
+       enddo
+    enddo
+
+  end subroutine lect_vitbil_Lebrocq
   !____________________________________________________________________________________
 
@@ -242 +252 @@
   subroutine init_dragging      ! Cette routine fait l'initialisation du dragging.
 
+    use module3D_phy, only: mstream_mx,mstream_my,drag_mx,drag_my,gzmx,gzmy,flgzmx,flgzmy,flotmx,flotmy
+    
     implicit none
 
@@ -265 +277 @@
 
   subroutine force_balance_vel
-
+    
+    use module3D_phy, only: uxbar, uybar,debug_3D
+    use runparam, only : itracebug
     if (itracebug.eq.1)  call tracebug(' Subroutine force_balance_vel')
 
@@ -281 +295 @@
   !<
   subroutine calc_coef_vitbil
-
+    
+    use module3D_phy, only: ddbx,ddby,gzmx,gzmy,flgzmx,flgzmy,fleuvemx,fleuvemy,flot,flotmx,flotmy,&
+         coefmxbmelt,coefmybmelt,sdx,sdy,uxbar,uybar,uxdef,uydef,ubx,uby,debug_3D,iglen
+    use deformation_mod_2lois, only: n1poly,n2poly,sa_mx,sa_my,s2a_mx,s2a_my
+    use runparam, only : itracebug
+    
+    integer :: i,j,k
+    
     ddbx(:,:)=0.
     ddby(:,:)=0.
@@ -290 +311 @@
     fleuvemx(:,:)=.false.
     fleuvemy(:,:)=.false.
-    
+
     do j=2,ny
-      do i=2,nx
+       do i=2,nx
           if (flot(i,j).and.(flot(i-1,j))) then
-            flotmx(i,j)=.true.
-            flgzmx(i,j)=.true.
-            gzmx(i,j)=.false.
+             flotmx(i,j)=.true.
+             flgzmx(i,j)=.true.
+             gzmx(i,j)=.false.
           end if
           if (flot(i,j).and.(flot(i,j-1))) then
-            flotmy(i,j)=.true.
-            flgzmy(i,j)=.true.
-            gzmy(i,j)=.false.
+             flotmy(i,j)=.true.
+             flgzmy(i,j)=.true.
+             gzmy(i,j)=.false.
           end if
-      end do
-    end do   
-    
+       end do
+    end do
+
     where (coefmxbmelt(:,:).le.0.)
-      gzmx(:,:)=.false.
+       gzmx(:,:)=.false.
     endwhere
     where (coefmybmelt(:,:).le.0.)
-      gzmy(:,:)=.false.
+       gzmy(:,:)=.false.
     endwhere
     flgzmx(:,:) = flgzmx(:,:) .or. gzmx(:,:)
@@ -316 +337 @@
     fleuvemx(:,:)=gzmx(:,:)
     fleuvemy(:,:)=gzmy(:,:)
-    
+
     if (itracebug.eq.1)  call tracebug(' Subroutine  calc_coef_vitbil')
     call slope_surf
@@ -341 +362 @@
              coef_defmx(i,j) = 1.
              Uxbar(i,j)      = Vcol_x(i,j)
-! dmr below is a cast from a real*4 to logical*4
-! dmr cannot be implicit in gfortran
-! dmr              flgzmx(i,j)     = Vcol_x(i,j)
+             ! dmr below is a cast from a real*4 to logical*4
+             ! dmr cannot be implicit in gfortran
+             ! dmr              flgzmx(i,j)     = Vcol_x(i,j)
              flgzmx(i,j)     = (nint(Vcol_x(i,j)).ne.0)
              uxdef(i,j)      = 0.
@@ -365 +386 @@
                    coef_defmx(i,j)=Vcol_x(i,j)/Ux_deformation(i,j)
                    ddbx(i,j) = 0
-                    Ubx(i,j)  = 0.
+                   Ubx(i,j)  = 0.
                 else
                    coef_defmx(i,j)=1.                                  ! no rescaling of deformation
@@ -385 +406 @@
              coef_defmy(i,j) = 1.
              Uybar(i,j)      = Vcol_y(i,j)
-! dmr below is a cast from a real*4 to logical*4
-! dmr cannot be implicit in gfortran
-! dmr             flgzmy(i,j)     = Vcol_y(i,j)
+             ! dmr below is a cast from a real*4 to logical*4
+             ! dmr cannot be implicit in gfortran
+             ! dmr             flgzmy(i,j)     = Vcol_y(i,j)
              flgzmy(i,j)     = (nint(Vcol_y(i,j)).ne.0)
 
@@ -393 +414 @@
              Uby(i,j)        = Vcol_y(i,j)
 
- if ((itracebug.eq.1).and.(i.eq.313).and.(j.eq.143))  call tracebug(' test0')
+             if ((itracebug.eq.1).and.(i.eq.313).and.(j.eq.143))  call tracebug(' test0')
           else
              coldy:      if ((coefmybmelt(i,j).le.0.).and.(.not.gzmy(i,j))) then   !base froide 
@@ -402 +423 @@
                 if (abs(Uy_deformation(i,j)).gt.0.01) then             ! calcule par calc_ubar_def
                    coef_defmy(i,j)=Vcol_y(i,j)/Uy_deformation(i,j)
- if ((itracebug.eq.1).and.(i.eq.313).and.(j.eq.143))  call tracebug(' test1')
+                   if ((itracebug.eq.1).and.(i.eq.313).and.(j.eq.143))  call tracebug(' test1')
                 else
- if ((itracebug.eq.1).and.(i.eq.313).and.(j.eq.143))  call tracebug(' test2')
+                   if ((itracebug.eq.1).and.(i.eq.313).and.(j.eq.143))  call tracebug(' test2')
                    coef_defmy(i,j)=1.
                 end if
@@ -412 +433 @@
                    ddby(i,j) = 0.
                    Uby(i,j)  = 0.
- if ((itracebug.eq.1).and.(i.eq.313).and.(j.eq.143))  call tracebug(' test3')
+                   if ((itracebug.eq.1).and.(i.eq.313).and.(j.eq.143))  call tracebug(' test3')
                 else
                    coef_defmy(i,j)=1.                                  ! pas de rescaling de la deformation
                    ddby(i,j)=-(Vcol_y(i,j)-Uy_deformation(i,j))/sdy(i,j)        ! pb si directions opposees
                    Uby(i,j) =  Vcol_y(i,j)-Uy_deformation(i,j)
-if ((itracebug.eq.1).and.(i.eq.313).and.(j.eq.144))  call tracebug(' test4')
+                   if ((itracebug.eq.1).and.(i.eq.313).and.(j.eq.144))  call tracebug(' test4')
 
                 end if
@@ -426 +447 @@
        end do
     end do
- if (itracebug.eq.1)  call tracebug(' apres test floty')
+    if (itracebug.eq.1)  call tracebug(' apres test floty')
 
     ! remarque coef peut être <1  a cause de la pente locale, mais il joue toujours son role 
@@ -453 +474 @@
 
   subroutine limit_coef_vitbil
-
+    
+    use module3D_phy, only: sdx,sdy,flotmx,flotmy,coefmxbmelt,coefmybmelt,gzmx,gzmy,ddbx,ddby,debug_3D,&
+         iglen
+    use deformation_mod_2lois, only: n1poly,n2poly,sa_mx,sa_my,s2a_mx,s2a_my
+    use runparam, only : itracebug
+    
+    integer :: i,j,k
+    
     call slope_surf
 
@@ -583 +611 @@
   subroutine calc_ubar_def   ! calculate velocity due to deformation before calibration
 
+    use module3D_phy, only: iglen,flotmx,flotmy,hmx,hmy,sdx,sdy,slope2mx,slope2my
+    use deformation_mod_2lois, only: n1poly,n2poly,glen,s2a_mx,s2a_my,ddx,ddy
+    use geography, only: dx,dy
+    use param_phy_mod, only: rog
+    use runparam, only: itracebug
+    
     implicit none              ! extrait de diffusiv pour calculer la partie deformation
 
@@ -588 +622 @@
     real :: inv_4dx ! inverse de dx pour eviter les divisions =1/(4*dx)
     real :: inv_4dy ! inverse de dy pour eviter les divisions =1/(4*dy)
+    integer :: i,j
 
     inv_4dx=1./(4*dx)
@@ -646 +681 @@
 
   subroutine ajust_ghf
-
+    
+    use module3D_phy, only: debug_3D,flot,ibase,ghf,ghf0,secyear
+    
     implicit none
+    
     real,dimension(nx,ny) :: coefdef_maj     !< coefficient deformation
     real :: increment_ghf
     real :: ghf_min
+    integer :: i,j
+    
     increment_ghf=0.00001                    !< exprime en  W/m2   
     ghf_min=0.025