Changeset 3963


Ignore:
Timestamp:
2013-07-09T17:41:20+02:00 (7 years ago)
Author:
clem
Message:

bugs correction + creation of glob_max and glob_min in lib_fortran.F90, see ticket:#1116

Location:
branches/2013/dev_r3406_CNRS_LIM3/NEMOGCM/NEMO
Files:
17 edited

Legend:

Unmodified
Added
Removed
  • branches/2013/dev_r3406_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limdiahsb.F90

    r3938 r3963  
    127127      CALL iom_put( 'ibgfrcemp',frc_vol * 1.e-9 )           !vol - forcing : km3 (equivalent liquid)          
    128128      CALL iom_put( 'ibgfrcemps',frc_sal * 1.e-9 )          !sal - forcing : psu*km3    
    129       CALL iom_put( 'ibgemps',zbg_emps *86400)             
    130       CALL iom_put( 'ibgemp',zbg_emp *86400)               
    131       CALL iom_put( 'ibgfsbri',zbg_fsbri *86400)          
    132       CALL iom_put( 'ibgfseqv',zbg_fseqv *86400)          
    133       CALL iom_put( 'ibgfsaltres',zbg_fsalt_res *86400)          
    134       CALL iom_put( 'ibgfsaltrpo',zbg_fsalt_rpo *86400)          
     129      CALL iom_put( 'ibgemps',zbg_emps *86400.)             
     130      CALL iom_put( 'ibgemp',zbg_emp *86400.)               
     131      CALL iom_put( 'ibgfsbri',zbg_fsbri *86400.)          
     132      CALL iom_put( 'ibgfseqv',zbg_fseqv *86400.)          
     133      CALL iom_put( 'ibgfsaltres',zbg_fsalt_res *86400.)          
     134      CALL iom_put( 'ibgfsaltrpo',zbg_fsalt_rpo *86400.)          
    135135      CALL iom_put( 'ibggrpme',bg_grpme * rhoic / rau0 * 1.e-9 ) ! km3 (equivalent liquid)          
    136136      ! 
  • branches/2013/dev_r3406_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limdyn.F90

    r3938 r3963  
    6666      REAL(wp), POINTER, DIMENSION(:,:) ::   zu_io, zv_io   ! ice-ocean velocity 
    6767      REAL(wp) :: zchk_v_i, zchk_smv, zchk_fs, zchk_fw, zchk_v_i_b, zchk_smv_b, zchk_fs_b, zchk_fw_b ! Check conservation (C Rousset) 
     68      REAL(wp) :: zchk_vmin, zchk_amin, zchk_amax ! Check errors (C Rousset) 
    6869     !!--------------------------------------------------------------------- 
    6970 
     
    234235         zchk_smv = ( glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zchk_smv_b ) / rdt_ice + ( zchk_fs / rhoic ) 
    235236 
     237         zchk_vmin = glob_min(v_i) 
     238         zchk_amax = glob_max(SUM(a_i,dim=3)) 
     239         zchk_amin = glob_min(a_i) 
     240 
    236241         IF(lwp) THEN 
    237             IF (    ABS( zchk_v_i   ) >  1.e-5 ) WRITE(numout,*) 'violation volume [m3/day]     (limdyn) = ',(zchk_v_i * 86400.) 
    238             IF (    ABS( zchk_smv   ) >  1.e-4  ) WRITE(numout,*) 'violation saline [psu*m3/day] (limdyn) = ',(zchk_smv * 86400.) 
    239             IF ( MINVAL( v_i(:,:,:) ) <  0.    ) WRITE(numout,*) 'violation v_i<0  [mm]         (limdyn) = ',(MINVAL(v_i) * 1.e-3) 
    240             IF ( MAXVAL( SUM(a_i(:,:,:),dim=3) ) > amax+1.e-10 ) WRITE(numout,*) 'violation a_i>amax    (limdyn) = ',MAXVAL(SUM(a_i,dim=3)) 
     242            IF ( ABS( zchk_v_i   ) >  1.e-5 ) WRITE(numout,*) 'violation volume [m3/day]     (limdyn) = ',(zchk_v_i * 86400.) 
     243            IF ( ABS( zchk_smv   ) >  1.e-4 ) WRITE(numout,*) 'violation saline [psu*m3/day] (limdyn) = ',(zchk_smv * 86400.) 
     244            IF ( zchk_vmin <  0.            ) WRITE(numout,*) 'violation v_i<0  [mm]         (limdyn) = ',(zchk_vmin * 1.e-3) 
     245            !IF ( zchk_amax >  amax+1.e-10   ) WRITE(numout,*) 'violation a_i>amax            (limdyn) = ',zchk_amax 
     246            IF ( zchk_amin <  0.            ) WRITE(numout,*) 'violation a_i<0               (limdyn) = ',zchk_amin 
    241247         ENDIF 
    242248         !CALL ctl_opn( numhsb , cl_name , 'UNKNOWN' , 'FORMATTED' , 'SEQUENTIAL' , 1 , numout , lwp , 1 ) 
  • branches/2013/dev_r3406_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90

    r3938 r3963  
    142142      REAL(wp), POINTER, DIMENSION(:,:) ::   vt_i_init, vt_i_final  !  ice volume summed over categories 
    143143      REAL(wp) :: zchk_v_i, zchk_smv, zchk_fs, zchk_fw, zchk_v_i_b, zchk_smv_b, zchk_fs_b, zchk_fw_b ! Check conservation (C Rousset) 
     144      REAL(wp) :: zchk_vmin, zchk_amin, zchk_amax ! Check errors (C Rousset) 
    144145      ! mass and salt flux (clem) 
    145146      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zviold, zvsold, zsmvold   ! old ice volume... 
     
    473474         zchk_smv = ( glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zchk_smv_b ) / rdt_ice + ( zchk_fs / rhoic ) 
    474475 
     476         zchk_vmin = glob_min(v_i) 
     477         zchk_amax = glob_max(SUM(a_i,dim=3)) 
     478         zchk_amin = glob_min(a_i) 
     479        
    475480         IF(lwp) THEN 
    476             IF (    ABS( zchk_v_i   ) >  1.e-5 ) WRITE(numout,*) 'violation volume [m3/day]     (limitd_me) = ',(zchk_v_i * 86400.) 
    477             IF (    ABS( zchk_smv   ) >  1.e-4 ) WRITE(numout,*) 'violation saline [psu*m3/day] (limitd_me) = ',(zchk_smv * 86400.) 
    478             IF ( MINVAL( v_i(:,:,:) ) <  0.    ) WRITE(numout,*) 'violation v_i<0  [mm]         (limitd_me) = ',(MINVAL(v_i) * 1.e-3) 
    479             IF ( MAXVAL( SUM(a_i(:,:,:),dim=3) ) > amax+epsi10 ) WRITE(numout,*) 'violation a_i>amax            (limitd_me) = ',MAXVAL(SUM(a_i,dim=3)) 
    480             IF ( MINVAL( a_i(:,:,:) ) <  0.    ) WRITE(numout,*) 'violation a_i<0               (limitd_me) = ',MINVAL(a_i) 
     481            IF ( ABS( zchk_v_i   ) >  1.e-5 ) WRITE(numout,*) 'violation volume [m3/day]     (limitd_me) = ',(zchk_v_i * 86400.) 
     482            IF ( ABS( zchk_smv   ) >  1.e-4 ) WRITE(numout,*) 'violation saline [psu*m3/day] (limitd_me) = ',(zchk_smv * 86400.) 
     483            IF ( zchk_vmin <  0.            ) WRITE(numout,*) 'violation v_i<0  [mm]         (limitd_me) = ',(zchk_vmin * 1.e-3) 
     484            IF ( zchk_amax >  amax+epsi10   ) WRITE(numout,*) 'violation a_i>amax            (limitd_me) = ',zchk_amax 
     485            IF ( zchk_amin <  0.            ) WRITE(numout,*) 'violation a_i<0               (limitd_me) = ',zchk_amin 
    481486         ENDIF 
    482487      ENDIF 
  • branches/2013/dev_r3406_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limitd_th.F90

    r3938 r3963  
    6868      INTEGER ::   jl, ja, jm, jbnd1, jbnd2   ! ice types    dummy loop index          
    6969      REAL(wp) :: zchk_v_i, zchk_smv, zchk_fs, zchk_fw, zchk_v_i_b, zchk_smv_b, zchk_fs_b, zchk_fw_b ! Check conservation (C Rousset) 
     70      REAL(wp) :: zchk_vmin, zchk_amin, zchk_amax ! Check errors (C Rousset) 
    7071      !!------------------------------------------------------------------ 
    7172      ! ------------------------------- 
     
    165166         zchk_smv = ( glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zchk_smv_b ) / rdt_ice + ( zchk_fs / rhoic ) 
    166167 
     168         zchk_vmin = glob_min(v_i) 
     169         zchk_amax = glob_max(SUM(a_i,dim=3)) 
     170         zchk_amin = glob_min(a_i) 
     171 
    167172         IF(lwp) THEN 
    168             IF (    ABS( zchk_v_i   ) >  1.e-5 ) WRITE(numout,*) 'violation volume [m3/day]     (limitd_th) = ',(zchk_v_i * 86400.) 
    169             IF (    ABS( zchk_smv   ) >  1.e-4 ) WRITE(numout,*) 'violation saline [psu*m3/day] (limitd_th) = ',(zchk_smv * 86400.) 
    170             IF ( MINVAL( v_i(:,:,:) ) <  0.    ) WRITE(numout,*) 'violation v_i<0  [mm]         (limitd_th) = ',(MINVAL(v_i) * 1.e-3) 
    171             IF ( MAXVAL( SUM(a_i(:,:,:),dim=3) ) >  amax+epsi10 ) WRITE(numout,*) 'violation a_i>amax    (limitd_th) = ',MAXVAL(SUM(a_i,dim=3)) 
     173            IF ( ABS( zchk_v_i   ) >  1.e-5 ) WRITE(numout,*) 'violation volume [m3/day]     (limitd_th) = ',(zchk_v_i * 86400.) 
     174            IF ( ABS( zchk_smv   ) >  1.e-4 ) WRITE(numout,*) 'violation saline [psu*m3/day] (limitd_th) = ',(zchk_smv * 86400.) 
     175            IF ( zchk_vmin <  0.            ) WRITE(numout,*) 'violation v_i<0  [mm]         (limitd_th) = ',(zchk_vmin * 1.e-3) 
     176            IF ( zchk_amax >  amax+epsi10   ) WRITE(numout,*) 'violation a_i>amax            (limitd_th) = ',zchk_amax 
     177            IF ( zchk_amin <  0.            ) WRITE(numout,*) 'violation a_i<0               (limitd_th) = ',zchk_amin 
    172178         ENDIF 
    173179       ENDIF 
     
    206212      INTEGER  ::   zji, zjj, nd   ! local integer 
    207213      REAL(wp) ::   zx1, zwk1, zdh0, zetamin, zdamax   ! local scalars 
    208       REAL(wp) ::   zx2, zwk2, zda0, zetamax, zhimin   !   -      - 
     214      REAL(wp) ::   zx2, zwk2, zda0, zetamax           !   -      - 
    209215      REAL(wp) ::   zx3,             zareamin, zindb   !   -      - 
    210216      CHARACTER (len = 15) :: fieldid 
     
    242248      CALL wrk_alloc( jpi,jpj, zhb0,zhb1,vt_i_init,vt_i_final,vt_s_init,vt_s_final,et_i_init,et_i_final,et_s_init,et_s_final ) 
    243249 
    244       zhimin   = 0.1      !minimum ice thickness tolerated by the model 
    245250      zareamin = epsi10   !minimum area in thickness categories tolerated by the conceptors of the model 
    246251 
     
    527532         zji = nind_i(ji) 
    528533         zjj = nind_j(ji) 
    529          IF ( ( zhimin .GT. 0.0 ) .AND. &  
    530             ( ( a_i(zji,zjj,1) .GT. epsi10 ) .AND. ( ht_i(zji,zjj,1) .LT. zhimin ) ) & 
    531             ) THEN 
    532             a_i(zji,zjj,1)  = a_i(zji,zjj,1) * ht_i(zji,zjj,1) / zhimin  
    533             ht_i(zji,zjj,1) = zhimin 
    534             v_i(zji,zjj,1)  = a_i(zji,zjj,1)*ht_i(zji,zjj,1) !clem@useless 
     534         IF ( ( a_i(zji,zjj,1) .GT. epsi10 ) .AND. ( ht_i(zji,zjj,1) .LT. hiclim ) ) THEN 
     535            a_i(zji,zjj,1)  = a_i(zji,zjj,1) * ht_i(zji,zjj,1) / hiclim  
     536            ht_i(zji,zjj,1) = hiclim 
     537            v_i(zji,zjj,1)  = a_i(zji,zjj,1) * ht_i(zji,zjj,1) !clem@useless 
    535538         ENDIF 
    536539      END DO !ji 
  • branches/2013/dev_r3406_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90

    r3938 r3963  
    9494      REAL(wp), POINTER, DIMENSION(:,:) ::   zqlbsbq   ! link with lead energy budget qldif 
    9595      REAL(wp) :: zchk_v_i, zchk_smv, zchk_fs, zchk_fw, zchk_v_i_b, zchk_smv_b, zchk_fs_b, zchk_fw_b ! Check conservation (C Rousset) 
     96      REAL(wp) :: zchk_vmin, zchk_amin, zchk_amax ! Check errors (C Rousset) 
    9697      !!------------------------------------------------------------------- 
    9798 
     
    476477         zchk_smv = ( glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zchk_smv_b ) / rdt_ice + ( zchk_fs / rhoic ) 
    477478 
     479         zchk_vmin = glob_min(v_i) 
     480         zchk_amax = glob_max(SUM(a_i,dim=3)) 
     481         zchk_amin = glob_min(a_i) 
     482 
    478483         IF(lwp) THEN 
    479             IF (    ABS( zchk_v_i   ) >  1.e-5 ) WRITE(numout,*) 'violation volume [m3/day]     (limthd) = ',(zchk_v_i * 86400.) 
    480             IF (    ABS( zchk_smv   ) >  1.e-4 ) WRITE(numout,*) 'violation saline [psu*m3/day] (limthd) = ',(zchk_smv * 86400.) 
    481             IF ( MINVAL( v_i(:,:,:) ) <  0.    ) WRITE(numout,*) 'violation v_i<0  [mm]         (limthd) = ',(MINVAL(v_i) * 1.e-3) 
    482             IF ( MAXVAL( SUM(a_i(:,:,:),dim=3) ) >  amax+epsi10 ) WRITE(numout,*) 'violation a_i>amax    (limthd) = ',MAXVAL(SUM(a_i,dim=3)) 
     484            IF ( ABS( zchk_v_i   ) >  1.e-5 ) WRITE(numout,*) 'violation volume [m3/day]     (limthd) = ',(zchk_v_i * 86400.) 
     485            IF ( ABS( zchk_smv   ) >  1.e-4 ) WRITE(numout,*) 'violation saline [psu*m3/day] (limthd) = ',(zchk_smv * 86400.) 
     486            IF ( zchk_vmin <  0.            ) WRITE(numout,*) 'violation v_i<0  [mm]         (limthd) = ',(zchk_vmin * 1.e-3) 
     487            IF ( zchk_amax >  amax+epsi10   ) WRITE(numout,*) 'violation a_i>amax            (limthd) = ',zchk_amax 
     488            IF ( zchk_amin <  0.            ) WRITE(numout,*) 'violation a_i<0               (limthd) = ',zchk_amin 
    483489         ENDIF 
    484490      ENDIF 
  • branches/2013/dev_r3406_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90

    r3938 r3963  
    121121 
    122122      ! mass and salt flux (clem) 
    123       REAL(wp) :: zdvres 
     123      REAL(wp) :: zdvres, zdvsur, zdvbot 
    124124      REAL(wp), POINTER, DIMENSION(:) ::   zviold, zvsold   ! old ice volume... 
    125125 
     
    255255         zhn            =  1.0 - MAX( zzero , SIGN( zone , - zhsnew ) ) 
    256256         ht_s_b(ji)     =  MAX( zzero , zhsnew ) 
     257         ! we recompute dh_s_tot (clem)  
     258         dh_s_tot (ji)  =  ht_s_b(ji) - zhsold(ji) 
    257259         ! Volume and mass variations of snow 
    258260         ! dvsbq_1d  (ji) =  a_i_b(ji) * ( ht_s_b(ji) - zhsold(ji) - zdh_s_mel(ji) ) 
     
    283285            zdq_i    (ji)    = zdq_i(ji) + zdeltah(ji,jk) * q_i_b(ji,jk) / rdt_ice 
    284286            ! 
    285             ! IOVINO 
    286             !zfsalt_melt(ji) = zfsalt_melt(ji) - sm_i_b(ji) * a_i_b(ji)    & 
    287             !   &                              * MIN( zdeltah(ji,jk) , 0.e0 ) * rhoic / rdt_ice 
     287            ! clem 
    288288            fseqv_1d(ji) = fseqv_1d(ji) - sm_i_b(ji) * a_i_b(ji)    & 
    289289               &                              * MIN( zdeltah(ji,jk) , 0._wp ) * rhoic / rdt_ice 
     290!            fseqv_1d(ji) = fseqv_1d(ji) - sm_i_b(ji) * a_i_b(ji) * zdeltah(ji,jk) * rhoic / rdt_ice 
    290291         END DO 
    291292      END DO 
     
    449450               ! Salinity update 
    450451               ! entrapment during bottom growth 
    451                dsm_i_se_1d(ji) = ( s_i_new(ji) * dh_i_bott(ji) + sm_i_b(ji) * ht_i_b(ji) )    & 
    452                   &            / MAX( ht_i_b(ji) + dh_i_bott(ji) ,epsi13 ) - sm_i_b(ji) 
     452               !clem:movedown dsm_i_se_1d(ji) = ( s_i_new(ji) * dh_i_bott(ji) + sm_i_b(ji) * ht_i_b(ji) )    & 
     453               !clem:movedown   &            / MAX( ht_i_b(ji) + dh_i_bott(ji) ,epsi13 ) - sm_i_b(ji) 
    453454               fseqv_1d(ji) = fseqv_1d(ji) - s_i_new(ji) * a_i_b(ji) * dh_i_bott(ji) * rhoic / rdt_ice 
    454455            ENDIF ! heat budget 
     
    489490                  zdq_i(ji)       = zdq_i(ji) + zdeltah(ji,jk) * q_i_b(ji,jk) / rdt_ice 
    490491               ENDIF 
    491                ! IOVINO contribution to salt flux 
    492                !zfsalt_melt(ji) = zfsalt_melt(ji) - sm_i_b(ji) * a_i_b(ji)   & 
    493                !     &                              * MIN( zdeltah(ji,jk) , 0._wp ) * rhoic / rdt_ice 
     492               ! clem 
    494493               fseqv_1d(ji) = fseqv_1d(ji) - sm_i_b(ji) * a_i_b(ji)    & 
    495494                    &                              * MIN( zdeltah(ji,jk) , 0._wp ) * rhoic / rdt_ice 
     495!               fseqv_1d(ji) = fseqv_1d(ji) - sm_i_b(ji) * a_i_b(ji) * zdeltah(ji,jk) * rhoic / rdt_ice 
    496496            ENDIF 
    497497         END DO ! ji 
     
    549549         !                     ! excessive energy is sent to lateral ablation 
    550550         fsup     (ji) =  rhoic * lfus * at_i_b(ji) / MAX( 1.0 - at_i_b(ji) , epsi13 ) * zdvres / rdt_ice 
    551          !                     !since ice volume is only used for outputs, we keep it global for all categories 
    552          dvbbq_1d (ji) = a_i_b(ji) * dh_i_bott(ji) 
    553          !                     !new ice thickness 
    554          zhgnew   (ji) = ht_i_b(ji) + dh_i_surf(ji) + dh_i_bott(ji) 
    555551      END DO 
    556552 
     
    564560         ! Adapt the remaining energy if too much ice melts 
    565561         !-------------------------------------------------- 
    566          zdvres     = MAX( 0._wp, - zhgnew(ji) ) 
    567          zihgnew    =  1.0 - MAX( zzero , SIGN( zone , - zhgnew(ji) ) ) !1 if ice 
    568          dh_i_bott (ji) = dh_i_bott(ji) + zdvres ! clem@bug 
    569          zhgnew    (ji) =         zihgnew   * zhgnew(ji)      ! ice thickness is put to 0 
    570          ! remaining heat 
     562         zdvres     = MAX( 0._wp, - ht_i_b(ji) - dh_i_surf(ji) - dh_i_bott(ji) ) 
     563         zdvsur     = MIN( 0._wp, dh_i_surf(ji) + zdvres ) - dh_i_surf(ji) ! fill the surface first 
     564         zdvbot     = MAX( 0._wp, zdvres - zdvsur ) ! then the bottom 
     565         dh_i_surf (ji) = dh_i_surf(ji) + zdvsur ! clem 
     566         dh_i_bott (ji) = dh_i_bott(ji) + zdvbot ! clem 
     567 
     568         ! new ice thickness (clem) 
     569         zhgnew(ji) = ht_i_b(ji) + dh_i_surf(ji) + dh_i_bott(ji) 
     570         zihgnew    = 1.0 - MAX( zzero , SIGN( zone , - zhgnew(ji) ) ) !1 if ice 
     571         zhgnew(ji) = zihgnew * zhgnew(ji)      ! ice thickness is put to 0 
     572  
     573         !                     !since ice volume is only used for outputs, we keep it global for all categories 
     574         dvbbq_1d (ji) = a_i_b(ji) * dh_i_bott(ji) 
     575 
     576        ! remaining heat 
    571577         zfdt_final(ji) = ( 1.0 - zihgnew ) * ( zqfont_su(ji) +  zqfont_bo(ji) )  
    572578 
     
    582588         ht_s_b(ji)     =  MAX( zzero , zhnfi ) 
    583589         zqt_s(ji)      =  zqt_s(ji) * ht_s_b(ji) 
     590         ! we recompute dh_s_tot (clem) 
     591         dh_s_tot (ji)  =  ht_s_b(ji) - zhsold(ji) 
    584592 
    585593         ! Mass variations of ice and snow 
     
    598606         !--------------------------------- 
    599607         focea(ji)  = - zfdt_final(ji) / rdt_ice         ! focea is in W.m-2 * dt 
    600          ! salt flux 
    601          fseqv_1d(ji)  = fseqv_1d(ji) - sm_i_b(ji) * a_i_b(ji) * zdvres * rhoic / rdt_ice 
     608 
     609         ! residual salt flux (clem) 
     610         !-------------------------- 
     611         ! surface 
     612         fseqv_1d(ji)  = fseqv_1d(ji) - sm_i_b(ji) * a_i_b(ji) * zdvsur * rhoic / rdt_ice 
     613         ! bottom 
     614         IF     ( ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) >= 0._wp ) THEN ! melting 
     615            fseqv_1d(ji)  = fseqv_1d(ji) - sm_i_b(ji) * a_i_b(ji) * zdvbot * rhoic / rdt_ice 
     616         ELSE ! growth 
     617            fseqv_1d(ji)  = fseqv_1d(ji) - s_i_new(ji) * a_i_b(ji) * zdvbot * rhoic / rdt_ice 
     618         ENDIF 
     619         ! 
    602620         !                     ! diagnostic ( bottom ice growth ) 
    603621         zji = MOD( npb(ji) - 1, jpi ) + 1 
     
    700718         !ENDIF 
    701719         ! entrapment during snow ice formation 
    702          i_ice_switch = 1.0 - MAX( 0.e0 , SIGN( 1.0 , - ht_i_b(ji) + epsi13 ) ) 
    703          isnowic      = 1.0 - MAX( 0.e0 , SIGN( 1.0 , - dh_snowice(ji) ) ) * i_ice_switch 
    704720 
    705721         !clem IF(  num_sal == 2  .OR.  num_sal == 4  )   & 
     
    707723         !clem   &               + sm_i_b(ji) * ht_i_b(ji) / MAX( ht_i_b(ji) + dh_snowice(ji), epsi13)   & 
    708724         !clem   &               - sm_i_b(ji) ) * isnowic 
    709          IF (  num_sal == 2  .OR.  num_sal == 4  )   & 
    710             & dsm_i_si_1d(ji) = ( ( zsm_snowice * dh_snowice(ji) + sm_i_b(ji) * ht_i_b(ji) ) & 
    711             &                    / MAX( ht_i_b(ji) + dh_snowice(ji), epsi13 ) - sm_i_b(ji) ) * isnowic      
     725 
     726         ! clem: new salinity difference stored (to be used in limthd_ent.F90) 
     727         IF (  num_sal == 2  .OR.  num_sal == 4  ) THEN 
     728            i_ice_switch = 1.0 - MAX( 0.e0 , SIGN( 1.0 , - zhgnew(ji) + epsi13 ) ) 
     729            ! salinity dif due to snow-ice formation 
     730            dsm_i_si_1d(ji) = ( zsm_snowice - sm_i_b(ji) ) * dh_snowice(ji) / MAX( zhgnew(ji), epsi13 ) * i_ice_switch      
     731            ! salinity dif due to bottom growth  
     732            IF ( ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) < 0._wp ) THEN 
     733               dsm_i_se_1d(ji) = ( s_i_new(ji) - sm_i_b(ji) ) * dh_i_bott(ji) / MAX( zhgnew(ji), epsi13 ) * i_ice_switch 
     734            ENDIF 
     735         ENDIF 
    712736 
    713737         !  Actualize new snow and ice thickness. 
  • branches/2013/dev_r3406_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90

    r3938 r3963  
    8181      INTEGER ::   layer, nbpac     ! local integers  
    8282      INTEGER ::   zji, zjj, iter   !   -       - 
    83       REAL(wp)  ::   ztmelts, zdv, zqold, zfrazb, zweight, zalphai, zindb, zde  ! local scalars 
     83      REAL(wp)  ::   ztmelts, zdv, zqold, zfrazb, zweight, zalphai, zindb, zinda, zde  ! local scalars 
    8484      REAL(wp) ::   zgamafr, zvfrx, zvgx, ztaux, ztwogp, zf , zhicol_new        !   -      - 
    8585      REAL(wp) ::   ztenagm, zvfry, zvgy, ztauy, zvrel2, zfp, zsqcd , zhicrit   !   -      - 
     
    454454            zji                  = MOD( npac(ji) - 1, jpi ) + 1 
    455455            zjj                  = ( npac(ji) - 1 ) / jpi + 1 
    456             diag_lat_gr(zji,zjj) = zv_newice(ji) / rdt_ice 
     456            diag_lat_gr(zji,zjj) = diag_lat_gr(zji,zjj) + zv_newice(ji) / rdt_ice ! clem 
    457457          END DO !ji 
    458458 
     
    547547            DO ji = 1, nbpac 
    548548               zindb = MAX( 0._wp, SIGN( 1._wp , zdv_res(ji) ) ) 
    549                zv_i_ac(ji,jl) = zv_i_ac(ji,jl) + zindb * zdv_res(ji) * za_i_ac(ji,jl) / MAX( zat_i_lev(ji) , epsi06 ) 
     549               zinda = MAX( 0._wp, SIGN( 1._wp , zat_i_lev(ji) - epsi06 ) )  ! clem 
     550               zv_i_ac(ji,jl) = zv_i_ac(ji,jl) + zindb * zinda * zdv_res(ji) * za_i_ac(ji,jl) / MAX( zat_i_lev(ji) , epsi06 ) 
    550551            END DO 
    551552         END DO 
     
    623624         ! Update salinity 
    624625         !----------------- 
    625          IF(  num_sal == 2  .OR.  num_sal == 4  ) THEN 
     626         !clem IF(  num_sal == 2  .OR.  num_sal == 4  ) THEN 
    626627            DO jl = 1, jpl 
    627628               DO ji = 1, nbpac 
    628629                  zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - zv_i_ac(ji,jl) + epsi10 ) )  ! 0 if no ice and 1 if yes 
    629630                  zdv   = zv_i_ac(ji,jl) - zv_old(ji,jl) 
    630                   zsmv_i_ac(ji,jl) = ( zsmv_i_ac(ji,jl) + zdv * zs_newice(ji) ) * zindb 
     631                  zsmv_i_ac(ji,jl) = zsmv_i_ac(ji,jl) + zdv * zs_newice(ji) * zindb ! clem modif 
    631632               END DO 
    632633            END DO    
    633          ENDIF 
     634         !clem ENDIF 
    634635 
    635636         !-------------------------------- 
     
    640641               zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - zv_i_ac(ji,jl) + epsi10 ) )  ! 0 if no ice and 1 if yes 
    641642               zdv   = zv_i_ac(ji,jl) - zv_old(ji,jl) 
    642                rdmicif_1d(ji) = rdmicif_1d(ji) + zdv * rhoic !* zindb 
     643               rdmicif_1d(ji) = rdmicif_1d(ji) + zdv * rhoic * zindb 
    643644               fseqv_1d(ji)   =   fseqv_1d(ji) - zdv * rhoic * zs_newice(ji) / rdt_ice * zindb 
    644645           END DO 
     
    652653            CALL tab_1d_2d( nbpac, v_i (:,:,jl), npac(1:nbpac), zv_i_ac (1:nbpac,jl), jpi, jpj ) 
    653654            CALL tab_1d_2d( nbpac, oa_i(:,:,jl), npac(1:nbpac), zoa_i_ac(1:nbpac,jl), jpi, jpj ) 
    654             IF (  num_sal == 2  .OR.  num_sal == 4  )   & 
     655            !clem IF (  num_sal == 2  .OR.  num_sal == 4  )   & 
    655656               CALL tab_1d_2d( nbpac, smv_i (:,:,jl), npac(1:nbpac), zsmv_i_ac(1:nbpac,jl) , jpi, jpj ) 
    656657            DO jk = 1, nlay_i 
  • branches/2013/dev_r3406_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limthd_sal.F90

    r3938 r3963  
    127127            i_ice_switch = 1._wp - MAX ( 0._wp, SIGN( 1._wp , - ht_i_b(ji) ) ) 
    128128            sm_i_b(ji)   = i_ice_switch * sm_i_b(ji) + s_i_min * ( 1._wp - i_ice_switch ) 
    129          END DO ! ji 
     129 
     130            !---------------------------- 
     131            ! Heat flux - brine drainage 
     132            !---------------------------- 
     133            fhbri_1d(ji) = 0._wp 
     134 
     135            !---------------------------- 
     136            ! Salt flux - brine drainage 
     137            !---------------------------- 
     138            fsbri_1d(ji) = fsbri_1d(ji) - i_ice_switch * rhoic * a_i_b(ji) * ht_i_b(ji) * ( sm_i_b(ji) - zsiold(ji) ) / rdt_ice 
     139            IF( num_sal == 4 ) fsbri_1d(ji) = 0._wp 
     140 
     141        END DO ! ji 
    130142 
    131143         ! Salinity profile 
     
    136148         !---------------------------- 
    137149 
    138          DO ji = kideb, kiut 
     150!clem:move         DO ji = kideb, kiut 
    139151!!gm useless 
    140152            ! iflush  : 1 if summer  
    141             iflush  =  MAX( 0._wp , SIGN ( 1._wp , t_su_b(ji) - rtt ) )  
     153!clem            iflush  =  MAX( 0._wp , SIGN ( 1._wp , t_su_b(ji) - rtt ) )  
    142154            ! igravdr : 1 if t_su lt t_bo 
    143             igravdr =  MAX( 0._wp , SIGN ( 1._wp , t_bo_b(ji) - t_su_b(ji) ) )  
     155!clem            igravdr =  MAX( 0._wp , SIGN ( 1._wp , t_bo_b(ji) - t_su_b(ji) ) )  
    144156            ! iaccrbo : 1 if bottom accretion 
    145             iaccrbo =  MAX( 0._wp , SIGN ( 1._wp , dh_i_bott(ji) ) ) 
     157!clem            iaccrbo =  MAX( 0._wp , SIGN ( 1._wp , dh_i_bott(ji) ) ) 
    146158!!gm end useless 
    147159            ! 
    148             fhbri_1d(ji) = 0._wp 
    149          END DO ! ji 
     160!clem:move            fhbri_1d(ji) = 0._wp 
     161!clem:move         END DO ! ji 
    150162 
    151163         !---------------------------- 
    152164         ! Salt flux - brine drainage 
    153165         !---------------------------- 
    154           DO ji = kideb, kiut 
    155             i_ice_switch = 1._wp - MAX ( 0._wp, SIGN( 1._wp , - ht_i_b(ji) ) ) 
    156             fsbri_1d(ji) = fsbri_1d(ji) - i_ice_switch * rhoic * a_i_b(ji) * ht_i_b(ji) * ( sm_i_b(ji) - zsiold(ji) ) / rdt_ice 
     166!clem:move          DO ji = kideb, kiut 
     167!clem:move            i_ice_switch = 1._wp - MAX ( 0._wp, SIGN( 1._wp , - ht_i_b(ji) ) ) 
     168!clem:move            fsbri_1d(ji) = fsbri_1d(ji) - i_ice_switch * rhoic * a_i_b(ji) * ht_i_b(ji) * ( sm_i_b(ji) - zsiold(ji) ) / rdt_ice 
    157169            !i_ice_switch = 1._wp - MAX ( 0._wp, SIGN( 1._wp , - zhiold(ji) ) ) 
    158170            !fsbri_1d(ji) = fsbri_1d(ji) - i_ice_switch * rhoic * a_i_b(ji) * zhiold(ji) * ( sm_i_b(ji) - zsiold(ji) ) / rdt_ice 
    159171            !clem fsbri_1d(ji) = fsbri_1d(ji) - i_ice_switch * rhoic * a_i_b(ji) * ht_i_b(ji)         & 
    160172            !clem     &         * ( MAX(dsm_i_gd_1d(ji) + dsm_i_fl_1d(ji), sm_i_b(ji) - zsiold(ji) ) ) / rdt_ice 
    161             IF( num_sal == 4 ) fsbri_1d(ji) = 0._wp 
    162           END DO ! ji 
     173!clem:move            IF( num_sal == 4 ) fsbri_1d(ji) = 0._wp 
     174!clem:move          END DO ! ji 
    163175 
    164176         ! Only necessary for conservation check since salinity is modified 
  • branches/2013/dev_r3406_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90

    r3938 r3963  
    8181      REAL(wp), POINTER, DIMENSION(:,:,:,:)  ::   zs0e 
    8282      REAL(wp) :: zchk_v_i, zchk_smv, zchk_fs, zchk_fw, zchk_v_i_b, zchk_smv_b, zchk_fs_b, zchk_fw_b ! Check conservation (C Rousset) 
     83      REAL(wp) :: zchk_vmin, zchk_amin, zchk_amax ! Check errors (C Rousset) 
    8384      ! mass and salt flux (clem) 
    8485      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zviold   ! old ice volume... 
     
    170171         zusnit = 1.0 / REAL( initad )  
    171172         IF( zcfl > 0.5 .AND. lwp )   & 
    172             WRITE(numout,*) 'lim_trp_2 : CFL violation at day ', nday, ', cfl = ', zcfl,   & 
     173            WRITE(numout,*) 'lim_trp  : CFL violation at day ', nday, ', cfl = ', zcfl,   & 
    173174               &                        ': the ice time stepping is split in two' 
    174175 
     
    430431               DO ji = 1, jpi 
    431432 
    432                   IF ( v_i(ji,jj,jl) > 0. ) THEN 
     433                  IF ( v_i(ji,jj,jl) > 0._wp ) THEN 
    433434                     zvi = v_i(ji,jj,jl) 
    434435                     zvs = v_s(ji,jj,jl) 
     
    437438                      
    438439                     IF ( ( zdv > 0.0 .AND. ht_i(ji,jj,jl) > zhimax(ji,jj,jl) .AND. SUM( zaiold(ji,jj,1:jpl) ) < 0.80 ) .OR. & 
    439                           ( zdv < 0.0 .AND. ht_i(ji,jj,jl) > zhimax(ji,jj,jl) ) ) THEN                                           
     440                        & ( zdv < 0.0 .AND. ht_i(ji,jj,jl) > zhimax(ji,jj,jl) ) ) THEN                                           
    440441                        ht_i(ji,jj,jl) = MIN( zhimax(ji,jj,jl), hi_max(jl) ) 
    441442                        zindh   =  MAX( rzero, SIGN( rone, ht_i(ji,jj,jl) - epsi10 ) ) 
     
    447448                     ENDIF 
    448449 
    449                      !zindh   =  MAX( rzero, SIGN( rone, ht_i(ji,jj,jl) - epsi10 ) ) 
    450                      v_i(ji,jj,jl) = a_i(ji,jj,jl) * ht_i(ji,jj,jl) 
    451                      v_s(ji,jj,jl) = a_i(ji,jj,jl) * ht_s(ji,jj,jl) 
    452  
    453                      ! Update mass fluxes (clem) 
     450                     ! small correction due to *zindh for a_i 
     451                     v_i(ji,jj,jl) = zindh * v_i(ji,jj,jl) 
     452                     v_s(ji,jj,jl) = zindh * v_s(ji,jj,jl) 
     453 
     454                     ! Update mass fluxes 
    454455                     rdmicif(ji,jj) = rdmicif(ji,jj) + ( v_i(ji,jj,jl) - zvi ) * rhoic 
    455456                     rdmsnif(ji,jj) = rdmsnif(ji,jj) + ( v_s(ji,jj,jl) - zvs ) * rhosn 
     
    587588         zchk_smv = ( glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zchk_smv_b ) / rdt_ice + ( zchk_fs / rhoic ) 
    588589 
     590         zchk_vmin = glob_min(v_i) 
     591         zchk_amax = glob_max(SUM(a_i,dim=3)) 
     592         zchk_amin = glob_min(a_i) 
     593 
    589594         IF(lwp) THEN 
    590             IF (    ABS( zchk_v_i   ) >  1.e-5 ) WRITE(numout,*) 'violation volume [m3/day]     (limtrp) = ',(zchk_v_i * 86400.) 
    591             IF (    ABS( zchk_smv   ) >  1.e-4 ) WRITE(numout,*) 'violation saline [psu*m3/day] (limtrp) = ',(zchk_smv * 86400.) 
    592             IF ( MINVAL( v_i(:,:,:) ) <  0.    ) WRITE(numout,*) 'violation v_i<0  [mm]         (limtrp) = ',(MINVAL(v_i) * 1.e-3) 
    593             IF ( MINVAL( a_i(:,:,:) ) <  0.    ) WRITE(numout,*) 'violation a_i<0               (limtrp) = ',MINVAL(a_i) 
     595            IF ( ABS( zchk_v_i   ) >  1.e-5 ) WRITE(numout,*) 'violation volume [m3/day]     (limtrp) = ',(zchk_v_i * 86400.) 
     596            IF ( ABS( zchk_smv   ) >  1.e-4 ) WRITE(numout,*) 'violation saline [psu*m3/day] (limtrp) = ',(zchk_smv * 86400.) 
     597            IF ( zchk_vmin <  0.            ) WRITE(numout,*) 'violation v_i<0  [mm]         (limtrp) = ',(zchk_vmin * 1.e-3) 
     598            IF ( zchk_amin <  0.            ) WRITE(numout,*) 'violation a_i<0               (limtrp) = ',zchk_amin 
    594599         ENDIF 
    595600      ENDIF 
  • branches/2013/dev_r3406_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limupdate1.F90

    r3938 r3963  
    8686      INTEGER ::   i_ice_switch 
    8787      INTEGER ::   ind_im, layer      ! indices for internal melt 
    88       REAL(wp) ::   zweight, zesum, z_da_i 
    89       REAL(wp) ::   zindb, zindsn, zindic 
     88      REAL(wp) ::   zweight, zesum, z_da_i, zhimax 
     89      REAL(wp) ::   zinda, zindb, zindsn, zindic 
    9090      REAL(wp) ::   zindg, zh, zdvres, zviold2 
    9191      REAL(wp) ::   zbigvalue, zvsold2, z_da_ex 
     
    9494      REAL(wp), POINTER, DIMENSION(:) ::   zthick0, zqm0      ! thickness of the layers and heat contents for 
    9595      REAL(wp) :: zchk_v_i, zchk_smv, zchk_fs, zchk_fw, zchk_v_i_b, zchk_smv_b, zchk_fs_b, zchk_fw_b ! Check conservation (C Rousset) 
     96      REAL(wp) :: zchk_vmin, zchk_amin, zchk_amax ! Check errors (C Rousset) 
    9697      ! mass and salt flux (clem) 
    9798      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zviold, zvsold, zsmvold   ! old ice volume... 
     
    182183      ! 2. Review of all pathological cases 
    183184      !-------------------------------------- 
     185 
     186      !------------------------------------------- 
     187      ! 2.1) Advection of ice in an ice-free cell 
     188      !------------------------------------------- 
     189      ! should be removed since it is treated after dynamics now 
     190 
     191!      !IF( ln_nicep ) THEN   
     192!         WRITE(numout,*) ' limupdate1 - before h correction ' 
     193!         WRITE(numout,*) ' a_i  : ', a_i(12, 44, 1:jpl) 
     194!         WRITE(numout,*) ' at_i : ', at_i(12,44) 
     195!         WRITE(numout,*) ' v_i  : ', v_i(12, 44, 1:jpl) 
     196!         WRITE(numout,*) ' ht_i : ', ht_i(12, 44, 1:jpl) 
     197!      !ENDIF 
     198! 
     199      zhimax = 1._wp 
     200      ! first category 
     201      DO jj = 1, jpj 
     202         DO ji = 1, jpi 
     203            !--- the thickness of such an ice is often out of bounds 
     204            !--- thus we recompute a new area while conserving ice volume 
     205            zat_i_old = SUM( old_a_i(ji,jj,:) ) 
     206            zindb     = MAX( 0._wp, SIGN( 1._wp, ABS( d_a_i_trp(ji,jj,1)) - epsi10 ) )  
     207            IF ( ( ABS(d_v_i_trp(ji,jj,1))/MAX(ABS(d_a_i_trp(ji,jj,1)),epsi10)*zindb .GT. zhimax) & 
     208                 .AND.( ( v_i(ji,jj,1)/MAX(a_i(ji,jj,1),epsi10)*zindb) .GT. zhimax ) & 
     209                 .AND.( zat_i_old .LT. 1.e-6 ) )  THEN ! new line 
     210               ht_i(ji,jj,1) = hi_max(1) / 2.0 
     211               a_i (ji,jj,1) = v_i(ji,jj,1) / ht_i(ji,jj,1) 
     212            ENDIF 
     213         END DO 
     214      END DO 
     215 
     216 
     217!      !IF( ln_nicep ) THEN   
     218!         at_i(:,:) = 0._wp 
     219!         DO jl = 1, jpl 
     220!            at_i(:,:) = a_i(:,:,jl) + at_i(:,:) 
     221!         END DO 
     222!         WRITE(numout,*) ' limupdate1 - after h correction 1 ' 
     223!         WRITE(numout,*) ' a_i  : ', a_i(12, 44, 1:jpl) 
     224!         WRITE(numout,*) ' at_i : ', at_i(12,44) 
     225!         WRITE(numout,*) ' v_i  : ', v_i(12, 44, 1:jpl) 
     226!         WRITE(numout,*) ' ht_i : ', ht_i(12, 44, 1:jpl) 
     227!      !ENDIF 
     228 
     229      zhimax = 10._wp 
     230      ! other categories 
     231      DO jl = 2, jpl 
     232         jm = ice_types(jl) 
     233         DO jj = 1, jpj 
     234            DO ji = 1, jpi 
     235               zindb =  MAX( rzero, SIGN( rone, ABS(d_a_i_trp(ji,jj,jl)) - epsi10 ) )  
     236               ! this correction is very tricky... sometimes, advection gets wrong i don't know why 
     237               ! it makes problems when the advected volume and concentration do not seem to be  
     238               ! related with each other 
     239               ! the new thickness is sometimes very big! 
     240               ! and sometimes d_a_i_trp and d_v_i_trp have different sign 
     241               ! which of course is plausible 
     242               ! but fuck! it fucks everything up :) 
     243               IF ( (ABS(d_v_i_trp(ji,jj,jl))/MAX(ABS(d_a_i_trp(ji,jj,jl)),epsi10)*zindb .GT. zhimax) & 
     244                    .AND.(v_i(ji,jj,jl)/MAX(a_i(ji,jj,jl),epsi10)*zindb) .GT. zhimax ) THEN 
     245                  ht_i(ji,jj,jl) = ( hi_max_typ(jl-ice_cat_bounds(jm,1),jm) + hi_max_typ(jl-ice_cat_bounds(jm,1)+1,jm) ) / 2.0 
     246                  a_i (ji,jj,jl) = v_i(ji,jj,jl) / ht_i(ji,jj,jl) 
     247               ENDIF 
     248            END DO ! ji 
     249         END DO !jj 
     250      END DO !jl 
    184251      
    185252      at_i(:,:) = 0._wp 
     
    187254         at_i(:,:) = a_i(:,:,jl) + at_i(:,:) 
    188255      END DO 
     256 
     257!      !IF( ln_nicep ) THEN   
     258!         WRITE(numout,*) ' limupdate1 - after h correction 2 ' 
     259!         WRITE(numout,*) ' a_i  : ', a_i(12, 44, 1:jpl) 
     260!         WRITE(numout,*) ' at_i : ', at_i(12,44) 
     261!         WRITE(numout,*) ' v_i  : ', v_i(12, 44, 1:jpl) 
     262!         WRITE(numout,*) ' ht_i : ', ht_i(12, 44, 1:jpl) 
     263!      !ENDIF 
    189264 
    190265      !---------------------------------------------------- 
     
    247322               !---------------------------------------------------------------------------- 
    248323               zindg          = tms(ji,jj) *  MAX( 0._wp, SIGN( 1._wp, -v_i(ji,jj,jl) ) ) 
    249                v_i(ji,jj,jl)  = v_i(ji,jj,jl) + zindg * rhosn * v_s(ji,jj,jl) / rau0 
    250                v_s(ji,jj,jl)  = ( 1._wp - zindg ) * v_s(ji,jj,jl) + zindg * v_i(ji,jj,jl) * ( rau0 - rhoic ) / rhosn 
     324               zdvres         = zindg * rhosn * v_s(ji,jj,jl) / rau0 
     325               v_i(ji,jj,jl)  = v_i(ji,jj,jl)  + zdvres 
     326 
     327               zdvres         = zindsn*zindb * ( - zindg * v_s(ji,jj,jl) + zindg * v_i(ji,jj,jl) * ( rau0 - rhoic ) / rhosn ) 
     328               v_s(ji,jj,jl)  = v_s(ji,jj,jl)  + zdvres 
    251329 
    252330               !--- 2.7 Correction to ice concentrations  
     
    286364 
    287365      !--- 2.13 ice concentration should not exceed amax  
     366      !         (it should not be the case)  
    288367      !----------------------------------------------------- 
    289368      DO jj = 1, jpj 
    290369         DO ji = 1, jpi 
    291370            z_da_ex =  MAX( at_i(ji,jj) - amax , 0.0 )         
     371            zindb   =  MAX( rzero, SIGN( rone, at_i(ji,jj) - epsi06 ) )  
    292372            DO jl  = 1, jpl 
    293                zindb   =  MAX( rzero, SIGN( rone, v_i(ji,jj,jl) - epsi04 ) )  
    294373               z_da_i = a_i(ji,jj,jl) * z_da_ex / MAX( at_i(ji,jj), epsi06 ) * zindb 
    295                a_i(ji,jj,jl) = a_i(ji,jj,jl) - z_da_i 
    296            END DO 
    297          END DO  
    298       END DO  
     374               a_i(ji,jj,jl) = MAX( 0._wp, a_i(ji,jj,jl) - z_da_i ) 
     375               ! 
     376               zinda   =  MAX( rzero, SIGN( rone, a_i(ji,jj,jl) - epsi06 ) )  
     377               ht_i(ji,jj,jl) = v_i(ji,jj,jl) / MAX( a_i(ji,jj,jl), epsi06 ) * zinda 
     378               !v_i(ji,jj,jl) = ht_i(ji,jj,jl) * a_i(ji,jj,jl) ! makes ice shrinken but should not be used 
     379            END DO 
     380         END DO 
     381      END DO 
    299382      at_i(:,:) = a_i(:,:,1) 
    300383      DO jl = 2, jpl 
    301384         at_i(:,:) = a_i(:,:,jl) + at_i(:,:) 
    302385      END DO 
     386 
    303387 
    304388      ! Final thickness distribution rebinning 
     
    311395         ENDIF 
    312396      END DO 
     397 
    313398 
    314399      !--------------------- 
     
    356441      !- check conservation (C Rousset) 
    357442      IF (ln_limdiahsb) THEN 
     443 
    358444         zchk_fs  = glob_sum( ( fsbri(:,:) + fseqv(:,:) + fsalt_res(:,:) + fsalt_rpo(:,:) ) * area(:,:) * tms(:,:) ) - zchk_fs_b 
    359445         zchk_fw  = glob_sum( rdmicif(:,:) * area(:,:) * tms(:,:) ) - zchk_fw_b 
     
    362448         zchk_smv = ( glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zchk_smv_b ) / rdt_ice + ( zchk_fs / rhoic ) 
    363449 
     450         zchk_vmin = glob_min(v_i) 
     451         zchk_amax = glob_max(SUM(a_i,dim=3)) 
     452         zchk_amin = glob_min(a_i) 
     453 
    364454         IF(lwp) THEN 
    365             IF (    ABS( zchk_v_i   ) >  1.e-5 ) WRITE(numout,*) 'violation volume [m3/day]     (limupdate1) = ',(zchk_v_i * 86400.) 
    366             IF (    ABS( zchk_smv   ) >  1.e-4 ) WRITE(numout,*) 'violation saline [psu*m3/day] (limupdate1) = ',(zchk_smv * 86400.) 
    367             IF ( MINVAL( v_i(:,:,:) ) <  0.    ) WRITE(numout,*) 'violation v_i<0  [mm]         (limupdate1) = ',(MINVAL(v_i) * 1.e-3) 
    368             IF ( MAXVAL( SUM(a_i(:,:,:),dim=3) ) >  amax+epsi10 ) WRITE(numout,*) 'violation a_i>amax    (limupdate1) = ',MAXVAL(SUM(a_i,dim=3)) 
    369             IF ( MINVAL( a_i(:,:,:) ) <  0.    ) WRITE(numout,*) 'violation a_i<0               (limupdate1) = ',MINVAL(a_i) 
     455            IF ( ABS( zchk_v_i   ) >  1.e-5 ) WRITE(numout,*) 'violation volume [m3/day]     (limupdate1) = ',(zchk_v_i * 86400.) 
     456            IF ( ABS( zchk_smv   ) >  1.e-4 ) WRITE(numout,*) 'violation saline [psu*m3/day] (limupdate1) = ',(zchk_smv * 86400.) 
     457            IF ( zchk_vmin <  0.            ) WRITE(numout,*) 'violation v_i<0  [mm]         (limupdate1) = ',(zchk_vmin * 1.e-3) 
     458            IF ( zchk_amax >  amax+epsi10   ) WRITE(numout,*) 'violation a_i>amax            (limupdate1) = ',zchk_amax 
     459            IF ( zchk_amin <  0.            ) WRITE(numout,*) 'violation a_i<0               (limupdate1) = ',zchk_amin 
    370460         ENDIF 
    371461      ENDIF 
  • branches/2013/dev_r3406_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limupdate2.F90

    r3938 r3963  
    8383      INTEGER ::   ind_im, layer      ! indices for internal melt 
    8484      REAL(wp) ::   zweight, zesum, zhimax, z_da_i 
    85       REAL(wp) ::   zindb, zindsn, zindic 
     85      REAL(wp) ::   zinda, zindb, zindsn, zindic 
    8686      REAL(wp) ::   zindg, zh, zdvres, zviold2 
    8787      REAL(wp) ::   zbigvalue, zvsold2, z_da_ex 
     
    9191      REAL(wp), POINTER, DIMENSION(:) ::   zthick0, zqm0      ! thickness of the layers and heat contents for 
    9292      REAL(wp) :: zchk_v_i, zchk_smv, zchk_fs, zchk_fw, zchk_v_i_b, zchk_smv_b, zchk_fs_b, zchk_fw_b ! Check conservation (C Rousset) 
    93  
     93      REAL(wp) :: zchk_vmin, zchk_amin, zchk_amax ! Check errors (C Rousset) 
    9494      ! mass and salt flux (clem) 
    9595      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zviold, zvsold, zsmvold   ! old ice volume... 
     
    177177      END DO 
    178178 
    179  
    180179      !-------------------------------------- 
    181180      ! 2. Review of all pathological cases 
    182181      !-------------------------------------- 
    183  
     182      !------------------------------------------- 
     183      ! 2.1) Advection of ice in an ice-free cell 
     184      !------------------------------------------- 
     185      ! should be removed since it is treated after dynamics now 
     186 
     187!      !IF( ln_nicep ) THEN   
     188!         WRITE(numout,*) ' limupdate2 - before h correction ' 
     189!         WRITE(numout,*) ' a_i  : ', a_i(12, 44, 1:jpl) 
     190!         WRITE(numout,*) ' at_i : ', at_i(12,44) 
     191!         WRITE(numout,*) ' v_i  : ', v_i(12, 44, 1:jpl) 
     192!         WRITE(numout,*) ' ht_i : ', ht_i(12, 44, 1:jpl) 
     193!      !ENDIF 
     194 
     195      zhimax = 1._wp 
     196      ! first category 
     197      DO jj = 1, jpj 
     198         DO ji = 1, jpi 
     199            !--- the thickness of such an ice is often out of bounds 
     200            !--- thus we recompute a new area while conserving ice volume 
     201            zat_i_old = SUM( old_a_i(ji,jj,:) ) 
     202            zindb     = MAX( 0._wp, SIGN( 1._wp, ABS( d_a_i_thd(ji,jj,1)) - epsi10 ) )  
     203            IF ( ( ABS(d_v_i_thd(ji,jj,1))/MAX(ABS(d_a_i_thd(ji,jj,1)),epsi10)*zindb .GT. zhimax) & 
     204                 .AND.( ( v_i(ji,jj,1)/MAX(a_i(ji,jj,1),epsi10)*zindb) .GT. zhimax ) & 
     205                 .AND.( zat_i_old .LT. 1.e-6 ) )  THEN ! new line 
     206               ht_i(ji,jj,1) = hi_max(1) / 2.0 
     207               a_i (ji,jj,1) = v_i(ji,jj,1) / ht_i(ji,jj,1) 
     208            ENDIF 
     209         END DO 
     210      END DO 
     211 
     212!      !IF( ln_nicep ) THEN   
     213!         at_i(:,:) = 0._wp 
     214!         DO jl = 1, jpl 
     215!            at_i(:,:) = a_i(:,:,jl) + at_i(:,:) 
     216!         END DO 
     217!         WRITE(numout,*) ' limupdate2 - after h correction 1 ' 
     218!         WRITE(numout,*) ' a_i  : ', a_i(12, 44, 1:jpl) 
     219!         WRITE(numout,*) ' at_i : ', at_i(12,44) 
     220!         WRITE(numout,*) ' v_i  : ', v_i(12, 44, 1:jpl) 
     221!         WRITE(numout,*) ' ht_i : ', ht_i(12, 44, 1:jpl) 
     222!      !ENDIF 
     223!  
     224      zhimax = 10._wp 
     225      ! other categories 
     226      DO jl = 2, jpl 
     227         jm = ice_types(jl) 
     228         DO jj = 1, jpj 
     229            DO ji = 1, jpi 
     230               zindb =  MAX( rzero, SIGN( rone, ABS(d_a_i_thd(ji,jj,jl)) - epsi10 ) )  
     231              ! this correction is very tricky... sometimes, advection gets wrong i don't know why 
     232               ! it makes problems when the advected volume and concentration do not seem to be  
     233               ! related with each other 
     234               ! the new thickness is sometimes very big! 
     235               ! and sometimes d_a_i_trp and d_v_i_trp have different sign 
     236               ! which of course is plausible 
     237               ! but fuck! it fucks everything up :) 
     238               IF ( (ABS(d_v_i_thd(ji,jj,jl))/MAX(ABS(d_a_i_thd(ji,jj,jl)),epsi10)*zindb .GT. zhimax) & 
     239                    .AND.(v_i(ji,jj,jl)/MAX(a_i(ji,jj,jl),epsi10)*zindb) .GT. zhimax ) THEN 
     240                  ht_i(ji,jj,jl) = ( hi_max_typ(jl-ice_cat_bounds(jm,1),jm) + hi_max_typ(jl-ice_cat_bounds(jm,1)+1,jm) ) / 2.0 
     241                  a_i (ji,jj,jl) = v_i(ji,jj,jl) / ht_i(ji,jj,jl) 
     242               ENDIF 
     243            END DO ! ji 
     244         END DO !jj 
     245      END DO !jl 
     246      
     247      at_i(:,:) = 0._wp 
     248      DO jl = 1, jpl 
     249         at_i(:,:) = a_i(:,:,jl) + at_i(:,:) 
     250      END DO 
     251 
     252!      !IF( ln_nicep ) THEN   
     253!         WRITE(numout,*) ' limupdate2 - after h correction 2 ' 
     254!         WRITE(numout,*) ' a_i  : ', a_i(12, 44, 1:jpl) 
     255!         WRITE(numout,*) ' at_i : ', at_i(12,44) 
     256!         WRITE(numout,*) ' v_i  : ', v_i(12, 44, 1:jpl) 
     257!         WRITE(numout,*) ' ht_i : ', ht_i(12, 44, 1:jpl) 
     258!      !ENDIF 
    184259 
    185260      !---------------------------------------------------- 
     
    363438               v_i(ji,jj,jl)  = v_i(ji,jj,jl)  + zdvres 
    364439 
    365                zdvres         = - zindg * v_s(ji,jj,jl) + zindg * v_i(ji,jj,jl) * ( rau0 - rhoic ) / rhosn 
     440               zdvres         = zindsn*zindb * ( - zindg * v_s(ji,jj,jl) + zindg * v_i(ji,jj,jl) * ( rau0 - rhoic ) / rhosn ) 
    366441               v_s(ji,jj,jl)  = v_s(ji,jj,jl)  + zdvres 
    367442 
     
    421496       
    422497      !--- 2.13 ice concentration should not exceed amax  
     498      !         (it should not be the case)  
    423499      !----------------------------------------------------- 
    424500      DO jj = 1, jpj 
    425501         DO ji = 1, jpi 
    426502            z_da_ex =  MAX( at_i(ji,jj) - amax , 0.0 )         
     503            zindb   =  MAX( rzero, SIGN( rone, at_i(ji,jj) - epsi06 ) )  
    427504            DO jl  = 1, jpl 
    428                zindb   =  MAX( rzero, SIGN( rone, v_i(ji,jj,jl) - epsi04 ) )  
    429505               z_da_i = a_i(ji,jj,jl) * z_da_ex / MAX( at_i(ji,jj), epsi06 ) * zindb 
    430                a_i(ji,jj,jl) = a_i(ji,jj,jl) - z_da_i 
     506               a_i(ji,jj,jl) = MAX( 0._wp, a_i(ji,jj,jl) - z_da_i ) 
     507               ! 
     508               zinda   =  MAX( rzero, SIGN( rone, a_i(ji,jj,jl) - epsi06 ) )  
     509               ht_i(ji,jj,jl) = v_i(ji,jj,jl) / MAX( a_i(ji,jj,jl), epsi06 ) * zinda 
     510               !v_i(ji,jj,jl) = ht_i(ji,jj,jl) * a_i(ji,jj,jl) ! makes ice shrinken but should not be used 
    431511            END DO 
    432512         END DO 
     
    517597         zchk_smv = ( glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zchk_smv_b ) / rdt_ice + ( zchk_fs / rhoic ) 
    518598 
     599         zchk_vmin = glob_min(v_i) 
     600         zchk_amax = glob_max(SUM(a_i,dim=3)) 
     601         zchk_amin = glob_min(a_i) 
     602 
    519603         IF(lwp) THEN 
    520             IF (    ABS( zchk_v_i   ) >  1.e-5 ) WRITE(numout,*) 'violation volume [m3/day]     (limupdate2) = ',(zchk_v_i * 86400.) 
    521             IF (    ABS( zchk_smv   ) >  1.e-4 ) WRITE(numout,*) 'violation saline [psu*m3/day] (limupdate2) = ',(zchk_smv * 86400.) 
    522             IF ( MINVAL( v_i(:,:,:) ) <  0.    ) WRITE(numout,*) 'violation v_i<0  [mm]         (limupdate2) = ',(MINVAL(v_i) * 1.e-3) 
    523             IF ( MAXVAL( SUM(a_i(:,:,:),dim=3) ) > amax+epsi10 ) WRITE(numout,*) 'violation a_i>amax    (limupdate2) = ',MAXVAL(SUM(a_i,dim=3)) 
    524             IF ( MINVAL( a_i(:,:,:) ) <  0.    ) WRITE(numout,*) 'violation a_i<0               (limupdate2) = ',MINVAL(a_i) 
     604            IF ( ABS( zchk_v_i   ) >  1.e-5 ) WRITE(numout,*) 'violation volume [m3/day]     (limupdate2) = ',(zchk_v_i * 86400.) 
     605            IF ( ABS( zchk_smv   ) >  1.e-4 ) WRITE(numout,*) 'violation saline [psu*m3/day] (limupdate2) = ',(zchk_smv * 86400.) 
     606            IF ( zchk_vmin <  0.            ) WRITE(numout,*) 'violation v_i<0  [mm]         (limupdate2) = ',(zchk_vmin * 1.e-3) 
     607            IF ( zchk_amax >  amax+epsi10   ) WRITE(numout,*) 'violation a_i>amax            (limupdate2) = ',zchk_amax 
     608            IF ( zchk_amin <  0.            ) WRITE(numout,*) 'violation a_i<0               (limupdate2) = ',zchk_amin 
    525609         ENDIF 
    526610      ENDIF 
  • branches/2013/dev_r3406_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limvar.F90

    r3938 r3963  
    180180               ht_s(ji,jj,jl) = v_s (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi10 ) * zindb 
    181181               o_i(ji,jj,jl)  = oa_i(ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi10 ) * zindb 
     182               a_i(ji,jj,jl) = a_i (ji,jj,jl) * zindb ! clem correction 
    182183            END DO 
    183184         END DO 
  • branches/2013/dev_r3406_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limwri.F90

    r3938 r3963  
    122122         ! Normal file 
    123123         !------------- 
    124  
    125          zsto     = rdt_ice 
    126          IF( ln_mskland )   THEN   ;   clop = "ave(only(x))"   ! put 1.e+20 on land (very expensive!!) 
    127          ELSE                      ;   clop = "ave(x)"         ! no use of the mask value (require less cpu time) 
    128          ENDIF 
    129          zout     = nwrite * rdt_ice / nn_fsbc 
    130124         niter    = ( nit000 - 1 ) / nn_fsbc 
    131          zdept(1) = 0. 
    132  
    133125         CALL ymds2ju ( nyear, nmonth, nday, rdt, zjulian ) 
    134126         zjulian = zjulian - adatrj   !   set calendar origin to the beginning of the experiment 
    135          CALL dia_nam ( clhstnam, nwrite, 'icemod_old' ) 
    136          CALL histbeg ( clhstnam, jpi, glamt, jpj, gphit, 1, jpi, 1, jpj, niter, zjulian, rdt_ice,   & 
    137             &           nhorid, nice, domain_id=nidom, snc4chunks=snc4set ) 
    138          CALL histvert( nice, "deptht", "Vertical T levels", "m", 1, zdept, ndepid, "down") 
    139          CALL wheneq  ( jpij , tmask(:,:,1), 1, 1., ndex51, ndim) 
    140  
    141          DO jf = 1 , noumef 
    142             IF(lwp) WRITE(numout,*) 'jf', jf 
    143             IF ( nc(jf) == 1 ) THEN 
    144                CALL histdef( nice, nam(jf), titn(jf), uni(jf), jpi, jpj & 
    145                   , nhorid, 1, 1, 1, -99, 32, clop, zsto, zout ) 
    146                IF(lwp) WRITE(numout,*) 'nice, nam(jf), titn(jf), uni(jf), nhorid, clop, zsto, zout' 
    147                IF(lwp) WRITE(numout,*)  nice, nam(jf), titn(jf), uni(jf), nhorid, clop, zsto, zout  
    148             ENDIF 
    149          END DO 
    150  
    151          CALL histend(nice, snc4set) 
    152  
     127!clem 
     128!         zsto     = rdt_ice 
     129!         IF( ln_mskland )   THEN   ;   clop = "ave(only(x))"   ! put 1.e+20 on land (very expensive!!) 
     130!         ELSE                      ;   clop = "ave(x)"         ! no use of the mask value (require less cpu time) 
     131!         ENDIF 
     132!         zout     = nwrite * rdt_ice / nn_fsbc 
     133!         zdept(1) = 0. 
     134! 
     135!         CALL dia_nam ( clhstnam, nwrite, 'icemod_old' ) 
     136!         CALL histbeg ( clhstnam, jpi, glamt, jpj, gphit, 1, jpi, 1, jpj, niter, zjulian, rdt_ice,   & 
     137!            &           nhorid, nice, domain_id=nidom, snc4chunks=snc4set ) 
     138!         CALL histvert( nice, "deptht", "Vertical T levels", "m", 1, zdept, ndepid, "down") 
     139!         CALL wheneq  ( jpij , tmask(:,:,1), 1, 1., ndex51, ndim) 
     140! 
     141!         DO jf = 1 , noumef 
     142!            IF(lwp) WRITE(numout,*) 'jf', jf 
     143!            IF ( nc(jf) == 1 ) THEN 
     144!               CALL histdef( nice, nam(jf), titn(jf), uni(jf), jpi, jpj & 
     145!                  , nhorid, 1, 1, 1, -99, 32, clop, zsto, zout ) 
     146!               IF(lwp) WRITE(numout,*) 'nice, nam(jf), titn(jf), uni(jf), nhorid, clop, zsto, zout' 
     147!               IF(lwp) WRITE(numout,*)  nice, nam(jf), titn(jf), uni(jf), nhorid, clop, zsto, zout  
     148!            ENDIF 
     149!         END DO 
     150! 
     151!         CALL histend(nice, snc4set) 
     152!clem 
     153         ! 
    153154         !----------------- 
    154155         ! ITD file output 
     
    280281      ! ecriture d'un fichier netcdf 
    281282      ! 
     283  
    282284      niter = niter + 1 
    283       DO jf = 1 , noumef 
    284          ! 
    285          zfield(:,:) = zcmo(:,:,jf) * cmulti(jf) + cadd(jf) 
    286          ! 
    287          IF( jf == 7  .OR. jf == 8  .OR. jf == 15 .OR. jf == 16 ) THEN   ;   CALL lbc_lnk( zfield, 'T', -1. ) 
    288          ELSE                                                            ;   CALL lbc_lnk( zfield, 'T',  1. ) 
    289          ENDIF 
    290          ! 
    291          IF( ln_nicep ) THEN  
    292             WRITE(numout,*) 
    293             WRITE(numout,*) 'nc(jf), nice, nam(jf), niter, ndim' 
    294             WRITE(numout,*) nc(jf), nice, nam(jf), niter, ndim 
    295          ENDIF 
    296          IF( nc(jf) == 1 ) CALL histwrite( nice, nam(jf), niter, zfield, ndim, ndex51 ) 
    297          ! 
    298       END DO 
    299  
    300       IF( ( nn_fsbc * niter ) >= nitend .OR. kindic < 0 ) THEN 
    301          IF( lwp) WRITE(numout,*) ' Closing the icemod file ' 
    302          CALL histclo( nice ) 
    303       ENDIF 
    304  
     285!clem 
     286!      DO jf = 1 , noumef 
     287!         ! 
     288!         zfield(:,:) = zcmo(:,:,jf) * cmulti(jf) + cadd(jf) 
     289!         ! 
     290!         IF( jf == 7  .OR. jf == 8  .OR. jf == 15 .OR. jf == 16 ) THEN   ;   CALL lbc_lnk( zfield, 'T', -1. ) 
     291!         ELSE                                                            ;   CALL lbc_lnk( zfield, 'T',  1. ) 
     292!         ENDIF 
     293!         ! 
     294!         IF( ln_nicep ) THEN  
     295!            WRITE(numout,*) 
     296!            WRITE(numout,*) 'nc(jf), nice, nam(jf), niter, ndim' 
     297!            WRITE(numout,*) nc(jf), nice, nam(jf), niter, ndim 
     298!         ENDIF 
     299!         IF( nc(jf) == 1 ) CALL histwrite( nice, nam(jf), niter, zfield, ndim, ndex51 ) 
     300!         ! 
     301!      END DO 
     302! 
     303!      IF( ( nn_fsbc * niter ) >= nitend .OR. kindic < 0 ) THEN 
     304!         IF( lwp) WRITE(numout,*) ' Closing the icemod file ' 
     305!         CALL histclo( nice ) 
     306!      ENDIF 
     307!clem 
     308      ! 
    305309       CALL iom_put ('iceconc', zcmo(:,:,1) )          ! field1: ice concentration 
    306310       CALL iom_put ('icethic_cea', zcmo(:,:,2) )      ! field2: ice thickness (i.e. icethi(:,:)) 
     
    328332       CALL iom_put ('isnowhc', zcmo(:,:,26) )         ! field 26: snow total heat content 
    329333       CALL iom_put ('icest', zcmo(:,:,27) )           ! field 27: ice surface temperature 
    330        CALL iom_put ('fsbri', zcmo(:,:,28) * 86400 )           ! field 28: brine salt flux 
    331        CALL iom_put ('fseqv', zcmo(:,:,29) * 86400 )           ! field 29: equivalent FW salt flux 
     334       CALL iom_put ('fsbri', zcmo(:,:,28) * 86400. )           ! field 28: brine salt flux 
     335       CALL iom_put ('fseqv', zcmo(:,:,29) * 86400. )           ! field 29: equivalent FW salt flux 
    332336       CALL iom_put ('ibrinv', zcmo(:,:,30) *100 )     ! field 30: brine volume 
    333337       CALL iom_put ('icecolf', zcmo(:,:,31) )         ! field 31: frazil ice collection thickness 
     
    341345       CALL iom_put ('icevolu', zcmo(:,:,39) ) ! field 39: ice volume 
    342346       CALL iom_put ('snowvol', zcmo(:,:,40) ) ! field 40: snow volume 
    343        CALL iom_put ('fsrpo', zcmo(:,:,41) * 86400 )           ! field 41: salt flux from ridging rafting 
    344        CALL iom_put ('fsres', zcmo(:,:,42) * 86400 )           ! field 42: salt flux from limupdate (resultant) 
     347       CALL iom_put ('fsrpo', zcmo(:,:,41) * 86400. )           ! field 41: salt flux from ridging rafting 
     348       CALL iom_put ('fsres', zcmo(:,:,42) * 86400. )           ! field 42: salt flux from limupdate (resultant) 
    345349       CALL iom_put ('icetrp', zcmo(:,:,43) )    ! field 43: ice volume transport 
    346350 
  • branches/2013/dev_r3406_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limwri_dimg.h90

    r3938 r3963  
    1414   !!  modif : 03/06/98 
    1515   !!------------------------------------------------------------------- 
     16   USE lib_fortran     ! to use key_nosignedzero 
    1617   USE  diawri, ONLY : dia_wri_dimg 
    1718   REAL(wp),DIMENSION(1) ::   zdept 
  • branches/2013/dev_r3406_CNRS_LIM3/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90

    r3938 r3963  
    4141   USE sbc_ice         ! Surface boundary condition: ice fields 
    4242#endif 
     43   USE lib_fortran     ! to use key_nosignedzero 
    4344 
    4445   IMPLICIT NONE 
     
    530531            DO ji = fs_2, fs_jpim1   ! vect. opt. 
    531532               p_taui(ji,jj) = zcoef_wnorm2 * ( z_wnds_t(ji+1,jj  ) + z_wnds_t(ji,jj) )                          & 
    532                   &          * ( 0.5 * (sf(jp_wndi)%fnow(ji+1,jj,1) + sf(jp_wndi)%fnow(ji,jj,1) ) - pui(ji,jj) ) 
     533                  &          * ( 0.5 * (sf(jp_wndi)%fnow(ji+1,jj,1) + sf(jp_wndi)%fnow(ji,jj,1) ) - rn_vfac * pui(ji,jj) ) 
    533534               p_tauj(ji,jj) = zcoef_wnorm2 * ( z_wnds_t(ji,jj+1  ) + z_wnds_t(ji,jj) )                          & 
    534                   &          * ( 0.5 * (sf(jp_wndj)%fnow(ji,jj+1,1) + sf(jp_wndj)%fnow(ji,jj,1) ) - pvi(ji,jj) ) 
     535                  &          * ( 0.5 * (sf(jp_wndj)%fnow(ji,jj+1,1) + sf(jp_wndj)%fnow(ji,jj,1) ) - rn_vfac * pvi(ji,jj) ) 
    535536            END DO 
    536537         END DO 
     
    560561               ! iovino 
    561562               IF( ff(ji,jj) .GT. 0._wp ) THEN 
    562                   z_qlw(ji,jj,jl) = ( 0.95 * sf(jp_qlw)%fnow(ji,jj,1) - Stef * pst(ji,jj,jl) * zst3 ) * tmask(ji,jj,1) ! MV test 2 garde juste cette ligne ci 
     563                  z_qlw(ji,jj,jl) = ( 0.95 * sf(jp_qlw)%fnow(ji,jj,1) - Stef * pst(ji,jj,jl) * zst3 ) * tmask(ji,jj,1) 
    563564               ELSE 
    564                   z_qlw(ji,jj,jl) = 0.95 * ( sf(jp_qlw)%fnow(ji,jj,1) - Stef * pst(ji,jj,jl) * zst3 ) * tmask(ji,jj,1) ! MV test 1 garde juste cette ligne ci 
     565                  z_qlw(ji,jj,jl) = 0.95 * ( sf(jp_qlw)%fnow(ji,jj,1) - Stef * pst(ji,jj,jl) * zst3 ) * tmask(ji,jj,1) 
    565566               ENDIF 
    566567               ! lw sensitivity 
     
    610611!CDIR COLLAPSE 
    611612      p_spr(:,:) = sf(jp_snow)%fnow(:,:,1) * rn_pfac      ! solid precipitation [kg/m2/s] 
    612       CALL iom_put( 'snowpre', p_spr * 86400 )                  ! Snow precipitation  
    613       CALL iom_put( 'precip', p_tpr * 86400 )                   ! Total precipitation  
     613      CALL iom_put( 'snowpre', p_spr * 86400. )                  ! Snow precipitation  
     614      CALL iom_put( 'precip', p_tpr * 86400. )                   ! Total precipitation  
    614615      ! 
    615616      IF(ln_ctl) THEN 
  • branches/2013/dev_r3406_CNRS_LIM3/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90

    r3938 r3963  
    9797 
    9898      !                          ! overwrite namelist parameter using CPP key information 
    99       !IF( Agrif_Root() ) THEN                ! AGRIF zoom 
    100       !  IF( lk_lim2 )   nn_ice      = 2 
    101       !  IF( lk_lim3 )   nn_ice      = 3 
    102       !  IF( lk_cice )   nn_ice      = 4 
    103       !ENDIF 
     99      IF( Agrif_Root() ) THEN                ! AGRIF zoom 
     100        IF( lk_lim2 )   nn_ice      = 2 
     101        IF( lk_lim3 )   nn_ice      = 3 
     102        IF( lk_cice )   nn_ice      = 4 
     103      ENDIF 
    104104      IF( cp_cfg == 'gyre' ) THEN            ! GYRE configuration 
    105105          ln_ana      = .TRUE.    
  • branches/2013/dev_r3406_CNRS_LIM3/NEMOGCM/NEMO/OPA_SRC/lib_fortran.F90

    r3294 r3963  
    55   !!====================================================================== 
    66   !! History :  3.2  !  2010-05  (M. Dunphy, R. Benshila)  Original code 
     7   !!            3.4  !  2013-06  (C. Rousset)  add glob_min, glob_max  
     8   !!                                           + 3d dim. of input is fexible (jpk, jpl...)  
    79   !!---------------------------------------------------------------------- 
    810 
     
    2325 
    2426   PUBLIC glob_sum 
     27   PUBLIC glob_min, glob_max 
    2528#if defined key_nosignedzero 
    2629   PUBLIC SIGN 
     
    2932   INTERFACE glob_sum 
    3033      MODULE PROCEDURE glob_sum_2d, glob_sum_3d,glob_sum_2d_a, glob_sum_3d_a  
     34   END INTERFACE 
     35   INTERFACE glob_min 
     36      MODULE PROCEDURE glob_min_2d, glob_min_3d,glob_min_2d_a, glob_min_3d_a  
     37   END INTERFACE 
     38   INTERFACE glob_max 
     39      MODULE PROCEDURE glob_max_2d, glob_max_3d,glob_max_2d_a, glob_max_3d_a  
    3140   END INTERFACE 
    3241 
     
    4756 
    4857#if ! defined key_mpp_rep 
     58 
     59   ! --- SUM --- 
    4960   FUNCTION glob_sum_2d( ptab )  
    5061      !!----------------------------------------------------------------------- 
     
    6172      ! 
    6273   END FUNCTION glob_sum_2d 
    63     
    64     
     74   
    6575   FUNCTION glob_sum_3d( ptab )  
    6676      !!----------------------------------------------------------------------- 
     
    7383      !! 
    7484      INTEGER :: jk 
    75       !!----------------------------------------------------------------------- 
     85      INTEGER :: zjpk ! local variable: size of the 3d dimension of ptab 
     86      !!----------------------------------------------------------------------- 
     87      ! 
     88      zjpk = SIZE(ptab,3) 
    7689      ! 
    7790      glob_sum_3d = 0.e0 
    78       DO jk = 1, jpk 
     91      DO jk = 1, zjpk 
    7992         glob_sum_3d = glob_sum_3d + SUM( ptab(:,:,jk)*tmask_i(:,:) ) 
    8093      END DO 
     
    111124      !! 
    112125      INTEGER :: jk 
    113       !!----------------------------------------------------------------------- 
     126      INTEGER :: zjpk ! local variable: size of the 3d dimension of ptab 
     127      !!----------------------------------------------------------------------- 
     128      ! 
     129      zjpk = SIZE(ptab1,3) 
    114130      ! 
    115131      glob_sum_3d_a(:) = 0.e0 
    116       DO jk = 1, jpk 
     132      DO jk = 1, zjpk 
    117133         glob_sum_3d_a(1) = glob_sum_3d_a(1) + SUM( ptab1(:,:,jk)*tmask_i(:,:) ) 
    118134         glob_sum_3d_a(2) = glob_sum_3d_a(2) + SUM( ptab2(:,:,jk)*tmask_i(:,:) ) 
     
    121137      ! 
    122138   END FUNCTION glob_sum_3d_a 
     139  
     140 
     141   ! --- MIN --- 
     142   FUNCTION glob_min_2d( ptab )  
     143      !!----------------------------------------------------------------------- 
     144      !!                  ***  FUNCTION  glob_min_2D  *** 
     145      !! 
     146      !! ** Purpose : perform a masked min on the inner global domain of a 2D array 
     147      !!----------------------------------------------------------------------- 
     148      REAL(wp), INTENT(in), DIMENSION(:,:) ::   ptab          ! input 2D array 
     149      REAL(wp)                             ::   glob_min_2d   ! global masked min 
     150      !!----------------------------------------------------------------------- 
     151      ! 
     152      glob_min_2d = MINVAL( ptab(:,:)*tmask_i(:,:) ) 
     153      IF( lk_mpp )   CALL mpp_min( glob_min_2d ) 
     154      ! 
     155   END FUNCTION glob_min_2d 
     156  
     157   FUNCTION glob_min_3d( ptab )  
     158      !!----------------------------------------------------------------------- 
     159      !!                  ***  FUNCTION  glob_min_3D  *** 
     160      !! 
     161      !! ** Purpose : perform a masked min on the inner global domain of a 3D array 
     162      !!----------------------------------------------------------------------- 
     163      REAL(wp), INTENT(in), DIMENSION(:,:,:) ::   ptab          ! input 3D array 
     164      REAL(wp)                               ::   glob_min_3d   ! global masked min 
     165      !! 
     166      INTEGER :: jk 
     167      INTEGER :: zjpk ! local variable: size of the 3d dimension of ptab 
     168      !!----------------------------------------------------------------------- 
     169      ! 
     170      zjpk = SIZE(ptab,3) 
     171      ! 
     172      glob_min_3d = 0.e0 
     173      DO jk = 1, zjpk 
     174         glob_min_3d = glob_min_3d + MINVAL( ptab(:,:,jk)*tmask_i(:,:) ) 
     175      END DO 
     176      IF( lk_mpp )   CALL mpp_min( glob_min_3d ) 
     177      ! 
     178   END FUNCTION glob_min_3d 
     179 
     180 
     181   FUNCTION glob_min_2d_a( ptab1, ptab2 )  
     182      !!----------------------------------------------------------------------- 
     183      !!                  ***  FUNCTION  glob_min_2D _a *** 
     184      !! 
     185      !! ** Purpose : perform a masked min on the inner global domain of two 2D array 
     186      !!----------------------------------------------------------------------- 
     187      REAL(wp), INTENT(in), DIMENSION(:,:) ::   ptab1, ptab2    ! input 2D array 
     188      REAL(wp)            , DIMENSION(2)   ::   glob_min_2d_a   ! global masked min 
     189      !!----------------------------------------------------------------------- 
     190      !              
     191      glob_min_2d_a(1) = MINVAL( ptab1(:,:)*tmask_i(:,:) ) 
     192      glob_min_2d_a(2) = MINVAL( ptab2(:,:)*tmask_i(:,:) ) 
     193      IF( lk_mpp )   CALL mpp_min( glob_min_2d_a, 2 ) 
     194      ! 
     195   END FUNCTION glob_min_2d_a 
     196  
     197  
     198   FUNCTION glob_min_3d_a( ptab1, ptab2 )  
     199      !!----------------------------------------------------------------------- 
     200      !!                  ***  FUNCTION  glob_min_3D_a *** 
     201      !! 
     202      !! ** Purpose : perform a masked min on the inner global domain of two 3D array 
     203      !!----------------------------------------------------------------------- 
     204      REAL(wp), INTENT(in), DIMENSION(:,:,:) ::   ptab1, ptab2    ! input 3D array 
     205      REAL(wp)            , DIMENSION(2)     ::   glob_min_3d_a   ! global masked min 
     206      !! 
     207      INTEGER :: jk 
     208      INTEGER :: zjpk ! local variable: size of the 3d dimension of ptab 
     209      !!----------------------------------------------------------------------- 
     210      ! 
     211      zjpk = SIZE(ptab1,3) 
     212      ! 
     213      glob_min_3d_a(:) = 0.e0 
     214      DO jk = 1, zjpk 
     215         glob_min_3d_a(1) = glob_min_3d_a(1) + MINVAL( ptab1(:,:,jk)*tmask_i(:,:) ) 
     216         glob_min_3d_a(2) = glob_min_3d_a(2) + MINVAL( ptab2(:,:,jk)*tmask_i(:,:) ) 
     217      END DO 
     218      IF( lk_mpp )   CALL mpp_min( glob_min_3d_a, 2 ) 
     219      ! 
     220   END FUNCTION glob_min_3d_a 
     221 
     222   ! --- MAX --- 
     223   FUNCTION glob_max_2d( ptab )  
     224      !!----------------------------------------------------------------------- 
     225      !!                  ***  FUNCTION  glob_max_2D  *** 
     226      !! 
     227      !! ** Purpose : perform a masked max on the inner global domain of a 2D array 
     228      !!----------------------------------------------------------------------- 
     229      REAL(wp), INTENT(in), DIMENSION(:,:) ::   ptab          ! input 2D array 
     230      REAL(wp)                             ::   glob_max_2d   ! global masked max 
     231      !!----------------------------------------------------------------------- 
     232      ! 
     233      glob_max_2d = MAXVAL( ptab(:,:)*tmask_i(:,:) ) 
     234      IF( lk_mpp )   CALL mpp_max( glob_max_2d ) 
     235      ! 
     236   END FUNCTION glob_max_2d 
     237  
     238   FUNCTION glob_max_3d( ptab )  
     239      !!----------------------------------------------------------------------- 
     240      !!                  ***  FUNCTION  glob_max_3D  *** 
     241      !! 
     242      !! ** Purpose : perform a masked max on the inner global domain of a 3D array 
     243      !!----------------------------------------------------------------------- 
     244      REAL(wp), INTENT(in), DIMENSION(:,:,:) ::   ptab          ! input 3D array 
     245      REAL(wp)                               ::   glob_max_3d   ! global masked max 
     246      !! 
     247      INTEGER :: jk 
     248      INTEGER :: zjpk ! local variable: size of the 3d dimension of ptab 
     249      !!----------------------------------------------------------------------- 
     250      ! 
     251      zjpk = SIZE(ptab,3) 
     252      ! 
     253      glob_max_3d = 0.e0 
     254      DO jk = 1, zjpk 
     255         glob_max_3d = glob_max_3d + MAXVAL( ptab(:,:,jk)*tmask_i(:,:) ) 
     256      END DO 
     257      IF( lk_mpp )   CALL mpp_max( glob_max_3d ) 
     258      ! 
     259   END FUNCTION glob_max_3d 
     260 
     261 
     262   FUNCTION glob_max_2d_a( ptab1, ptab2 )  
     263      !!----------------------------------------------------------------------- 
     264      !!                  ***  FUNCTION  glob_max_2D _a *** 
     265      !! 
     266      !! ** Purpose : perform a masked max on the inner global domain of two 2D array 
     267      !!----------------------------------------------------------------------- 
     268      REAL(wp), INTENT(in), DIMENSION(:,:) ::   ptab1, ptab2    ! input 2D array 
     269      REAL(wp)            , DIMENSION(2)   ::   glob_max_2d_a   ! global masked max 
     270      !!----------------------------------------------------------------------- 
     271      !              
     272      glob_max_2d_a(1) = MAXVAL( ptab1(:,:)*tmask_i(:,:) ) 
     273      glob_max_2d_a(2) = MAXVAL( ptab2(:,:)*tmask_i(:,:) ) 
     274      IF( lk_mpp )   CALL mpp_max( glob_max_2d_a, 2 ) 
     275      ! 
     276   END FUNCTION glob_max_2d_a 
     277  
     278  
     279   FUNCTION glob_max_3d_a( ptab1, ptab2 )  
     280      !!----------------------------------------------------------------------- 
     281      !!                  ***  FUNCTION  glob_max_3D_a *** 
     282      !! 
     283      !! ** Purpose : perform a masked max on the inner global domain of two 3D array 
     284      !!----------------------------------------------------------------------- 
     285      REAL(wp), INTENT(in), DIMENSION(:,:,:) ::   ptab1, ptab2    ! input 3D array 
     286      REAL(wp)            , DIMENSION(2)     ::   glob_max_3d_a   ! global masked max 
     287      !! 
     288      INTEGER :: jk 
     289      INTEGER :: zjpk ! local variable: size of the 3d dimension of ptab 
     290      !!----------------------------------------------------------------------- 
     291      ! 
     292      zjpk = SIZE(ptab1,3) 
     293      ! 
     294      glob_max_3d_a(:) = 0.e0 
     295      DO jk = 1, zjpk 
     296         glob_max_3d_a(1) = glob_max_3d_a(1) + MAXVAL( ptab1(:,:,jk)*tmask_i(:,:) ) 
     297         glob_max_3d_a(2) = glob_max_3d_a(2) + MAXVAL( ptab2(:,:,jk)*tmask_i(:,:) ) 
     298      END DO 
     299      IF( lk_mpp )   CALL mpp_max( glob_max_3d_a, 2 ) 
     300      ! 
     301   END FUNCTION glob_max_3d_a 
     302 
    123303 
    124304#else   
     
    127307   !!---------------------------------------------------------------------- 
    128308    
     309   ! --- SUM --- 
    129310   FUNCTION glob_sum_2d( ptab )  
    130311      !!---------------------------------------------------------------------- 
     
    133314      !! ** Purpose : perform a sum in calling DDPDD routine 
    134315      !!---------------------------------------------------------------------- 
    135       REAL(wp), INTENT(in), DIMENSION(jpi,jpj) ::   ptab 
    136       REAL(wp)                                 ::   glob_sum_2d   ! global masked sum 
     316      REAL(wp), INTENT(in), DIMENSION(:,:) ::   ptab 
     317      REAL(wp)                             ::   glob_sum_2d   ! global masked sum 
    137318      !! 
    138319      COMPLEX(wp)::   ctmp 
    139320      REAL(wp)   ::   ztmp 
    140321      INTEGER    ::   ji, jj   ! dummy loop indices 
    141       !!----------------------------------------------------------------------- 
    142       ! 
    143       ztmp = 0.e0 
    144       ctmp = CMPLX( 0.e0, 0.e0, wp ) 
    145       DO jj = 1, jpj 
    146          DO ji =1, jpi 
     322      INTEGER    ::   zjpi, zjpj ! local variables: size of ptab 
     323      !!----------------------------------------------------------------------- 
     324      zjpi = SIZE(ptab,1) 
     325      zjpj = SIZE(ptab,2) 
     326      ! 
     327      ztmp = 0.e0 
     328      ctmp = CMPLX( 0.e0, 0.e0, wp ) 
     329      DO jj = 1, zjpj 
     330         DO ji =1, zjpi 
    147331         ztmp =  ptab(ji,jj) * tmask_i(ji,jj) 
    148332         CALL DDPDD( CMPLX( ztmp, 0.e0, wp ), ctmp ) 
     
    161345      !! ** Purpose : perform a sum on a 3D array in calling DDPDD routine 
    162346      !!---------------------------------------------------------------------- 
    163       REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk) ::   ptab 
    164       REAL(wp)                                     ::   glob_sum_3d   ! global masked sum 
     347      REAL(wp), INTENT(in), DIMENSION(:,:,:) ::   ptab 
     348      REAL(wp)                               ::   glob_sum_3d   ! global masked sum 
    165349      !! 
    166350      COMPLEX(wp)::   ctmp 
    167351      REAL(wp)   ::   ztmp 
    168352      INTEGER    ::   ji, jj, jk   ! dummy loop indices 
    169       !!----------------------------------------------------------------------- 
    170       ! 
    171       ztmp = 0.e0 
    172       ctmp = CMPLX( 0.e0, 0.e0, wp ) 
    173       DO jk = 1, jpk 
    174          DO jj = 1, jpj 
    175             DO ji =1, jpi 
     353      INTEGER    ::   zjpi, zjpj, zjpk ! local variables: size of ptab 
     354      !!----------------------------------------------------------------------- 
     355      ! 
     356      zjpi = SIZE(ptab,1) 
     357      zjpj = SIZE(ptab,2) 
     358      zjpk = SIZE(ptab,3) 
     359      ! 
     360      ztmp = 0.e0 
     361      ctmp = CMPLX( 0.e0, 0.e0, wp ) 
     362      DO jk = 1, zjpk 
     363         DO jj = 1, zjpj 
     364            DO ji =1, zjpi 
    176365            ztmp =  ptab(ji,jj,jk) * tmask_i(ji,jj) 
    177366            CALL DDPDD( CMPLX( ztmp, 0.e0, wp ), ctmp ) 
     
    191380      !! ** Purpose : perform a sum on two 2D arrays in calling DDPDD routine 
    192381      !!---------------------------------------------------------------------- 
    193       REAL(wp), INTENT(in), DIMENSION(jpi,jpj) ::   ptab1, ptab2 
    194       REAL(wp)                                 ::   glob_sum_2d_a   ! global masked sum 
     382      REAL(wp), INTENT(in), DIMENSION(:,:) ::   ptab1, ptab2 
     383      REAL(wp)                             ::   glob_sum_2d_a   ! global masked sum 
    195384      !! 
    196385      COMPLEX(wp)::   ctmp 
    197386      REAL(wp)   ::   ztmp 
    198387      INTEGER    ::   ji, jj   ! dummy loop indices 
    199       !!----------------------------------------------------------------------- 
    200       ! 
    201       ztmp = 0.e0 
    202       ctmp = CMPLX( 0.e0, 0.e0, wp ) 
    203       DO jj = 1, jpj 
    204          DO ji =1, jpi 
     388      INTEGER    ::   zjpi, zjpj ! local variables: size of ptab 
     389      !!----------------------------------------------------------------------- 
     390      ! 
     391      zjpi = SIZE(ptab1,1) 
     392      zjpj = SIZE(ptab1,2) 
     393      ! 
     394      ztmp = 0.e0 
     395      ctmp = CMPLX( 0.e0, 0.e0, wp ) 
     396      DO jj = 1, zjpj 
     397         DO ji =1, zjpi 
    205398         ztmp =  ptab1(ji,jj) * tmask_i(ji,jj) 
    206399         CALL DDPDD( CMPLX( ztmp, 0.e0, wp ), ctmp ) 
     
    221414      !! ** Purpose : perform a sum on two 3D array in calling DDPDD routine 
    222415      !!---------------------------------------------------------------------- 
    223       REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk) ::   ptab1, ptab2 
    224       REAL(wp)                                     ::   glob_sum_3d_a   ! global masked sum 
     416      REAL(wp), INTENT(in), DIMENSION(:,:,:) ::   ptab1, ptab2 
     417      REAL(wp)                               ::   glob_sum_3d_a   ! global masked sum 
    225418      !! 
    226419      COMPLEX(wp)::   ctmp 
    227420      REAL(wp)   ::   ztmp 
    228421      INTEGER    ::   ji, jj, jk   ! dummy loop indices 
    229       !!----------------------------------------------------------------------- 
    230       ! 
    231       ztmp = 0.e0 
    232       ctmp = CMPLX( 0.e0, 0.e0, wp ) 
    233       DO jk = 1, jpk 
    234          DO jj = 1, jpj 
    235             DO ji =1, jpi 
     422      INTEGER    ::   zjpi, zjpj, zjpk ! local variables: size of ptab 
     423      !!----------------------------------------------------------------------- 
     424      ! 
     425      zjpi = SIZE(ptab1,1) 
     426      zjpj = SIZE(ptab1,2) 
     427      zjpk = SIZE(ptab1,3) 
     428      ! 
     429      ztmp = 0.e0 
     430      ctmp = CMPLX( 0.e0, 0.e0, wp ) 
     431      DO jk = 1, zjpk 
     432         DO jj = 1, zjpj 
     433            DO ji =1, zjpi 
    236434            ztmp =  ptab1(ji,jj,jk) * tmask_i(ji,jj) 
    237435            CALL DDPDD( CMPLX( ztmp, 0.e0, wp ), ctmp ) 
     
    247445 
    248446 
     447   ! --- MIN --- 
     448   FUNCTION glob_min_2d( ptab )  
     449      !!---------------------------------------------------------------------- 
     450      !!                  ***  FUNCTION  glob_min_2d *** 
     451      !! 
     452      !! ** Purpose : perform a min in calling DDPDD routine 
     453      !!---------------------------------------------------------------------- 
     454      REAL(wp), INTENT(in), DIMENSION(:,:) ::   ptab 
     455      REAL(wp)                             ::   glob_min_2d   ! global masked min 
     456      !! 
     457      COMPLEX(wp)::   ctmp 
     458      REAL(wp)   ::   ztmp 
     459      INTEGER    ::   ji, jj   ! dummy loop indices 
     460      INTEGER    ::   zjpi, zjpj ! local variables: size of ptab 
     461      !!----------------------------------------------------------------------- 
     462      zjpi = SIZE(ptab,1) 
     463      zjpj = SIZE(ptab,2) 
     464      ! 
     465      ztmp = 0.e0 
     466      ctmp = CMPLX( 0.e0, 0.e0, wp ) 
     467      DO jj = 1, zjpj 
     468         DO ji =1, zjpi 
     469         ztmp =  ptab(ji,jj) * tmask_i(ji,jj) 
     470         CALL DDPDD( CMPLX( ztmp, 0.e0, wp ), ctmp ) 
     471         END DO 
     472      END DO 
     473      IF( lk_mpp )   CALL mpp_min( ctmp )   ! min over the global domain 
     474      glob_min_2d = REAL(ctmp,wp) 
     475      ! 
     476   END FUNCTION glob_min_2d    
     477 
     478 
     479   FUNCTION glob_min_3d( ptab )  
     480      !!---------------------------------------------------------------------- 
     481      !!                  ***  FUNCTION  glob_min_3d *** 
     482      !! 
     483      !! ** Purpose : perform a min on a 3D array in calling DDPDD routine 
     484      !!---------------------------------------------------------------------- 
     485      REAL(wp), INTENT(in), DIMENSION(:,:,:) ::   ptab 
     486      REAL(wp)                               ::   glob_min_3d   ! global masked min 
     487      !! 
     488      COMPLEX(wp)::   ctmp 
     489      REAL(wp)   ::   ztmp 
     490      INTEGER    ::   ji, jj, jk   ! dummy loop indices 
     491      INTEGER    ::   zjpi, zjpj, zjpk ! local variables: size of ptab 
     492      !!----------------------------------------------------------------------- 
     493      ! 
     494      zjpi = SIZE(ptab,1) 
     495      zjpj = SIZE(ptab,2) 
     496      zjpk = SIZE(ptab,3) 
     497      ! 
     498      ztmp = 0.e0 
     499      ctmp = CMPLX( 0.e0, 0.e0, wp ) 
     500      DO jk = 1, zjpk 
     501         DO jj = 1, zjpj 
     502            DO ji =1, zjpi 
     503            ztmp =  ptab(ji,jj,jk) * tmask_i(ji,jj) 
     504            CALL DDPDD( CMPLX( ztmp, 0.e0, wp ), ctmp ) 
     505            END DO 
     506         END DO     
     507      END DO 
     508      IF( lk_mpp )   CALL mpp_min( ctmp )   ! min over the global domain 
     509      glob_min_3d = REAL(ctmp,wp) 
     510      ! 
     511   END FUNCTION glob_min_3d    
     512 
     513 
     514   FUNCTION glob_min_2d_a( ptab1, ptab2 )  
     515      !!---------------------------------------------------------------------- 
     516      !!                  ***  FUNCTION  glob_min_2d_a *** 
     517      !! 
     518      !! ** Purpose : perform a min on two 2D arrays in calling DDPDD routine 
     519      !!---------------------------------------------------------------------- 
     520      REAL(wp), INTENT(in), DIMENSION(:,:) ::   ptab1, ptab2 
     521      REAL(wp)                             ::   glob_min_2d_a   ! global masked min 
     522      !! 
     523      COMPLEX(wp)::   ctmp 
     524      REAL(wp)   ::   ztmp 
     525      INTEGER    ::   ji, jj   ! dummy loop indices 
     526      INTEGER    ::   zjpi, zjpj ! local variables: size of ptab 
     527      !!----------------------------------------------------------------------- 
     528      ! 
     529      zjpi = SIZE(ptab1,1) 
     530      zjpj = SIZE(ptab1,2) 
     531      ! 
     532      ztmp = 0.e0 
     533      ctmp = CMPLX( 0.e0, 0.e0, wp ) 
     534      DO jj = 1, zjpj 
     535         DO ji =1, zjpi 
     536         ztmp =  ptab1(ji,jj) * tmask_i(ji,jj) 
     537         CALL DDPDD( CMPLX( ztmp, 0.e0, wp ), ctmp ) 
     538         ztmp =  ptab2(ji,jj) * tmask_i(ji,jj) 
     539         CALL DDPDD( CMPLX( ztmp, 0.e0, wp ), ctmp ) 
     540         END DO 
     541      END DO 
     542      IF( lk_mpp )   CALL mpp_min( ctmp )   ! min over the global domain 
     543      glob_min_2d_a = REAL(ctmp,wp) 
     544      ! 
     545   END FUNCTION glob_min_2d_a    
     546 
     547 
     548   FUNCTION glob_min_3d_a( ptab1, ptab2 )  
     549      !!---------------------------------------------------------------------- 
     550      !!                  ***  FUNCTION  glob_min_3d_a *** 
     551      !! 
     552      !! ** Purpose : perform a min on two 3D array in calling DDPDD routine 
     553      !!---------------------------------------------------------------------- 
     554      REAL(wp), INTENT(in), DIMENSION(:,:,:) ::   ptab1, ptab2 
     555      REAL(wp)                               ::   glob_min_3d_a   ! global masked min 
     556      !! 
     557      COMPLEX(wp)::   ctmp 
     558      REAL(wp)   ::   ztmp 
     559      INTEGER    ::   ji, jj, jk   ! dummy loop indices 
     560      INTEGER    ::   zjpi, zjpj, zjpk ! local variables: size of ptab 
     561      !!----------------------------------------------------------------------- 
     562      ! 
     563      zjpi = SIZE(ptab1,1) 
     564      zjpj = SIZE(ptab1,2) 
     565      zjpk = SIZE(ptab1,3) 
     566      ! 
     567      ztmp = 0.e0 
     568      ctmp = CMPLX( 0.e0, 0.e0, wp ) 
     569      DO jk = 1, zjpk 
     570         DO jj = 1, zjpj 
     571            DO ji =1, zjpi 
     572            ztmp =  ptab1(ji,jj,jk) * tmask_i(ji,jj) 
     573            CALL DDPDD( CMPLX( ztmp, 0.e0, wp ), ctmp ) 
     574            ztmp =  ptab2(ji,jj,jk) * tmask_i(ji,jj) 
     575            CALL DDPDD( CMPLX( ztmp, 0.e0, wp ), ctmp ) 
     576            END DO 
     577         END DO     
     578      END DO 
     579      IF( lk_mpp )   CALL mpp_min( ctmp )   ! min over the global domain 
     580      glob_min_3d_a = REAL(ctmp,wp) 
     581      ! 
     582   END FUNCTION glob_min_3d_a    
     583 
     584  
     585   ! --- MAX --- 
     586   FUNCTION glob_max_2d( ptab )  
     587      !!---------------------------------------------------------------------- 
     588      !!                  ***  FUNCTION  glob_max_2d *** 
     589      !! 
     590      !! ** Purpose : perform a max in calling DDPDD routine 
     591      !!---------------------------------------------------------------------- 
     592      REAL(wp), INTENT(in), DIMENSION(:,:) ::   ptab 
     593      REAL(wp)                             ::   glob_max_2d   ! global masked max 
     594      !! 
     595      COMPLEX(wp)::   ctmp 
     596      REAL(wp)   ::   ztmp 
     597      INTEGER    ::   ji, jj   ! dummy loop indices 
     598      INTEGER    ::   zjpi, zjpj ! local variables: size of ptab 
     599      !!----------------------------------------------------------------------- 
     600      zjpi = SIZE(ptab,1) 
     601      zjpj = SIZE(ptab,2) 
     602      ! 
     603      ztmp = 0.e0 
     604      ctmp = CMPLX( 0.e0, 0.e0, wp ) 
     605      DO jj = 1, zjpj 
     606         DO ji =1, zjpi 
     607         ztmp =  ptab(ji,jj) * tmask_i(ji,jj) 
     608         CALL DDPDD( CMPLX( ztmp, 0.e0, wp ), ctmp ) 
     609         END DO 
     610      END DO 
     611      IF( lk_mpp )   CALL mpp_max( ctmp )   ! max over the global domain 
     612      glob_max_2d = REAL(ctmp,wp) 
     613      ! 
     614   END FUNCTION glob_max_2d    
     615 
     616 
     617   FUNCTION glob_max_3d( ptab )  
     618      !!---------------------------------------------------------------------- 
     619      !!                  ***  FUNCTION  glob_max_3d *** 
     620      !! 
     621      !! ** Purpose : perform a max on a 3D array in calling DDPDD routine 
     622      !!---------------------------------------------------------------------- 
     623      REAL(wp), INTENT(in), DIMENSION(:,:,:) ::   ptab 
     624      REAL(wp)                               ::   glob_max_3d   ! global masked max 
     625      !! 
     626      COMPLEX(wp)::   ctmp 
     627      REAL(wp)   ::   ztmp 
     628      INTEGER    ::   ji, jj, jk   ! dummy loop indices 
     629      INTEGER    ::   zjpi, zjpj, zjpk ! local variables: size of ptab 
     630      !!----------------------------------------------------------------------- 
     631      ! 
     632      zjpi = SIZE(ptab,1) 
     633      zjpj = SIZE(ptab,2) 
     634      zjpk = SIZE(ptab,3) 
     635      ! 
     636      ztmp = 0.e0 
     637      ctmp = CMPLX( 0.e0, 0.e0, wp ) 
     638      DO jk = 1, zjpk 
     639         DO jj = 1, zjpj 
     640            DO ji =1, zjpi 
     641            ztmp =  ptab(ji,jj,jk) * tmask_i(ji,jj) 
     642            CALL DDPDD( CMPLX( ztmp, 0.e0, wp ), ctmp ) 
     643            END DO 
     644         END DO     
     645      END DO 
     646      IF( lk_mpp )   CALL mpp_max( ctmp )   ! max over the global domain 
     647      glob_max_3d = REAL(ctmp,wp) 
     648      ! 
     649   END FUNCTION glob_max_3d    
     650 
     651 
     652   FUNCTION glob_max_2d_a( ptab1, ptab2 )  
     653      !!---------------------------------------------------------------------- 
     654      !!                  ***  FUNCTION  glob_max_2d_a *** 
     655      !! 
     656      !! ** Purpose : perform a max on two 2D arrays in calling DDPDD routine 
     657      !!---------------------------------------------------------------------- 
     658      REAL(wp), INTENT(in), DIMENSION(:,:) ::   ptab1, ptab2 
     659      REAL(wp)                             ::   glob_max_2d_a   ! global masked max 
     660      !! 
     661      COMPLEX(wp)::   ctmp 
     662      REAL(wp)   ::   ztmp 
     663      INTEGER    ::   ji, jj   ! dummy loop indices 
     664      INTEGER    ::   zjpi, zjpj ! local variables: size of ptab 
     665      !!----------------------------------------------------------------------- 
     666      ! 
     667      zjpi = SIZE(ptab1,1) 
     668      zjpj = SIZE(ptab1,2) 
     669      ! 
     670      ztmp = 0.e0 
     671      ctmp = CMPLX( 0.e0, 0.e0, wp ) 
     672      DO jj = 1, zjpj 
     673         DO ji =1, zjpi 
     674         ztmp =  ptab1(ji,jj) * tmask_i(ji,jj) 
     675         CALL DDPDD( CMPLX( ztmp, 0.e0, wp ), ctmp ) 
     676         ztmp =  ptab2(ji,jj) * tmask_i(ji,jj) 
     677         CALL DDPDD( CMPLX( ztmp, 0.e0, wp ), ctmp ) 
     678         END DO 
     679      END DO 
     680      IF( lk_mpp )   CALL mpp_max( ctmp )   ! max over the global domain 
     681      glob_max_2d_a = REAL(ctmp,wp) 
     682      ! 
     683   END FUNCTION glob_max_2d_a    
     684 
     685 
     686   FUNCTION glob_max_3d_a( ptab1, ptab2 )  
     687      !!---------------------------------------------------------------------- 
     688      !!                  ***  FUNCTION  glob_max_3d_a *** 
     689      !! 
     690      !! ** Purpose : perform a max on two 3D array in calling DDPDD routine 
     691      !!---------------------------------------------------------------------- 
     692      REAL(wp), INTENT(in), DIMENSION(:,:,:) ::   ptab1, ptab2 
     693      REAL(wp)                               ::   glob_max_3d_a   ! global masked max 
     694      !! 
     695      COMPLEX(wp)::   ctmp 
     696      REAL(wp)   ::   ztmp 
     697      INTEGER    ::   ji, jj, jk   ! dummy loop indices 
     698      INTEGER    ::   zjpi, zjpj, zjpk ! local variables: size of ptab 
     699      !!----------------------------------------------------------------------- 
     700      ! 
     701      zjpi = SIZE(ptab1,1) 
     702      zjpj = SIZE(ptab1,2) 
     703      zjpk = SIZE(ptab1,3) 
     704      ! 
     705      ztmp = 0.e0 
     706      ctmp = CMPLX( 0.e0, 0.e0, wp ) 
     707      DO jk = 1, zjpk 
     708         DO jj = 1, zjpj 
     709            DO ji =1, zjpi 
     710            ztmp =  ptab1(ji,jj,jk) * tmask_i(ji,jj) 
     711            CALL DDPDD( CMPLX( ztmp, 0.e0, wp ), ctmp ) 
     712            ztmp =  ptab2(ji,jj,jk) * tmask_i(ji,jj) 
     713            CALL DDPDD( CMPLX( ztmp, 0.e0, wp ), ctmp ) 
     714            END DO 
     715         END DO     
     716      END DO 
     717      IF( lk_mpp )   CALL mpp_max( ctmp )   ! max over the global domain 
     718      glob_max_3d_a = REAL(ctmp,wp) 
     719      ! 
     720   END FUNCTION glob_max_3d_a    
     721 
     722 
    249723   SUBROUTINE DDPDD( ydda, yddb ) 
    250724      !!---------------------------------------------------------------------- 
Note: See TracChangeset for help on using the changeset viewer.