Changeset 8291 for branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO

Timestamp:

2017-07-06T11:49:56+02:00 (7 years ago)

Author:

clem

Message:

correct last commit changes. This is the start version for the modifications of the ice model. All sette tests passed except restartability for agrif and sas (though sas_biper is restartable). Configurations CREG, SPITZ, FER & SASBIPER are running.

Location:

branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO

Files:

• LIM_SRC_3/ice.F90 (modified) (1 diff)
• LIM_SRC_3/limistate.F90 (modified) (1 diff)
• LIM_SRC_3/limwri.F90 (modified) (4 diffs)
• OPA_SRC/SBC/sbcice_lim.F90 (modified) (1 diff)
• SAS_SRC/sbcssm.F90 (modified) (1 diff)

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

@@ -640 +640 @@
                 diag_intstrx(jpi,jpj) , diag_intstry(jpi,jpj) ,    &
                 diag_sig1(jpi,jpj)    , diag_sig2(jpi,jpj)    ,    &
-                STAT = ierr(ii) )
+                diag_shear(jpi,jpj)   , STAT = ierr(ii) )
 
       ice_alloc = MAXVAL( ierr(:) )

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

@@ -531 +531 @@
       INTEGER :: ifpr, ierror
       !
-      CHARACTER(len=100) ::  cn_dir          ! Root directory for location of ice files
+      CHARACTER(len=256) ::  cn_dir          ! Root directory for location of ice files
       TYPE(FLD_N)                    ::   sn_hti, sn_hts, sn_ati, sn_tsu, sn_tmi, sn_smi
       TYPE(FLD_N), DIMENSION(jpfldi) ::   slf_i                 ! array of namelist informations on the fields to read

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

@@ -54 +54 @@
       INTEGER  ::  ji, jj, jk, jl  ! dummy loop indices
       REAL(wp) ::  z2da, z2db, ztmp, zrho1, zrho2, zmiss_val
+      REAL(wp) ::  zs12, zshear
       REAL(wp), POINTER, DIMENSION(:,:,:) ::  zswi2, zmiss2
       REAL(wp), POINTER, DIMENSION(:,:)   ::  z2d, zswi, zmiss !  2D workspace
       REAL(wp), POINTER, DIMENSION(:,:)   ::  zfb              ! ice freeboard
       REAL(wp), POINTER, DIMENSION(:,:)   ::  zamask, zamask15 ! 15% concentration mask
+      REAL(wp), POINTER, DIMENSION(:,:)   ::  zsig1, zsig2, zsig3
 
       ! Global ice diagnostics (SIMIP)
@@ -71 +73 @@
       IF( nn_timing == 1 )  CALL timing_start('limwri')
 
-      CALL wrk_alloc( jpi, jpj, jpl, zswi2, zmiss2 )
-      CALL wrk_alloc( jpi, jpj     , z2d, zswi, zmiss )
-      CALL wrk_alloc( jpi, jpj     , zfb, zamask, zamask15 )
+      CALL wrk_alloc( jpi,jpj, jpl, zswi2, zmiss2 )
+      CALL wrk_alloc( jpi,jpj     , z2d, zswi, zmiss )
+      CALL wrk_alloc( jpi,jpj     , zfb, zamask, zamask15 )
+      CALL wrk_alloc( jpi,jpj     , zsig1, zsig2, zsig3 )
 
       !----------------------------------------
@@ -226 +229 @@
       CALL iom_put ('hfxspr'     , hfx_spr(:,:)         )   ! Heat content of snow precip 
 
+      ! specific outputs for EVP rheology
+      IF( iom_use( "isig1" ) .OR. iom_use( "isig2" ) .OR. iom_use( "isig3" ) ) THEN
+         zsig1(:,:) = 0._wp; zsig2(:,:) = 0._wp; zsig3(:,:) = 0._wp;
+         DO jj = 2, jpjm1
+            DO ji = 2, jpim1
+               zs12 = ( zswi(ji-1,jj) * stress12_i(ji-1,jj) + zswi(ji  ,jj-1) * stress12_i(ji  ,jj-1) +  &  ! stress12_i at T-point
+                  &     zswi(ji  ,jj) * stress12_i(ji  ,jj) + zswi(ji-1,jj-1) * stress12_i(ji-1,jj-1) )  &
+                  &   / MAX( 1._wp, zswi(ji-1,jj) + zswi(ji,jj-1) + zswi(ji,jj) + zswi(ji-1,jj-1) )
+
+               zshear = SQRT( stress2_i(ji,jj) * stress2_i(ji,jj) + 4._wp * zs12 * zs12 ) ! shear stress  
+
+               z2da = zswi(ji,jj) / MAX( 1._wp, strength(ji,jj) )
+
+!!               zsig1(ji,jj) = 0.5_wp * z2da * ( stress1_i(ji,jj) + zshear ) ! principal stress (y-direction, see Hunke & Dukowicz 2002)
+!!               zsig2(ji,jj) = 0.5_wp * z2da * ( stress1_i(ji,jj) - zshear ) ! principal stress (x-direction, see Hunke & Dukowicz 2002)
+!!               zsig3(ji,jj) = z2da**2 * ( ( stress1_i(ji,jj) + strength(ji,jj) )**2 + ( rn_ecc * zshear )**2 ) ! quadratic relation linking compressive stress to shear stress
+!!                                                                                                               ! (scheme converges if this value is ~1, see Bouillon et al 2009 (eq. 11))
+               zsig1(ji,jj) = 0.5_wp * z2da * ( stress1_i(ji,jj) )          ! compressive stress, see Bouillon et al. 2015
+               zsig2(ji,jj) = 0.5_wp * z2da * ( zshear )                    ! shear stress
+               zsig3(ji,jj) = z2da**2 * ( ( stress1_i(ji,jj) + strength(ji,jj) )**2 + ( rn_ecc * zshear )**2 )
+            END DO
+         END DO
+         CALL lbc_lnk_multi( zsig1, 'T', 1., zsig2, 'T', 1., zsig3, 'T', 1. )
+         CALL iom_put( "isig1" , zsig1 )
+         CALL iom_put( "isig2" , zsig2 )
+         CALL iom_put( "isig3" , zsig3 )
+      ENDIF
+
       ! MV MP 2016
       IF ( ln_pnd ) THEN
@@ -400 +431 @@
       CALL wrk_dealloc( jpi, jpj     , z2d, zswi, zmiss )
       CALL wrk_dealloc( jpi, jpj     , zfb, zamask, zamask15 )
+      CALL wrk_dealloc( jpi, jpj     , zsig1, zsig2, zsig3 )
 
       IF( nn_timing == 1 )  CALL timing_stop('limwri')

branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90

@@ -179 +179 @@
                u_ice(:,:) = rn_uice * umask(:,:,1)             !     or prescribed velocity
                v_ice(:,:) = rn_vice * vmask(:,:,1)
+               !!CALL RANDOM_NUMBER(u_ice(:,:))
+               !!CALL RANDOM_NUMBER(v_ice(:,:))
             ENDIF
                                       CALL lim_trp( kt )       ! -- Ice transport (Advection/diffusion)

branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/SAS_SRC/sbcssm.F90

@@ -115 +115 @@
       ELSE
          sst_m(:,:) = 0._wp
-         sss_m(:,:) = 0._wp
+         sss_m(:,:) = 35._wp
          ssu_m(:,:) = 0._wp
          ssv_m(:,:) = 0._wp