Changeset 8504 for branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceadv_umx.F90

Timestamp:

2017-09-06T17:47:28+02:00 (7 years ago)

Author:

clem

Message:

changes in style - part4 - clarify ice advection routines

File:

• branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceadv_umx.F90 (modified) (2 diffs)

branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceadv_umx.F90

@@ -43 +43 @@
 CONTAINS
 
-   SUBROUTINE ice_adv_umx( kt, pdt, puc, pvc, pubox, pvbox, ptc )
+   SUBROUTINE ice_adv_umx( kt, pu_ice, pv_ice,  &
+      &                    pato_i, pv_i, pv_s, psmv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE ice_adv_umx  ***
+      !! 
+      !! **  Purpose :   Compute the now trend due to total advection of 
+      !!                 tracers and add it to the general trend of tracer equations
+      !!                 using an "Ultimate-Macho" scheme
+      !!
+      !! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74. 
+      !!----------------------------------------------------------------------
+      INTEGER                     , INTENT(in   ) ::   kt         ! time step
+      REAL(wp), DIMENSION(:,:)    , INTENT(in   ) ::   pu_ice     ! ice i-velocity
+      REAL(wp), DIMENSION(:,:)    , INTENT(in   ) ::   pv_ice     ! ice j-velocity
+      REAL(wp), DIMENSION(:,:)    , INTENT(inout) ::   pato_i     ! open water area
+      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_i       ! ice volume
+      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_s       ! snw volume
+      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   psmv_i     ! salt content
+      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   poa_i      ! age content
+      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pa_i       ! ice concentration
+      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pa_ip      ! melt pond fraction
+      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_ip      ! melt pond volume
+      REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_s       ! snw heat content
+      REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_i       ! ice heat content
+      !
+      INTEGER  ::   ji, jj, jk, jl, jt      ! dummy loop indices
+      INTEGER  ::   initad                  ! number of sub-timestep for the advection
+      REAL(wp) ::   zcfl , zusnit, zdt      !   -      -
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zudy, zvdx, zcu_box, zcv_box
+      !!----------------------------------------------------------------------
+      !
+      IF( kt == nit000 .AND. lwp ) THEN
+         WRITE(numout,*) 
+         WRITE(numout,*) 'ice_adv_umx : Ultimate-MACHO advection scheme'
+         WRITE(numout,*) '~~~~~~~~~~~'
+      ENDIF
+      !
+      ALLOCATE( zudy(jpi,jpj) , zvdx(jpi,jpj) , zcu_box(jpi,jpj) , zcv_box(jpi,jpj) )
+      !
+      ! --- If ice drift field is too fast, use an appropriate time step for advection (CFL test for stability) --- !        
+      zcfl  =            MAXVAL( ABS( pu_ice(:,:) ) * rdt_ice * r1_e1u(:,:) )
+      zcfl  = MAX( zcfl, MAXVAL( ABS( pv_ice(:,:) ) * rdt_ice * r1_e2v(:,:) ) )
+      IF( lk_mpp )   CALL mpp_max( zcfl )
+
+      IF( zcfl > 0.5 ) THEN   ;   initad = 2   ;   zusnit = 0.5_wp
+      ELSE                    ;   initad = 1   ;   zusnit = 1.0_wp
+      ENDIF
+
+      zdt = rdt_ice / REAL(initad)
+
+      ! --- transport --- !
+      zudy(:,:) = pu_ice(:,:) * e2u(:,:)
+      zvdx(:,:) = pv_ice(:,:) * e1v(:,:)
+
+      ! --- define velocity for advection: u*grad(H) --- !
+      DO jj = 2, jpjm1
+         DO ji = fs_2, fs_jpim1
+            IF    ( pu_ice(ji,jj) * pu_ice(ji-1,jj) <= 0._wp ) THEN   ;   zcu_box(ji,jj) = 0._wp
+            ELSEIF( pu_ice(ji,jj)                   >  0._wp ) THEN   ;   zcu_box(ji,jj) = pu_ice(ji-1,jj)
+            ELSE                                                      ;   zcu_box(ji,jj) = pu_ice(ji  ,jj)
+            ENDIF
+
+            IF    ( pv_ice(ji,jj) * pv_ice(ji,jj-1) <= 0._wp ) THEN   ;   zcv_box(ji,jj) = 0._wp
+            ELSEIF( pv_ice(ji,jj)                   >  0._wp ) THEN   ;   zcv_box(ji,jj) = pv_ice(ji,jj-1)
+            ELSE                                                      ;   zcv_box(ji,jj) = pv_ice(ji,jj  )
+            ENDIF
+         END DO
+      END DO
+
+      !---------------!
+      !== advection ==!
+      !---------------!
+      DO jt = 1, initad
+         CALL adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, pato_i(:,:) )             ! Open water area 
+         DO jl = 1, jpl
+            CALL adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, pa_i(:,:,jl) )         ! Ice area
+            CALL adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, pv_i(:,:,jl) )         ! Ice  volume
+            CALL adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, psmv_i(:,:,jl) )       ! Salt content
+            CALL adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, poa_i (:,:,jl) )       ! Age content
+            DO jk = 1, nlay_i
+               CALL adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, pe_i(:,:,jk,jl) )   ! Ice  heat content
+            END DO
+            CALL adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, pv_s(:,:,jl) )         ! Snow volume
+            CALL adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, pe_s(:,:,1,jl) )       ! Snow heat content
+            IF ( nn_pnd_scheme > 0 ) THEN
+               CALL adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, pa_ip(:,:,jl) )     ! Melt pond fraction
+               CALL adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, pv_ip(:,:,jl) )     ! Melt pond volume
+            ENDIF
+         END DO
+      END DO
+      !
+      DEALLOCATE( zudy, zvdx, zcu_box, zcv_box )
+      !
+   END SUBROUTINE ice_adv_umx
+   
+   SUBROUTINE adv_umx( kt, pdt, puc, pvc, pubox, pvbox, ptc )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE adv_umx  ***
       !! 
       !! **  Purpose :   Compute the now trend due to total advection of 
@@ -146 +241 @@
       IF( nn_timing == 1 )  CALL timing_stop('ice_adv_umx')
       !
-   END SUBROUTINE ice_adv_umx
+   END SUBROUTINE adv_umx