New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Changeset 4177 for branches/2013/dev_CMCC_2013/NEMOGCM – NEMO

Ignore:
Timestamp:
2013-11-11T12:15:42+01:00 (11 years ago)
Author:
vichi
Message:

ticket #1173 step 5: Add in changes from the trunk between revisions 3996 and 4119

Location:
branches/2013/dev_CMCC_2013/NEMOGCM
Files:
1 deleted
19 edited
1 copied

Legend:

Unmodified
Added
Removed
  • branches/2013/dev_CMCC_2013/NEMOGCM/NEMO/NST_SRC/agrif_opa_sponge.F90

    r3918 r4177  
    185185      INTEGER  :: ji,jj,jk 
    186186      INTEGER  :: ispongearea, ilci, ilcj 
    187       REAL(wp) :: z1spongearea 
    188       REAL(wp), POINTER, DIMENSION(:,:) :: zlocalviscsponge 
     187      LOGICAL  :: ll_spdone 
     188      REAL(wp) :: z1spongearea, zramp 
     189      REAL(wp), POINTER, DIMENSION(:,:) :: ztabramp 
    189190 
    190191#if defined SPONGE || defined SPONGE_TOP 
    191  
    192       CALL wrk_alloc( jpi, jpj, zlocalviscsponge ) 
    193  
    194       ispongearea  = 2 + 2 * Agrif_irhox() 
    195       ilci = nlci - ispongearea 
    196       ilcj = nlcj - ispongearea  
    197       z1spongearea = 1._wp / REAL( ispongearea - 2 ) 
    198       spbtr2(:,:) = 1. / ( e1t(:,:) * e2t(:,:) ) 
     192      ll_spdone=.TRUE. 
     193      IF (( .NOT. spongedoneT ).OR.( .NOT. spongedoneU )) THEN 
     194         ! Define ramp from boundaries towards domain interior 
     195         ! at T-points 
     196         ! Store it in ztabramp 
     197         ll_spdone=.FALSE. 
     198 
     199         CALL wrk_alloc( jpi, jpj, ztabramp ) 
     200 
     201         ispongearea  = 2 + 2 * Agrif_irhox() 
     202         ilci = nlci - ispongearea 
     203         ilcj = nlcj - ispongearea  
     204         z1spongearea = 1._wp / REAL( ispongearea - 2 ) 
     205         spbtr2(:,:) = 1. / ( e1t(:,:) * e2t(:,:) ) 
     206 
     207         ztabramp(:,:) = 0. 
     208 
     209         IF( (nbondi == -1) .OR. (nbondi == 2) ) THEN 
     210            DO jj = 1, jpj 
     211               IF ( umask(2,jj,1) == 1._wp ) THEN 
     212                 DO ji = 2, ispongearea                   
     213                    ztabramp(ji,jj) = ( ispongearea-ji ) * z1spongearea 
     214                 END DO 
     215               ENDIF 
     216            ENDDO 
     217         ENDIF 
     218 
     219         IF( (nbondi == 1) .OR. (nbondi == 2) ) THEN 
     220            DO jj = 1, jpj 
     221               IF ( umask(nlci-2,jj,1) == 1._wp ) THEN 
     222                  DO ji = ilci+1,nlci-1 
     223                     zramp = (ji - (ilci+1) ) * z1spongearea 
     224                     ztabramp(ji,jj) = MAX( ztabramp(ji,jj), zramp ) 
     225                  ENDDO 
     226               ENDIF 
     227            ENDDO 
     228         ENDIF 
     229 
     230         IF( (nbondj == -1) .OR. (nbondj == 2) ) THEN 
     231            DO ji = 1, jpi 
     232               IF ( vmask(ji,2,1) == 1._wp ) THEN 
     233                  DO jj = 2, ispongearea 
     234                     zramp = ( ispongearea-jj ) * z1spongearea 
     235                     ztabramp(ji,jj) = MAX( ztabramp(ji,jj), zramp ) 
     236                  END DO 
     237               ENDIF 
     238            ENDDO 
     239         ENDIF 
     240 
     241         IF( (nbondj == 1) .OR. (nbondj == 2) ) THEN 
     242            DO ji = 1, jpi 
     243               IF ( vmask(ji,nlcj-2,1) == 1._wp ) THEN 
     244                  DO jj = ilcj+1,nlcj-1 
     245                     zramp = (jj - (ilcj+1) ) * z1spongearea 
     246                     ztabramp(ji,jj) = MAX( ztabramp(ji,jj), zramp ) 
     247                  END DO 
     248               ENDIF 
     249            ENDDO 
     250         ENDIF 
     251 
     252      ENDIF 
    199253 
    200254      ! Tracers 
    201255      IF( .NOT. spongedoneT ) THEN 
    202          zlocalviscsponge(:,:) = 0. 
    203256         spe1ur(:,:) = 0. 
    204257         spe2vr(:,:) = 0. 
    205258 
    206259         IF( (nbondi == -1) .OR. (nbondi == 2) ) THEN 
    207             DO ji = 2, ispongearea 
    208                zlocalviscsponge(ji,:) = visc_tra * ( ispongearea-ji ) * z1spongearea 
    209             ENDDO 
    210             spe1ur(2:ispongearea-1,:      ) = 0.5 * ( zlocalviscsponge(2:ispongearea-1,:      )   & 
    211                &                         +            zlocalviscsponge(3:ispongearea  ,:      ) ) & 
    212                &                         * e2u(2:ispongearea-1,:      ) / e1u(2:ispongearea-1,:      ) 
    213             spe2vr(2:ispongearea  ,1:jpjm1) = 0.5 * ( zlocalviscsponge(2:ispongearea  ,1:jpjm1)   & 
    214                &                         +            zlocalviscsponge(2:ispongearea,2  :jpj  ) ) & 
    215                &                         * e1v(2:ispongearea  ,1:jpjm1) / e2v(2:ispongearea  ,1:jpjm1) 
     260            spe1ur(2:ispongearea-1,:       ) = visc_tra                                        & 
     261               &                             *    0.5 * (  ztabramp(2:ispongearea-1,:      )   & 
     262               &                                         + ztabramp(3:ispongearea  ,:      ) ) & 
     263               &                             * e2u(2:ispongearea-1,:) / e1u(2:ispongearea-1,:) 
     264 
     265            spe2vr(2:ispongearea  ,1:jpjm1 ) = visc_tra                                        & 
     266               &                             *    0.5 * (  ztabramp(2:ispongearea  ,1:jpjm1)   & 
     267               &                                         + ztabramp(2:ispongearea,2  :jpj  ) ) & 
     268               &                             * e1v(2:ispongearea,1:jpjm1) / e2v(2:ispongearea,1:jpjm1) 
    216269         ENDIF 
    217270 
    218271         IF( (nbondi == 1) .OR. (nbondi == 2) ) THEN 
    219             DO ji = ilci+1,nlci-1 
    220                zlocalviscsponge(ji,:) = visc_tra * (ji - (ilci+1) ) * z1spongearea 
    221             ENDDO 
    222    
    223             spe1ur(ilci+1:nlci-2,:      ) = 0.5 * (  zlocalviscsponge(ilci+1:nlci-2,:)    &  
    224                &                          +          zlocalviscsponge(ilci+2:nlci-1,:) )  & 
    225                &                          * e2u(ilci+1:nlci-2,:) / e1u(ilci+1:nlci-2,:) 
    226  
    227             spe2vr(ilci+1:nlci-1,1:jpjm1) = 0.5 * (  zlocalviscsponge(ilci+1:nlci-1,1:jpjm1)    &  
    228                &                            +        zlocalviscsponge(ilci+1:nlci-1,2:jpj  )  ) &  
    229                &                                   * e1v(ilci+1:nlci-1,1:jpjm1) / e2v(ilci+1:nlci-1,1:jpjm1) 
     272            spe1ur(ilci+1:nlci-2,:        ) = visc_tra                                   & 
     273               &                            * 0.5 * (  ztabramp(ilci+1:nlci-2,:      )   &  
     274               &                                     + ztabramp(ilci+2:nlci-1,:      ) ) & 
     275               &                            * e2u(ilci+1:nlci-2,:) / e1u(ilci+1:nlci-2,:) 
     276 
     277            spe2vr(ilci+1:nlci-1,1:jpjm1  )  = visc_tra                                  & 
     278               &                            * 0.5 * (  ztabramp(ilci+1:nlci-1,1:jpjm1)   &  
     279               &                                     + ztabramp(ilci+1:nlci-1,2:jpj  ) ) &  
     280               &                            * e1v(ilci+1:nlci-1,1:jpjm1) / e2v(ilci+1:nlci-1,1:jpjm1) 
    230281         ENDIF 
    231282 
    232283         IF( (nbondj == -1) .OR. (nbondj == 2) ) THEN 
    233             DO jj = 2, ispongearea 
    234                zlocalviscsponge(:,jj) = visc_tra * ( ispongearea-jj ) * z1spongearea 
    235             ENDDO 
    236             spe1ur(1:jpim1,2:ispongearea  ) = 0.5 * ( zlocalviscsponge(1:jpim1,2:ispongearea  ) &  
    237                &                            +         zlocalviscsponge(2:jpi  ,2:ispongearea) ) & 
     284            spe1ur(1:jpim1,2:ispongearea  ) = visc_tra                                     & 
     285               &                            * 0.5 * (  ztabramp(1:jpim1,2:ispongearea  )   &  
     286               &                                     + ztabramp(2:jpi  ,2:ispongearea  ) ) & 
    238287               &                            * e2u(1:jpim1,2:ispongearea) / e1u(1:jpim1,2:ispongearea) 
    239288    
    240             spe2vr(:      ,2:ispongearea-1) = 0.5 * ( zlocalviscsponge(:,2:ispongearea-1)       & 
    241                &                            +         zlocalviscsponge(:,3:ispongearea  )     ) & 
     289            spe2vr(:      ,2:ispongearea-1) = visc_tra                                     & 
     290               &                            * 0.5 * (  ztabramp(:      ,2:ispongearea-1)   & 
     291               &                                     + ztabramp(:      ,3:ispongearea  ) ) & 
    242292               &                            * e1v(:,2:ispongearea-1) / e2v(:,2:ispongearea-1) 
    243293         ENDIF 
    244294 
    245295         IF( (nbondj == 1) .OR. (nbondj == 2) ) THEN 
    246             DO jj = ilcj+1,nlcj-1 
    247                zlocalviscsponge(:,jj) = visc_tra * (jj - (ilcj+1) ) * z1spongearea 
    248             ENDDO 
    249             spe1ur(1:jpim1,ilcj+1:nlcj-1) = 0.5 * ( zlocalviscsponge(1:jpim1,ilcj+1:nlcj-1)   & 
    250                &                          +         zlocalviscsponge(2:jpi  ,ilcj+1:nlcj-1) ) & 
     296            spe1ur(1:jpim1,ilcj+1:nlcj-1) = visc_tra                                   & 
     297               &                          * 0.5 * (  ztabramp(1:jpim1,ilcj+1:nlcj-1)   & 
     298               &                                   + ztabramp(2:jpi  ,ilcj+1:nlcj-1) ) & 
    251299               &                                * e2u(1:jpim1,ilcj+1:nlcj-1) / e1u(1:jpim1,ilcj+1:nlcj-1) 
    252             spe2vr(:      ,ilcj+1:nlcj-2) = 0.5 * ( zlocalviscsponge(:,ilcj+1:nlcj-2      )   & 
    253                &                          +         zlocalviscsponge(:,ilcj+2:nlcj-1)     )   & 
     300 
     301            spe2vr(:      ,ilcj+1:nlcj-2) = visc_tra                                   & 
     302               &                          * 0.5 * (  ztabramp(:      ,ilcj+1:nlcj-2)   & 
     303               &                                   + ztabramp(:      ,ilcj+2:nlcj-1) ) & 
    254304               &                                * e1v(:,ilcj+1:nlcj-2) / e2v(:,ilcj+1:nlcj-2) 
    255305         ENDIF 
     
    259309      ! Dynamics 
    260310      IF( .NOT. spongedoneU ) THEN 
    261          zlocalviscsponge(:,:) = 0. 
    262311         spe1ur2(:,:) = 0. 
    263312         spe2vr2(:,:) = 0. 
    264313 
    265314         IF( (nbondi == -1) .OR. (nbondi == 2) ) THEN 
    266             DO ji = 2, ispongearea 
    267                zlocalviscsponge(ji,:) = visc_dyn * ( ispongearea-ji ) * z1spongearea 
    268             ENDDO 
    269             spe1ur2(2:ispongearea-1,:      ) = 0.5 * ( zlocalviscsponge(2:ispongearea-1,:      ) & 
    270                                              &     +   zlocalviscsponge(3:ispongearea,:    ) ) 
    271             spe2vr2(2:ispongearea  ,1:jpjm1) = 0.5 * ( zlocalviscsponge(2:ispongearea  ,1:jpjm1) & 
    272                                              &     +   zlocalviscsponge(2:ispongearea,2:jpj) )  
     315            spe1ur2(2:ispongearea-1,:      ) = visc_dyn                                   & 
     316               &                             * 0.5 * (  ztabramp(2:ispongearea-1,:      ) & 
     317               &                                      + ztabramp(3:ispongearea  ,:      ) ) 
     318            spe2vr2(2:ispongearea  ,1:jpjm1) = visc_dyn                                   & 
     319               &                             * 0.5 * (  ztabramp(2:ispongearea  ,1:jpjm1) & 
     320               &                                      + ztabramp(2:ispongearea  ,2:jpj  ) )  
    273321         ENDIF 
    274322 
    275323         IF( (nbondi == 1) .OR. (nbondi == 2) ) THEN 
    276             DO ji = ilci+1,nlci-1 
    277                zlocalviscsponge(ji,:) = visc_dyn * (ji - (ilci+1) ) * z1spongearea 
    278             ENDDO 
    279             spe1ur2(ilci+1:nlci-2,:      ) = 0.5 * (  zlocalviscsponge(ilci+1:nlci-2,:) & 
    280                                            &        + zlocalviscsponge(ilci+2:nlci-1,:) )   
    281             spe2vr2(ilci+1:nlci-1,1:jpjm1) = 0.5 * (  zlocalviscsponge(ilci+1:nlci-1,1:jpjm1) & 
    282                                            &        + zlocalviscsponge(ilci+1:nlci-1,2:jpj  )  )  
     324            spe1ur2(ilci+1:nlci-2  ,:      ) = visc_dyn                                   & 
     325               &                             * 0.5 * (  ztabramp(ilci+1:nlci-2, :       ) & 
     326               &                                      + ztabramp(ilci+2:nlci-1, :       ) )                       
     327            spe2vr2(ilci+1:nlci-1  ,1:jpjm1) = visc_dyn                                   & 
     328               &                             * 0.5 * (  ztabramp(ilci+1:nlci-1,1:jpjm1  ) & 
     329               &                                      + ztabramp(ilci+1:nlci-1,2:jpj    ) )  
    283330         ENDIF 
    284331 
    285332         IF( (nbondj == -1) .OR. (nbondj == 2) ) THEN 
    286             DO jj = 2, ispongearea 
    287                zlocalviscsponge(:,jj) = visc_dyn * ( ispongearea-jj ) * z1spongearea 
    288             ENDDO 
    289             spe1ur2(1:jpim1,2:ispongearea  ) = 0.5 * ( zlocalviscsponge(1:jpim1,2:ispongearea) & 
    290                                              &      + zlocalviscsponge(2:jpi,2:ispongearea) )  
    291             spe2vr2(:      ,2:ispongearea-1) = 0.5 * ( zlocalviscsponge(:,2:ispongearea-1)     & 
    292                                              &      + zlocalviscsponge(:,3:ispongearea)     ) 
     333            spe1ur2(1:jpim1,2:ispongearea  ) = visc_dyn                                   &   
     334               &                             * 0.5 * (  ztabramp(1:jpim1,2:ispongearea  ) & 
     335               &                                      + ztabramp(2:jpi  ,2:ispongearea  ) )  
     336            spe2vr2(:      ,2:ispongearea-1) = visc_dyn                                   & 
     337               &                             * 0.5 * (  ztabramp(:      ,2:ispongearea-1) & 
     338               &                                      + ztabramp(:      ,3:ispongearea  ) ) 
    293339         ENDIF 
    294340 
    295341         IF( (nbondj == 1) .OR. (nbondj == 2) ) THEN 
    296             DO jj = ilcj+1,nlcj-1 
    297                zlocalviscsponge(:,jj) = visc_dyn * (jj - (ilcj+1) ) * z1spongearea 
    298             ENDDO 
    299             spe1ur2(1:jpim1,ilcj+1:nlcj-1) = 0.5 * ( zlocalviscsponge(1:jpim1,ilcj+1:nlcj-1) & 
    300                                            &         + zlocalviscsponge(2:jpi,ilcj+1:nlcj-1) )  
    301             spe2vr2(:      ,ilcj+1:nlcj-2) = 0.5 * ( zlocalviscsponge(:,ilcj+1:nlcj-2      ) & 
    302                                            &         + zlocalviscsponge(:,ilcj+2:nlcj-1)     ) 
     342            spe1ur2(1:jpim1,ilcj+1:nlcj-1  ) = visc_dyn                                   & 
     343               &                             * 0.5 * (  ztabramp(1:jpim1,ilcj+1:nlcj-1  ) & 
     344               &                                      + ztabramp(2:jpi  ,ilcj+1:nlcj-1  ) )  
     345            spe2vr2(:      ,ilcj+1:nlcj-2  ) = visc_dyn                                   & 
     346               &                             * 0.5 * (  ztabramp(:      ,ilcj+1:nlcj-2  ) & 
     347               &                                      + ztabramp(:      ,ilcj+2:nlcj-1  ) ) 
    303348         ENDIF 
    304349         spongedoneU = .TRUE. 
     
    306351      ENDIF 
    307352      ! 
    308       CALL wrk_dealloc( jpi, jpj, zlocalviscsponge ) 
     353      IF (.NOT.ll_spdone) CALL wrk_dealloc( jpi, jpj, ztabramp ) 
    309354      ! 
    310355#endif 
  • branches/2013/dev_CMCC_2013/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90

    r3785 r4177  
    682682      ! used to prevent the applied increments taking the temperature below the local freezing point  
    683683 
    684 #if defined key_cice  
    685         fzptnz(:,:,:) = -1.8_wp 
    686 #else  
    687         DO jk = 1, jpk 
    688            DO jj = 1, jpj 
    689               DO ji = 1, jpk 
    690                  fzptnz (ji,jj,jk) = ( -0.0575_wp + 1.710523e-3_wp * SQRT( tsn(ji,jj,jk,jp_sal) )                   &  
    691                                                   - 2.154996e-4_wp *       tsn(ji,jj,jk,jp_sal)   ) * tsn(ji,jj,jk,jp_sal)  &  
    692                                                   - 7.53e-4_wp * fsdepw(ji,jj,jk)       ! (pressure in dbar)  
    693               END DO 
    694            END DO 
    695         END DO 
    696 #endif  
     684      DO jk=1, jpkm1 
     685         fzptnz (:,:,jk) = tfreez( tsn(:,:,jk,jp_sal), fsdept(:,:,jk) ) 
     686      ENDDO 
    697687 
    698688      IF ( ln_asmiau ) THEN 
  • branches/2013/dev_CMCC_2013/NEMOGCM/NEMO/OPA_SRC/DIA/diahsb.F90

    r3625 r4177  
    2121   USE bdy_par         ! (for lk_bdy) 
    2222   USE timing          ! preformance summary 
     23   USE lib_fortran 
     24   USE sbcrnf 
    2325 
    2426   IMPLICIT NONE 
     
    3335   REAL(dp)                                ::   surf_tot   , vol_tot             ! 
    3436   REAL(dp)                                ::   frc_t      , frc_s     , frc_v   ! global forcing trends 
     37   REAL(dp)                                ::   frc_wn_t      , frc_wn_s ! global forcing trends 
    3538   REAL(dp)                                ::   fact1                            ! conversion factors 
    3639   REAL(dp)                                ::   fact21    , fact22               !     -         - 
     
    3841   REAL(dp), DIMENSION(:,:)  , ALLOCATABLE ::   surf      , ssh_ini              ! 
    3942   REAL(dp), DIMENSION(:,:,:), ALLOCATABLE ::   hc_loc_ini, sc_loc_ini, e3t_ini  ! 
     43   REAL(dp), DIMENSION(:,:)  , ALLOCATABLE ::   ssh_hc_loc_ini, ssh_sc_loc_ini 
    4044 
    4145   !! * Substitutions 
     
    6771      INTEGER    ::   jk                          ! dummy loop indice 
    6872      REAL(dp)   ::   zdiff_hc    , zdiff_sc      ! heat and salt content variations 
     73      REAL(dp)   ::   zdiff_hc1   , zdiff_sc1     ! heat and salt content variations of ssh 
    6974      REAL(dp)   ::   zdiff_v1    , zdiff_v2      ! volume variation 
     75      REAL(dp)   ::   zerr_hc1    , zerr_sc1      ! Non conservation due to free surface 
    7076      REAL(dp)   ::   z1_rau0                     ! local scalars 
    7177      REAL(dp)   ::   zdeltat                     !    -     - 
    7278      REAL(dp)   ::   z_frc_trd_t , z_frc_trd_s   !    -     - 
    7379      REAL(dp)   ::   z_frc_trd_v                 !    -     - 
     80      REAL(dp)   ::   z_wn_trd_t , z_wn_trd_s   !    -     - 
     81      REAL(dp)   ::   z_ssh_hc , z_ssh_sc   !    -     - 
    7482      !!--------------------------------------------------------------------------- 
    7583      IF( nn_timing == 1 )   CALL timing_start('dia_hsb') 
     
    7987      ! ------------------------- ! 
    8088      z1_rau0 = 1.e0 / rau0 
    81       z_frc_trd_v = z1_rau0 * SUM( - ( emp(:,:) - rnf(:,:) ) * surf(:,:) )     ! volume fluxes 
    82       z_frc_trd_t =           SUM( sbc_tsc(:,:,jp_tem) * surf(:,:) )     ! heat fluxes 
    83       z_frc_trd_s =           SUM( sbc_tsc(:,:,jp_sal) * surf(:,:) )     ! salt fluxes 
     89      z_frc_trd_v = z1_rau0 * glob_sum( - ( emp(:,:) - rnf(:,:) ) * surf(:,:) )     ! volume fluxes 
     90      z_frc_trd_t =           glob_sum( sbc_tsc(:,:,jp_tem) * surf(:,:) )     ! heat fluxes 
     91      z_frc_trd_s =           glob_sum( sbc_tsc(:,:,jp_sal) * surf(:,:) )     ! salt fluxes 
     92      ! Add runoff heat & salt input 
     93      IF( ln_rnf    )   z_frc_trd_t = z_frc_trd_t + glob_sum( rnf_tsc(:,:,jp_tem) * surf(:,:) ) 
     94      IF( ln_rnf_sal)   z_frc_trd_s = z_frc_trd_s + glob_sum( rnf_tsc(:,:,jp_sal) * surf(:,:) ) 
    8495      ! Add penetrative solar radiation 
    85       IF( ln_traqsr )   z_frc_trd_t = z_frc_trd_t + r1_rau0_rcp * SUM( qsr     (:,:) * surf(:,:) ) 
     96      IF( ln_traqsr )   z_frc_trd_t = z_frc_trd_t + r1_rau0_rcp * glob_sum( qsr     (:,:) * surf(:,:) ) 
    8697      ! Add geothermal heat flux 
    87       IF( ln_trabbc )   z_frc_trd_t = z_frc_trd_t + r1_rau0_rcp * SUM( qgh_trd0(:,:) * surf(:,:) ) 
    88       IF( lk_mpp ) THEN 
    89          CALL mpp_sum( z_frc_trd_v ) 
    90          CALL mpp_sum( z_frc_trd_t ) 
    91       ENDIF 
     98      IF( ln_trabbc )   z_frc_trd_t = z_frc_trd_t +  glob_sum( qgh_trd0(:,:) * surf(:,:) ) 
     99      IF( .NOT. lk_vvl ) THEN 
     100         z_wn_trd_t = - glob_sum( surf(:,:) * wn(:,:,1) * tsb(:,:,1,jp_tem) ) 
     101         z_wn_trd_s = - glob_sum( surf(:,:) * wn(:,:,1) * tsb(:,:,1,jp_sal) ) 
     102      ENDIF 
     103 
    92104      frc_v = frc_v + z_frc_trd_v * rdt 
    93105      frc_t = frc_t + z_frc_trd_t * rdt 
    94106      frc_s = frc_s + z_frc_trd_s * rdt 
     107      !                                          ! Advection flux through fixed surface (z=0) 
     108      IF( .NOT. lk_vvl ) THEN 
     109         frc_wn_t = frc_wn_t + z_wn_trd_t * rdt 
     110         frc_wn_s = frc_wn_s + z_wn_trd_s * rdt 
     111      ENDIF 
    95112 
    96113      ! ----------------------- ! 
     
    100117      zdiff_hc = 0.d0 
    101118      zdiff_sc = 0.d0 
     119 
    102120      ! volume variation (calculated with ssh) 
    103       zdiff_v1 = SUM( surf(:,:) * tmask(:,:,1) * ( sshn(:,:) - ssh_ini(:,:) ) ) 
     121      zdiff_v1 = glob_sum( surf(:,:) * ( sshn(:,:) - ssh_ini(:,:) ) ) 
     122 
     123      ! heat & salt content variation (associated with ssh) 
     124      IF( .NOT. lk_vvl ) THEN 
     125         z_ssh_hc = glob_sum( surf(:,:) * ( tsn(:,:,1,jp_tem) * sshn(:,:) - ssh_hc_loc_ini(:,:) ) ) 
     126         z_ssh_sc = glob_sum( surf(:,:) * ( tsn(:,:,1,jp_sal) * sshn(:,:) - ssh_sc_loc_ini(:,:) ) ) 
     127      ENDIF 
     128 
    104129      DO jk = 1, jpkm1 
    105          ! volume variation (calculated with scale factors) 
    106          zdiff_v2 = zdiff_v2 + SUM( surf(:,:) * tmask(:,:,jk)   & 
     130        ! volume variation (calculated with scale factors) 
     131         zdiff_v2 = zdiff_v2 + glob_sum( surf(:,:) * tmask(:,:,jk)   & 
    107132            &                       * ( fse3t_n(:,:,jk)         & 
    108133            &                           - e3t_ini(:,:,jk) ) ) 
    109134         ! heat content variation 
    110          zdiff_hc = zdiff_hc + SUM( surf(:,:) * tmask(:,:,jk)          & 
     135         zdiff_hc = zdiff_hc + glob_sum( surf(:,:) * tmask(:,:,jk)          & 
    111136            &                       * ( fse3t_n(:,:,jk) * tsn(:,:,jk,jp_tem)   & 
    112137            &                           - hc_loc_ini(:,:,jk) ) ) 
    113138         ! salt content variation 
    114          zdiff_sc = zdiff_sc + SUM( surf(:,:) * tmask(:,:,jk)          & 
     139         zdiff_sc = zdiff_sc + glob_sum( surf(:,:) * tmask(:,:,jk)          & 
    115140            &                       * ( fse3t_n(:,:,jk) * tsn(:,:,jk,jp_sal)   & 
    116141            &                           - sc_loc_ini(:,:,jk) ) ) 
    117142      ENDDO 
    118143 
    119       IF( lk_mpp ) THEN 
    120          CALL mpp_sum( zdiff_hc ) 
    121          CALL mpp_sum( zdiff_sc ) 
    122          CALL mpp_sum( zdiff_v1 ) 
    123          CALL mpp_sum( zdiff_v2 ) 
    124       ENDIF 
    125  
    126144      ! Substract forcing from heat content, salt content and volume variations 
    127145      zdiff_v1 = zdiff_v1 - frc_v 
    128       zdiff_v2 = zdiff_v2 - frc_v 
     146      IF( lk_vvl )   zdiff_v2 = zdiff_v2 - frc_v 
    129147      zdiff_hc = zdiff_hc - frc_t 
    130148      zdiff_sc = zdiff_sc - frc_s 
     149      IF( .NOT. lk_vvl ) THEN 
     150         zdiff_hc1 = zdiff_hc + z_ssh_hc  
     151         zdiff_sc1 = zdiff_sc + z_ssh_sc 
     152         zerr_hc1  = z_ssh_hc - frc_wn_t 
     153         zerr_sc1  = z_ssh_sc - frc_wn_s 
     154      ENDIF 
    131155       
    132156      ! ----------------------- ! 
     
    134158      ! ----------------------- ! 
    135159      zdeltat  = 1.e0 / ( ( kt - nit000 + 1 ) * rdt ) 
    136       WRITE(numhsb , 9020) kt , zdiff_hc / vol_tot , zdiff_hc * fact1  * zdeltat,                                & 
    137          &                      zdiff_sc / vol_tot , zdiff_sc * fact21 * zdeltat, zdiff_sc * fact22 * zdeltat,   & 
    138          &                      zdiff_v1           , zdiff_v1 * fact31 * zdeltat, zdiff_v1 * fact32 * zdeltat,   & 
    139          &                      zdiff_v2           , zdiff_v2 * fact31 * zdeltat, zdiff_v2 * fact32 * zdeltat 
     160      IF( lk_vvl ) THEN 
     161         WRITE(numhsb , 9020) kt , zdiff_hc / vol_tot , zdiff_hc * fact1  * zdeltat,                                & 
     162            &                      zdiff_sc / vol_tot , zdiff_sc * fact21 * zdeltat, zdiff_sc * fact22 * zdeltat,   & 
     163            &                      zdiff_v1           , zdiff_v1 * fact31 * zdeltat, zdiff_v1 * fact32 * zdeltat,   & 
     164            &                      zdiff_v2           , zdiff_v2 * fact31 * zdeltat, zdiff_v2 * fact32 * zdeltat 
     165      ELSE 
     166         WRITE(numhsb , 9030) kt , zdiff_hc1 / vol_tot , zdiff_hc1 * fact1  * zdeltat,                                & 
     167            &                      zdiff_sc1 / vol_tot , zdiff_sc1 * fact21 * zdeltat, zdiff_sc1 * fact22 * zdeltat,   & 
     168            &                      zdiff_v1            , zdiff_v1  * fact31 * zdeltat, zdiff_v1  * fact32 * zdeltat,   & 
     169            &                      zerr_hc1 / vol_tot  , zerr_sc1 / vol_tot 
     170      ENDIF 
    140171 
    141172      IF ( kt == nitend ) CLOSE( numhsb ) 
     
    144175 
    1451769020  FORMAT(I5,11D15.7) 
     1779030  FORMAT(I5,10D15.7) 
    146178      ! 
    147179   END SUBROUTINE dia_hsb 
     
    179211 
    180212      IF( .NOT. ln_diahsb )   RETURN 
     213      IF( .NOT. lk_mpp_rep ) & 
     214        CALL ctl_stop (' Your global mpp_sum if performed in single precision - 64 bits -', & 
     215             &         ' whereas the global sum to be precise must be done in double precision ',& 
     216             &         ' please add key_mpp_rep') 
    181217 
    182218      ! ------------------- ! 
    183219      ! 1 - Allocate memory ! 
    184220      ! ------------------- ! 
    185       ALLOCATE( hc_loc_ini(jpi,jpj,jpk), STAT=ierror ) 
     221      ALLOCATE( hc_loc_ini(jpi,jpj,jpk), sc_loc_ini(jpi,jpj,jpk), & 
     222         &      ssh_hc_loc_ini(jpi,jpj), ssh_sc_loc_ini(jpi,jpj), & 
     223         &      e3t_ini(jpi,jpj,jpk)                            , & 
     224         &      surf(jpi,jpj),  ssh_ini(jpi,jpj), STAT=ierror ) 
    186225      IF( ierror > 0 ) THEN 
    187226         CALL ctl_stop( 'dia_hsb: unable to allocate hc_loc_ini' )   ;   RETURN 
    188       ENDIF 
    189       ALLOCATE( sc_loc_ini(jpi,jpj,jpk), STAT=ierror ) 
    190       IF( ierror > 0 ) THEN 
    191          CALL ctl_stop( 'dia_hsb: unable to allocate sc_loc_ini' )   ;   RETURN 
    192       ENDIF 
    193       ALLOCATE( e3t_ini(jpi,jpj,jpk)   , STAT=ierror ) 
    194       IF( ierror > 0 ) THEN 
    195          CALL ctl_stop( 'dia_hsb: unable to allocate e3t_ini' )      ;   RETURN 
    196       ENDIF 
    197       ALLOCATE( surf(jpi,jpj)          , STAT=ierror ) 
    198       IF( ierror > 0 ) THEN 
    199          CALL ctl_stop( 'dia_hsb: unable to allocate surf' )         ;   RETURN 
    200       ENDIF 
    201       ALLOCATE( ssh_ini(jpi,jpj)       , STAT=ierror ) 
    202       IF( ierror > 0 ) THEN 
    203          CALL ctl_stop( 'dia_hsb: unable to allocate ssh_ini' )      ;   RETURN 
    204227      ENDIF 
    205228 
     
    214237      cl_name    = 'heat_salt_volume_budgets.txt'                         ! name of output file 
    215238      surf(:,:) = e1t(:,:) * e2t(:,:) * tmask(:,:,1) * tmask_i(:,:)      ! masked surface grid cell area 
    216       surf_tot  = SUM( surf(:,:) )                                       ! total ocean surface area 
     239      surf_tot  = glob_sum( surf(:,:) )                                       ! total ocean surface area 
    217240      vol_tot   = 0.d0                                                   ! total ocean volume 
    218241      DO jk = 1, jpkm1 
    219          vol_tot  = vol_tot + SUM( surf(:,:) * tmask(:,:,jk)     & 
    220             &                      * fse3t_n(:,:,jk)         ) 
     242         vol_tot  = vol_tot + glob_sum( surf(:,:) * tmask(:,:,jk)     & 
     243            &                         * fse3t_n(:,:,jk)         ) 
    221244      END DO 
    222       IF( lk_mpp ) THEN  
    223          CALL mpp_sum( vol_tot ) 
    224          CALL mpp_sum( surf_tot ) 
    225       ENDIF 
    226245 
    227246      CALL ctl_opn( numhsb , cl_name , 'UNKNOWN' , 'FORMATTED' , 'SEQUENTIAL' , 1 , numout , lwp , 1 ) 
    228       !                   12345678901234567890123456789012345678901234567890123456789012345678901234567890 -> 80 
    229       WRITE( numhsb, 9010 ) "kt   |     heat content budget     |            salt content budget             ",   & 
    230          !                                                   123456789012345678901234567890123456789012345 -> 45 
    231          &                                                  "|            volume budget (ssh)             ",   & 
    232          !                                                   678901234567890123456789012345678901234567890 -> 45 
    233          &                                                  "|            volume budget (e3t)             " 
    234       WRITE( numhsb, 9010 ) "     |      [C]         [W/m2]     |     [psu]        [mmm/s]          [SV]     ",   & 
    235          &                                                  "|     [m3]         [mmm/s]          [SV]     ",   & 
    236          &                                                  "|     [m3]         [mmm/s]          [SV]     " 
    237  
     247      IF( lk_vvl ) THEN 
     248         !                   12345678901234567890123456789012345678901234567890123456789012345678901234567890 -> 80 
     249         WRITE( numhsb, 9010 ) "kt   |     heat content budget     |            salt content budget             ",   & 
     250            !                                                   123456789012345678901234567890123456789012345 -> 45 
     251            &                                                  "|            volume budget (ssh)             ",   & 
     252            !                                                   678901234567890123456789012345678901234567890 -> 45 
     253            &                                                  "|            volume budget (e3t)             " 
     254         WRITE( numhsb, 9010 ) "     |      [C]         [W/m2]     |     [psu]        [mmm/s]          [SV]     ",   & 
     255            &                                                  "|     [m3]         [mmm/s]          [SV]     ",   & 
     256            &                                                  "|     [m3]         [mmm/s]          [SV]     " 
     257      ELSE 
     258         !                   12345678901234567890123456789012345678901234567890123456789012345678901234567890 -> 80 
     259         WRITE( numhsb, 9011 ) "kt   |     heat content budget     |            salt content budget             ",   & 
     260            !                                                   123456789012345678901234567890123456789012345 -> 45 
     261            &                                                  "|            volume budget (ssh)             ",   & 
     262            !                                                   678901234567890123456789012345678901234567890 -> 45 
     263            &                                                  "|  Non conservation due to free surface      " 
     264         WRITE( numhsb, 9011 ) "     |      [C]         [W/m2]     |     [psu]        [mmm/s]          [SV]     ",   & 
     265            &                                                  "|     [m3]         [mmm/s]          [SV]     ",   & 
     266            &                                                  "|  [heat - C]     [salt - psu]                " 
     267      ENDIF 
    238268      ! --------------- ! 
    239269      ! 3 - Conversions ! (factors will be multiplied by duration afterwards) 
     
    261291      frc_t = 0.d0                                           ! heat content   -    -   -    -    
    262292      frc_s = 0.d0                                           ! salt content   -    -   -    -          
     293      IF( .NOT. lk_vvl ) THEN 
     294         ssh_hc_loc_ini(:,:) = tsn(:,:,1,jp_tem) * ssh_ini(:,:)   ! initial heat content associated with ssh 
     295         ssh_sc_loc_ini(:,:) = tsn(:,:,1,jp_sal) * ssh_ini(:,:)   ! initial salt content associated with ssh 
     296         frc_wn_t = 0.d0 
     297         frc_wn_s = 0.d0 
     298      ENDIF 
    263299      ! 
    2643009010  FORMAT(A80,A45,A45) 
     3019011  FORMAT(A80,A45,A45) 
    265302      ! 
    266303   END SUBROUTINE dia_hsb_init 
  • branches/2013/dev_CMCC_2013/NEMOGCM/NEMO/OPA_SRC/DOM/closea.F90

    r3632 r4177  
    108108            ncsi1(2)   =  97  ;  ncsj1(2)   = 107 
    109109            ncsi2(2)   = 103  ;  ncsj2(2)   = 111 
    110             ncsir(2,1) = 110  ;  ncsjr(2,1) = 111 
    111             !                                            ! Black Sea 1 : west part of the Black Sea  
    112             ncsnr(3)   = 1    ; ncstt(3)   =   2            !            (ie west of the cyclic b.c.) 
    113             ncsi1(3)   = 174  ; ncsj1(3)   = 107            ! put in Med Sea 
    114             ncsi2(3)   = 181  ; ncsj2(3)   = 112 
    115             ncsir(3,1) = 171  ; ncsjr(3,1) = 106  
    116             !                                            ! Black Sea 2 : est part of the Black Sea  
    117             ncsnr(4)   =   1  ;  ncstt(4)   =   2           !               (ie est of the cyclic b.c.) 
    118             ncsi1(4)   =   2  ;  ncsj1(4)   = 107           ! put in Med Sea 
    119             ncsi2(4)   =   6  ;  ncsj2(4)   = 112 
    120             ncsir(4,1) = 171  ;  ncsjr(4,1) = 106  
     110            ncsir(2,1) = 110  ;  ncsjr(2,1) = 111            
     111            !                                            ! Black Sea (crossed by the cyclic boundary condition) 
     112            ncsnr(3:4) =   4  ;  ncstt(3:4) =   2           ! put in Med Sea (north of Aegean Sea) 
     113            ncsir(3:4,1) = 171;  ncsjr(3:4,1) = 106         ! 
     114            ncsir(3:4,2) = 170;  ncsjr(3:4,2) = 106  
     115            ncsir(3:4,3) = 171;  ncsjr(3:4,3) = 105  
     116            ncsir(3:4,4) = 170;  ncsjr(3:4,4) = 105  
     117            ncsi1(3)   = 174  ;  ncsj1(3)   = 107           ! 1 : west part of the Black Sea       
     118            ncsi2(3)   = 181  ;  ncsj2(3)   = 112           !            (ie west of the cyclic b.c.) 
     119            ncsi1(4)   =   2  ;  ncsj1(4)   = 107           ! 2 : east part of the Black Sea  
     120            ncsi2(4)   =   6  ;  ncsj2(4)   = 112           !           (ie east of the cyclic b.c.) 
     121              
     122           
     123 
    121124            !                                        ! ======================= 
    122125         CASE ( 4 )                                  !  ORCA_R4 configuration 
     
    372375      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   p_rnfmsk   ! river runoff mask (rnfmsk array) 
    373376      ! 
    374       INTEGER  ::   jc, jn      ! dummy loop indices 
    375       INTEGER  ::   ii, ij      ! temporary integer 
     377      INTEGER  ::   jc, jn, ji, jj      ! dummy loop indices 
    376378      !!---------------------------------------------------------------------- 
    377379      ! 
     
    379381         IF( ncstt(jc) >= 1 ) THEN            ! runoff mask set to 1 at closed sea outflows 
    380382             DO jn = 1, 4 
    381                ii = mi0( ncsir(jc,jn) ) 
    382                ij = mj0( ncsjr(jc,jn) ) 
    383                p_rnfmsk(ii,ij) = MAX( p_rnfmsk(ii,ij), 1.0_wp ) 
     383                DO jj =    mj0( ncsjr(jc,jn) ), mj1( ncsjr(jc,jn) ) 
     384                   DO ji = mi0( ncsir(jc,jn) ), mi1( ncsir(jc,jn) ) 
     385                      p_rnfmsk(ji,jj) = MAX( p_rnfmsk(ji,jj), 1.0_wp ) 
     386                   END DO 
     387                END DO 
    384388            END DO  
    385389         ENDIF  
  • branches/2013/dev_CMCC_2013/NEMOGCM/NEMO/OPA_SRC/DOM/daymod.F90

    r3851 r4177  
    238238               nday_year = 1 
    239239               nsec_year = ndt05 
     240               IF( nsec1jan000 >= 2 * (2**30 - nsecd * nyear_len(1) / 2 ) ) THEN   ! test integer 4 max value 
     241                  CALL ctl_stop( 'The number of seconds between Jan. 1st 00h of nit000 year and Jan. 1st 00h ',   & 
     242                     &           'of the current year is exceeding the INTEGER 4 max VALUE: 2^31-1 -> 68.09 years in seconds', & 
     243                     & 'You must do a restart at higher frequency (or remove this STOP and recompile everything in I8)' ) 
     244               ENDIF 
    240245               nsec1jan000 = nsec1jan000 + nsecd * nyear_len(1) 
    241246               IF( nleapy == 1 )   CALL day_mth 
  • branches/2013/dev_CMCC_2013/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90

    r4175 r4177  
    11021102      INTEGER  ::   iip1, ijp1, iim1, ijm1   ! temporary integers 
    11031103      REAL(wp) ::   zrmax, ztaper   ! temporary scalars 
    1104       REAL(wp) ::   zrfact   ! temporary scalars 
    1105       REAL(wp), POINTER, DIMENSION(:,:  ) :: ztmpi1, ztmpi2, ztmpj1, ztmpj2 
    1106  
    1107       ! 
    1108       REAL(wp), POINTER, DIMENSION(:,:  ) :: zenv, zri, zrj, zhbat 
     1104      ! 
     1105      REAL(wp), POINTER, DIMENSION(:,:  ) :: zenv, ztmp, zmsk, zri, zrj, zhbat 
    11091106 
    11101107      NAMELIST/namzgr_sco/ln_s_sh94, ln_s_sf12, ln_sigcrit, rn_sbot_min, rn_sbot_max, rn_hc, rn_rmax,rn_theta, & 
     
    11141111      IF( nn_timing == 1 )  CALL timing_start('zgr_sco') 
    11151112      ! 
    1116       CALL wrk_alloc( jpi, jpj,      ztmpi1, ztmpi2, ztmpj1, ztmpj2         ) 
    1117       CALL wrk_alloc( jpi, jpj,      zenv, zri, zrj, zhbat     ) 
    1118      ! 
     1113      CALL wrk_alloc( jpi, jpj,      zenv, ztmp, zmsk, zri, zrj, zhbat                           ) 
     1114      ! 
    11191115      REWIND( numnam )                       ! Read Namelist namzgr_sco : sigma-stretching parameters 
    11201116      READ  ( numnam, namzgr_sco ) 
     
    11631159      !                                        ! ============================= 
    11641160      ! use r-value to create hybrid coordinates 
    1165 !     DO jj = 1, jpj 
    1166 !        DO ji = 1, jpi 
    1167 !           zenv(ji,jj) = MAX( bathy(ji,jj), 0._wp ) 
    1168 !        END DO 
    1169 !     END DO 
    1170 !     CALL lbc_lnk( zenv, 'T', 1._wp ) 
    1171       zenv(:,:) = bathy(:,:) 
     1161      DO jj = 1, jpj 
     1162         DO ji = 1, jpi 
     1163            zenv(ji,jj) = MAX( bathy(ji,jj), rn_sbot_min ) 
     1164         END DO 
     1165      END DO 
    11721166      !  
    11731167      ! Smooth the bathymetry (if required) 
     
    11771171      jl = 0 
    11781172      zrmax = 1._wp 
    1179       !      
    1180       ! set scaling factor used in reducing vertical gradients 
    1181       zrfact = ( 1._wp - rn_rmax ) / ( 1._wp + rn_rmax )  
    1182       ! 
    1183       ! initialise temporary evelope depth arrays 
    1184       ztmpi1(:,:) = zenv(:,:) 
    1185       ztmpi2(:,:) = zenv(:,:) 
    1186       ztmpj1(:,:) = zenv(:,:) 
    1187       ztmpj2(:,:) = zenv(:,:) 
    1188       ! 
    1189       ! initialise temporary r-value arrays 
    1190       zri(:,:) = 1._wp 
    1191       zrj(:,:) = 1._wp 
    1192       !                                                            ! ================ ! 
    1193       DO WHILE( jl <= 10000 .AND. ( zrmax - rn_rmax ) > 1.e-8_wp ) !  Iterative loop  ! 
    1194          !                                                         ! ================ ! 
     1173      !                                                     ! ================ ! 
     1174      DO WHILE( jl <= 10000 .AND. zrmax > rn_rmax )         !  Iterative loop  ! 
     1175         !                                                  ! ================ ! 
    11951176         jl = jl + 1 
    11961177         zrmax = 0._wp 
    1197          ! we set zrmax from previous r-values (zri abd zrj) first 
    1198          ! if set after current r-value calculation (as previously) 
    1199          ! we could exit DO WHILE prematurely before checking r-value 
    1200          ! of current zenv 
    1201          DO jj = 1, nlcj 
    1202             DO ji = 1, nlci 
    1203                zrmax = MAX( zrmax, ABS(zri(ji,jj)), ABS(zrj(ji,jj)) ) 
    1204             END DO 
    1205          END DO 
    1206          zri(:,:) = 0._wp 
    1207          zrj(:,:) = 0._wp 
     1178         zmsk(:,:) = 0._wp 
    12081179         DO jj = 1, nlcj 
    12091180            DO ji = 1, nlci 
    12101181               iip1 = MIN( ji+1, nlci )      ! force zri = 0 on last line (ji=ncli+1 to jpi) 
    12111182               ijp1 = MIN( jj+1, nlcj )      ! force zrj = 0 on last raw  (jj=nclj+1 to jpj) 
    1212                IF( (zenv(ji,jj) > 0._wp) .AND. (zenv(iip1,jj) > 0._wp)) THEN 
    1213                   zri(ji,jj) = ( zenv(iip1,jj  ) - zenv(ji,jj) ) / ( zenv(iip1,jj  ) + zenv(ji,jj) ) 
    1214                END IF 
    1215                IF( (zenv(ji,jj) > 0._wp) .AND. (zenv(ji,ijp1) > 0._wp)) THEN 
    1216                   zrj(ji,jj) = ( zenv(ji  ,ijp1) - zenv(ji,jj) ) / ( zenv(ji  ,ijp1) + zenv(ji,jj) ) 
    1217                END IF 
    1218                IF( zri(ji,jj) >  rn_rmax )   ztmpi1(ji  ,jj  ) = zenv(iip1,jj  ) * zrfact 
    1219                IF( zri(ji,jj) < -rn_rmax )   ztmpi2(iip1,jj  ) = zenv(ji  ,jj  ) * zrfact  
    1220                IF( zrj(ji,jj) >  rn_rmax )   ztmpj1(ji  ,jj  ) = zenv(ji  ,ijp1) * zrfact 
    1221                IF( zrj(ji,jj) < -rn_rmax )   ztmpj2(ji  ,ijp1) = zenv(ji  ,jj  ) * zrfact 
     1183               zri(ji,jj) = ABS( zenv(iip1,jj  ) - zenv(ji,jj) ) / ( zenv(iip1,jj  ) + zenv(ji,jj) ) 
     1184               zrj(ji,jj) = ABS( zenv(ji  ,ijp1) - zenv(ji,jj) ) / ( zenv(ji  ,ijp1) + zenv(ji,jj) ) 
     1185               zrmax = MAX( zrmax, zri(ji,jj), zrj(ji,jj) ) 
     1186               IF( zri(ji,jj) > rn_rmax )   zmsk(ji  ,jj  ) = 1._wp 
     1187               IF( zri(ji,jj) > rn_rmax )   zmsk(iip1,jj  ) = 1._wp 
     1188               IF( zrj(ji,jj) > rn_rmax )   zmsk(ji  ,jj  ) = 1._wp 
     1189               IF( zrj(ji,jj) > rn_rmax )   zmsk(ji  ,ijp1) = 1._wp 
    12221190            END DO 
    12231191         END DO 
    12241192         IF( lk_mpp )   CALL mpp_max( zrmax )   ! max over the global domain 
     1193         ! lateral boundary condition on zmsk: keep 1 along closed boundary (use of MAX) 
     1194         ztmp(:,:) = zmsk(:,:)   ;   CALL lbc_lnk( zmsk, 'T', 1._wp ) 
     1195         DO jj = 1, nlcj 
     1196            DO ji = 1, nlci 
     1197                zmsk(ji,jj) = MAX( zmsk(ji,jj), ztmp(ji,jj) ) 
     1198            END DO 
     1199         END DO 
    12251200         ! 
    1226          IF(lwp)WRITE(numout,*) 'zgr_sco :   iter= ',jl, ' rmax= ', zrmax 
     1201         IF(lwp)WRITE(numout,*) 'zgr_sco :   iter= ',jl, ' rmax= ', zrmax, ' nb of pt= ', INT( SUM(zmsk(:,:) ) ) 
    12271202         ! 
    12281203         DO jj = 1, nlcj 
    12291204            DO ji = 1, nlci 
    1230                zenv(ji,jj) = MAX(zenv(ji,jj), ztmpi1(ji,jj), ztmpi2(ji,jj), ztmpj1(ji,jj), ztmpj2(ji,jj) ) 
     1205               iip1 = MIN( ji+1, nlci )     ! last  line (ji=nlci) 
     1206               ijp1 = MIN( jj+1, nlcj )     ! last  raw  (jj=nlcj) 
     1207               iim1 = MAX( ji-1,  1  )      ! first line (ji=nlci) 
     1208               ijm1 = MAX( jj-1,  1  )      ! first raw  (jj=nlcj) 
     1209               IF( zmsk(ji,jj) == 1._wp ) THEN 
     1210                  ztmp(ji,jj) =   (                                                                                   & 
     1211             &      zenv(iim1,ijp1)*zmsk(iim1,ijp1) + zenv(ji,ijp1)*zmsk(ji,ijp1) + zenv(iip1,ijp1)*zmsk(iip1,ijp1)   & 
     1212             &    + zenv(iim1,jj  )*zmsk(iim1,jj  ) + zenv(ji,jj  )*    2._wp     + zenv(iip1,jj  )*zmsk(iip1,jj  )   & 
     1213             &    + zenv(iim1,ijm1)*zmsk(iim1,ijm1) + zenv(ji,ijm1)*zmsk(ji,ijm1) + zenv(iip1,ijm1)*zmsk(iip1,ijm1)   & 
     1214             &                    ) / (                                                                               & 
     1215             &                      zmsk(iim1,ijp1) +               zmsk(ji,ijp1) +                 zmsk(iip1,ijp1)   & 
     1216             &    +                 zmsk(iim1,jj  ) +                   2._wp     +                 zmsk(iip1,jj  )   & 
     1217             &    +                 zmsk(iim1,ijm1) +               zmsk(ji,ijm1) +                 zmsk(iip1,ijm1)   & 
     1218             &                        ) 
     1219               ENDIF 
    12311220            END DO 
    12321221         END DO 
    12331222         ! 
    1234          CALL lbc_lnk( zenv, 'T', 1._wp ) 
     1223         DO jj = 1, nlcj 
     1224            DO ji = 1, nlci 
     1225               IF( zmsk(ji,jj) == 1._wp )   zenv(ji,jj) = MAX( ztmp(ji,jj), bathy(ji,jj) ) 
     1226            END DO 
     1227         END DO 
     1228         ! 
     1229         ! Apply lateral boundary condition   CAUTION: keep the value when the lbc field is zero 
     1230         ztmp(:,:) = zenv(:,:)   ;   CALL lbc_lnk( zenv, 'T', 1._wp ) 
     1231         DO jj = 1, nlcj 
     1232            DO ji = 1, nlci 
     1233               IF( zenv(ji,jj) == 0._wp )   zenv(ji,jj) = ztmp(ji,jj) 
     1234            END DO 
     1235         END DO 
    12351236         !                                                  ! ================ ! 
    12361237      END DO                                                !     End loop     ! 
    12371238      !                                                     ! ================ ! 
    12381239      ! 
    1239 !     DO jj = 1, jpj 
    1240 !        DO ji = 1, jpi 
    1241 !           zenv(ji,jj) = MAX( zenv(ji,jj), rn_sbot_min ) ! set all points to avoid undefined scale values 
    1242 !        END DO 
    1243 !     END DO 
     1240      ! Fill ghost rows with appropriate values to avoid undefined e3 values with some mpp decompositions 
     1241      DO ji = nlci+1, jpi  
     1242         zenv(ji,1:nlcj) = zenv(nlci,1:nlcj) 
     1243      END DO 
     1244      ! 
     1245      DO jj = nlcj+1, jpj 
     1246         zenv(:,jj) = zenv(:,nlcj) 
     1247      END DO 
    12441248      ! 
    12451249      ! Envelope bathymetry saved in hbatt 
    12461250      hbatt(:,:) = zenv(:,:)  
    1247  
    12481251      IF( MINVAL( gphit(:,:) ) * MAXVAL( gphit(:,:) ) <= 0._wp ) THEN 
    12491252         CALL ctl_warn( ' s-coordinates are tapered in vicinity of the Equator' ) 
    12501253         DO jj = 1, jpj 
    12511254            DO ji = 1, jpi 
    1252                ztaper = EXP( -(gphit(ji,jj)/8._wp)**2 ) 
     1255               ztaper = EXP( -(gphit(ji,jj)/8._wp)**2._wp ) 
    12531256               hbatt(ji,jj) = rn_sbot_max * ztaper + hbatt(ji,jj) * ( 1._wp - ztaper ) 
    12541257            END DO 
     
    13651368      fsde3w(:,:,:) = gdep3w(:,:,:) 
    13661369      ! 
    1367       where (e3t   (:,:,:).eq.0.0)  e3t(:,:,:) = 1.0 
    1368       where (e3u   (:,:,:).eq.0.0)  e3u(:,:,:) = 1.0 
    1369       where (e3v   (:,:,:).eq.0.0)  e3v(:,:,:) = 1.0 
    1370       where (e3f   (:,:,:).eq.0.0)  e3f(:,:,:) = 1.0 
    1371       where (e3w   (:,:,:).eq.0.0)  e3w(:,:,:) = 1.0 
    1372       where (e3uw  (:,:,:).eq.0.0)  e3uw(:,:,:) = 1.0 
    1373       where (e3vw  (:,:,:).eq.0.0)  e3vw(:,:,:) = 1.0 
     1370      where (e3t   (:,:,:).eq.0.0)  e3t(:,:,:) = 1._wp 
     1371      where (e3u   (:,:,:).eq.0.0)  e3u(:,:,:) = 1._wp 
     1372      where (e3v   (:,:,:).eq.0.0)  e3v(:,:,:) = 1._wp 
     1373      where (e3f   (:,:,:).eq.0.0)  e3f(:,:,:) = 1._wp 
     1374      where (e3w   (:,:,:).eq.0.0)  e3w(:,:,:) = 1._wp 
     1375      where (e3uw  (:,:,:).eq.0.0)  e3uw(:,:,:) = 1._wp 
     1376      where (e3vw  (:,:,:).eq.0.0)  e3vw(:,:,:) = 1._wp 
    13741377 
    13751378#if defined key_agrif 
     
    15191522      END DO 
    15201523      ! 
    1521       CALL wrk_dealloc( jpi, jpj,      zenv, ztmpi1, ztmpi2, ztmpj1, ztmpj2, zri, zrj, zhbat                           )      ! 
     1524      CALL wrk_dealloc( jpi, jpj,      zenv, ztmp, zmsk, zri, zrj, zhbat                           ) 
     1525      ! 
    15221526      IF( nn_timing == 1 )  CALL timing_stop('zgr_sco') 
    15231527      ! 
     
    17481752      ENDDO 
    17491753      ! 
    1750       CALL lbc_lnk(e3t ,'T',1.) ; CALL lbc_lnk(e3u ,'T',1.) 
    1751       CALL lbc_lnk(e3v ,'T',1.) ; CALL lbc_lnk(e3f ,'T',1.) 
    1752       CALL lbc_lnk(e3w ,'T',1.) 
    1753       CALL lbc_lnk(e3uw,'T',1.) ; CALL lbc_lnk(e3vw,'T',1.) 
    1754       ! 
    17551754      !                                               ! ============= 
    17561755 
     
    18491848      !!---------------------------------------------------------------------- 
    18501849      ! 
    1851       pf =   (   TANH( rn_theta * ( -(pk-0.5_wp) / REAL(jpkm1) + rn_thetb )  )   & 
     1850      pf =   (   TANH( rn_theta * ( -(pk-0.5_wp) / REAL(jpkm1,wp) + rn_thetb )  )   & 
    18521851         &     - TANH( rn_thetb * rn_theta                                )  )   & 
    18531852         & * (   COSH( rn_theta                           )                      & 
     
    18751874      ! 
    18761875      IF ( rn_theta == 0 ) then      ! uniform sigma 
    1877          pf1 = - ( pk1 - 0.5_wp ) / REAL( jpkm1 ) 
     1876         pf1 = - ( pk1 - 0.5_wp ) / REAL( jpkm1,wp ) 
    18781877      ELSE                        ! stretched sigma 
    1879          pf1 =   ( 1._wp - pbb ) * ( SINH( rn_theta*(-(pk1-0.5_wp)/REAL(jpkm1)) ) ) / SINH( rn_theta )              & 
    1880             &  + pbb * (  (TANH( rn_theta*( (-(pk1-0.5_wp)/REAL(jpkm1)) + 0.5_wp) ) - TANH( 0.5_wp * rn_theta )  )  & 
     1878         pf1 =   ( 1._wp - pbb ) * ( SINH( rn_theta*(-(pk1-0.5_wp)/REAL(jpkm1,wp)) ) ) / SINH( rn_theta )              & 
     1879            &  + pbb * (  (TANH( rn_theta*( (-(pk1-0.5_wp)/REAL(jpkm1,wp)) + 0.5_wp) ) - TANH( 0.5_wp * rn_theta )  )  & 
    18811880            &        / ( 2._wp * TANH( 0.5_wp * rn_theta ) )  ) 
    18821881      ENDIF 
  • branches/2013/dev_CMCC_2013/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_flt.F90

    r3765 r4177  
    109109      INTEGER  ::   ji, jj, jk   ! dummy loop indices 
    110110      REAL(wp) ::   z2dt, z2dtg, zgcb, zbtd, ztdgu, ztdgv   ! local scalars 
    111       REAL(wp), POINTER, DIMENSION(:,:,:) ::  zub, zvb 
    112111      !!---------------------------------------------------------------------- 
    113112      ! 
    114113      IF( nn_timing == 1 )  CALL timing_start('dyn_spg_flt') 
    115114      ! 
    116       CALL wrk_alloc( jpi,jpj,jpk, zub, zvb ) 
    117115      ! 
    118116      IF( kt == nit000 ) THEN 
     
    213211         DO jk = 1, jpkm1 
    214212            DO ji = 1, jpij 
    215                spgu(ji,1) = spgu(ji,1) + fse3u(ji,1,jk) * ua(ji,1,jk) 
    216                spgv(ji,1) = spgv(ji,1) + fse3v(ji,1,jk) * va(ji,1,jk) 
     213               spgu(ji,1) = spgu(ji,1) + fse3u_a(ji,1,jk) * ua(ji,1,jk) 
     214               spgv(ji,1) = spgv(ji,1) + fse3v_a(ji,1,jk) * va(ji,1,jk) 
    217215            END DO 
    218216         END DO 
     
    221219            DO jj = 2, jpjm1 
    222220               DO ji = 2, jpim1 
    223                   spgu(ji,jj) = spgu(ji,jj) + fse3u(ji,jj,jk) * ua(ji,jj,jk) 
    224                   spgv(ji,jj) = spgv(ji,jj) + fse3v(ji,jj,jk) * va(ji,jj,jk) 
     221                  spgu(ji,jj) = spgu(ji,jj) + fse3u_a(ji,jj,jk) * ua(ji,jj,jk) 
     222                  spgv(ji,jj) = spgv(ji,jj) + fse3v_a(ji,jj,jk) * va(ji,jj,jk) 
    225223               END DO 
    226224            END DO 
     
    360358      IF( lrst_oce ) CALL flt_rst( kt, 'WRITE' ) 
    361359      ! 
    362       CALL wrk_dealloc( jpi,jpj,jpk, zub, zvb ) 
    363360      ! 
    364361      IF( nn_timing == 1 )  CALL timing_stop('dyn_spg_flt') 
  • branches/2013/dev_CMCC_2013/NEMOGCM/NEMO/OPA_SRC/LBC/lib_mpp.F90

    r4174 r4177  
    21812181!!gm Remark : this is very time consumming!!! 
    21822182      !                                         ! ------------------------ ! 
    2183             IF( ijpt0 > ijpt1 .OR. iipt0 > iipt1 ) THEN 
     2183        IF(((nbondi .ne. 0) .AND. (ktype .eq. 2)) .OR. ((nbondj .ne. 0) .AND. (ktype .eq. 1))) THEN 
    21842184            ! there is nothing to be migrated 
    2185                lmigr = .FALSE. 
     2185              lmigr = .TRUE. 
    21862186            ELSE 
    2187               lmigr = .TRUE. 
     2187              lmigr = .FALSE. 
    21882188            ENDIF 
    21892189 
  • branches/2013/dev_CMCC_2013/NEMOGCM/NEMO/OPA_SRC/SBC/geo2ocean.F90

    r2715 r4177  
    187187         &      gsinf(jpi,jpj), gcosf(jpi,jpj), STAT=ierr ) 
    188188      IF(lk_mpp)   CALL mpp_sum( ierr ) 
    189       IF( ierr /= 0 )   CALL ctl_stop('STOP', 'angle_msh_geo: unable to allocate arrays' ) 
     189      IF( ierr /= 0 )   CALL ctl_stop('angle: unable to allocate arrays' ) 
    190190 
    191191      ! ============================= ! 
     
    361361            &      gsinlat(jpi,jpj,4) , gcoslat(jpi,jpj,4) , STAT=ierr ) 
    362362         IF( lk_mpp    )   CALL mpp_sum( ierr ) 
    363          IF( ierr /= 0 )   CALL ctl_stop('STOP', 'angle_msh_geo: unable to allocate arrays' ) 
     363         IF( ierr /= 0 )   CALL ctl_stop('geo2oce: unable to allocate arrays' ) 
    364364      ENDIF 
    365365 
     
    438438      !!---------------------------------------------------------------------- 
    439439 
    440       IF( ALLOCATED( gsinlon ) ) THEN 
     440      IF( .NOT. ALLOCATED( gsinlon ) ) THEN 
    441441         ALLOCATE( gsinlon(jpi,jpj,4) , gcoslon(jpi,jpj,4) ,   & 
    442442            &      gsinlat(jpi,jpj,4) , gcoslat(jpi,jpj,4) , STAT=ierr ) 
    443443         IF( lk_mpp    )   CALL mpp_sum( ierr ) 
    444          IF( ierr /= 0 )   CALL ctl_stop('STOP', 'angle_msh_geo: unable to allocate arrays' ) 
     444         IF( ierr /= 0 )   CALL ctl_stop('oce2geo: unable to allocate arrays' ) 
    445445      ENDIF 
    446446 
  • branches/2013/dev_CMCC_2013/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90

    r3914 r4177  
    388388      ! 
    389389      IF( TRIM( sn_rcv_tau%cldes ) /= 'oce and ice' ) THEN        ! 'oce and ice' case ocean stress on ocean mesh used 
    390          srcv(jpr_itz1:jpr_itz2)%laction = .FALSE.    ! ice components not received (itx1 and ity1 used later) 
     390         srcv(jpr_itx1:jpr_itz2)%laction = .FALSE.    ! ice components not received 
    391391         srcv(jpr_itx1)%clgrid = 'U'                  ! ocean stress used after its transformation 
    392392         srcv(jpr_ity1)%clgrid = 'V'                  ! i.e. it is always at U- & V-points for i- & j-comp. resp. 
     
    407407      SELECT CASE( TRIM( sn_rcv_emp%cldes ) ) 
    408408      CASE( 'oce only'      )   ;   srcv(                                 jpr_oemp   )%laction = .TRUE.  
    409       CASE( 'conservative'  )   ;   srcv( (/jpr_rain, jpr_snow, jpr_ievp, jpr_tevp/) )%laction = .TRUE. 
     409      CASE( 'conservative'  ) 
     410         srcv( (/jpr_rain, jpr_snow, jpr_ievp, jpr_tevp/) )%laction = .TRUE. 
     411         IF ( k_ice <= 1 )  srcv(jpr_ivep)%laction = .FALSE. 
    410412      CASE( 'oce and ice'   )   ;   srcv( (/jpr_ievp, jpr_sbpr, jpr_semp, jpr_oemp/) )%laction = .TRUE. 
    411413      CASE default              ;   CALL ctl_stop( 'sbc_cpl_init: wrong definition of sn_rcv_emp%cldes' ) 
     
    465467         CALL ctl_stop( 'sbc_cpl_init: namsbc_cpl namelist mismatch between sn_rcv_qns%cldes and sn_rcv_dqnsdt%cldes' ) 
    466468      !                                                      ! ------------------------- ! 
    467       !                                                      !    Ice Qsr penetration    !    
    468       !                                                      ! ------------------------- ! 
    469       ! fraction of net shortwave radiation which is not absorbed in the thin surface layer  
    470       ! and penetrates inside the ice cover ( Maykut and Untersteiner, 1971 ; Elbert anbd Curry, 1993 ) 
    471       ! Coupled case: since cloud cover is not received from atmosphere  
    472       !               ===> defined as constant value -> definition done in sbc_cpl_init 
    473       fr1_i0(:,:) = 0.18 
    474       fr2_i0(:,:) = 0.82 
    475       !                                                      ! ------------------------- ! 
    476469      !                                                      !      10m wind module      !    
    477470      !                                                      ! ------------------------- ! 
     
    508501      ! Allocate taum part of frcv which is used even when not received as coupling field 
    509502      IF ( .NOT. srcv(jpr_taum)%laction ) ALLOCATE( frcv(jpr_taum)%z3(jpi,jpj,srcv(jn)%nct) ) 
     503      ! Allocate itx1 and ity1 as they are used in sbc_cpl_ice_tau even if srcv(jpr_itx1)%laction = .FALSE. 
     504      IF( k_ice /= 0 ) THEN 
     505         IF ( .NOT. srcv(jpr_itx1)%laction ) ALLOCATE( frcv(jpr_itx1)%z3(jpi,jpj,srcv(jn)%nct) ) 
     506         IF ( .NOT. srcv(jpr_ity1)%laction ) ALLOCATE( frcv(jpr_ity1)%z3(jpi,jpj,srcv(jn)%nct) ) 
     507      END IF 
    510508 
    511509      ! ================================ ! 
     
    13291327      END SELECT 
    13301328 
     1329      !    Ice Qsr penetration used (only?)in lim2 or lim3  
     1330      ! fraction of net shortwave radiation which is not absorbed in the thin surface layer  
     1331      ! and penetrates inside the ice cover ( Maykut and Untersteiner, 1971 ; Elbert anbd Curry, 1993 ) 
     1332      ! Coupled case: since cloud cover is not received from atmosphere  
     1333      !               ===> defined as constant value -> definition done in sbc_cpl_init 
     1334      fr1_i0(:,:) = 0.18 
     1335      fr2_i0(:,:) = 0.82 
     1336 
     1337 
    13311338      CALL wrk_dealloc( jpi,jpj, zcptn, ztmp, zicefr ) 
    13321339      ! 
  • branches/2013/dev_CMCC_2013/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90

    r3905 r4177  
    221221      ENDIF 
    222222      ! 
     223                          CALL sbc_ssm_init               ! Sea-surface mean fields initialisation 
     224      ! 
    223225      IF( ln_ssr      )   CALL sbc_ssr_init               ! Sea-Surface Restoring initialisation 
    224226      ! 
  • branches/2013/dev_CMCC_2013/NEMOGCM/NEMO/OPA_SRC/TRA/eosbn2.F90

    r3625 r4177  
    675675 
    676676 
    677    FUNCTION tfreez( psal ) RESULT( ptf ) 
     677   FUNCTION tfreez( psal, pdep ) RESULT( ptf ) 
    678678      !!---------------------------------------------------------------------- 
    679679      !!                 ***  ROUTINE eos_init  *** 
     
    688688      !!---------------------------------------------------------------------- 
    689689      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   psal   ! salinity             [psu] 
     690      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ), OPTIONAL ::   pdep   ! depth      [decibars] 
    690691      ! Leave result array automatic rather than making explicitly allocated 
    691692      REAL(wp), DIMENSION(jpi,jpj)                ::   ptf    ! freezing temperature [Celcius] 
     
    694695      ptf(:,:) = ( - 0.0575_wp + 1.710523e-3_wp * SQRT( psal(:,:) )   & 
    695696         &                     - 2.154996e-4_wp *       psal(:,:)   ) * psal(:,:) 
     697      IF ( PRESENT( pdep ) ) THEN    
     698         ptf(:,:) = ptf(:,:) - 7.53e-4_wp * pdep(:,:) 
     699      ENDIF 
    696700      ! 
    697701   END FUNCTION tfreez 
  • branches/2013/dev_CMCC_2013/NEMOGCM/NEMO/SAS_SRC/daymod.F90

    r3851 r4177  
    246246               nday_year = 1 
    247247               nsec_year = ndt05 
     248               IF( nsec1jan000 >= 2 * (2**30 - nsecd * nyear_len(1) / 2 ) ) THEN   ! test integer 4 max value 
     249                  CALL ctl_stop( 'The number of seconds between Jan. 1st 00h of nit000 year and Jan. 1st 00h ',   & 
     250                     &           'of the current year is exceeding the INTEGER 4 max VALUE: 2^31-1 -> 68.09 years in seconds', & 
     251                     & 'You must do a restart at higher frequency (or remove this STOP and recompile everything in I8)' ) 
     252               ENDIF 
    248253               nsec1jan000 = nsec1jan000 + nsecd * nyear_len(1) 
    249254               IF( nleapy == 1 )   CALL day_mth 
  • branches/2013/dev_CMCC_2013/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsed.F90

    r3905 r4177  
    8282      IF( nn_timing == 1 )  CALL timing_start('p4z_sed') 
    8383      ! 
    84       IF( kt == nit000 .AND. jnt == 1 )  THEN 
     84      IF( kt == nittrc000 .AND. jnt == 1 )  THEN 
    8585         ryyss    = nyear_len(1) * rday    ! number of seconds per year and per month 
    8686         rmtss    = ryyss / raamo 
  • branches/2013/dev_CMCC_2013/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsms.F90

    r4175 r4177  
    7676      ENDIF 
    7777      ! 
    78       IF( ln_rsttr .AND. kt == nittrc000 )                         CALL p4z_rst( nittrc000, 'READ' )  !* read or initialize all required fields  
     78      IF( kt == nittrc000 ) THEN 
     79        ! 
     80        CALL p4z_che                              ! initialize the chemical constants 
     81        ! 
     82        IF( .NOT. ln_rsttr ) THEN  ;   CALL p4z_ph_ini   !  set PH at kt=nit000  
     83        ELSE                       ;   CALL p4z_rst( nittrc000, 'READ' )  !* read or initialize all required fields  
     84        ENDIF 
     85        ! 
     86      ENDIF 
     87 
    7988      IF( ln_pisdmp .AND. MOD( kt - nn_dttrc, nn_pisdmp ) == 0 )   CALL p4z_dmp( kt )      ! Relaxation of some tracers 
    8089      ! 
     
    238247   END SUBROUTINE p4z_sms_init 
    239248 
     249   SUBROUTINE p4z_ph_ini 
     250      !!--------------------------------------------------------------------- 
     251      !!                   ***  ROUTINE p4z_ini_ph  *** 
     252      !! 
     253      !!  ** Purpose : Initialization of chemical variables of the carbon cycle 
     254      !!--------------------------------------------------------------------- 
     255      INTEGER  ::  ji, jj, jk 
     256      REAL(wp) ::  zcaralk, zbicarb, zco3 
     257      REAL(wp) ::  ztmas, ztmas1 
     258      !!--------------------------------------------------------------------- 
     259 
     260      ! Set PH from  total alkalinity, borat (???), akb3 (???) and ak23 (???) 
     261      ! -------------------------------------------------------- 
     262      DO jk = 1, jpk 
     263         DO jj = 1, jpj 
     264            DO ji = 1, jpi 
     265               ztmas   = tmask(ji,jj,jk) 
     266               ztmas1  = 1. - tmask(ji,jj,jk) 
     267               zcaralk = trn(ji,jj,jk,jptal) - borat(ji,jj,jk) / (  1. + 1.E-8 / ( rtrn + akb3(ji,jj,jk) )  ) 
     268               zco3    = ( zcaralk - trn(ji,jj,jk,jpdic) ) * ztmas + 0.5e-3 * ztmas1 
     269               zbicarb = ( 2. * trn(ji,jj,jk,jpdic) - zcaralk ) 
     270               hi(ji,jj,jk) = ( ak23(ji,jj,jk) * zbicarb / zco3 ) * ztmas + 1.e-9 * ztmas1 
     271            END DO 
     272         END DO 
     273     END DO 
     274     ! 
     275   END SUBROUTINE p4z_ph_ini 
     276 
    240277   SUBROUTINE p4z_rst( kt, cdrw ) 
    241278      !!--------------------------------------------------------------------- 
     
    266303         ELSE 
    267304!            hi(:,:,:) = 1.e-9  
    268             ! Set PH from  total alkalinity, borat (???), akb3 (???) and ak23 (???) 
    269             ! -------------------------------------------------------- 
    270             DO jk = 1, jpk 
    271                DO jj = 1, jpj 
    272                   DO ji = 1, jpi 
    273                      ztmas   = tmask(ji,jj,jk) 
    274                      ztmas1  = 1. - tmask(ji,jj,jk) 
    275                      zcaralk = trn(ji,jj,jk,jptal) - borat(ji,jj,jk) / (  1. + 1.E-8 / ( rtrn + akb3(ji,jj,jk) )  ) 
    276                      zco3    = ( zcaralk - trn(ji,jj,jk,jpdic) ) * ztmas + 0.5e-3 * ztmas1 
    277                      zbicarb = ( 2. * trn(ji,jj,jk,jpdic) - zcaralk ) 
    278                      hi(ji,jj,jk) = ( ak23(ji,jj,jk) * zbicarb / zco3 ) * ztmas + 1.e-9 * ztmas1 
    279                   END DO 
    280                END DO 
    281             END DO 
     305            CALL p4z_ph_ini 
    282306         ENDIF 
    283307         CALL iom_get( numrtr, jpdom_autoglo, 'Silicalim', xksi(:,:) ) 
  • branches/2013/dev_CMCC_2013/NEMOGCM/NEMO/TOP_SRC/PISCES/trcini_pisces.F90

    r3757 r4177  
    122122      rdenita =   3._wp /  5._wp 
    123123      o2ut    = 131._wp / 122._wp 
    124  
    125       CALL p4z_che        ! initialize the chemical constants 
    126124 
    127125      ! Initialization of tracer concentration in case of  no restart  
     
    162160         xksi(:,:)    = 2.e-6 
    163161         xksimax(:,:) = xksi(:,:) 
    164  
    165          ! Initialization of chemical variables of the carbon cycle 
    166          ! -------------------------------------------------------- 
    167          DO jk = 1, jpk 
    168             DO jj = 1, jpj 
    169                DO ji = 1, jpi 
    170                   ztmas   = tmask(ji,jj,jk) 
    171                   ztmas1  = 1. - tmask(ji,jj,jk) 
    172                   zcaralk = trn(ji,jj,jk,jptal) - borat(ji,jj,jk) / (  1. + 1.E-8 / ( rtrn + akb3(ji,jj,jk) )  ) 
    173                   zco3    = ( zcaralk - trn(ji,jj,jk,jpdic) ) * ztmas + 0.5e-3 * ztmas1 
    174                   zbicarb = ( 2. * trn(ji,jj,jk,jpdic) - zcaralk ) 
    175                   hi(ji,jj,jk) = ( ak23(ji,jj,jk) * zbicarb / zco3 ) * ztmas + 1.e-9 * ztmas1 
    176                END DO 
    177             END DO 
    178          END DO 
    179          ! 
     162        ! 
    180163      END IF 
    181164 
  • branches/2013/dev_CMCC_2013/NEMOGCM/TOOLS/COMPILE/Fcheck_archfile.sh

    r3925 r4177  
    4040# :: 
    4141# 
    42 #  $ ./Fcheck_archfile.sh ARCHFILE COMPILER 
     42#  $ ./Fcheck_archfile.sh ARCHFILE CPPFILE COMPILER 
    4343# 
    4444# 
     
    9494   else 
    9595       if [ -f ${COMPIL_DIR}/$1 ]; then 
    96       # has the cpp keys file been changed since we copied the arch file in ${COMPIL_DIR}? 
    97       mycpp=$( ls -l ${COMPIL_DIR}/$2 | sed -e "s/.* -> //" ) 
    98       if [ "$mycpp" != "$( cat ${COMPIL_DIR}/cpp.history )" ]; then 
    99           echo $mycpp > ${COMPIL_DIR}/cpp.history 
    100           cpeval ${myarch} ${COMPIL_DIR}/$1 
     96      if [ "$2" != "nocpp" ]  
     97      then 
     98          # has the cpp keys file been changed since we copied the arch file in ${COMPIL_DIR}? 
     99          mycpp=$( ls -l ${COMPIL_DIR}/$2 | sed -e "s/.* -> //" ) 
     100          if [ "$mycpp" != "$( cat ${COMPIL_DIR}/cpp.history )" ]; then 
     101         echo $mycpp > ${COMPIL_DIR}/cpp.history 
     102         cpeval ${myarch} ${COMPIL_DIR}/$1 
     103          fi 
     104          # has the cpp keys file been updated since we copied the arch file in ${COMPIL_DIR}? 
     105          mycpp=$( find -L ${COMPIL_DIR} -cnewer ${COMPIL_DIR}/$1 -name $2 -print ) 
     106          [ ${#mycpp} -ne 0 ] && cpeval ${myarch} ${COMPIL_DIR}/$1 
    101107      fi 
    102       # has the cpp keys file been updated since we copied the arch file in ${COMPIL_DIR}? 
    103       mycpp=$( find -L ${COMPIL_DIR} -cnewer ${COMPIL_DIR}/$1 -name $2 -print ) 
    104       [ ${#mycpp} -ne 0 ] && cpeval ${myarch} ${COMPIL_DIR}/$1 
    105108      # has myarch file been updated since we copied it in ${COMPIL_DIR}? 
    106109      myarchdir=$( dirname ${myarch} ) 
     
    134137    if [ "$myarch" == "$( cat ${COMPIL_DIR}/arch.history )" ]; then  
    135138   if [ -f ${COMPIL_DIR}/$1 ]; then 
    136        # has the cpp keys file been changed since we copied the arch file in ${COMPIL_DIR}? 
    137        mycpp=$( ls -l ${COMPIL_DIR}/$2 | sed -e "s/.* -> //" ) 
    138        if [ "$mycpp" != "$( cat ${COMPIL_DIR}/cpp.history )" ]; then 
    139       echo $mycpp > ${COMPIL_DIR}/cpp.history 
    140       cpeval ${myarch} ${COMPIL_DIR}/$1 
     139       if [ "$2" != "nocpp" ]  
     140       then 
     141      # has the cpp keys file been changed since we copied the arch file in ${COMPIL_DIR}? 
     142      mycpp=$( ls -l ${COMPIL_DIR}/$2 | sed -e "s/.* -> //" ) 
     143      if [ "$mycpp" != "$( cat ${COMPIL_DIR}/cpp.history )" ]; then 
     144          echo $mycpp > ${COMPIL_DIR}/cpp.history 
     145          cpeval ${myarch} ${COMPIL_DIR}/$1 
     146      fi 
     147      # has the cpp keys file been updated since we copied the arch file in ${COMPIL_DIR}? 
     148      mycpp=$( find -L ${COMPIL_DIR} -cnewer ${COMPIL_DIR}/$1 -name $2 -print ) 
     149      [ ${#mycpp} -ne 0 ] && cpeval ${myarch} ${COMPIL_DIR}/$1 
    141150       fi 
    142        # has the cpp keys file been updated since we copied the arch file in ${COMPIL_DIR}? 
    143        mycpp=$( find -L ${COMPIL_DIR} -cnewer ${COMPIL_DIR}/$1 -name $2 -print ) 
    144        [ ${#mycpp} -ne 0 ] && cpeval ${myarch} ${COMPIL_DIR}/$1 
    145151       # has myarch file been updated since we copied it in ${COMPIL_DIR}? 
    146152       myarch=$( find -L ${MAIN_DIR}/ARCH -cnewer ${COMPIL_DIR}/$1 -name arch-${3}.fcm -print ) 
     
    150156   fi 
    151157    else 
    152    ls -l ${COMPIL_DIR}/$2 | sed -e "s/.* -> //" > ${COMPIL_DIR}/cpp.history 
     158   if [ "$2" != "nocpp" ]  
     159   then 
     160       ls -l ${COMPIL_DIR}/$2 | sed -e "s/.* -> //" > ${COMPIL_DIR}/cpp.history 
     161   fi 
    153162   echo ${myarch} > ${COMPIL_DIR}/arch.history 
    154163   cpeval ${myarch} ${COMPIL_DIR}/$1 
     
    157166 
    158167#- do we need xios library? 
    159 use_iom=$( sed -e "s/#.*$//" ${COMPIL_DIR}/$2 | grep -c key_iomput ) 
     168if [ "$2" != "nocpp" ]  
     169then 
     170    use_iom=$( sed -e "s/#.*$//" ${COMPIL_DIR}/$2 | grep -c key_iomput ) 
     171else 
     172    use_iom=0 
     173fi 
    160174have_lxios=$( sed -e "s/#.*$//" ${COMPIL_DIR}/$1 | grep -c "\-lxios" ) 
    161175if [[ ( $use_iom -eq 0 ) && ( $have_lxios -ge 1 ) ]] 
     
    166180 
    167181#- do we need oasis libraries? 
    168 use_oasis=$( sed -e "s/#.*$//" ${COMPIL_DIR}/$2 | grep -c key_oasis3 ) 
     182if [ "$2" != "nocpp" ]  
     183then 
     184    use_oasis=$( sed -e "s/#.*$//" ${COMPIL_DIR}/$2 | grep -c key_oasis3 ) 
     185else 
     186    use_oasis=0 
     187fi 
    169188for liboa in psmile.MPI1 mct mpeu scrip mpp_io 
    170189do 
  • branches/2013/dev_CMCC_2013/NEMOGCM/TOOLS/MISCELLANEOUS/chk_iomput.sh

    r4175 r4177  
    5959#------------------------------------------------ 
    6060# 
    61 external=$( grep -c "<field_definition.* src=" $xmlfile ) 
     61external=$( grep -c "<field_definition  *\([^ ].* \)*src=" $xmlfile ) 
    6262if [ $external -eq 1 ] 
    6363then 
    64     xmlfield_def=$( grep "<field_definition.* src=" $xmlfile | sed -e 's/.*src="\([^"]*\)".*/\1/' ) 
     64    xmlfield_def=$( grep "<field_definition  *\([^ ].* \)*src=" $xmlfile | sed -e 's/.*src="\([^"]*\)".*/\1/' ) 
    6565    xmlfield_def=$( dirname $xmlfile )/$xmlfield_def    
    6666else 
    6767    xmlfield_def=$xmlfile 
    6868fi 
    69 [ $inxml -eq 1 ] && grep "< *field * id *=" $xmlfield_def 
     69[ $inxml -eq 1 ] && grep "< *field  *\([^ ].* \)*id *=" $xmlfield_def 
    7070[ $insrc -eq 1 ] && find $srcdir -name "*.[Ffh]90" -exec grep -iH "^[^\!]*call  *iom_put *(" {} \; 
    7171[ $(( $insrc + $inxml )) -ge 1 ] && exit 
     
    9191# list of variables defined in the xml file 
    9292# 
    93 varlistxml=$( grep "< *field.* id *=" $xmlfield_def  | sed -e "s/^.*< *field.* id *= *[\"\']\([^\"\']*\)[\"\'].*/\1/" | sort -d ) 
     93varlistxml=$( grep "< *field  *\([^ ].* \)*id *=" $xmlfield_def  | sed -e "s/^.*< *field .*id *= *[\"\']\([^\"\']*\)[\"\'].*/\1/" | sort -d ) 
    9494# 
    9595# list of variables to be outputed in the xml file 
    9696# 
    97 varlistout=$( grep "< *field.* field_ref *=" $xmlfile  | sed -e "s/^.*< *field.* field_ref *= *[\"\']\([^\"\']*\)[\"\'].*/\1/" | sort -d ) 
     97varlistout=$( grep "< *field  *\([^ ].* \)*field_ref *=" $xmlfile  | sed -e "s/^.*< *field .*field_ref *= *[\"\']\([^\"\']*\)[\"\'].*/\1/" | sort -d ) 
    9898# 
    9999echo "--------------------------------------------------" 
  • branches/2013/dev_CMCC_2013/NEMOGCM/TOOLS/maketools

    r3294 r4177  
    146146 
    147147#- When used for the first time, choose a compiler --- 
    148 . ${COMPIL_DIR}/Fcheck_archfile.sh arch_tools.fcm ${CMP_NAM} || exit 
     148. ${COMPIL_DIR}/Fcheck_archfile.sh arch_tools.fcm nocpp ${CMP_NAM} || exit 
    149149 
    150150#- Choose a default tool if needed --- 
Note: See TracChangeset for help on using the changeset viewer.