Changeset 13544


Ignore:
Timestamp:
2020-09-30T09:03:49+02:00 (4 months ago)
Author:
vancop
Message:

now VP rheology compiles

Location:
NEMO/branches/2020/SI3-03_VP_rheology/src/ICE
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2020/SI3-03_VP_rheology/src/ICE/icedyn_rhg.F90

    r13536 r13544  
    1717   USE ice            ! sea-ice: variables 
    1818   USE icedyn_rhg_evp ! sea-ice: EVP rheology 
     19   USE icedyn_rhg_vp  ! sea-ice: VP  rheology 
    1920   USE icectl         ! sea-ice: control prints 
    2021   ! 
     
    8283         !                             !---------------------------------! 
    8384         CALL ice_dyn_rhg_evp( kt, stress1_i, stress2_i, stress12_i, shear_i, divu_i, delta_i ) 
    84          !          
    85       END SELECT 
    86  
    8785      !                                !---------------------------------------------! 
    8886      CASE( np_rhgVP  )                ! Viscous-Plastic rheology                    ! 
     
    166164      ! 
    167165      IF( ln_rhg_EVP  )   CALL rhg_evp_rst( 'READ' )  !* read or initialize all required files 
    168       IF( ln_rhg_VP   )   CALL rhg_lsr_rst( 'READ' )  !* read or initialize all required files 
    169166      ! 
    170167   END SUBROUTINE ice_dyn_rhg_init 
  • NEMO/branches/2020/SI3-03_VP_rheology/src/ICE/icedyn_rhg_evp.F90

    r13536 r13544  
    169169      REAL(wp), DIMENSION(jpi,jpj) ::   zu_ice, zv_ice 
    170170      !! --- diags 
    171       REAL(wp) ::   zsig1, zsig2, zsig12 
     171      REAL(wp) ::   zsig1, zsig2, zsig12, zfac, z1_strength 
    172172      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zsig_I, zsig_II, zsig1_p, zsig2_p, zten_i          
    173173      !! --- SIMIP diags 
  • NEMO/branches/2020/SI3-03_VP_rheology/src/ICE/icedyn_rhg_vp.F90

    r13536 r13544  
    1 ! TODOLOIST 
    2 ! 
    3 ! Commit 
    4 ! Compile 
    5 ! 
    6 ! - Re-calculate internal force diagnostic (currently incorrect, see end of the routine) 
    7 ! - QC:compare code to mitGCM 
    8 ! - QC: comparing EVP and VP codes  
    9 ! 
    10 ! Clarify 
    11 ! --- Boundary conditions --> how to enforce them - is the fmask strategy sufficient ? 
    12 ! --- Swap of independent term in the U-V solvers ? -> JF Lemieux says maybe 
    13 ! --- Is stress tensor for restart needed ? No!!! 
    14 ! 
    15 ! Test 
    16 ! --- Try ADI for u-v solver 
    17 ! 
    18 ! Add 
    19 ! - Charge ellipse as an output (good one from Lemieux and Dupont) 
    20 ! - Bottom drag (landfast ice) 
    21 ! - Tensile strength 
    22 ! - agrif 
    23 ! - bdy 
    24 ! 
    25 ! Write documentation 
    26 ! 
    27 ! Optimize 
    28 ! - subroutinize common parts (diagnostics) 
    29 ! - namelist: rationalize common parameters EVP/VP 
    30  
    31  
    321MODULE icedyn_rhg_vp 
    332   !!====================================================================== 
     
    7342 
    7443   !! for convergence tests 
    75    INTEGER ::   ncvgid   ! netcdf file id 
     44   INTEGER ::   ncvgid        ! netcdf file id 
    7645   INTEGER ::   nvarid_ucvg   ! netcdf variable id 
    7746   INTEGER ::   nvarid_ures 
     
    8251   INTEGER ::   nvarid_mke 
    8352   REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zmsk00, zmsk15 
     53 
    8454   !!---------------------------------------------------------------------- 
    8555   !! NEMO/ICE 4.0 , NEMO Consortium (2018) 
     
    157127      REAL(wp), DIMENSION(:,:), INTENT(  out) ::   pshear_i  , pdivu_i   , pdelta_i      ! 
    158128      !! 
    159       LOGICAL ::   ll_u_iterate, ll_viterate   ! continue iteration or not 
     129      LOGICAL ::   ll_u_iterate, ll_v_iterate   ! continue iteration or not 
    160130      ! 
    161131      INTEGER ::   ji, jj, jn          ! dummy loop indices 
     
    171141      REAL(wp) ::   zm2, zm3, zmassU, zmassV                            ! ice/snow mass and volume 
    172142      REAL(wp) ::   zdeltat, zds2, zdt, zdt2, zdiv, zdiv2               ! temporary scalars 
    173       REAL(wp) ::   zp_deltastar_f 
     143      REAL(wp) ::   zp_deltastar_f                                      !  
     144      REAL(wp) ::   zu_cV, zv_cU                                        !  
    174145      REAL(wp) ::   zfac, zfac1, zfac2, zfac3 
    175146      REAL(wp) ::   zt12U, zt11U, zt22U, zt21U, zt122U, zt121U 
     
    178149      ! 
    179150      REAL(wp), DIMENSION(jpi,jpj) ::   zfmask                          ! mask at F points for the ice 
    180       REAL(wp), DIMENSION(jpi,jpj) ::   zaU  , zaV                      ! ice fraction on U/V points 
     151      REAL(wp), DIMENSION(jpi,jpj) ::   za_iU  , za_iV                      ! ice fraction on U/V points 
    181152      REAL(wp), DIMENSION(jpi,jpj) ::   zmU_t, zmV_t                    ! Acceleration term contribution to RHS 
    182153      REAL(wp), DIMENSION(jpi,jpj) ::   zmassU_t, zmassV_t              ! Mass per unit area divided by time step 
     
    191162      REAL(wp), DIMENSION(jpi,jpj) ::   v_oceU, u_oceV, v_iceU, u_iceV  ! ocean/ice u/v component on V/U points 
    192163      REAL(wp), DIMENSION(jpi,jpj) ::   zu_c, zv_c                      ! "current" ice velocity (m/s), average of previous two OL iterates 
    193       REAL(wp), DIMENSION(jpi,jpj) ::   zu_cV, zv_cU                    ! "current" u (v) ice velocity at V (U) point (m/s) 
    194164      REAL(wp), DIMENSION(jpi,jpj) ::   zu_b, zv_b                      ! velocity at previous sub-iterate 
    195165      REAL(wp), DIMENSION(jpi,jpj) ::   zuerr, zverr                    ! absolute U/Vvelocity difference between current and previous sub-iterates 
     
    209179      REAL(wp), DIMENSION(jpi,jpj) ::   zrhsu, zrhsv                    ! U and V RHS  
    210180      REAL(wp), DIMENSION(jpi,jpj) ::   zAU, zBU, zCU, zDU, zEU         ! Linear system coefficients, U equation 
    211       REAL(wp), DIMENSION(jpi,jpj) ::   zAV, zBV, zCV, zDV, zEV, zFV    ! Linear system coefficients, V equation 
     181      REAL(wp), DIMENSION(jpi,jpj) ::   zAV, zBV, zCV, zDV, zEV         ! Linear system coefficients, V equation 
    212182      REAL(wp), DIMENSION(jpi,jpj) ::   zFU, zFU_prime, zBU_prime       ! Rearranged linear system coefficients, U equation 
    213183      REAL(wp), DIMENSION(jpi,jpj) ::   zFV, zFV_prime, zBV_prime       ! Rearranged linear system coefficients, V equation 
     
    222192      REAL(wp), PARAMETER          ::   zamin  = 0.001_wp               ! ice concentration below which ice velocity becomes very small 
    223193      !! --- diags 
    224       REAL(wp) ::   zsig1, zsig2, zsig12 
     194      REAL(wp) ::   zsig1, zsig2, zsig12, zdelta, z1_strength, zfac_x, zfac_y 
    225195      REAL(wp), DIMENSION(jpi,jpj) ::   zs1, zs2, zs12, zs12f, zten_i   ! stress tensor components 
    226196      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zsig_I, zsig_II, zsig1_p, zsig2_p 
    227       REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_xmtrp_ice ! X-component of ice mass transport (kg/s, SIMIP) 
    228       REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_ymtrp_ice ! Y-component of ice mass transport (kg/s, SIMIP) 
    229       REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_xmtrp_snw ! X-component of snow mass transport (kg/s, SIMIP) 
    230       REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_ymtrp_snw ! Y-component of snow mass transport (kg/s, SIMIP) 
    231       REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_xatrp     ! X-component of area transport (m2/s, SIMIP) 
    232       REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_yatrp     ! Y-component of area transport (m2/s, SIMIP)       
     197      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   ztaux_oi, ztauy_oi 
     198      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_xmtrp_ice, zdiag_ymtrp_ice ! X/Y-component of ice mass transport (kg/s, SIMIP) 
     199      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_xmtrp_snw, zdiag_ymtrp_snw ! X/Y-component of snow mass transport (kg/s, SIMIP) 
     200      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_xatrp, zdiag_yatrp         ! X/Y-component of area transport (m2/s, SIMIP) 
    233201                         
    234202      !!---------------------------------------------------------------------------------------------------------------------- 
     
    251219      END DO 
    252220       
    253       IF ( ln_zebra_vp ) THEN; nn_nzebra_vp = 2; ELSE; nn_nzebra_vp = 1; ENDIF  
    254        
    255       nn_nvp = nn_out_vp * nn_inn_vp ! maximum number of iterations 
     221      IF ( ln_zebra_vp ) THEN; nn_zebra_vp = 2; ELSE; nn_zebra_vp = 1; ENDIF  
     222       
     223      nn_nvp = nn_nout_vp * nn_ninn_vp ! maximum number of iterations 
    256224       
    257225      zrhoco = rau0 * rn_cio  
     
    341309 
    342310            ! Ice fraction at U-V points 
    343             zaU(ji,jj)      = 0.5_wp * ( at_i(ji,jj) * e1e2t(ji,jj) + at_i(ji+1,jj) * e1e2t(ji+1,jj) ) * r1_e1e2u(ji,jj) * umask(ji,jj,1) 
    344             zaV(ji,jj)      = 0.5_wp * ( at_i(ji,jj) * e1e2t(ji,jj) + at_i(ji,jj+1) * e1e2t(ji,jj+1) ) * r1_e1e2v(ji,jj) * vmask(ji,jj,1) 
     311            za_iU(ji,jj)    = 0.5_wp * ( at_i(ji,jj) * e1e2t(ji,jj) + at_i(ji+1,jj) * e1e2t(ji+1,jj) ) * r1_e1e2u(ji,jj) * umask(ji,jj,1) 
     312            za_iV(ji,jj)    = 0.5_wp * ( at_i(ji,jj) * e1e2t(ji,jj) + at_i(ji,jj+1) * e1e2t(ji,jj+1) ) * r1_e1e2v(ji,jj) * vmask(ji,jj,1) 
    345313 
    346314            ! Snow and ice mass at U-V points 
     
    368336             
    369337            ! Wind stress 
    370             ztaux_ai(ji,jj) = zaU(ji,jj) * utau_ice(ji,jj) 
    371             ztauy_ai(ji,jj) = zaV(ji,jj) * vtau_ice(ji,jj) 
     338            ztaux_ai(ji,jj) = za_iU(ji,jj) * utau_ice(ji,jj) 
     339            ztauy_ai(ji,jj) = za_iV(ji,jj) * vtau_ice(ji,jj) 
    372340 
    373341            ! Force due to sea surface tilt(- m*g*GRAD(ssh)) 
     
    380348 
    381349            ! Mask for lots of ice (1) or little ice (0) 
    382             IF( zmassU <= zmmin .AND. zaU(ji,jj) <= zamin ) THEN   ;   zmsk01x(ji,jj) = 0._wp 
     350            IF( zmassU <= zmmin .AND. za_iU(ji,jj) <= zamin ) THEN   ;   zmsk01x(ji,jj) = 0._wp 
    383351            ELSE                                                   ;   zmsk01x(ji,jj) = 1._wp   ;   ENDIF 
    384             IF( zmassV <= zmmin .AND. zaV(ji,jj) <= zamin ) THEN   ;   zmsk01y(ji,jj) = 0._wp 
     352            IF( zmassV <= zmmin .AND. za_iV(ji,jj) <= zamin ) THEN   ;   zmsk01y(ji,jj) = 0._wp 
    385353            ELSE                                                   ;   zmsk01y(ji,jj) = 1._wp   ;   ENDIF               
    386354 
     
    499467                 
    500468                !--- non-linear drag coefficients (need to be updated at each outer loop, see Lemieux and Tremblay JGR09, p.3, beginning of Section 3) 
    501                 zCwU(ji,jj)          = zaU(ji,jj) * zrhoco * SQRT( ( zu_c (ji,jj) - u_oce (ji,jj) ) * ( zu_c (ji,jj) - u_oce (ji,jj) )  & 
    502                   &                                              + ( zv_cU - v_oceU(ji,jj) ) * ( zv_cU - v_oceU(ji,jj) ) ) 
    503                 zCwV(ji,jj)          = zaV(ji,jj) * zrhoco * SQRT( ( zv_c (ji,jj) - v_oce (ji,jj) ) * ( zv_c (ji,jj) - v_oce (ji,jj) )  & 
    504                   &                                              + ( zu_cV - u_oceV(ji,jj) ) * ( zu_cV - u_oceV(ji,jj) ) ) 
     469                zCwU(ji,jj)          = za_iU(ji,jj) * zrhoco * SQRT( ( zu_c (ji,jj) - u_oce (ji,jj) ) * ( zu_c (ji,jj) - u_oce (ji,jj) )  & 
     470                  &                                                + ( zv_cU - v_oceU(ji,jj) ) * ( zv_cU - v_oceU(ji,jj) ) ) 
     471                zCwV(ji,jj)          = za_iV(ji,jj) * zrhoco * SQRT( ( zv_c (ji,jj) - v_oce (ji,jj) ) * ( zv_c (ji,jj) - v_oce (ji,jj) )  & 
     472                  &                                                + ( zu_cV - u_oceV(ji,jj) ) * ( zu_cV - u_oceV(ji,jj) ) ) 
    505473                  
    506474                !--- Ocean-ice drag contributions to RHS  
     
    577545                              (    e2u(ji,jj)    * ( zs1_rhsu(ji+1,jj) - zs1_rhsu(ji,jj) )                                                                 & 
    578546                  &      +         r1_e2u(ji,jj) * ( e2t(ji+1,jj) * e2t(ji+1,jj) * zs2_rhsu(ji+1,jj) - e2t(ji,jj)   * e2t(ji,jj)   * zs2_rhsu(ji,jj) )     & 
    579                   &      + 2._wp * r1_e1u(ji,jj) * ( e1f(ji,jj)   * e1f(ji,jj)   * zs12_rhsu(ji,jj)  - e1f(ji,jj-1) * e1f(ji,jj-1) * zs12_rhsu(ji,jj-1) )  
     547                  &      + 2._wp * r1_e1u(ji,jj) * ( e1f(ji,jj)   * e1f(ji,jj)   * zs12_rhsu(ji,jj)  - e1f(ji,jj-1) * e1f(ji,jj-1) * zs12_rhsu(ji,jj-1) ) ) 
    580548                   
    581549               ! --- V component of internal force contribution to RHS at V points 
     
    583551                  &           (    e1v(ji,jj)    * ( zs1_rhsv(ji,jj+1) - zs1_rhsv(ji,jj) )                                                                 & 
    584552                  &      +         r1_e1v(ji,jj) * ( e1t(ji,jj+1) * e1t(ji,jj+1) * zs2_rhsv(ji,jj+1) - e1t(ji,jj) * e1t(ji,jj)     * zs2_rhsv(ji,jj) )     & 
    585                   &      + 2._wp * r1_e2v(ji,jj) * ( e2f(ji,jj)   * e2f(ji,jj)   * zs12_rhsv(ji,jj)  - e2f(ji-1,jj) * e2f(ji-1,jj) * zs12_rhsv(ji-1,jj) ) 
     553                  &      + 2._wp * r1_e2v(ji,jj) * ( e2f(ji,jj)   * e2f(ji,jj)   * zs12_rhsv(ji,jj)  - e2f(ji-1,jj) * e2f(ji-1,jj) * zs12_rhsv(ji-1,jj) ) ) 
    586554                   
    587555            END DO 
     
    723691                  ! what follows could be subroutinized... 
    724692                   
    725                   DO jn = 1, nn_nzebra_vp ! "zebra" loop (! red-black-sor!!! ) 
     693                  DO jn = 1, nn_zebra_vp ! "zebra" loop (! red-black-sor!!! ) 
    726694                   
    727695                     ! OPT: could be even better optimized with a true red-black SOR 
     
    735703                     zBU_prime(:,:) = 0._wp 
    736704          
    737                      DO jj = jj_min, jpj - 1, nn_nzebra_vp 
     705                     DO jj = jj_min, jpj - 1, nn_zebra_vp 
    738706 
    739707                        !------------------------------------------- 
     
    756724      
    757725                           ! right hand side 
    758                            zFU(ji,jj) = ( zrhsu(ji,jj) &                                 ! right-hand side terms 
    759                                &      +   zAA3                                           ! boundary condition translation 
    760                                &      +   zDU(ji,jj) * u_ice(ji,jj-1)                    ! internal force, j-1 
    761                                &      +   zEU(ji,jj) * u_ice(ji,jj+1) ) * umask(ji,jj,1) ! internal force, j+1 
     726                           zFU(ji,jj) = ( zrhsu(ji,jj) &                                   ! right-hand side terms 
     727                               &      +   zAA3                                           & ! boundary condition translation 
     728                               &      +   zDU(ji,jj) * u_ice(ji,jj-1)                    & ! internal force, j-1 
     729                               &      +   zEU(ji,jj) * u_ice(ji,jj+1) ) * umask(ji,jj,1)   ! internal force, j+1 
    762730                     
    763731                        END DO 
     
    806774                  ! MV OPT: what follows could be subroutinized... 
    807775                   
    808                   DO jn = 1, nn_nzebra_vp ! "zebra" loop 
     776                  DO jn = 1, nn_zebra_vp ! "zebra" loop 
    809777       
    810778                     IF ( jn == 1 ) THEN   ;   ji_min = 2  
     
    816784                     zBV_prime(:,:) = 0._wp 
    817785 
    818                      DO ji = ji_min, jpi - 1, nn_nzebra_vp  
     786                     DO ji = ji_min, jpi - 1, nn_zebra_vp  
    819787                         
    820788                        !!! It is intentional to have a ji then jj loop for V-velocity 
     
    836804      
    837805                           ! right hand side 
    838                            zFV(ji,jj) = ( zrhsv(ji,jj) &              ! right-hand side terms 
    839                                &        + zAA3                        ! boundary condition translation 
    840                                &        + zDV(ji,jj) * v_ice(ji-1,jj) ! internal force, j-1 
    841                                &        + zEV(ji,jj) * v_ice(ji+1,jj) ) * vmask(ji,jj,1) ! internal force, j+1 
     806                           zFV(ji,jj) = ( zrhsv(ji,jj) &                                   ! right-hand side terms 
     807                               &        + zAA3                                           & ! boundary condition translation 
     808                               &        + zDV(ji,jj) * v_ice(ji-1,jj)                    & ! internal force, j-1 
     809                               &        + zEV(ji,jj) * v_ice(ji+1,jj) ) * vmask(ji,jj,1)   ! internal force, j+1 
    842810                     
    843811                        END DO 
     
    855823                        v_ice(ji,jpj-1)  = zFV_prime(ji,jpj-1) / zBV_prime(ji,jpj-1) * vmask(ji,jj,jpj-1)  ! last row 
    856824                        DO jj = jpj-2, 2, -1 ! can be turned into forward row by substituting jj if needed 
    857                            v_ice(ji,jj)   = zFV_prime(ji,jj) - zCV(ji,jj) * v_ice(ji,jj+1) * vmask(ji,jj) / zBV_prime(ji,jj) 
     825                           v_ice(ji,jj)   = zFV_prime(ji,jj) - zCV(ji,jj) * v_ice(ji,jj+1) * vmask(ji,jj,1) / zBV_prime(ji,jj) 
    858826                        END DO             
    859827             
     
    902870                   
    903871                  ! - Stop if error is too large ("safeguard against bad forcing" of original Zhang routine) 
    904                   IF ( i_inn > 1 && zuerr_max > rn_uerr_max_vp ) THEN 
    905                       IF ( lwp ) " VP rheology error was too large : ", zu_err_max, " in outer U-iteration ", i_out, " after ", i_inn, " iterations, we stopped " 
     872                  IF ( i_inn > 1 .AND. zuerr_max > rn_uerr_max_vp ) THEN 
     873                      IF ( lwp ) WRITE(numout,*) " VP rheology error was too large : ", zuerr_max, " in outer U-iteration ", i_out, " after ", i_inn, " iterations, we stopped " 
    906874                      ll_u_iterate = .FALSE. 
    907875                  ENDIF 
     
    909877                  ! - Stop if error small enough 
    910878                  IF ( zuerr_max < rn_uerr_min_vp ) THEN                                         
    911                       IF ( lwp ) " VP rheology nicely done in outer U-iteration ", i_out, " after ", i_inn, " iterations, finished! " 
     879                      IF ( lwp ) WRITE(numout,*) " VP rheology nicely done in outer U-iteration ", i_out, " after ", i_inn, " iterations, finished! " 
    912880                      ll_u_iterate = .FALSE. 
    913881                  ENDIF 
     
    933901                   
    934902                  ! - Stop if error is too large ("safeguard against bad forcing" of original Zhang routine) 
    935                   IF ( i_inn > 1 && zverr_max > rn_uerr_max_vp ) THEN 
    936                       IF ( lwp ) " VP rheology error was too large : ", zv_err_max, " in outer V-iteration ", i_out, " after ", i_inn, " iterations, we stopped " 
     903                  IF ( i_inn > 1 .AND. zverr_max > rn_uerr_max_vp ) THEN 
     904                      IF ( lwp ) WRITE(numout,*) " VP rheology error was too large : ", zverr_max, " in outer V-iteration ", i_out, " after ", i_inn, " iterations, we stopped " 
    937905                      ll_v_iterate = .FALSE. 
    938906                  ENDIF 
    939907                   
    940908                  ! - Stop if error small enough 
    941                   IF ( zverr_max < rn_verr_min ) THEN                                         
    942                       IF ( lwp ) " VP rheology nicely done in outer V-iteration ", i_out, " after ", i_inn, " iterations, finished! " 
     909                  IF ( zverr_max < rn_uerr_min_vp ) THEN                                         
     910                      IF ( lwp ) WRITE(numout,*) " VP rheology nicely done in outer V-iteration ", i_out, " after ", i_inn, " iterations, finished! " 
    943911                      ll_v_iterate = .FALSE. 
    944912                  ENDIF 
     
    965933      !------------------------------------------------------------------------------! 
    966934 
    967       END ! End of outer loop (i_out) ============================================================================================= 
     935      END DO ! End of outer loop (i_out) ============================================================================================= 
    968936 
    969937      !------------------------------------------------------------------------------! 
     
    10771045      IF(  iom_use('utau_oi') .OR. iom_use('vtau_oi') .OR. iom_use('utau_ai') .OR. iom_use('vtau_ai') .OR. & 
    10781046         & iom_use('utau_bi') .OR. iom_use('vtau_bi') ) THEN 
    1079          ! 
    1080          CALL lbc_lnk_multi( 'icedyn_rhg_vp', ztaux_oi, 'U', -1., ztauy_oi, 'V', -1., ztaux_ai, 'U', -1., ztauy_ai, 'V', -1., & 
    1081             &                                 ztaux_bi, 'U', -1., ztauy_bi, 'V', -1. ) 
     1047 
     1048         ALLOCATE( ztaux_oi(jpi,jpj) , ztauy_oi(jpi,jpj) ) 
     1049 
     1050         !--- Recalculate oceanic stress at last inner iteration 
     1051         DO jj = 2, jpj - 1 
     1052            DO ji = 2, jpi - 1 
     1053 
     1054                !--- ice u-velocity @V points, v-velocity @U points (for non-linear drag computation) 
     1055                zu_cV            = 0.25_wp * ( u_ice(ji,jj) + u_ice(ji-1,jj) + u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * vmask(ji,jj,1) 
     1056                zv_cU            = 0.25_wp * ( v_ice(ji,jj) + v_ice(ji,jj-1) + v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * umask(ji,jj,1) 
     1057                 
     1058                !--- non-linear drag coefficients (need to be updated at each outer loop, see Lemieux and Tremblay JGR09, p.3, beginning of Section 3) 
     1059                zCwU(ji,jj)          = za_iU(ji,jj) * zrhoco * SQRT( ( u_ice(ji,jj) - u_oce (ji,jj) ) * ( u_ice(ji,jj) - u_oce (ji,jj) )  & 
     1060                  &                                                + ( zv_cU - v_oceU(ji,jj) ) * ( zv_cU - v_oceU(ji,jj) ) ) 
     1061                zCwV(ji,jj)          = za_iV(ji,jj) * zrhoco * SQRT( ( v_ice(ji,jj) - v_oce (ji,jj) ) * ( v_ice(ji,jj) - v_oce (ji,jj) )  & 
     1062                  &                                                + ( zu_cV - u_oceV(ji,jj) ) * ( zu_cV - u_oceV(ji,jj) ) ) 
     1063                  
     1064                !--- Ocean-ice stress 
     1065                ztaux_oi(ji,jj) = zCwU(ji,jj) * ( u_oce(ji,jj) - u_ice(ji,jj) ) 
     1066                ztauy_oi(ji,jj) = zCwV(ji,jj) * ( v_oce(ji,jj) - v_ice(ji,jj) ) 
     1067                 
     1068            END DO 
     1069         END DO 
     1070          
     1071         ! 
     1072         CALL lbc_lnk_multi( 'icedyn_rhg_vp', ztaux_oi, 'U', -1., ztauy_oi, 'V', -1., ztaux_ai, 'U', -1., ztauy_ai, 'V', -1. ) !, & 
     1073!            &                                 ztaux_bi, 'U', -1., ztauy_bi, 'V', -1. ) 
    10821074         ! 
    10831075         CALL iom_put( 'utau_oi' , ztaux_oi * zmsk00 ) 
     
    10851077         CALL iom_put( 'utau_ai' , ztaux_ai * zmsk00 ) 
    10861078         CALL iom_put( 'vtau_ai' , ztauy_ai * zmsk00 ) 
    1087          CALL iom_put( 'utau_bi' , ztaux_bi * zmsk00 ) 
    1088          CALL iom_put( 'vtau_bi' , ztauy_bi * zmsk00 ) 
     1079!        CALL iom_put( 'utau_bi' , ztaux_bi * zmsk00 ) 
     1080!        CALL iom_put( 'vtau_bi' , ztauy_bi * zmsk00 ) 
     1081 
     1082         DEALLOCATE( ztaux_oi , ztauy_oi ) 
    10891083 
    10901084      ENDIF 
     
    11671161      IF(  iom_use('dssh_dx') .OR. iom_use('dssh_dy') .OR. & 
    11681162         & iom_use('corstrx') .OR. iom_use('corstry') ) THEN 
     1163 
     1164         ! --- Recalculate Coriolis stress at last inner iteration 
     1165         DO jj = 2, jpj - 1 
     1166            DO ji = 2, jpi - 1 
     1167                ! --- U-component  
     1168                zCorU(ji,jj)         =   0.25_wp * r1_e1u(ji,jj) *  & 
     1169                           &             ( zmf(ji  ,jj) * ( e1v(ji  ,jj) * v_ice(ji  ,jj) + e1v(ji  ,jj-1) * v_ice(ji  ,jj-1) )  & 
     1170                           &             + zmf(ji+1,jj) * ( e1v(ji+1,jj) * v_ice(ji+1,jj) + e1v(ji+1,jj-1) * v_ice(ji+1,jj-1) ) ) 
     1171                zCorV(ji,jj)         = - 0.25_wp * r1_e2v(ji,jj) *  & 
     1172                           &             ( zmf(ji,jj  ) * ( e2u(ji,jj  ) * u_ice(ji,jj  ) + e2u(ji-1,jj  ) * u_ice(ji-1,jj  ) )  & 
     1173                           &             + zmf(ji,jj+1) * ( e2u(ji,jj+1) * u_ice(ji,jj+1) + e2u(ji-1,jj+1) * u_ice(ji-1,jj+1) ) ) 
     1174            END DO 
     1175         END DO 
    11691176         ! 
    11701177         CALL lbc_lnk_multi( 'icedyn_rhg_vp', zspgU, 'U', -1., zspgV, 'V', -1., & 
    11711178            &                                 zCorU, 'U', -1., zCorV, 'V', -1. ) 
    1172  
     1179         ! 
    11731180         CALL iom_put( 'dssh_dx' , zspgU * zmsk00 )   ! Sea-surface tilt term in force balance (x) 
    11741181         CALL iom_put( 'dssh_dy' , zspgV * zmsk00 )   ! Sea-surface tilt term in force balance (y) 
    11751182         CALL iom_put( 'corstrx' , zCorU * zmsk00 )   ! Coriolis force term in force balance (x) 
    11761183         CALL iom_put( 'corstry' , zCorV * zmsk00 )   ! Coriolis force term in force balance (y) 
     1184 
    11771185      ENDIF 
    11781186       
    11791187      IF ( iom_use('intstrx') .OR. iom_use('intstry') ) THEN 
    11801188 
    1181  
    1182          ! Recalculate internal forces (divergence of stress tensor) 
     1189         ! Recalculate internal forces (divergence of stress tensor) at last inner iteration 
    11831190         DO jj = 2, jpj - 1 
    11841191            DO ji = 2, jpi - 1 
     
    12911298      INTEGER           ::   it, idtime, istatus 
    12921299      INTEGER           ::   ji, jj          ! dummy loop indices 
    1293       REAL(wp)          ::   zveldif, zu_res_mean, zv_res_mean, zmke, zu, zv ! local scalars 
    1294       REAL(wp), DIMENSION(jpi,jpj)   ::   zu_res(:,:), zv_res(:,:), zvel2(:,:) ! local arrays 
     1300      REAL(wp)          ::   zveldif, zu_res_mean, zv_res_mean, zvelres, zmke, zu, zv ! local scalars 
     1301      REAL(wp), DIMENSION(jpi,jpj) ::   zu_res, zv_res, zvel2                        ! local arrays 
    12951302                                                                              
    12961303      CHARACTER(len=20) ::   clname 
    1297       ::   zres           ! check convergence 
    12981304      !!---------------------------------------------------------------------- 
    12991305 
     
    13521358      DO jj = 2, jpj - 1 
    13531359         DO ji = 2, jpi - 1                                       
    1354             zu_res(ji,jj)  = zrhsu(ji,jj) + zDU(ji,jj) * pu(ji,jj-1) + zEU(ji,jj) * pu(ji,jj+1) & 
    1355                &           - zAU(ji,jj) * pu(ji-1,jj) - zBU(ji,jj) * pu(ji,jj) - zCU(ji,jj) * pu(ji+1,jj) 
     1360            zu_res(ji,jj)  = prhsu(ji,jj) + pDU(ji,jj) * pu(ji,jj-1) + pEU(ji,jj) * pu(ji,jj+1) & 
     1361               &           - pAU(ji,jj) * pu(ji-1,jj) - pBU(ji,jj) * pu(ji,jj) - pCU(ji,jj) * pu(ji+1,jj) 
    13561362                            
    1357             zv_res(ji,jj)  = zrhsv(ji,jj) + zDV(ji,jj) * pv(ji-1,jj) + zEV(ji,jj) * pv(ji+1,jj) & 
    1358                &           - zAV(ji,jj) * pv(ji,jj-1) - zBV(ji,jj) * pv(ji,jj) - zCV(ji,jj) * pv(ji,jj+1)                            
     1363            zv_res(ji,jj)  = prhsv(ji,jj) + pDV(ji,jj) * pv(ji-1,jj) + pEV(ji,jj) * pv(ji+1,jj) & 
     1364               &           - pAV(ji,jj) * pv(ji,jj-1) - pBV(ji,jj) * pv(ji,jj) - pCV(ji,jj) * pv(ji,jj+1)                            
    13591365          END DO 
    13601366       END DO                   
    1361        zu_res_mean = glob_sum( 'ice_rhg_vp', zu_res(:,:) * zu_res(:,:) * e1e2u(:,:) * umask(:,:) ) 
    1362        zv_res_mean = glob_sum( 'ice_rhg_vp', zv_res(:,:) * zv_res(:,:) * e1e2v(:,:) * vmask(:,:) ) 
    1363        zu_res_mean = SQRT( zu_resmean / pglob_area ) 
    1364        zv_res_mean = SQRT( zv_resmean / pglob_area ) 
    1365        zvelres     = MEAN( zu_res_mean, zv_res_mean ) 
     1367       zu_res_mean = glob_sum( 'ice_rhg_vp', zu_res(:,:) * zu_res(:,:) * e1e2u(:,:) * umask(:,:,1) ) 
     1368       zv_res_mean = glob_sum( 'ice_rhg_vp', zv_res(:,:) * zv_res(:,:) * e1e2v(:,:) * vmask(:,:,1) ) 
     1369       zu_res_mean = SQRT( zu_res_mean / pglob_area ) 
     1370       zv_res_mean = SQRT( zv_res_mean / pglob_area ) 
     1371       zvelres     = 0.5_wp * ( zu_res_mean + zv_res_mean ) 
    13661372                                           
    13671373       ! -- Global mean kinetic energy per unit area (J/m2)   
     
    13741380       END DO 
    13751381        
    1376        zmke = 0.5_wp * globsum( 'ice_rhg_vp', zmt(:,:) * e1e2t(:,:) * zvel2(:,:) ) 
     1382       zmke = 0.5_wp * glob_sum( 'ice_rhg_vp', pmt(:,:) * e1e2t(:,:) * zvel2(:,:) ) 
    13771383                   
    13781384         !                                                ! ==================== ! 
Note: See TracChangeset for help on using the changeset viewer.