Changeset 187 for branches/iLoveclim/SOURCES/furst_schoof_mod.f90

Timestamp:

03/15/18 14:18:24 (6 years ago)

Author:

aquiquet

Message:

Grisli-iloveclim branch merged to trunk at revision 185

Location:

branches/iLoveclim

Files:

• . (modified) (1 prop)
• SOURCES/furst_schoof_mod.f90 (modified) (13 diffs)

branches/iLoveclim/SOURCES/furst_schoof_mod.f90

@@ -14 +14 @@
 real,dimension(nx,ny) :: back_force_y
 real, parameter :: m_weert = 1.0 ! afq: is this defined elsewhere in the model???
+real, parameter :: inv_mweert = 1./m_weert
 
 contains
@@ -54 +55 @@
   integer,dimension(nx,ny) :: countx, county   ! how often do we modify ux/uy
   real :: phi_prescr
-  real :: archim
   real :: toutpetit = 1e-6
   real :: denom, prodscal
+
+  real :: bfx, bfy
 
   !debug
@@ -102 +104 @@
               xpos_tab(i,j)=xpos
               Hglx_tab(i,j)=Hgl
-
+              
+              ! afq: the back force is on the staggered grid, the GL can be either West or East to this point.
+              if (xpos .lt. -dx/2.) then
+                 bfx = back_force_x(i,j)
+              else
+                 bfx = back_force_x(i+1,j)
+              endif
               if (gr_select.eq.1) then ! flux de Tsai
-                 call flux_Tsai4Schoof (Hgl,Abar(i,j),frot_coef,alpha_flot,glen(1),back_force_x(i-1,j),phi_prescr)
+                 call flux_Tsai4Schoof (Hgl,Abar(i,j),frot_coef,alpha_flot,glen(1),bfx,phi_prescr)
               else if (gr_select.eq.2) then ! flux de Schoof
                  ! afq: the dragging coef. is on the staggered grid, the GL can be either West or East to this point.
@@ -112 +120 @@
                     frot_coef = betamx(i+1,j)
                  endif
-                 call flux_Schoof4Schoof (Hgl,Abar(i,j),frot_coef,alpha_flot,glen(1),m_weert,back_force_x(i-1,j),phi_prescr)
+                 call flux_Schoof4Schoof (Hgl,Abar(i,j),frot_coef,alpha_flot,glen(1),m_weert,bfx,phi_prescr)
               else
                  print*,'ATTENTION FLUX AUTRE QUE TSAI OU SCHOOF NON IMPLEMENTE'
@@ -190 +198 @@
               xpos_tab(i,j)=xpos
               Hglx_tab(i,j)=Hgl
+              
+              ! afq: the back force is on the staggered grid, the GL can be either West or East to this point.
+              if (xpos .lt. dx/2.) then
+                 bfx = back_force_x(i,j)
+              else
+                 bfx = back_force_x(i+1,j)
+              endif
               if (gr_select.eq.1) then ! flux de Tsai
-                 call flux_Tsai4Schoof (Hgl,Abar(i,j),frot_coef,alpha_flot,glen(1),back_force_x(i+1,j),phi_prescr)
+                 call flux_Tsai4Schoof (Hgl,Abar(i,j),frot_coef,alpha_flot,glen(1),bfx,phi_prescr)
               else if (gr_select.eq.2) then ! flux de Schoof
                  ! afq: the dragging coef. is on the staggered grid, the GL can be either West or East to this point.
-                 if (xpos .lt. -dx/2.) then
+                 if (xpos .lt. dx/2.) then
                     frot_coef = betamx(i,j)
                  else
                     frot_coef = betamx(i+1,j)
                  endif
-                 call flux_Schoof4Schoof (Hgl,Abar(i,j),frot_coef,alpha_flot,glen(1),m_weert,back_force_x(i+1,j),phi_prescr)
+                 call flux_Schoof4Schoof (Hgl,Abar(i,j),frot_coef,alpha_flot,glen(1),m_weert,bfx,phi_prescr)
               else
                  print*,'ATTENTION FLUX AUTRE QUE TSAI OU SCHOOF NON IMPLEMENTE'
@@ -205 +220 @@
               endif
               phi_prescr_tabx(i,j)=phi_prescr
-              if (xpos .lt. -dx/2.) then  !  GL a west du i staggered,                          o centre, x stag
+              if (xpos .lt. dx/2.) then  !  GL a west du i staggered,                          o centre, x stag
 
                  !  grille    ! .....x......o......x......o......x......O..G...x......o......x......o
@@ -274 +289 @@
               ypos_tab(i,j)=ypos
               Hgly_tab(i,j)=Hgl
+              
+              ! afq: the back force is on the staggered grid, the GL can be either South or North to this point.
+              if (ypos .lt. -dy/2.) then
+                 bfy = back_force_y(i,j)
+              else
+                 bfy = back_force_y(i,j+1)
+              endif
               if (gr_select.eq.1) then ! flux de Tsai
-                 call flux_Tsai4Schoof (Hgl,Abar(i,j),frot_coef,alpha_flot,glen(1),back_force_y(i,j-1),phi_prescr)
+                 call flux_Tsai4Schoof (Hgl,Abar(i,j),frot_coef,alpha_flot,glen(1),bfy,phi_prescr)
               else if (gr_select.eq.2) then ! flux de Schoof
                  ! afq: the dragging coef. is on the staggered grid, the GL can be either South or North to this point.
@@ -283 +305 @@
                     frot_coef = betamy(i,j+1)
                  endif
-                 call flux_Schoof4Schoof (Hgl,Abar(i,j),frot_coef,alpha_flot,glen(1),m_weert,back_force_y(i,j-1),phi_prescr)
+                 call flux_Schoof4Schoof (Hgl,Abar(i,j),frot_coef,alpha_flot,glen(1),m_weert,bfy,phi_prescr)
               else
                  print*,'ATTENTION FLUX AUTRE QUE TSAI OU SCHOOF NON IMPLEMENTE'
@@ -355 +377 @@
               ypos_tab(i,j)=ypos
               Hgly_tab(i,j)=Hgl
+              
+              ! afq: the back force is on the staggered grid, the GL can be either South or North to this point.
+              if (ypos .lt. dy/2.) then
+                 bfy = back_force_y(i,j)
+              else
+                 bfy = back_force_y(i,j+1)
+              endif
               if (gr_select.eq.1) then ! flux de Tsai
-                 call flux_Tsai4Schoof (Hgl,Abar(i,j),frot_coef,alpha_flot,glen(1),back_force_y(i,j+1),phi_prescr)
+                 call flux_Tsai4Schoof (Hgl,Abar(i,j),frot_coef,alpha_flot,glen(1),bfy,phi_prescr)
               else if (gr_select.eq.2) then ! flux de Schoof
                  ! afq: the dragging coef. is on the staggered grid, the GL can be either South or North to this point.
-                 if (ypos .lt. -dy/2.) then
+                 if (ypos .lt. dy/2.) then
                     frot_coef = betamy(i,j)
                  else
                     frot_coef = betamy(i,j+1)
                  endif
-                 call flux_Schoof4Schoof (Hgl,Abar(i,j),frot_coef,alpha_flot,glen(1),m_weert,back_force_y(i,j+1),phi_prescr)
+                 call flux_Schoof4Schoof (Hgl,Abar(i,j),frot_coef,alpha_flot,glen(1),m_weert,bfy,phi_prescr)
               else
                  print*,'ATTENTION FLUX AUTRE QUE TSAI OU SCHOOF NON IMPLEMENTE'
@@ -370 +399 @@
               endif
               phi_prescr_taby(i,j)=phi_prescr
-              if (ypos .lt. -dy/2.) then  !  GL au sud du j staggered,                          o centre, x stag
+              if (ypos .lt. dy/2.) then  !  GL au sud du j staggered,                          o centre, x stag
 
                  !  grille    ! .....x......o......x......o......x......O..G...x......o......x......o
@@ -471 +500 @@
 !!$     enddo
 !!$  enddo
+
+  debug_3D(:,:,66) = phi_prescr_tabx(:,:)
+  debug_3D(:,:,67) = phi_prescr_taby(:,:)
+  
 end subroutine interpol_glflux
   
@@ -561 +594 @@
 ! afq: in standard GRISLI (19/05/17), tob= -beta ub => m_weert = 1 and C_frot = beta
 
-
-phi_schoof = (((A_mean * rog**(n_glen+1) * (1 - alpha)**n_glen) / (4**n_glen *C_frot)) **(1./(1.+1./m_weert)))*back_force_coef
-phi_schoof = phi_schoof * Hgl**((n_glen+3.+1/m_weert)/(1.+1./m_weert))
+!phi_schoof = (((A_mean * rog**(n_glen+1) * (1 - alpha)**n_glen) / (4**n_glen *C_frot)) **(1./(1.+1./m_weert)))*back_force_coef**(n_glen/(1+1./m_weert))
+!phi_schoof = phi_schoof * Hgl**((n_glen+3.+1/m_weert)/(1.+1./m_weert))
+phi_schoof = (((A_mean * rog**(n_glen+1.) * (1. - alpha)**n_glen) / (4.**n_glen*C_frot)) **(1./(1.+inv_mweert)))*back_force_coef**(n_glen/(1.+inv_mweert))
+phi_schoof = phi_schoof * Hgl**((n_glen+3.+inv_mweert)/(1.+inv_mweert))
 
 
@@ -590 +624 @@
 
 
-phi_Tsai = Q0 * ((8 * A_mean * rog**(n_glen) * (1 - alpha)**(n_glen-1.)) / (4**n_glen *f_frot))*back_force_coef
-phi_Tsai = phi_Tsai * Hgl**(n_glen+2)
+phi_Tsai = Q0 * ((8. * A_mean * rog**(n_glen) * (1. - alpha)**(n_glen-1.)) / (4.**n_glen *f_frot))*back_force_coef**(n_glen-1.)
+phi_Tsai = phi_Tsai * Hgl**(n_glen+2.)
 
 end subroutine flux_Tsai4Schoof