Ignore:
Timestamp:
2015-07-15T17:46:12+02:00 (5 years ago)
Author:
andrewryan
Message:

merged in latest version of trunk alongside changes to SAO_SRC to be compatible with latest OBS

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/2014/dev_r4650_UKMO14.12_STAND_ALONE_OBSOPER/NEMOGCM/NEMO/LIM_SRC_3/limdyn.F90

    r5034 r5600  
    66   !! history :  1.0  ! 2002-08  (C. Ethe, G. Madec)  original VP code  
    77   !!            3.0  ! 2007-03  (MA Morales Maqueda, S. Bouillon, M. Vancoppenolle)  LIM3: EVP-Cgrid 
    8    !!            4.0  ! 2011-02  (G. Madec) dynamical allocation 
     8   !!            3.5  ! 2011-02  (G. Madec) dynamical allocation 
    99   !!---------------------------------------------------------------------- 
    1010#if defined key_lim3 
     
    2020   USE sbc_ice          ! Surface boundary condition: ice   fields 
    2121   USE ice              ! LIM-3 variables 
    22    USE par_ice          ! LIM-3 parameters 
    2322   USE dom_ice          ! LIM-3 domain 
    2423   USE limrhg           ! LIM-3 rheology 
     
    3130   USE timing          ! Timing 
    3231   USE limcons        ! conservation tests 
     32   USE limvar 
    3333 
    3434   IMPLICIT NONE 
     
    7676      CALL wrk_alloc( jpj, zswitch, zmsk ) 
    7777 
     78      CALL lim_var_agg(1)             ! aggregate ice categories 
     79 
    7880      IF( kt == nit000 )   CALL lim_dyn_init   ! Initialization (first time-step only) 
    7981 
     
    8385         IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limdyn', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 
    8486 
    85          u_ice_b(:,:) = u_ice(:,:) * tmu(:,:) 
    86          v_ice_b(:,:) = v_ice(:,:) * tmv(:,:) 
     87         u_ice_b(:,:) = u_ice(:,:) * umask(:,:,1) 
     88         v_ice_b(:,:) = v_ice(:,:) * vmask(:,:,1) 
    8789 
    8890         ! Rheology (ice dynamics) 
     
    101103            DO jj = 1, jpj 
    102104               zswitch(jj) = SUM( 1.0 - at_i(:,jj) )   ! = REAL(jpj) if ocean everywhere on a j-line 
    103                zmsk(jj) = SUM( tmask(:,jj,1)    )   ! = 0         if land  everywhere on a j-line 
     105               zmsk   (jj) = SUM( tmask(:,jj,1)    )   ! = 0         if land  everywhere on a j-line 
    104106            END DO 
    105107 
     
    157159         zv_io(:,:) = v_ice(:,:) - ssv_m(:,:) 
    158160         ! frictional velocity at T-point 
    159          zcoef = 0.5_wp * cw 
     161         zcoef = 0.5_wp * rn_cio 
    160162         DO jj = 2, jpjm1  
    161163            DO ji = fs_2, fs_jpim1   ! vector opt. 
    162164               ust2s(ji,jj) = zcoef * (  zu_io(ji,jj) * zu_io(ji,jj) + zu_io(ji-1,jj) * zu_io(ji-1,jj)   & 
    163                   &                    + zv_io(ji,jj) * zv_io(ji,jj) + zv_io(ji,jj-1) * zv_io(ji,jj-1)   ) * tms(ji,jj) 
     165                  &                    + zv_io(ji,jj) * zv_io(ji,jj) + zv_io(ji,jj-1) * zv_io(ji,jj-1) ) * tmask(ji,jj,1) 
    164166            END DO 
    165167         END DO 
     
    170172      ELSE      ! no ice dynamics : transmit directly the atmospheric stress to the ocean 
    171173         ! 
    172          zcoef = SQRT( 0.5_wp ) / rau0 
     174         zcoef = SQRT( 0.5_wp ) * r1_rau0 
    173175         DO jj = 2, jpjm1 
    174176            DO ji = fs_2, fs_jpim1   ! vector opt. 
    175177               ust2s(ji,jj) = zcoef * SQRT(  utau(ji,jj) * utau(ji,jj) + utau(ji-1,jj) * utau(ji-1,jj)   & 
    176                   &                        + vtau(ji,jj) * vtau(ji,jj) + vtau(ji,jj-1) * vtau(ji,jj-1)   ) * tms(ji,jj) 
     178                  &                        + vtau(ji,jj) * vtau(ji,jj) + vtau(ji,jj-1) * vtau(ji,jj-1) ) * tmask(ji,jj,1) 
    177179            END DO 
    178180         END DO 
     
    189191         CALL prt_ctl(tab2d_1=delta_i   , clinfo1=' lim_dyn  : delta_i   :') 
    190192         CALL prt_ctl(tab2d_1=strength  , clinfo1=' lim_dyn  : strength  :') 
    191          CALL prt_ctl(tab2d_1=area      , clinfo1=' lim_dyn  : cell area :') 
     193         CALL prt_ctl(tab2d_1=e12t      , clinfo1=' lim_dyn  : cell area :') 
    192194         CALL prt_ctl(tab2d_1=at_i      , clinfo1=' lim_dyn  : at_i      :') 
    193195         CALL prt_ctl(tab2d_1=vt_i      , clinfo1=' lim_dyn  : vt_i      :') 
     
    241243      !!------------------------------------------------------------------- 
    242244      INTEGER  ::   ios                 ! Local integer output status for namelist read 
    243       NAMELIST/namicedyn/ epsd, om, cw, pstar,   & 
    244          &                c_rhg, creepl, ecc, ahi0,     & 
    245          &                nevp, relast, alphaevp, hminrhg 
     245      NAMELIST/namicedyn/ nn_icestr, ln_icestr_bvf, rn_pe_rdg, rn_pstar, rn_crhg, rn_cio,  rn_creepl, rn_ecc, & 
     246         &                nn_nevp, rn_relast, nn_ahi0, rn_ahi0_ref 
     247      INTEGER  ::   ji, jj 
     248      REAL(wp) ::   za00, zd_max 
    246249      !!------------------------------------------------------------------- 
    247250 
     
    259262         WRITE(numout,*) 'lim_dyn_init : ice parameters for ice dynamics ' 
    260263         WRITE(numout,*) '~~~~~~~~~~~~' 
    261          WRITE(numout,*) '   tolerance parameter                              epsd   = ', epsd 
    262          WRITE(numout,*) '   relaxation constant                              om     = ', om 
    263          WRITE(numout,*) '   drag coefficient for oceanic stress              cw     = ', cw 
    264          WRITE(numout,*) '   first bulk-rheology parameter                    pstar  = ', pstar 
    265          WRITE(numout,*) '   second bulk-rhelogy parameter                    c_rhg  = ', c_rhg 
    266          WRITE(numout,*) '   creep limit                                      creepl = ', creepl 
    267          WRITE(numout,*) '   eccentricity of the elliptical yield curve       ecc    = ', ecc 
    268          WRITE(numout,*) '   horizontal diffusivity coeff. for sea-ice        ahi0   = ', ahi0 
    269          WRITE(numout,*) '   number of iterations for subcycling              nevp   = ', nevp 
    270          WRITE(numout,*) '   ratio of elastic timescale over ice time step    relast = ', relast 
    271          WRITE(numout,*) '   coefficient for the solution of int. stresses  alphaevp = ', alphaevp 
    272          WRITE(numout,*) '   min ice thickness for rheology calculations     hminrhg = ', hminrhg 
     264         WRITE(numout,*)'    ice strength parameterization (0=Hibler 1=Rothrock)  nn_icestr     = ', nn_icestr  
     265         WRITE(numout,*)'    Including brine volume in ice strength comp.         ln_icestr_bvf = ', ln_icestr_bvf 
     266         WRITE(numout,*)'   Ratio of ridging work to PotEner change in ridging    rn_pe_rdg     = ', rn_pe_rdg  
     267         WRITE(numout,*) '   drag coefficient for oceanic stress                  rn_cio        = ', rn_cio 
     268         WRITE(numout,*) '   first bulk-rheology parameter                        rn_pstar      = ', rn_pstar 
     269         WRITE(numout,*) '   second bulk-rhelogy parameter                        rn_crhg       = ', rn_crhg 
     270         WRITE(numout,*) '   creep limit                                          rn_creepl     = ', rn_creepl 
     271         WRITE(numout,*) '   eccentricity of the elliptical yield curve           rn_ecc        = ', rn_ecc 
     272         WRITE(numout,*) '   number of iterations for subcycling                  nn_nevp       = ', nn_nevp 
     273         WRITE(numout,*) '   ratio of elastic timescale over ice time step        rn_relast     = ', rn_relast 
     274         WRITE(numout,*) '   horizontal diffusivity calculation                   nn_ahi0       = ', nn_ahi0 
     275         WRITE(numout,*) '   horizontal diffusivity coeff. (orca2 grid)           rn_ahi0_ref   = ', rn_ahi0_ref 
    273276      ENDIF 
    274277      ! 
    275       usecc2 = 1._wp / ( ecc * ecc ) 
    276       rhoco  = rau0  * cw 
    277  
    278       ! elastic damping 
    279       telast = relast * rdt_ice 
    280  
    281       !  Diffusion coefficients. 
    282       ahiu(:,:) = ahi0 * umask(:,:,1) 
    283       ahiv(:,:) = ahi0 * vmask(:,:,1) 
    284       ! 
     278      usecc2 = 1._wp / ( rn_ecc * rn_ecc ) 
     279      rhoco  = rau0  * rn_cio 
     280      ! 
     281      !  Diffusion coefficients 
     282      SELECT CASE( nn_ahi0 ) 
     283 
     284      CASE( 0 ) 
     285         ahiu(:,:) = rn_ahi0_ref 
     286         ahiv(:,:) = rn_ahi0_ref 
     287 
     288         IF(lwp) WRITE(numout,*) '' 
     289         IF(lwp) WRITE(numout,*) '   laplacian operator: ahim constant = rn_ahi0_ref' 
     290 
     291      CASE( 1 )  
     292 
     293         zd_max = MAX( MAXVAL( e1t(:,:) ), MAXVAL( e2t(:,:) ) ) 
     294         IF( lk_mpp )   CALL mpp_max( zd_max )          ! max over the global domain 
     295          
     296         ahiu(:,:) = rn_ahi0_ref * zd_max * 1.e-05_wp   ! 1.e05 = 100km = max grid space at 60° latitude in orca2 
     297                                                        !                    (60° = min latitude for ice cover)   
     298         ahiv(:,:) = rn_ahi0_ref * zd_max * 1.e-05_wp 
     299 
     300         IF(lwp) WRITE(numout,*) '' 
     301         IF(lwp) WRITE(numout,*) '   laplacian operator: ahim proportional to max of e1 e2 over the domain (', zd_max, ')' 
     302         IF(lwp) WRITE(numout,*) '   value for ahim = ', rn_ahi0_ref * zd_max * 1.e-05_wp  
     303          
     304      CASE( 2 )  
     305 
     306         zd_max = MAX( MAXVAL( e1t(:,:) ), MAXVAL( e2t(:,:) ) ) 
     307         IF( lk_mpp )   CALL mpp_max( zd_max )   ! max over the global domain 
     308          
     309         za00 = rn_ahi0_ref * 1.e-05_wp          ! 1.e05 = 100km = max grid space at 60° latitude in orca2 
     310                                                 !                    (60° = min latitude for ice cover)   
     311         DO jj = 1, jpj 
     312            DO ji = 1, jpi 
     313               ahiu(ji,jj) = za00 * MAX( e1t(ji,jj), e2t(ji,jj) ) * umask(ji,jj,1) 
     314               ahiv(ji,jj) = za00 * MAX( e1f(ji,jj), e2f(ji,jj) ) * vmask(ji,jj,1) 
     315            END DO 
     316         END DO 
     317         ! 
     318         IF(lwp) WRITE(numout,*) '' 
     319         IF(lwp) WRITE(numout,*) '   laplacian operator: ahim proportional to e1' 
     320         IF(lwp) WRITE(numout,*) '   maximum grid-spacing = ', zd_max, ' maximum value for ahim = ', za00*zd_max 
     321          
     322      END SELECT 
     323 
    285324   END SUBROUTINE lim_dyn_init 
    286325 
Note: See TracChangeset for help on using the changeset viewer.