Changeset 888 for trunk/NEMO/LIM_SRC_3/limdyn.F90

Timestamp:

2008-04-11T19:05:03+02:00 (16 years ago)

Author:

ctlod

Message:

merge dev_001_SBC branche with the trunk to include the New Surface Module package, see ticket: #113

File:

• trunk/NEMO/LIM_SRC_3/limdyn.F90 (modified) (7 diffs)

trunk/NEMO/LIM_SRC_3/limdyn.F90

@@ -16 +16 @@
    USE dom_ice
    USE dom_oce         ! ocean space and time domain
-   USE taumod
    USE ice
    USE par_ice
+   USE sbc_ice         ! Surface boundary condition: ice fields
    USE ice_oce
    USE iceini
@@ -41 +41 @@
    !!----------------------------------------------------------------------
    !!   LIM 3.0,  UCL-ASTR-LOCEAN-IPSL (2008)
-   !! $Header: /home/opalod/NEMOCVSROOT/NEMO/LIM_SRC/limdyn.F90,v 1.5 2005/03/27 18:34:41 opalod Exp $
-   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt
+   !! $ Id: $
+   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
 
@@ -90 +90 @@
       IF ( ln_limdyn ) THEN
 
-         ! ocean velocity
-         u_oce(:,:)  = u_io(:,:) * tmu(:,:)
-         v_oce(:,:)  = v_io(:,:) * tmv(:,:)
-         
          old_u_ice(:,:) = u_ice(:,:) * tmu(:,:)
          old_v_ice(:,:) = v_ice(:,:) * tmv(:,:)
@@ -162 +158 @@
          ENDIF
 
-         ! Ice-Ocean stress
-         ! ================
-         DO jj = 2, jpjm1
-            zsang  = SIGN(1.e0, gphif(1,jj-1) ) * sangvg
-
-            DO ji = fs_2, fs_jpim1
-               ! computation of wind stress over ocean in X and Y direction
-#if defined key_coupled && defined key_lim_cp1
-!              ztairx =  ( 1.0 - at_i(ji-1,jj)   ) * gtaux(ji-1,jj)   + &
-!                        ( 1.0 - at_i(ji,jj)     ) * gtaux(ji,jj  )   + &
-!                        ( 1.0 - at_i(ji-1,jj-1) ) * gtaux(ji-1,jj-1) + & 
-!                        ( 1.0 - at_i(ji,jj-1)   ) * gtaux(ji,jj-1)
-
-!              ztairy =  ( 1.0 - at_i(ji-1,jj)   ) * gtauy(ji-1,jj  ) + &
-!                        ( 1.0 - at_i(ji,jj  )   ) * gtauy(ji,jj    ) + &
-!                        ( 1.0 - at_i(ji-1,jj-1) ) * gtauy(ji-1,jj-1) + & 
-!                        ( 1.0 - at_i(ji,jj-1)   ) * gtauy(ji,jj-1)
-#else
-               ztairx =  ( 2.0 - at_i(ji,jj) - at_i(ji+1,jj) ) * gtaux(ji,jj) / cai * cao
-               ztairy =  ( 2.0 - at_i(ji,jj) - at_i(ji,jj+1) ) * gtauy(ji,jj) / cai * cao
-
-               zsfrldmx2 = at_i(ji,jj) + at_i(ji+1,jj)
-               zsfrldmy2 = at_i(ji,jj) + at_i(ji,jj+1)
-
-#endif
-               zu_ice   = u_ice(ji,jj) - u_oce(ji,jj)
-               zv_ice   = v_ice(ji,jj) - v_oce(ji,jj)
-               zmod     = SQRT( zu_ice * zu_ice + zv_ice * zv_ice ) 
-
-               ! quadratic drag formulation
-               ztglx   = zsfrldmx2 * rhoco * zmod * ( cangvg * zu_ice - zsang * zv_ice ) 
-               ztgly   = zsfrldmy2 * rhoco * zmod * ( cangvg * zv_ice + zsang * zu_ice ) 
-!
-!              ! IMPORTANT
-!              ! these lignes are bound to prevent numerical oscillations
-!              ! in the ice-ocean stress
-!              ! They are physically ill-based. There is a cleaner solution
-!              ! to try (remember discussion in Paris Gurvan)
-!
-               ztglx   = ztglx * exp( - zmod / 0.5 )
-               ztgly   = ztglx * exp( - zmod / 0.5 )
-
-               tio_u(ji,jj) = - ( ztairx + 1.0 * ztglx ) / ( 2. * rau0 )
-               tio_v(ji,jj) = - ( ztairy + 1.0 * ztgly ) / ( 2. * rau0 )
-            END DO
-         END DO
-         
          ! computation of friction velocity
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1
 
-               zu_ice   = u_ice(ji,jj) - u_io(ji,jj)
+               zu_ice   = u_ice(ji,jj) - u_oce(ji,jj)
                zt11  = rhoco * zu_ice * zu_ice
 
-               zu_ice   = u_ice(ji-1,jj) - u_io(ji-1,jj)
+               zu_ice   = u_ice(ji-1,jj) - u_oce(ji-1,jj)
                zt12  = rhoco * zu_ice * zu_ice
 
-               zv_ice   = v_ice(ji,jj) - v_io(ji,jj)
+               zv_ice   = v_ice(ji,jj) - v_oce(ji,jj)
                zt21  = rhoco * zv_ice * zv_ice
 
-               zv_ice   = v_ice(ji,jj-1) - v_io(ji,jj-1)
+               zv_ice   = v_ice(ji,jj-1) - v_oce(ji,jj-1)
                zt22  = rhoco * zv_ice * zv_ice
-               ztair2 = ( ( gtaux(ji,jj) + gtaux(ji-1,jj) ) / 2. )**2 + &
-                        ( ( gtauy(ji,jj) + gtauy(ji,jj-1) ) / 2. )**2
+               ztair2 = ( ( utaui_ice(ji,jj) + utaui_ice(ji-1,jj) ) / 2. )**2 + &
+                        ( ( vtaui_ice(ji,jj) + vtaui_ice(ji,jj-1) ) / 2. )**2
 
                ! should not be weighted
@@ -241 +190 @@
           DO jj = 2, jpjm1
              DO ji = fs_2, fs_jpim1
-#if defined key_coupled && defined key_lim_cp1
-                tio_u(ji,jj) = - (  gtaux(ji  ,jj  ) + gtaux(ji-1,jj  )       &
-                   &              + gtaux(ji-1,jj-1) + gtaux(ji  ,jj-1) ) / ( 4 * rau0 )
-
-                tio_v(ji,jj) = - (  gtauy(ji  ,jj )  + gtauy(ji-1,jj  )       &
-                   &              + gtauy(ji-1,jj-1) + gtauy(ji  ,jj-1) ) / ( 4 * rau0 )
-#else
-                tio_u(ji,jj) = - gtaux(ji,jj) / cai * cao / rau0
-                tio_v(ji,jj) = - gtauy(ji,jj) / cai * cao / rau0 
-#endif
-                ztair2 = ( ( gtaux(ji,jj) + gtaux(ji-1,jj) ) / 2. )**2 + &
-                         ( ( gtauy(ji,jj) + gtauy(ji,jj-1) ) / 2. )**2
+                ztair2 = ( ( utaui_ice(ji,jj) + utaui_ice(ji-1,jj) ) / 2. )**2 + &
+                         ( ( vtaui_ice(ji,jj) + vtaui_ice(ji,jj-1) ) / 2. )**2
                 zustm        = SQRT( ztair2  )
 
@@ -262 +201 @@
 
       CALL lbc_lnk( ust2s, 'T',  1. )   ! T-point
-      CALL lbc_lnk( tio_u, 'U', -1. )   ! I-point (i.e. ice U-V point)
-      CALL lbc_lnk( tio_v, 'V', -1. )   ! I-point (i.e. ice U-V point)
 
       IF(ln_ctl) THEN   ! Control print
@@ -269 +206 @@
          CALL prt_ctl_info(' - Cell values : ')
          CALL prt_ctl_info('   ~~~~~~~~~~~~~ ')
-         CALL prt_ctl(tab2d_1=tio_u     , clinfo1=' lim_dyn  : tio_u     :', tab2d_2=tio_v , clinfo2=' tio_v :')
          CALL prt_ctl(tab2d_1=ust2s     , clinfo1=' lim_dyn  : ust2s     :')
          CALL prt_ctl(tab2d_1=divu_i    , clinfo1=' lim_dyn  : divu_i    :')