Changeset 7597


Ignore:
Timestamp:
2017-01-20T19:40:06+01:00 (5 years ago)
Author:
clem
Message:

add the possibility to run without sea ice diffusion (by default)

Location:
branches/2016/v3_6_CMIP6_ice_diagnostics/NEMOGCM
Files:
4 edited

Legend:

Unmodified
Added
Removed
  • branches/2016/v3_6_CMIP6_ice_diagnostics/NEMOGCM/CONFIG/SHARED/namelist_ice_lim3_ref

    r7596 r7597  
    6969   rn_relast      =    0.333       !  ratio of elastic timescale to ice time step: Telast = dt_ice * rn_relast  
    7070                                   !     advised value: 1/3 (rn_nevp=120) or 1/9 (rn_nevp=300) 
    71    nn_ahi0        =    2           !  horizontal diffusivity computation 
    72                                    !     0: use rn_ahi0_ref 
    73                                    !     1: use rn_ahi0_ref x mean grid cell length / ( 2deg mean grid cell length ) 
    74                                    !     2: use rn_ahi0_ref x grid cell length      / ( 2deg mean grid cell length ) 
    75    rn_ahi0_ref    = 350.0          !  horizontal sea ice diffusivity (m2/s)  
    76                                    !     if nn_ahi0 > 0, rn_ahi0_ref is the reference value at a nominal 2 deg resolution 
    7771/ 
    7872!------------------------------------------------------------------------------ 
    7973&namicehdf     !   Ice horizontal diffusion 
    8074!------------------------------------------------------------------------------ 
     75   nn_ahi0        =    -1          !  horizontal diffusivity computation 
     76                                   !    -1: no diffusion (bypass limhdf) 
     77                                   !     0: use rn_ahi0_ref 
     78                                   !     1: use rn_ahi0_ref x mean grid cell length / ( 2deg mean grid cell length ) 
     79                                   !     2: use rn_ahi0_ref x grid cell length      / ( 2deg mean grid cell length ) 
     80   rn_ahi0_ref    = 350.0          !  horizontal sea ice diffusivity (m2/s) 
     81                                   !     if nn_ahi0 > 0, rn_ahi0_ref is the reference value at a nominal 2 deg resolution 
    8182   nn_convfrq     = 5              !  convergence check frequency of the Crant-Nicholson scheme (perf. optimization) 
    8283/ 
  • branches/2016/v3_6_CMIP6_ice_diagnostics/NEMOGCM/NEMO/LIM_SRC_3/limdyn.F90

    r5123 r7597  
    244244      INTEGER  ::   ios                 ! Local integer output status for namelist read 
    245245      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 
     246         &                nn_nevp, rn_relast 
    249247      !!------------------------------------------------------------------- 
    250248 
     
    272270         WRITE(numout,*) '   number of iterations for subcycling                  nn_nevp       = ', nn_nevp 
    273271         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 
    276272      ENDIF 
    277273      ! 
     
    279275      rhoco  = rau0  * rn_cio 
    280276      ! 
    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  
    324277   END SUBROUTINE lim_dyn_init 
    325278 
  • branches/2016/v3_6_CMIP6_ice_diagnostics/NEMOGCM/NEMO/LIM_SRC_3/limhdf.F90

    r6476 r7597  
    233233      !!------------------------------------------------------------------- 
    234234      INTEGER  ::   ios                 ! Local integer output status for namelist read 
    235       NAMELIST/namicehdf/ nn_convfrq  
    236       !!------------------------------------------------------------------- 
    237       ! 
    238       IF(lwp) THEN 
    239          WRITE(numout,*) 
    240          WRITE(numout,*) 'lim_hdf : Ice horizontal diffusion' 
    241          WRITE(numout,*) '~~~~~~~' 
    242       ENDIF 
     235      NAMELIST/namicehdf/  nn_ahi0, rn_ahi0_ref, nn_convfrq  
     236      INTEGER  ::   ji, jj 
     237      REAL(wp) ::   za00, zd_max 
     238      !!------------------------------------------------------------------- 
    243239      ! 
    244240      REWIND( numnam_ice_ref )              ! Namelist namicehdf in reference namelist : Ice horizontal diffusion 
     
    253249      IF(lwp) THEN                          ! control print 
    254250         WRITE(numout,*) 
    255          WRITE(numout,*)'   Namelist of ice parameters for ice horizontal diffusion computation ' 
    256          WRITE(numout,*)'      convergence check frequency of the Crant-Nicholson scheme   nn_convfrq   = ', nn_convfrq 
     251         WRITE(numout,*) 'lim_hdf_init : Ice horizontal diffusion' 
     252         WRITE(numout,*) '~~~~~~~~~~~' 
     253         WRITE(numout,*) '   horizontal diffusivity calculation                          nn_ahi0      = ', nn_ahi0 
     254         WRITE(numout,*) '   horizontal diffusivity coeff. (orca2 grid)                  rn_ahi0_ref  = ', rn_ahi0_ref 
     255         WRITE(numout,*) '   convergence check frequency of the Crant-Nicholson scheme   nn_convfrq   = ', nn_convfrq 
    257256      ENDIF 
     257      ! 
     258      !  Diffusion coefficients 
     259      SELECT CASE( nn_ahi0 ) 
     260 
     261      CASE( -1 ) 
     262         ahiu(:,:) = 0._wp 
     263         ahiv(:,:) = 0._wp 
     264 
     265         IF(lwp) WRITE(numout,*) '' 
     266         IF(lwp) WRITE(numout,*) '   No sea-ice diffusion applied' 
     267 
     268      CASE( 0 ) 
     269         ahiu(:,:) = rn_ahi0_ref 
     270         ahiv(:,:) = rn_ahi0_ref 
     271 
     272         IF(lwp) WRITE(numout,*) '' 
     273         IF(lwp) WRITE(numout,*) '   laplacian operator: ahim constant = rn_ahi0_ref' 
     274 
     275      CASE( 1 )  
     276 
     277         zd_max = MAX( MAXVAL( e1t(:,:) ), MAXVAL( e2t(:,:) ) ) 
     278         IF( lk_mpp )   CALL mpp_max( zd_max )          ! max over the global domain 
     279          
     280         ahiu(:,:) = rn_ahi0_ref * zd_max * 1.e-05_wp   ! 1.e05 = 100km = max grid space at 60deg latitude in orca2 
     281                                                        !                    (60deg = min latitude for ice cover)   
     282         ahiv(:,:) = rn_ahi0_ref * zd_max * 1.e-05_wp 
     283 
     284         IF(lwp) WRITE(numout,*) '' 
     285         IF(lwp) WRITE(numout,*) '   laplacian operator: ahim proportional to max of e1 e2 over the domain (', zd_max, ')' 
     286         IF(lwp) WRITE(numout,*) '   value for ahim = ', rn_ahi0_ref * zd_max * 1.e-05_wp  
     287          
     288      CASE( 2 )  
     289 
     290         zd_max = MAX( MAXVAL( e1t(:,:) ), MAXVAL( e2t(:,:) ) ) 
     291         IF( lk_mpp )   CALL mpp_max( zd_max )   ! max over the global domain 
     292          
     293         za00 = rn_ahi0_ref * 1.e-05_wp          ! 1.e05 = 100km = max grid space at 60deg latitude in orca2 
     294                                                 !                    (60deg = min latitude for ice cover)   
     295         DO jj = 1, jpj 
     296            DO ji = 1, jpi 
     297               ahiu(ji,jj) = za00 * MAX( e1t(ji,jj), e2t(ji,jj) ) * umask(ji,jj,1) 
     298               ahiv(ji,jj) = za00 * MAX( e1f(ji,jj), e2f(ji,jj) ) * vmask(ji,jj,1) 
     299            END DO 
     300         END DO 
     301         ! 
     302         IF(lwp) WRITE(numout,*) '' 
     303         IF(lwp) WRITE(numout,*) '   laplacian operator: ahim proportional to e1' 
     304         IF(lwp) WRITE(numout,*) '   maximum grid-spacing = ', zd_max, ' maximum value for ahim = ', za00*zd_max 
     305          
     306      END SELECT 
    258307      ! 
    259308   END SUBROUTINE lim_hdf_init 
  • branches/2016/v3_6_CMIP6_ice_diagnostics/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90

    r7506 r7597  
    325325         ! Diffusion of Ice fields                   
    326326         !------------------------------------------------------------------------------! 
    327          !------------------------------------ 
    328          !  Diffusion of other ice variables 
    329          !------------------------------------ 
    330          jm=1 
    331          DO jl = 1, jpl 
    332          !                             ! Masked eddy diffusivity coefficient at ocean U- and V-points 
    333          !   DO jj = 1, jpjm1                 ! NB: has not to be defined on jpj line and jpi row 
    334          !      DO ji = 1 , fs_jpim1   ! vector opt. 
    335          !         pahu(ji,jj) = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -a_i(ji  ,jj,jl) ) ) )   & 
    336          !            &        * ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -a_i(ji+1,jj,jl) ) ) ) * ahiu(ji,jj) 
    337          !         pahv(ji,jj) = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -a_i(ji,jj  ,jl) ) ) )   & 
    338          !            &        * ( 1._wp - MAX( 0._wp, SIGN( 1._wp,- a_i(ji,jj+1,jl) ) ) ) * ahiv(ji,jj) 
    339          !      END DO 
    340          !   END DO 
    341             DO jj = 1, jpjm1                 ! NB: has not to be defined on jpj line and jpi row 
     327         IF( nn_ahi0 /= -1 ) THEN 
     328 
     329            ! --- Prepare diffusion for variables with categories --- ! 
     330            !     mask eddy diffusivity coefficient at ocean U- and V-points 
     331            jm=1 
     332            DO jl = 1, jpl 
     333               DO jj = 1, jpjm1                 ! NB: has not to be defined on jpj line and jpi row 
     334                  DO ji = 1 , fs_jpim1   ! vector opt. 
     335                     pahu3D(ji,jj,jl) = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -a_i(ji  ,jj,  jl ) ) ) )   & 
     336                        &                * ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -a_i(ji+1,jj,  jl ) ) ) ) * ahiu(ji,jj) 
     337                     pahv3D(ji,jj,jl) = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -a_i(ji,  jj,  jl ) ) ) )   & 
     338                        &                * ( 1._wp - MAX( 0._wp, SIGN( 1._wp,- a_i(ji,  jj+1,jl ) ) ) ) * ahiv(ji,jj) 
     339                  END DO 
     340               END DO 
     341                
     342               zhdfptab(:,:,jm)= a_i  (:,:,  jl); jm = jm + 1     
     343               zhdfptab(:,:,jm)= v_i  (:,:,  jl); jm = jm + 1 
     344               zhdfptab(:,:,jm)= v_s  (:,:,  jl); jm = jm + 1  
     345               zhdfptab(:,:,jm)= smv_i(:,:,  jl); jm = jm + 1 
     346               zhdfptab(:,:,jm)= oa_i (:,:,  jl); jm = jm + 1 
     347               zhdfptab(:,:,jm)= e_s  (:,:,1,jl); jm = jm + 1 
     348               ! Sample of adding more variables to apply lim_hdf (ihdf_vars must be increased) 
     349               !   zhdfptab(:,:,jm) = variable_1 (:,:,1,jl); jm = jm + 1   
     350               !   zhdfptab(:,:,jm) = variable_2 (:,:,1,jl); jm = jm + 1  
     351               DO jk = 1, nlay_i 
     352                  zhdfptab(:,:,jm)=e_i(:,:,jk,jl); jm= jm+1 
     353               END DO 
     354            END DO 
     355            ! 
     356            ! --- Prepare diffusion for open water area --- ! 
     357            !     mask eddy diffusivity coefficient at ocean U- and V-points 
     358            DO jj = 1, jpjm1                    ! NB: has not to be defined on jpj line and jpi row 
    342359               DO ji = 1 , fs_jpim1   ! vector opt. 
    343                   pahu3D(ji,jj,jl) = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -a_i(ji  ,jj,  jl ) ) ) )   & 
    344                   &                * ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -a_i(ji+1,jj,  jl ) ) ) ) * ahiu(ji,jj) 
    345                   pahv3D(ji,jj,jl) = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -a_i(ji,  jj,  jl ) ) ) )   & 
    346                   &                * ( 1._wp - MAX( 0._wp, SIGN( 1._wp,- a_i(ji,  jj+1,jl ) ) ) ) * ahiv(ji,jj) 
    347                END DO 
    348             END DO 
    349  
    350             zhdfptab(:,:,jm)= a_i  (:,:,  jl); jm = jm + 1     
    351             zhdfptab(:,:,jm)= v_i  (:,:,  jl); jm = jm + 1 
    352             zhdfptab(:,:,jm)= v_s  (:,:,  jl); jm = jm + 1  
    353             zhdfptab(:,:,jm)= smv_i(:,:,  jl); jm = jm + 1 
    354             zhdfptab(:,:,jm)= oa_i (:,:,  jl); jm = jm + 1 
    355             zhdfptab(:,:,jm)= e_s  (:,:,1,jl); jm = jm + 1 
    356          ! Sample of adding more variables to apply lim_hdf using lim_hdf optimization--- 
    357          !   zhdfptab(:,:,jm) = variable_1 (:,:,1,jl); jm = jm + 1   
    358          !   zhdfptab(:,:,jm) = variable_2 (:,:,1,jl); jm = jm + 1  
    359          ! 
    360          ! and in this example the parameter ihdf_vars musb be changed to 8 (necessary for allocation) 
    361          !---------------------------------------------------------------------------------------- 
    362             DO jk = 1, nlay_i 
    363               zhdfptab(:,:,jm)=e_i(:,:,jk,jl); jm= jm+1 
    364             END DO 
    365          END DO 
    366          ! 
    367          !-------------------------------- 
    368          !  diffusion of open water area 
    369          !-------------------------------- 
    370          !                             ! Masked eddy diffusivity coefficient at ocean U- and V-points 
    371          !DO jj = 1, jpjm1                    ! NB: has not to be defined on jpj line and jpi row 
    372          !   DO ji = 1 , fs_jpim1   ! vector opt. 
    373          !      pahu(ji,jj) = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -at_i(ji  ,jj) ) ) )   & 
    374          !         &        * ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -at_i(ji+1,jj) ) ) ) * ahiu(ji,jj) 
    375          !      pahv(ji,jj) = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -at_i(ji,jj  ) ) ) )   & 
    376          !         &        * ( 1._wp - MAX( 0._wp, SIGN( 1._wp,- at_i(ji,jj+1) ) ) ) * ahiv(ji,jj) 
    377          !   END DO 
    378          !END DO 
     360                  pahu3D(ji,jj,jpl+1) = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -at_i(ji  ,jj) ) ) )   & 
     361                     &                * ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -at_i(ji+1,jj) ) ) ) * ahiu(ji,jj) 
     362                  pahv3D(ji,jj,jpl+1) = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -at_i(ji,jj  ) ) ) )   & 
     363                     &                * ( 1._wp - MAX( 0._wp, SIGN( 1._wp,- at_i(ji,jj+1) ) ) ) * ahiv(ji,jj) 
     364               END DO 
     365            END DO 
     366            ! 
     367            zhdfptab(:,:,jm)= ato_i  (:,:); 
     368 
     369            ! --- Apply diffusion --- ! 
     370            CALL lim_hdf( zhdfptab, ihdf_vars, jpl, nlay_i)  
     371             
     372            ! --- Recover properties --- ! 
     373            jm=1 
     374            DO jl = 1, jpl 
     375               a_i  (:,:,  jl) = zhdfptab(:,:,jm); jm = jm + 1       
     376               v_i  (:,:,  jl) = zhdfptab(:,:,jm); jm = jm + 1  
     377               v_s  (:,:,  jl) = zhdfptab(:,:,jm); jm = jm + 1  
     378               smv_i(:,:,  jl) = zhdfptab(:,:,jm); jm = jm + 1  
     379               oa_i (:,:,  jl) = zhdfptab(:,:,jm); jm = jm + 1  
     380               e_s  (:,:,1,jl) = zhdfptab(:,:,jm); jm = jm + 1  
     381               ! Sample of adding more variables to apply lim_hdf--------- 
     382               !   variable_1  (:,:,1,jl) = zhdfptab(:,:, jm  ) ; jm + 1  
     383               !   variable_2  (:,:,1,jl) = zhdfptab(:,:, jm  ) ; jm + 1 
     384               DO jk = 1, nlay_i 
     385                  e_i(:,:,jk,jl) = zhdfptab(:,:,jm);jm= jm + 1  
     386               END DO 
     387            END DO 
     388            ato_i  (:,:) = zhdfptab(:,:,jm) 
     389 
     390         ENDIF 
    379391          
    380          DO jj = 1, jpjm1                    ! NB: has not to be defined on jpj line and jpi row 
    381             DO ji = 1 , fs_jpim1   ! vector opt. 
    382                pahu3D(ji,jj,jpl+1) = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -at_i(ji  ,jj) ) ) )   & 
    383                   &                * ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -at_i(ji+1,jj) ) ) ) * ahiu(ji,jj) 
    384                pahv3D(ji,jj,jpl+1) = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -at_i(ji,jj  ) ) ) )   & 
    385                   &                * ( 1._wp - MAX( 0._wp, SIGN( 1._wp,- at_i(ji,jj+1) ) ) ) * ahiv(ji,jj) 
    386             END DO 
    387          END DO 
    388          ! 
    389          zhdfptab(:,:,jm)= ato_i  (:,:); 
    390          CALL lim_hdf( zhdfptab, ihdf_vars, jpl, nlay_i)  
    391  
    392          jm=1 
    393          DO jl = 1, jpl 
    394             a_i  (:,:,  jl) = zhdfptab(:,:,jm); jm = jm + 1       
    395             v_i  (:,:,  jl) = zhdfptab(:,:,jm); jm = jm + 1  
    396             v_s  (:,:,  jl) = zhdfptab(:,:,jm); jm = jm + 1  
    397             smv_i(:,:,  jl) = zhdfptab(:,:,jm); jm = jm + 1  
    398             oa_i (:,:,  jl) = zhdfptab(:,:,jm); jm = jm + 1  
    399             e_s  (:,:,1,jl) = zhdfptab(:,:,jm); jm = jm + 1  
    400          ! Sample of adding more variables to apply lim_hdf--------- 
    401          !   variable_1  (:,:,1,jl) = zhdfptab(:,:, jm  ) ; jm + 1  
    402          !   variable_2  (:,:,1,jl) = zhdfptab(:,:, jm  ) ; jm + 1 
    403          !----------------------------------------------------------- 
    404             DO jk = 1, nlay_i 
    405                e_i(:,:,jk,jl) = zhdfptab(:,:,jm);jm= jm + 1  
    406             END DO 
    407          END DO 
    408  
    409          ato_i  (:,:) = zhdfptab(:,:,jm) 
    410  
    411          !------------------------------------------------------------------------------! 
    412          ! limit ice properties after transport                            
    413          !------------------------------------------------------------------------------! 
    414 !!gm & cr   :  MAX should not be active if adv scheme is positive ! 
    415          DO jl = 1, jpl 
    416             DO jj = 1, jpj 
    417                DO ji = 1, jpi 
    418                   v_s  (ji,jj,jl)   = MAX( 0._wp, v_s  (ji,jj,jl) ) 
    419                   v_i  (ji,jj,jl)   = MAX( 0._wp, v_i  (ji,jj,jl) ) 
    420                   smv_i(ji,jj,jl)   = MAX( 0._wp, smv_i(ji,jj,jl) ) 
    421                   oa_i (ji,jj,jl)   = MAX( 0._wp, oa_i (ji,jj,jl) ) 
    422                   a_i  (ji,jj,jl)   = MAX( 0._wp, a_i  (ji,jj,jl) ) 
    423                   e_s  (ji,jj,1,jl) = MAX( 0._wp, e_s  (ji,jj,1,jl) ) 
    424                END DO 
    425             END DO 
    426  
    427             DO jk = 1, nlay_i 
    428                DO jj = 1, jpj 
    429                   DO ji = 1, jpi 
    430                      e_i(ji,jj,jk,jl) = MAX( 0._wp, e_i(ji,jj,jk,jl) ) 
    431                   END DO 
    432                END DO 
    433             END DO 
    434          END DO 
    435 !!gm & cr  
    436  
    437392         ! --- diags --- 
    438393         DO jj = 1, jpj 
Note: See TracChangeset for help on using the changeset viewer.