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

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

Ignore:
Timestamp:
2015-05-12T12:37:15+02:00 (9 years ago)
Author:
deazer
Message:

Merged branch with Trunk at revision 5253.
Checked with SETTE, passes modified iodef.xml for AMM12 experiment

Location:
branches/2014/dev_r4650_UKMO10_Tidally_Meaned_Diagnostics/NEMOGCM/NEMO/OPA_SRC/DIA
Files:
9 edited

Legend:

Unmodified
Added
Removed

  • branches/2014/dev_r4650_UKMO10_Tidally_Meaned_Diagnostics/NEMOGCM/NEMO/OPA_SRC/DIA/dia25h.F90

    r4756 r5260  
    88   USE oce             ! ocean dynamics and tracers variables 
    99   USE dom_oce         ! ocean space and time domain 
    10    USE insitu_tem, ONLY: insitu_t, theta2t 
     10   USE diainsitutem, ONLY: insitu_t, theta2t 
    1111   USE in_out_manager  ! I/O units 
    1212   USE iom             ! I/0 library 

  • branches/2014/dev_r4650_UKMO10_Tidally_Meaned_Diagnostics/NEMOGCM/NEMO/OPA_SRC/DIA/diaar5.F90

    r4313 r5260  
    2121   USE timing         ! preformance summary 
    2222   USE wrk_nemo       ! working arrays 
     23   USE fldread        ! type FLD_N 
     24   USE phycst         ! physical constant 
     25   USE in_out_manager  ! I/O manager 
    2326 
    2427   IMPLICIT NONE 
    … …  
    8386      CALL wrk_alloc( jpi , jpj , jpk , jpts , ztsn                 ) 
    8487 
    85       CALL iom_put( 'cellthc', fse3t(:,:,:) ) 
    86  
    8788      zarea_ssh(:,:) = area(:,:) * sshn(:,:) 
    8889 
    … …  
    104105         zbotpres(:,:) = zbotpres(:,:) + fse3t(:,:,jk) * zrhd(:,:,jk) 
    105106      END DO 
    106       IF( .NOT.lk_vvl )   zbotpres(:,:) = zbotpres(:,:) + sshn(:,:) * zrhd(:,:,1) 
     107      IF( .NOT.lk_vvl ) THEN 
     108         IF ( ln_isfcav ) THEN 
     109            DO ji=1,jpi 
     110               DO jj=1,jpj 
     111                  zbotpres(ji,jj) = zbotpres(ji,jj) + sshn(ji,jj) * zrhd(ji,jj,mikt(ji,jj)) + riceload(ji,jj) 
     112               END DO 
     113            END DO 
     114         ELSE 
     115            zbotpres(:,:) = zbotpres(:,:) + sshn(:,:) * zrhd(:,:,1) 
     116         END IF 
     117      END IF 
    107118      !                                          
    108119      zarho = SUM( area(:,:) * zbotpres(:,:) )  
    … …  
    120131         zbotpres(:,:) = zbotpres(:,:) + fse3t(:,:,jk) * zrhd(:,:,jk) 
    121132      END DO 
    122       IF( .NOT.lk_vvl )   zbotpres(:,:) = zbotpres(:,:) + sshn(:,:) * zrhd(:,:,1) 
     133      IF( .NOT.lk_vvl ) THEN 
     134         IF ( ln_isfcav ) THEN 
     135            DO ji=1,jpi 
     136               DO jj=1,jpj 
     137                  zbotpres(ji,jj) = zbotpres(ji,jj) + sshn(ji,jj) * zrhd(ji,jj,mikt(ji,jj)) + riceload(ji,jj) 
     138               END DO 
     139            END DO 
     140         ELSE 
     141            zbotpres(:,:) = zbotpres(:,:) + sshn(:,:) * zrhd(:,:,1) 
     142         END IF 
     143      END IF 
    123144      !     
    124145      zarho = SUM( area(:,:) * zbotpres(:,:) )  
    … …  
    145166      END DO 
    146167      IF( .NOT.lk_vvl ) THEN 
    147          ztemp = ztemp + SUM( zarea_ssh(:,:) * tsn(:,:,1,jp_tem) ) 
    148          zsal  = zsal  + SUM( zarea_ssh(:,:) * tsn(:,:,1,jp_sal) ) 
     168         IF ( ln_isfcav ) THEN 
     169            DO ji=1,jpi 
     170               DO jj=1,jpj 
     171                  ztemp = ztemp + zarea_ssh(ji,jj) * tsn(ji,jj,mikt(ji,jj),jp_tem)  
     172                  zsal  = zsal  + zarea_ssh(ji,jj) * tsn(ji,jj,mikt(ji,jj),jp_sal)  
     173               END DO 
     174            END DO 
     175         ELSE 
     176            ztemp = ztemp + SUM( zarea_ssh(:,:) * tsn(:,:,1,jp_tem) ) 
     177            zsal  = zsal  + SUM( zarea_ssh(:,:) * tsn(:,:,1,jp_sal) ) 
     178         END IF 
    149179      ENDIF 
    150180      IF( lk_mpp ) THEN   
    … …  
    181211      REAL(wp) ::   zztmp   
    182212      REAL(wp), POINTER, DIMENSION(:,:,:,:) ::   zsaldta   ! Jan/Dec levitus salinity 
     213      ! reading initial file 
     214      LOGICAL  ::   ln_tsd_init      !: T & S data flag 
     215      LOGICAL  ::   ln_tsd_tradmp    !: internal damping toward input data flag 
     216      CHARACTER(len=100)            ::   cn_dir 
     217      TYPE(FLD_N)                   ::  sn_tem,sn_sal 
     218      INTEGER  ::   ios=0 
     219 
     220      NAMELIST/namtsd/ ln_tsd_init,ln_tsd_tradmp,cn_dir,sn_tem,sn_sal 
     221      ! 
     222 
     223      REWIND( numnam_ref )              ! Namelist namtsd in reference namelist : 
     224      READ  ( numnam_ref, namtsd, IOSTAT = ios, ERR = 901) 
     225901   IF( ios /= 0 ) CALL ctl_nam ( ios , ' namtsd in reference namelist for dia_ar5', lwp ) 
     226      REWIND( numnam_cfg )              ! Namelist namtsd in configuration namelist : Parameters of the run 
     227      READ  ( numnam_cfg, namtsd, IOSTAT = ios, ERR = 902 ) 
     228902   IF( ios /= 0 ) CALL ctl_nam ( ios , ' namtsd in configuration namelist for dia_ar5', lwp ) 
     229      IF(lwm) WRITE ( numond, namtsd ) 
     230      ! 
    183231      !!---------------------------------------------------------------------- 
    184232      ! 
    … …  
    200248      END DO 
    201249      IF( lk_mpp )   CALL mpp_sum( vol0 ) 
    202        
    203       CALL iom_open ( 'data_1m_salinity_nomask', inum ) 
    204       CALL iom_get  ( inum, jpdom_data, 'vosaline', zsaldta(:,:,:,1), 1  ) 
    205       CALL iom_get  ( inum, jpdom_data, 'vosaline', zsaldta(:,:,:,2), 12 ) 
     250 
     251      CALL iom_open ( TRIM( cn_dir )//TRIM(sn_sal%clname), inum ) 
     252      CALL iom_get  ( inum, jpdom_data, TRIM(sn_sal%clvar), zsaldta(:,:,:,1), 1  ) 
     253      CALL iom_get  ( inum, jpdom_data, TRIM(sn_sal%clvar), zsaldta(:,:,:,2), 12 ) 
    206254      CALL iom_close( inum ) 
    207255      sn0(:,:,:) = 0.5_wp * ( zsaldta(:,:,:,1) + zsaldta(:,:,:,2) )         

  • branches/2014/dev_r4650_UKMO10_Tidally_Meaned_Diagnostics/NEMOGCM/NEMO/OPA_SRC/DIA/diadct.F90

    • Property svn:keywords set to Id
    r4624 r5260  
    4242#endif 
    4343#if defined key_lim3 
    44   USE par_ice 
    4544  USE ice 
    4645#endif 
    … …  
    113112  REAL(wp), ALLOCATABLE, DIMENSION(:,:,:)   ::  transports_2d   
    114113 
     114   !! $Id$ 
    115115CONTAINS 
    116116 
    … …  
    12981298   LOGICAL, PUBLIC, PARAMETER ::   lk_diadct = .FALSE.    !: diamht flag 
    12991299   PUBLIC  
     1300   !! $Id$ 
    13001301CONTAINS 
    13011302 

  • branches/2014/dev_r4650_UKMO10_Tidally_Meaned_Diagnostics/NEMOGCM/NEMO/OPA_SRC/DIA/diafwb.F90

    r4147 r5260  
    77   !!            8.5  !  02-06  (G. Madec)  F90: Free form and module 
    88   !!            9.0  !  05-11  (V. Garnier) Surface pressure gradient organization 
    9    !!---------------------------------------------------------------------- 
    10 #if ! defined key_coupled 
    11   
     9   !!----------------------------------------------------------------------  
    1210   !!---------------------------------------------------------------------- 
    1311   !!   Only for ORCA2 ORCA1 and ORCA025 
    … …  
    2927 
    3028   PUBLIC dia_fwb    ! routine called by step.F90 
    31  
    32    LOGICAL, PUBLIC, PARAMETER ::   lk_diafwb = .TRUE.    !: fresh water budget flag 
    3329 
    3430   REAL(wp)               ::   a_fwf ,          & 
    … …  
    453449   END SUBROUTINE dia_fwb 
    454450 
    455 #else 
    456    !!---------------------------------------------------------------------- 
    457    !!   Default option :                                       Dummy Module 
    458    !!---------------------------------------------------------------------- 
    459    LOGICAL, PUBLIC, PARAMETER ::   lk_diafwb = .FALSE.    !: fresh water budget flag 
    460 CONTAINS 
    461    SUBROUTINE dia_fwb( kt )        ! Empty routine 
    462       WRITE(*,*) 'dia_fwb: : You should not have seen this print! error?', kt 
    463    END SUBROUTINE dia_fwb 
    464 #endif 
    465  
    466451   !!====================================================================== 
    467452END MODULE diafwb 

  • branches/2014/dev_r4650_UKMO10_Tidally_Meaned_Diagnostics/NEMOGCM/NEMO/OPA_SRC/DIA/diaharm.F90

    • Property svn:keywords set to Id
    r4624 r5260  
    1818   USE daymod 
    1919   USE tide_mod 
     20   ! 
    2021   USE in_out_manager  ! I/O units 
    2122   USE iom             ! I/0 library 
    … …  
    3435   INTEGER, PARAMETER :: jpdimsparse  = jpincomax*300*24 
    3536 
    36    !                            !!!namelist variables 
     37   !                         !!** namelist variables ** 
    3738   INTEGER ::   nit000_han    ! First time step used for harmonic analysis 
    3839   INTEGER ::   nitend_han    ! Last time step used for harmonic analysis 
    3940   INTEGER ::   nstep_han     ! Time step frequency for harmonic analysis 
    40    INTEGER ::   nb_ana           ! Number of harmonics to analyse 
     41   INTEGER ::   nb_ana        ! Number of harmonics to analyse 
    4142 
    4243   INTEGER , ALLOCATABLE, DIMENSION(:)       ::   name 
    … …  
    5960   !!---------------------------------------------------------------------- 
    6061   !! NEMO/OPA 3.5 , NEMO Consortium (2013) 
    61    !! $Id:$ 
     62   !! $Id$ 
    6263   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
    6364   !!---------------------------------------------------------------------- 
    … …  
    119120            ENDIF 
    120121         END DO 
    121       ENDDO 
     122      END DO 
    122123      ! 
    123124      IF(lwp) THEN 
    … …  
    158159      ! ---------------------------- 
    159160      ALLOCATE( ana_temp(jpi,jpj,2*nb_ana,3) ) 
    160       ana_temp(:,:,:,:) = 0.e0 
     161      ana_temp(:,:,:,:) = 0._wp 
    161162 
    162163   END SUBROUTINE dia_harm_init 
    … …  
    179180      IF( nn_timing == 1 )   CALL timing_start('dia_harm') 
    180181 
    181       IF ( kt == nit000 ) CALL dia_harm_init 
    182  
    183       IF ( ((kt.GE.nit000_han).AND.(kt.LE.nitend_han)).AND. & 
    184            (MOD(kt,nstep_han).EQ.0) ) THEN 
    185  
    186         ztime = (kt-nit000+1)*rdt  
     182      IF( kt == nit000 ) CALL dia_harm_init 
     183 
     184      IF( kt >= nit000_han .AND. kt <= nitend_han .AND. MOD(kt,nstep_han) == 0 ) THEN 
     185 
     186         ztime = (kt-nit000+1) * rdt  
    187187        
    188         nhc = 0 
    189         DO jh = 1,nb_ana 
    190           DO jc = 1,2 
    191             nhc = nhc+1 
    192             ztemp =(     MOD(jc,2) * ft(jh) *COS(ana_freq(jh)*ztime + vt(jh) + ut(jh))  & 
    193                     +(1.-MOD(jc,2))* ft(jh) *SIN(ana_freq(jh)*ztime + vt(jh) + ut(jh))) 
    194  
    195             DO jj = 1,jpj 
    196               DO ji = 1,jpi 
    197                 ! Elevation 
    198                 ana_temp(ji,jj,nhc,1) = ana_temp(ji,jj,nhc,1) + ztemp*sshn(ji,jj)           *tmask(ji,jj,1)         
     188         nhc = 0 
     189         DO jh = 1, nb_ana 
     190            DO jc = 1, 2 
     191               nhc = nhc+1 
     192               ztemp =(     MOD(jc,2) * ft(jh) *COS(ana_freq(jh)*ztime + vt(jh) + ut(jh))  & 
     193                  &    +(1.-MOD(jc,2))* ft(jh) *SIN(ana_freq(jh)*ztime + vt(jh) + ut(jh))) 
     194 
     195               DO jj = 1,jpj 
     196                  DO ji = 1,jpi 
     197                     ! Elevation 
     198                     ana_temp(ji,jj,nhc,1) = ana_temp(ji,jj,nhc,1) + ztemp*sshn(ji,jj)           *tmask_i(ji,jj)         
    199199#if defined key_dynspg_ts 
    200                 ana_temp(ji,jj,nhc,2) = ana_temp(ji,jj,nhc,2) + ztemp*un_b(ji,jj)*hur(ji,jj)*umask(ji,jj,1) 
    201                 ana_temp(ji,jj,nhc,3) = ana_temp(ji,jj,nhc,3) + ztemp*vn_b(ji,jj)*hvr(ji,jj)*vmask(ji,jj,1) 
    202 #endif 
    203               END DO 
    204             END DO 
    205  
    206           END DO 
    207         END DO 
    208         
     200                     ana_temp(ji,jj,nhc,2) = ana_temp(ji,jj,nhc,2) + ztemp*un_b(ji,jj)*hur(ji,jj)*umask_i(ji,jj) 
     201                     ana_temp(ji,jj,nhc,3) = ana_temp(ji,jj,nhc,3) + ztemp*vn_b(ji,jj)*hvr(ji,jj)*vmask_i(ji,jj) 
     202#endif 
     203                  END DO 
     204               END DO 
     205               ! 
     206            END DO 
     207         END DO 
     208         !        
    209209      END IF 
    210210 
    … …  
    249249         keq = keq + 1 
    250250         kun = 0 
    251          DO jh = 1,nb_ana 
    252             DO jc = 1,2 
     251         DO jh = 1, nb_ana 
     252            DO jc = 1, 2 
    253253               kun = kun + 1 
    254254               ksp = ksp + 1 
    … …  
    294294               X1 = ana_amp(ji,jj,jh,1) 
    295295               X2 =-ana_amp(ji,jj,jh,2) 
    296                out_eta(ji,jj,jh       ) = X1 * tmask(ji,jj,1) 
    297                out_eta(ji,jj,jh+nb_ana) = X2 * tmask(ji,jj,1) 
    298             ENDDO 
    299          ENDDO 
    300       ENDDO 
     296               out_eta(ji,jj,jh       ) = X1 * tmask_i(ji,jj) 
     297               out_eta(ji,jj,jh+nb_ana) = X2 * tmask_i(ji,jj) 
     298            END DO 
     299         END DO 
     300      END DO 
    301301 
    302302      ! ubar: 
    … …  
    309309                  kun = kun + 1 
    310310                  ztmp4(kun)=ana_temp(ji,jj,kun,2) 
    311                ENDDO 
    312             ENDDO 
     311               END DO 
     312            END DO 
    313313 
    314314            CALL SUR_DETERMINE(jj+1) 
    … …  
    316316            ! Fill output array 
    317317            DO jh = 1, nb_ana 
    318                ana_amp(ji,jj,jh,1)=ztmp7((jh-1)*2+1) 
    319                ana_amp(ji,jj,jh,2)=ztmp7((jh-1)*2+2) 
     318               ana_amp(ji,jj,jh,1) = ztmp7((jh-1)*2+1) 
     319               ana_amp(ji,jj,jh,2) = ztmp7((jh-1)*2+2) 
    320320            END DO 
    321321 
    … …  
    326326         DO ji = 1, jpi 
    327327            DO jh = 1, nb_ana  
    328                X1=ana_amp(ji,jj,jh,1) 
     328               X1= ana_amp(ji,jj,jh,1) 
    329329               X2=-ana_amp(ji,jj,jh,2) 
    330                out_u(ji,jj,jh) = X1 * umask(ji,jj,1) 
    331                out_u (ji,jj,nb_ana+jh) = X2 * umask(ji,jj,1) 
     330               out_u(ji,jj,       jh) = X1 * umask_i(ji,jj) 
     331               out_u(ji,jj,nb_ana+jh) = X2 * umask_i(ji,jj) 
    332332            ENDDO 
    333333         ENDDO 
    … …  
    343343                  kun = kun + 1 
    344344                  ztmp4(kun)=ana_temp(ji,jj,kun,3) 
    345                ENDDO 
    346             ENDDO 
     345               END DO 
     346            END DO 
    347347 
    348348            CALL SUR_DETERMINE(jj+1) 
    … …  
    362362               X1=ana_amp(ji,jj,jh,1) 
    363363               X2=-ana_amp(ji,jj,jh,2) 
    364                out_v(ji,jj,jh)=X1 * vmask(ji,jj,1) 
    365                out_v(ji,jj,nb_ana+jh)=X2 * vmask(ji,jj,1) 
    366             ENDDO 
    367          ENDDO 
    368       ENDDO 
     364               out_v(ji,jj,       jh)=X1 * vmask_i(ji,jj) 
     365               out_v(ji,jj,nb_ana+jh)=X2 * vmask_i(ji,jj) 
     366            END DO 
     367         END DO 
     368      END DO 
    369369 
    370370      CALL dia_wri_harm ! Write results in files 
    … …  
    437437#else 
    438438      DO jh = 1, nb_ana 
    439          CALL iom_put( TRIM(tname(jh))//'x_v', out_u(:,:,jh       ) ) 
    440          CALL iom_put( TRIM(tname(jh))//'y_v', out_u(:,:,jh+nb_ana) ) 
    441       END DO 
    442 #endif 
    443  
     439         CALL iom_put( TRIM(tname(jh))//'x_v', out_v(:,:,jh       ) ) 
     440         CALL iom_put( TRIM(tname(jh))//'y_v', out_v(:,:,jh+nb_ana) ) 
     441      END DO 
     442#endif 
     443      ! 
    444444   END SUBROUTINE dia_wri_harm 
    445445 
    446446 
    447447   SUBROUTINE SUR_DETERMINE(init) 
    448    !!--------------------------------------------------------------------------------- 
    449    !!                      *** ROUTINE SUR_DETERMINE *** 
    450    !!     
    451    !!     
    452    !!        
    453    !!--------------------------------------------------------------------------------- 
    454    INTEGER, INTENT(in) ::   init  
    455    ! 
    456    INTEGER                         :: ji_sd, jj_sd, ji1_sd, ji2_sd, jk1_sd, jk2_sd 
    457    REAL(wp)                        :: zval1, zval2, zx1 
    458    REAL(wp), POINTER, DIMENSION(:) :: ztmpx, zcol1, zcol2 
    459    INTEGER , POINTER, DIMENSION(:) :: ipos2, ipivot 
    460    !--------------------------------------------------------------------------------- 
    461    CALL wrk_alloc( jpincomax , ztmpx , zcol1 , zcol2 ) 
    462    CALL wrk_alloc( jpincomax , ipos2 , ipivot        ) 
     448      !!--------------------------------------------------------------------------------- 
     449      !!                      *** ROUTINE SUR_DETERMINE *** 
     450      !!     
     451      !!     
     452      !!        
     453      !!--------------------------------------------------------------------------------- 
     454      INTEGER, INTENT(in) ::   init  
     455      ! 
     456      INTEGER                         :: ji_sd, jj_sd, ji1_sd, ji2_sd, jk1_sd, jk2_sd 
     457      REAL(wp)                        :: zval1, zval2, zx1 
     458      REAL(wp), POINTER, DIMENSION(:) :: ztmpx, zcol1, zcol2 
     459      INTEGER , POINTER, DIMENSION(:) :: ipos2, ipivot 
     460      !--------------------------------------------------------------------------------- 
     461      CALL wrk_alloc( jpincomax , ztmpx , zcol1 , zcol2 ) 
     462      CALL wrk_alloc( jpincomax , ipos2 , ipivot        ) 
    463463             
    464    IF( init == 1 ) THEN 
    465       IF( nsparse > jpdimsparse )   CALL ctl_stop( 'STOP', 'SUR_DETERMINE : nsparse .GT. jpdimsparse') 
    466       IF( ninco   > jpincomax   )   CALL ctl_stop( 'STOP', 'SUR_DETERMINE : ninco .GT. jpincomax') 
    467       ! 
    468       ztmp3(:,:) = 0._wp 
    469       ! 
    470       DO jk1_sd = 1, nsparse 
    471          DO jk2_sd = 1, nsparse 
    472             nisparse(jk2_sd) = nisparse(jk2_sd) 
    473             njsparse(jk2_sd) = njsparse(jk2_sd) 
    474             IF( nisparse(jk2_sd) == nisparse(jk1_sd) ) THEN 
    475                ztmp3(njsparse(jk1_sd),njsparse(jk2_sd)) = ztmp3(njsparse(jk1_sd),njsparse(jk2_sd))  & 
    476                                                         + valuesparse(jk1_sd)*valuesparse(jk2_sd) 
    477             ENDIF 
    478          END DO 
    479       END DO 
    480  
    481       DO jj_sd = 1 ,ninco 
    482           ipos1(jj_sd) = jj_sd 
    483           ipos2(jj_sd) = jj_sd 
    484       ENDDO 
    485  
    486       DO ji_sd = 1 , ninco 
    487  
    488          !find greatest non-zero pivot: 
    489          zval1 = ABS(ztmp3(ji_sd,ji_sd)) 
    490  
    491          ipivot(ji_sd) = ji_sd 
    492          DO jj_sd = ji_sd, ninco 
    493             zval2 = ABS(ztmp3(ji_sd,jj_sd)) 
    494             IF( zval2.GE.zval1 )THEN 
    495                ipivot(ji_sd) = jj_sd 
    496                zval1         = zval2 
    497             ENDIF 
    498          ENDDO 
    499  
    500          DO ji1_sd = 1, ninco 
    501             zcol1(ji1_sd)               = ztmp3(ji1_sd,ji_sd) 
    502             zcol2(ji1_sd)               = ztmp3(ji1_sd,ipivot(ji_sd)) 
    503             ztmp3(ji1_sd,ji_sd)         = zcol2(ji1_sd) 
    504             ztmp3(ji1_sd,ipivot(ji_sd)) = zcol1(ji1_sd) 
    505          ENDDO 
    506  
    507          ipos2(ji_sd)         = ipos1(ipivot(ji_sd)) 
    508          ipos2(ipivot(ji_sd)) = ipos1(ji_sd) 
    509          ipos1(ji_sd)         = ipos2(ji_sd) 
    510          ipos1(ipivot(ji_sd)) = ipos2(ipivot(ji_sd)) 
    511          zpivot(ji_sd)        = ztmp3(ji_sd,ji_sd) 
    512          DO jj_sd = 1, ninco 
    513             ztmp3(ji_sd,jj_sd) = ztmp3(ji_sd,jj_sd) / zpivot(ji_sd) 
    514          ENDDO 
    515  
     464      IF( init == 1 ) THEN 
     465         IF( nsparse > jpdimsparse )   CALL ctl_stop( 'STOP', 'SUR_DETERMINE : nsparse .GT. jpdimsparse') 
     466         IF( ninco   > jpincomax   )   CALL ctl_stop( 'STOP', 'SUR_DETERMINE : ninco .GT. jpincomax') 
     467         ! 
     468         ztmp3(:,:) = 0._wp 
     469         ! 
     470         DO jk1_sd = 1, nsparse 
     471            DO jk2_sd = 1, nsparse 
     472               nisparse(jk2_sd) = nisparse(jk2_sd) 
     473               njsparse(jk2_sd) = njsparse(jk2_sd) 
     474               IF( nisparse(jk2_sd) == nisparse(jk1_sd) ) THEN 
     475                  ztmp3(njsparse(jk1_sd),njsparse(jk2_sd)) = ztmp3(njsparse(jk1_sd),njsparse(jk2_sd))  & 
     476                     &                                     + valuesparse(jk1_sd)*valuesparse(jk2_sd) 
     477               ENDIF 
     478            END DO 
     479         END DO 
     480         ! 
     481         DO jj_sd = 1 ,ninco 
     482            ipos1(jj_sd) = jj_sd 
     483            ipos2(jj_sd) = jj_sd 
     484         END DO 
     485         ! 
     486         DO ji_sd = 1 , ninco 
     487            ! 
     488            !find greatest non-zero pivot: 
     489            zval1 = ABS(ztmp3(ji_sd,ji_sd)) 
     490            ! 
     491            ipivot(ji_sd) = ji_sd 
     492            DO jj_sd = ji_sd, ninco 
     493               zval2 = ABS(ztmp3(ji_sd,jj_sd)) 
     494               IF( zval2.GE.zval1 )THEN 
     495                  ipivot(ji_sd) = jj_sd 
     496                  zval1         = zval2 
     497               ENDIF 
     498            END DO 
     499            ! 
     500            DO ji1_sd = 1, ninco 
     501               zcol1(ji1_sd)               = ztmp3(ji1_sd,ji_sd) 
     502               zcol2(ji1_sd)               = ztmp3(ji1_sd,ipivot(ji_sd)) 
     503               ztmp3(ji1_sd,ji_sd)         = zcol2(ji1_sd) 
     504               ztmp3(ji1_sd,ipivot(ji_sd)) = zcol1(ji1_sd) 
     505            END DO 
     506            ! 
     507            ipos2(ji_sd)         = ipos1(ipivot(ji_sd)) 
     508            ipos2(ipivot(ji_sd)) = ipos1(ji_sd) 
     509            ipos1(ji_sd)         = ipos2(ji_sd) 
     510            ipos1(ipivot(ji_sd)) = ipos2(ipivot(ji_sd)) 
     511            zpivot(ji_sd)        = ztmp3(ji_sd,ji_sd) 
     512            DO jj_sd = 1, ninco 
     513               ztmp3(ji_sd,jj_sd) = ztmp3(ji_sd,jj_sd) / zpivot(ji_sd) 
     514            END DO 
     515            ! 
     516            DO ji2_sd = ji_sd+1, ninco 
     517               zpilier(ji2_sd,ji_sd)=ztmp3(ji2_sd,ji_sd) 
     518               DO jj_sd=1,ninco 
     519                  ztmp3(ji2_sd,jj_sd)=  ztmp3(ji2_sd,jj_sd) - ztmp3(ji_sd,jj_sd) * zpilier(ji2_sd,ji_sd) 
     520               END DO 
     521            END DO 
     522            ! 
     523         END DO 
     524         ! 
     525      ENDIF ! End init==1 
     526 
     527      DO ji_sd = 1, ninco 
     528         ztmp4(ji_sd) = ztmp4(ji_sd) / zpivot(ji_sd) 
    516529         DO ji2_sd = ji_sd+1, ninco 
    517             zpilier(ji2_sd,ji_sd)=ztmp3(ji2_sd,ji_sd) 
    518             DO jj_sd=1,ninco 
    519                ztmp3(ji2_sd,jj_sd)=  ztmp3(ji2_sd,jj_sd) - ztmp3(ji_sd,jj_sd) * zpilier(ji2_sd,ji_sd) 
    520             ENDDO 
    521          ENDDO 
    522  
    523       ENDDO 
    524  
    525    ENDIF ! End init==1 
    526  
    527    DO ji_sd = 1, ninco 
    528       ztmp4(ji_sd) = ztmp4(ji_sd) / zpivot(ji_sd) 
    529       DO ji2_sd = ji_sd+1, ninco 
    530          ztmp4(ji2_sd) = ztmp4(ji2_sd) - ztmp4(ji_sd) * zpilier(ji2_sd,ji_sd) 
    531       ENDDO 
    532    ENDDO 
    533  
    534    !system solving:  
    535    ztmpx(ninco) = ztmp4(ninco) / ztmp3(ninco,ninco) 
    536    ji_sd = ninco 
    537    DO ji_sd = ninco-1, 1, -1 
    538       zx1=0. 
    539       DO jj_sd = ji_sd+1, ninco 
    540          zx1 = zx1 + ztmpx(jj_sd) * ztmp3(ji_sd,jj_sd) 
    541       ENDDO 
    542       ztmpx(ji_sd) = ztmp4(ji_sd)-zx1 
    543    ENDDO 
    544  
    545    DO jj_sd =1, ninco 
    546       ztmp7(ipos1(jj_sd))=ztmpx(jj_sd) 
    547    ENDDO 
    548  
    549    CALL wrk_dealloc( jpincomax , ztmpx , zcol1 , zcol2 ) 
    550    CALL wrk_dealloc( jpincomax , ipos2 , ipivot        ) 
    551  
    552   END SUBROUTINE SUR_DETERMINE 
     530            ztmp4(ji2_sd) = ztmp4(ji2_sd) - ztmp4(ji_sd) * zpilier(ji2_sd,ji_sd) 
     531         END DO 
     532      END DO 
     533 
     534      !system solving:  
     535      ztmpx(ninco) = ztmp4(ninco) / ztmp3(ninco,ninco) 
     536      ji_sd = ninco 
     537      DO ji_sd = ninco-1, 1, -1 
     538         zx1 = 0._wp 
     539         DO jj_sd = ji_sd+1, ninco 
     540            zx1 = zx1 + ztmpx(jj_sd) * ztmp3(ji_sd,jj_sd) 
     541         END DO 
     542         ztmpx(ji_sd) = ztmp4(ji_sd)-zx1 
     543      END DO 
     544 
     545      DO jj_sd =1, ninco 
     546         ztmp7(ipos1(jj_sd))=ztmpx(jj_sd) 
     547      END DO 
     548 
     549      CALL wrk_dealloc( jpincomax , ztmpx , zcol1 , zcol2 ) 
     550      CALL wrk_dealloc( jpincomax , ipos2 , ipivot        ) 
     551      ! 
     552   END SUBROUTINE SUR_DETERMINE 
    553553 
    554554#else 

  • branches/2014/dev_r4650_UKMO10_Tidally_Meaned_Diagnostics/NEMOGCM/NEMO/OPA_SRC/DIA/diahsb.F90

    r4624 r5260  
    99 
    1010   !!---------------------------------------------------------------------- 
     11   !!   dia_hsb       : Diagnose the conservation of ocean heat and salt contents, and volume 
     12   !!   dia_hsb_rst   : Read or write DIA file in restart file 
     13   !!   dia_hsb_init  : Initialization of the conservation diagnostic 
     14   !!---------------------------------------------------------------------- 
    1115   USE oce             ! ocean dynamics and tracers 
    1216   USE dom_oce         ! ocean space and time domain 
    1317   USE phycst          ! physical constants 
    1418   USE sbc_oce         ! surface thermohaline fluxes 
    15    USE in_out_manager  ! I/O manager 
     19   USE sbcrnf          ! river runoff 
     20   USE sbcisf          ! ice shelves 
    1621   USE domvvl          ! vertical scale factors 
    1722   USE traqsr          ! penetrative solar radiation 
    1823   USE trabbc          ! bottom boundary condition  
    19    USE lib_mpp         ! distributed memory computing library 
    2024   USE trabbc          ! bottom boundary condition 
    2125   USE bdy_par         ! (for lk_bdy) 
     26   USE restart         ! ocean restart 
     27   ! 
     28   USE iom             ! I/O manager 
     29   USE in_out_manager  ! I/O manager 
     30   USE lib_fortran     ! glob_sum 
     31   USE lib_mpp         ! distributed memory computing library 
    2232   USE timing          ! preformance summary 
    23    USE iom             ! I/O manager 
    24    USE lib_fortran     ! glob_sum 
    25    USE restart         ! ocean restart 
    26    USE wrk_nemo         ! work arrays 
    27    USE sbcrnf         ! river runoffd 
     33   USE wrk_nemo        ! work arrays 
    2834 
    2935   IMPLICIT NONE 
    … …  
    3642   LOGICAL, PUBLIC ::   ln_diahsb   !: check the heat and salt budgets 
    3743 
    38    REAL(dp)                                ::   surf_tot                ! 
    39    REAL(dp)                                ::   frc_t      , frc_s     , frc_v   ! global forcing trends 
    40    REAL(dp)                                ::   frc_wn_t      , frc_wn_s ! global forcing trends 
    41    REAL(dp), DIMENSION(:,:)  , ALLOCATABLE ::   surf      , ssh_ini              ! 
    42    REAL(dp), DIMENSION(:,:,:), ALLOCATABLE ::   hc_loc_ini, sc_loc_ini, e3t_ini  ! 
    43    REAL(dp), DIMENSION(:,:)  , ALLOCATABLE ::   ssh_hc_loc_ini, ssh_sc_loc_ini 
     44   REAL(wp) ::   surf_tot              ! ocean surface 
     45   REAL(wp) ::   frc_t, frc_s, frc_v   ! global forcing trends 
     46   REAL(wp) ::   frc_wn_t, frc_wn_s    ! global forcing trends 
     47   ! 
     48   REAL(wp), DIMENSION(:,:)  , ALLOCATABLE ::   surf          , ssh_ini          ! 
     49   REAL(wp), DIMENSION(:,:)  , ALLOCATABLE ::   ssh_hc_loc_ini, ssh_sc_loc_ini   ! 
     50   REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   hc_loc_ini, sc_loc_ini, e3t_ini  ! 
    4451 
    4552   !! * Substitutions 
    4653#  include "domzgr_substitute.h90" 
    4754#  include "vectopt_loop_substitute.h90" 
    48  
    4955   !!---------------------------------------------------------------------- 
    5056   !! NEMO/OPA 3.3 , NEMO Consortium (2010) 
    … …  
    5258   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
    5359   !!---------------------------------------------------------------------- 
    54  
    5560CONTAINS 
    5661 
    … …  
    6772      !!--------------------------------------------------------------------------- 
    6873      INTEGER, INTENT(in) ::   kt   ! ocean time-step index 
    69       !! 
    70       INTEGER    ::   jk                          ! dummy loop indice 
    71       REAL(dp)   ::   zdiff_hc    , zdiff_sc      ! heat and salt content variations 
    72       REAL(dp)   ::   zdiff_hc1   , zdiff_sc1     ! -   -   -   -   -   -   -   -  
    73       REAL(dp)   ::   zdiff_v1    , zdiff_v2      ! volume variation 
    74       REAL(dp)   ::   zerr_hc1    , zerr_sc1       ! heat and salt content misfit 
    75       REAL(dp)   ::   zvol_tot                    ! volume 
    76       REAL(dp)   ::   z_frc_trd_t , z_frc_trd_s   !    -     - 
    77       REAL(dp)   ::   z_frc_trd_v                 !    -     - 
    78       REAL(dp)   ::   z_wn_trd_t , z_wn_trd_s   !    -     - 
    79       REAL(dp)   ::   z_ssh_hc , z_ssh_sc   !    -     - 
     74      ! 
     75      INTEGER    ::   ji, jj, jk                  ! dummy loop indice 
     76      REAL(wp)   ::   zdiff_hc    , zdiff_sc      ! heat and salt content variations 
     77      REAL(wp)   ::   zdiff_hc1   , zdiff_sc1     !  -         -     -        -  
     78      REAL(wp)   ::   zdiff_v1    , zdiff_v2      ! volume variation 
     79      REAL(wp)   ::   zerr_hc1    , zerr_sc1      ! heat and salt content misfit 
     80      REAL(wp)   ::   zvol_tot                    ! volume 
     81      REAL(wp)   ::   z_frc_trd_t , z_frc_trd_s   !    -     - 
     82      REAL(wp)   ::   z_frc_trd_v                 !    -     - 
     83      REAL(wp)   ::   z_wn_trd_t , z_wn_trd_s     !    -     - 
     84      REAL(wp)   ::   z_ssh_hc , z_ssh_sc         !    -     - 
     85      REAL(wp), DIMENSION(:,:), POINTER ::   z2d0, z2d1 
    8086      !!--------------------------------------------------------------------------- 
    8187      IF( nn_timing == 1 )   CALL timing_start('dia_hsb')       
    82  
     88      CALL wrk_alloc( jpi,jpj,   z2d0, z2d1 ) 
     89      ! 
     90      tsn(:,:,:,1) = tsn(:,:,:,1) * tmask(:,:,:) ; tsb(:,:,:,1) = tsb(:,:,:,1) * tmask(:,:,:) ; 
     91      tsn(:,:,:,2) = tsn(:,:,:,2) * tmask(:,:,:) ; tsb(:,:,:,2) = tsb(:,:,:,2) * tmask(:,:,:) ; 
    8392      ! ------------------------- ! 
    8493      ! 1 - Trends due to forcing ! 
    8594      ! ------------------------- ! 
    86       z_frc_trd_v = r1_rau0 * glob_sum( - ( emp(:,:) - rnf(:,:) ) * surf(:,:) ) ! volume fluxes 
    87       z_frc_trd_t =           glob_sum( sbc_tsc(:,:,jp_tem) * surf(:,:) )       ! heat fluxes 
    88       z_frc_trd_s =           glob_sum( sbc_tsc(:,:,jp_sal) * surf(:,:) )       ! salt fluxes 
    89       ! Add runoff heat & salt input 
     95      z_frc_trd_v = r1_rau0 * glob_sum( - ( emp(:,:) - rnf(:,:) + rdivisf * fwfisf(:,:) ) * surf(:,:) ) ! volume fluxes 
     96      z_frc_trd_t =           glob_sum( sbc_tsc(:,:,jp_tem) * surf(:,:) )                               ! heat fluxes 
     97      z_frc_trd_s =           glob_sum( sbc_tsc(:,:,jp_sal) * surf(:,:) )                               ! salt fluxes 
     98      ! Add runoff    heat & salt input 
    9099      IF( ln_rnf    )   z_frc_trd_t = z_frc_trd_t + glob_sum( rnf_tsc(:,:,jp_tem) * surf(:,:) ) 
    91100      IF( ln_rnf_sal)   z_frc_trd_s = z_frc_trd_s + glob_sum( rnf_tsc(:,:,jp_sal) * surf(:,:) ) 
     101      ! Add ice shelf heat & salt input 
     102      IF( nn_isf .GE. 1 )  THEN 
     103          z_frc_trd_t = z_frc_trd_t & 
     104              &   + glob_sum( ( risf_tsc(:,:,jp_tem) - rdivisf * fwfisf(:,:) * (-1.9) * r1_rau0 ) * surf(:,:) ) 
     105          z_frc_trd_s = z_frc_trd_s + (1.0_wp - rdivisf) * glob_sum( risf_tsc(:,:,jp_sal) * surf(:,:) ) 
     106      ENDIF 
    92107 
    93108      ! Add penetrative solar radiation 
    … …  
    97112      ! 
    98113      IF( .NOT. lk_vvl ) THEN 
    99          z_wn_trd_t = - glob_sum( surf(:,:) * wn(:,:,1) * tsb(:,:,1,jp_tem) ) 
    100          z_wn_trd_s = - glob_sum( surf(:,:) * wn(:,:,1) * tsb(:,:,1,jp_sal) ) 
     114         IF ( ln_isfcav ) THEN 
     115            DO ji=1,jpi 
     116               DO jj=1,jpj 
     117                  z2d0(ji,jj) = surf(ji,jj) * wn(ji,jj,mikt(ji,jj)) * tsb(ji,jj,mikt(ji,jj),jp_tem) 
     118                  z2d1(ji,jj) = surf(ji,jj) * wn(ji,jj,mikt(ji,jj)) * tsb(ji,jj,mikt(ji,jj),jp_sal) 
     119               ENDDO 
     120            ENDDO 
     121         ELSE 
     122            z2d0(:,:) = surf(:,:) * wn(:,:,1) * tsb(:,:,1,jp_tem) 
     123            z2d1(:,:) = surf(:,:) * wn(:,:,1) * tsb(:,:,1,jp_sal) 
     124         END IF 
     125         z_wn_trd_t = - glob_sum( z2d0 )  
     126         z_wn_trd_s = - glob_sum( z2d1 ) 
    101127      ENDIF 
    102128 
    … …  
    113139      ! 2 -  Content variations ! 
    114140      ! ------------------------ ! 
    115       zdiff_v2 = 0.d0 
    116       zdiff_hc = 0.d0 
    117       zdiff_sc = 0.d0 
     141      zdiff_v2 = 0._wp 
     142      zdiff_hc = 0._wp 
     143      zdiff_sc = 0._wp 
    118144 
    119145      ! volume variation (calculated with ssh) 
    … …  
    122148      ! heat & salt content variation (associated with ssh) 
    123149      IF( .NOT. lk_vvl ) THEN 
    124          z_ssh_hc = glob_sum( surf(:,:) * ( tsn(:,:,1,jp_tem) * sshn(:,:) - ssh_hc_loc_ini(:,:) ) ) 
    125          z_ssh_sc = glob_sum( surf(:,:) * ( tsn(:,:,1,jp_sal) * sshn(:,:) - ssh_sc_loc_ini(:,:) ) ) 
     150         IF ( ln_isfcav ) THEN 
     151            DO ji = 1, jpi 
     152               DO jj = 1, jpj 
     153                  z2d0(ji,jj) = surf(ji,jj) * ( tsn(ji,jj,mikt(ji,jj),jp_tem) * sshn(ji,jj) - ssh_hc_loc_ini(ji,jj) )  
     154                  z2d1(ji,jj) = surf(ji,jj) * ( tsn(ji,jj,mikt(ji,jj),jp_sal) * sshn(ji,jj) - ssh_sc_loc_ini(ji,jj) )  
     155               END DO 
     156            END DO 
     157         ELSE 
     158            z2d0(:,:) = surf(:,:) * ( tsn(:,:,1,jp_tem) * sshn(:,:) - ssh_hc_loc_ini(:,:) )  
     159            z2d1(:,:) = surf(:,:) * ( tsn(:,:,1,jp_sal) * sshn(:,:) - ssh_sc_loc_ini(:,:) )  
     160         END IF 
     161         z_ssh_hc = glob_sum( z2d0 )  
     162         z_ssh_sc = glob_sum( z2d1 )  
    126163      ENDIF 
    127164 
    … …  
    153190      ! 3 - Diagnostics writing ! 
    154191      ! ----------------------- ! 
    155       zvol_tot   = 0.d0                                                   ! total ocean volume 
     192      zvol_tot = 0._wp                    ! total ocean volume (calculated with scale factors) 
    156193      DO jk = 1, jpkm1 
    157194         zvol_tot  = zvol_tot + glob_sum( surf(:,:) * tmask(:,:,jk) * fse3t_n(:,:,jk) ) 
    158195      END DO 
     196 
     197!!gm to be added ? 
     198!      IF( .NOT. lk_vvl ) THEN            ! fixed volume, add the ssh contribution 
     199!        zvol_tot = zvol_tot + glob_sum( surf(:,:) * sshn(:,:) ) 
     200!      ENDIF 
     201!!gm end 
     202 
    159203 
    160204      IF( lk_vvl ) THEN 
    … …  
    183227      IF( lrst_oce )   CALL dia_hsb_rst( kt, 'WRITE' ) 
    184228 
     229      CALL wrk_dealloc( jpi,jpj,   z2d0, z2d1 ) 
     230 
    185231      IF( nn_timing == 1 )   CALL timing_stop('dia_hsb') 
    186 ! 
     232      ! 
    187233   END SUBROUTINE dia_hsb 
    188234 
    189  
    190    SUBROUTINE dia_hsb_init 
    191       !!--------------------------------------------------------------------------- 
    192       !!                  ***  ROUTINE dia_hsb  *** 
    193       !!      
    194       !! ** Purpose: Initialization for the heat salt volume budgets 
    195       !!  
    196       !! ** Method : Compute initial heat content, salt content and volume 
    197       !! 
    198       !! ** Action : - Compute initial heat content, salt content and volume 
    199       !!             - Initialize forcing trends 
    200       !!             - Compute coefficients for conversion 
    201       !!--------------------------------------------------------------------------- 
    202       INTEGER            ::   jk       ! dummy loop indice 
    203       INTEGER            ::   ierror   ! local integer 
    204       !! 
    205       NAMELIST/namhsb/ ln_diahsb 
    206       ! 
    207       INTEGER  ::   ios 
    208       !!---------------------------------------------------------------------- 
    209  
    210       IF(lwp) THEN 
    211          WRITE(numout,*) 
    212          WRITE(numout,*) 'dia_hsb_init : check the heat and salt budgets' 
    213          WRITE(numout,*) '~~~~~~~~ ' 
    214       ENDIF 
    215  
    216       REWIND( numnam_ref )              ! Namelist namhsb in reference namelist 
    217       READ  ( numnam_ref, namhsb, IOSTAT = ios, ERR = 901) 
    218 901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namhsb in reference namelist', lwp ) 
    219  
    220       REWIND( numnam_cfg )              ! Namelist namhsb in configuration namelist 
    221       READ  ( numnam_cfg, namhsb, IOSTAT = ios, ERR = 902 ) 
    222 902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namhsb in configuration namelist', lwp ) 
    223       IF(lwm) WRITE ( numond, namhsb ) 
    224  
    225       ! 
    226       IF(lwp) THEN                   ! Control print 
    227          WRITE(numout,*) 
    228          WRITE(numout,*) 'dia_hsb_init : check the heat and salt budgets' 
    229          WRITE(numout,*) '~~~~~~~~~~~~' 
    230          WRITE(numout,*) '   Namelist namhsb : set hsb parameters' 
    231          WRITE(numout,*) '      Switch for hsb diagnostic (T) or not (F)  ln_diahsb  = ', ln_diahsb 
    232          WRITE(numout,*) 
    233       ENDIF 
    234  
    235       IF( .NOT. ln_diahsb )   RETURN 
    236 !      IF( .NOT. lk_mpp_rep ) & 
    237 !        CALL ctl_stop (' Your global mpp_sum if performed in single precision - 64 bits -', & 
    238 !             &         ' whereas the global sum to be precise must be done in double precision ',& 
    239 !             &         ' please add key_mpp_rep') 
    240  
    241       ! ------------------- ! 
    242       ! 1 - Allocate memory ! 
    243       ! ------------------- ! 
    244       ALLOCATE( hc_loc_ini(jpi,jpj,jpk), sc_loc_ini(jpi,jpj,jpk), & 
    245          &      e3t_ini(jpi,jpj,jpk), surf(jpi,jpj),  ssh_ini(jpi,jpj), STAT=ierror ) 
    246       IF( ierror > 0 ) THEN 
    247          CALL ctl_stop( 'dia_hsb: unable to allocate hc_loc_ini' )   ;   RETURN 
    248       ENDIF 
    249  
    250       IF(.NOT. lk_vvl ) ALLOCATE( ssh_hc_loc_ini(jpi,jpj), ssh_sc_loc_ini(jpi,jpj),STAT=ierror ) 
    251       IF( ierror > 0 ) THEN 
    252          CALL ctl_stop( 'dia_hsb: unable to allocate hc_loc_ini' )   ;   RETURN 
    253       ENDIF 
    254  
    255       ! ----------------------------------------------- ! 
    256       ! 2 - Time independant variables and file opening ! 
    257       ! ----------------------------------------------- ! 
    258       IF(lwp) WRITE(numout,*) "dia_hsb: heat salt volume budgets activated" 
    259       IF(lwp) WRITE(numout,*) '~~~~~~~' 
    260       surf(:,:) = e1t(:,:) * e2t(:,:) * tmask(:,:,1) * tmask_i(:,:)      ! masked surface grid cell area 
    261       surf_tot  = glob_sum( surf(:,:) )                                       ! total ocean surface area 
    262  
    263       IF( lk_bdy ) CALL ctl_warn( 'dia_hsb does not take open boundary fluxes into account' )          
    264       ! 
    265       ! ---------------------------------- ! 
    266       ! 4 - initial conservation variables ! 
    267       ! ---------------------------------- ! 
    268       CALL dia_hsb_rst( nit000, 'READ' )  !* read or initialize all required files 
    269       ! 
    270    END SUBROUTINE dia_hsb_init 
    271235 
    272236   SUBROUTINE dia_hsb_rst( kt, cdrw ) 
    … …  
    281245     CHARACTER(len=*), INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag 
    282246     ! 
    283      INTEGER ::   jk   !  
    284      INTEGER ::   id1   ! local integers 
     247     INTEGER ::   ji, jj, jk   ! dummy loop indices 
     248     INTEGER ::   id1          ! local integers 
    285249     !!---------------------------------------------------------------------- 
    286250     ! 
    … …  
    317281             sc_loc_ini(:,:,jk) = tsn(:,:,jk,jp_sal) * fse3t_n(:,:,jk)   ! initial salt content 
    318282          END DO 
    319           frc_v = 0.d0                                           ! volume       trend due to forcing 
    320           frc_t = 0.d0                                           ! heat content   -    -   -    -    
    321           frc_s = 0.d0                                           ! salt content   -    -   -    -         
     283          frc_v = 0._wp                                           ! volume       trend due to forcing 
     284          frc_t = 0._wp                                           ! heat content   -    -   -    -    
     285          frc_s = 0._wp                                           ! salt content   -    -   -    -         
    322286          IF( .NOT. lk_vvl ) THEN 
    323              ssh_hc_loc_ini(:,:) = tsn(:,:,1,jp_tem) * sshn(:,:)   ! initial heat content in ssh 
    324              ssh_sc_loc_ini(:,:) = tsn(:,:,1,jp_sal) * sshn(:,:)   ! initial salt content in ssh 
    325              frc_wn_t = 0.d0                                       ! initial heat content misfit due to free surface 
    326              frc_wn_s = 0.d0                                       ! initial salt content misfit due to free surface 
     287             IF ( ln_isfcav ) THEN 
     288                DO ji=1,jpi 
     289                   DO jj=1,jpj 
     290                      ssh_hc_loc_ini(ji,jj) = tsn(ji,jj,mikt(ji,jj),jp_tem) * sshn(ji,jj)   ! initial heat content in ssh 
     291                      ssh_sc_loc_ini(ji,jj) = tsn(ji,jj,mikt(ji,jj),jp_sal) * sshn(ji,jj)   ! initial salt content in ssh 
     292                   ENDDO 
     293                ENDDO 
     294             ELSE 
     295                ssh_hc_loc_ini(:,:) = tsn(:,:,1,jp_tem) * sshn(:,:)   ! initial heat content in ssh 
     296                ssh_sc_loc_ini(:,:) = tsn(:,:,1,jp_sal) * sshn(:,:)   ! initial salt content in ssh 
     297             END IF 
     298             frc_wn_t = 0._wp                                       ! initial heat content misfit due to free surface 
     299             frc_wn_s = 0._wp                                       ! initial salt content misfit due to free surface 
    327300          ENDIF 
    328301       ENDIF 
    … …  
    354327   END SUBROUTINE dia_hsb_rst 
    355328 
     329 
     330   SUBROUTINE dia_hsb_init 
     331      !!--------------------------------------------------------------------------- 
     332      !!                  ***  ROUTINE dia_hsb  *** 
     333      !!      
     334      !! ** Purpose: Initialization for the heat salt volume budgets 
     335      !!  
     336      !! ** Method : Compute initial heat content, salt content and volume 
     337      !! 
     338      !! ** Action : - Compute initial heat content, salt content and volume 
     339      !!             - Initialize forcing trends 
     340      !!             - Compute coefficients for conversion 
     341      !!--------------------------------------------------------------------------- 
     342      INTEGER ::   jk       ! dummy loop indice 
     343      INTEGER ::   ierror   ! local integer 
     344      INTEGER ::   ios 
     345      ! 
     346      NAMELIST/namhsb/ ln_diahsb 
     347      !!---------------------------------------------------------------------- 
     348 
     349      IF(lwp) THEN 
     350         WRITE(numout,*) 
     351         WRITE(numout,*) 'dia_hsb_init : check the heat and salt budgets' 
     352         WRITE(numout,*) '~~~~~~~~ ' 
     353      ENDIF 
     354 
     355      REWIND( numnam_ref )              ! Namelist namhsb in reference namelist 
     356      READ  ( numnam_ref, namhsb, IOSTAT = ios, ERR = 901) 
     357901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namhsb in reference namelist', lwp ) 
     358 
     359      REWIND( numnam_cfg )              ! Namelist namhsb in configuration namelist 
     360      READ  ( numnam_cfg, namhsb, IOSTAT = ios, ERR = 902 ) 
     361902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namhsb in configuration namelist', lwp ) 
     362      IF(lwm) WRITE ( numond, namhsb ) 
     363 
     364      ! 
     365      IF(lwp) THEN                   ! Control print 
     366         WRITE(numout,*) 
     367         WRITE(numout,*) 'dia_hsb_init : check the heat and salt budgets' 
     368         WRITE(numout,*) '~~~~~~~~~~~~' 
     369         WRITE(numout,*) '   Namelist namhsb : set hsb parameters' 
     370         WRITE(numout,*) '      Switch for hsb diagnostic (T) or not (F)  ln_diahsb  = ', ln_diahsb 
     371         WRITE(numout,*) 
     372      ENDIF 
     373 
     374      IF( .NOT. ln_diahsb )   RETURN 
     375         !      IF( .NOT. lk_mpp_rep ) & 
     376         !        CALL ctl_stop (' Your global mpp_sum if performed in single precision - 64 bits -', & 
     377         !             &         ' whereas the global sum to be precise must be done in double precision ',& 
     378         !             &         ' please add key_mpp_rep') 
     379 
     380      ! ------------------- ! 
     381      ! 1 - Allocate memory ! 
     382      ! ------------------- ! 
     383      ALLOCATE( hc_loc_ini(jpi,jpj,jpk), sc_loc_ini(jpi,jpj,jpk), & 
     384         &      e3t_ini(jpi,jpj,jpk), surf(jpi,jpj),  ssh_ini(jpi,jpj), STAT=ierror ) 
     385      IF( ierror > 0 ) THEN 
     386         CALL ctl_stop( 'dia_hsb: unable to allocate hc_loc_ini' )   ;   RETURN 
     387      ENDIF 
     388 
     389      IF(.NOT. lk_vvl ) ALLOCATE( ssh_hc_loc_ini(jpi,jpj), ssh_sc_loc_ini(jpi,jpj),STAT=ierror ) 
     390      IF( ierror > 0 ) THEN 
     391         CALL ctl_stop( 'dia_hsb: unable to allocate hc_loc_ini' )   ;   RETURN 
     392      ENDIF 
     393 
     394      ! ----------------------------------------------- ! 
     395      ! 2 - Time independant variables and file opening ! 
     396      ! ----------------------------------------------- ! 
     397      IF(lwp) WRITE(numout,*) "dia_hsb: heat salt volume budgets activated" 
     398      IF(lwp) WRITE(numout,*) '~~~~~~~' 
     399      surf(:,:) = e1t(:,:) * e2t(:,:) * tmask_i(:,:)      ! masked surface grid cell area 
     400      surf_tot  = glob_sum( surf(:,:) )                                       ! total ocean surface area 
     401 
     402      IF( lk_bdy ) CALL ctl_warn( 'dia_hsb does not take open boundary fluxes into account' )          
     403      ! 
     404      ! ---------------------------------- ! 
     405      ! 4 - initial conservation variables ! 
     406      ! ---------------------------------- ! 
     407      CALL dia_hsb_rst( nit000, 'READ' )  !* read or initialize all required files 
     408      ! 
     409   END SUBROUTINE dia_hsb_init 
     410 
    356411   !!====================================================================== 
    357412END MODULE diahsb 

  • branches/2014/dev_r4650_UKMO10_Tidally_Meaned_Diagnostics/NEMOGCM/NEMO/OPA_SRC/DIA/diainsitutem.F90

    r4756 r5260  
    1 MODULE insitu_tem 
     1MODULE diainsitutem 
    22   !!====================================================================== 
    33   !!                       ***  MODULE  diaharm  *** 
    … …  
    122122   END SUBROUTINE ATG 
    123123 
    124 END MODULE insitu_tem 
     124END MODULE diainsitutem 

  • branches/2014/dev_r4650_UKMO10_Tidally_Meaned_Diagnostics/NEMOGCM/NEMO/OPA_SRC/DIA/diaptr.F90

    r4624 r5260  
    88   !!            3.2  ! 2010-03  (O. Marti, S. Flavoni) Add fields 
    99   !!            3.3  ! 2010-10  (G. Madec)  dynamical allocation 
     10   !!            3.6  ! 2014-12  (C. Ethe) use of IOM 
    1011   !!---------------------------------------------------------------------- 
    1112 
    … …  
    1314   !!   dia_ptr      : Poleward Transport Diagnostics module 
    1415   !!   dia_ptr_init : Initialization, namelist read 
    15    !!   dia_ptr_wri  : Output of poleward fluxes 
    16    !!   ptr_vjk      : "zonal" sum computation of a "meridional" flux array 
    17    !!   ptr_tjk      : "zonal" mean computation of a tracer field 
    18    !!   ptr_vj       : "zonal" and vertical sum computation of a "meridional" flux array 
    19    !!                   (Generic interface to ptr_vj_3d, ptr_vj_2d) 
     16   !!   ptr_sjk      : "zonal" mean computation of a field - tracer or flux array 
     17   !!   ptr_sj       : "zonal" and vertical sum computation of a "meridional" flux array 
     18   !!                   (Generic interface to ptr_sj_3d, ptr_sj_2d) 
    2019   !!---------------------------------------------------------------------- 
    2120   USE oce              ! ocean dynamics and active tracers 
    2221   USE dom_oce          ! ocean space and time domain 
    2322   USE phycst           ! physical constants 
    24    USE ldftra_oce       ! ocean active tracers: lateral physics 
    25    USE dianam           ! 
     23   ! 
    2624   USE iom              ! IOM library 
    27    USE ioipsl           ! IO-IPSL library 
    2825   USE in_out_manager   ! I/O manager 
    2926   USE lib_mpp          ! MPP library 
    30    USE lbclnk           ! lateral boundary condition - processor exchanges 
    3127   USE timing           ! preformance summary 
    32    USE wrk_nemo         ! working arrays 
    3328 
    3429   IMPLICIT NONE 
    3530   PRIVATE 
    3631 
    37    INTERFACE ptr_vj 
    38       MODULE PROCEDURE ptr_vj_3d, ptr_vj_2d 
     32   INTERFACE ptr_sj 
     33      MODULE PROCEDURE ptr_sj_3d, ptr_sj_2d 
    3934   END INTERFACE 
    4035 
    41    PUBLIC   dia_ptr_init   ! call in opa module 
     36   PUBLIC   ptr_sj         ! call by tra_ldf & tra_adv routines 
     37   PUBLIC   ptr_sjk        !  
     38   PUBLIC   dia_ptr_init   ! call in step module 
    4239   PUBLIC   dia_ptr        ! call in step module 
    43    PUBLIC   ptr_vj         ! call by tra_ldf & tra_adv routines 
    44    PUBLIC   ptr_vjk        ! call by tra_ldf & tra_adv routines 
    4540 
    4641   !                                  !!** namelist  namptr  ** 
    47    LOGICAL , PUBLIC ::   ln_diaptr     !: Poleward transport flag (T) or not (F) 
    48    LOGICAL , PUBLIC ::   ln_subbas     !: Atlantic/Pacific/Indian basins calculation 
    49    LOGICAL , PUBLIC ::   ln_diaznl     !: Add zonal means and meridional stream functions 
    50    LOGICAL , PUBLIC ::   ln_ptrcomp    !: Add decomposition : overturning (and gyre, soon ...) 
    51    INTEGER , PUBLIC ::   nn_fptr       !: frequency of ptr computation  [time step] 
    52    INTEGER , PUBLIC ::   nn_fwri       !: frequency of ptr outputs      [time step] 
    53  
    54    REAL(wp), ALLOCATABLE, SAVE, PUBLIC, DIMENSION(:) ::   htr_adv, htr_ldf, htr_ove   !: Heat TRansports (adv, diff, overturn.) 
    55    REAL(wp), ALLOCATABLE, SAVE, PUBLIC, DIMENSION(:) ::   str_adv, str_ldf, str_ove   !: Salt TRansports (adv, diff, overturn.) 
     42   REAL(wp), ALLOCATABLE, SAVE, PUBLIC, DIMENSION(:) ::   htr_adv, htr_ldf   !: Heat TRansports (adv, diff, overturn.) 
     43   REAL(wp), ALLOCATABLE, SAVE, PUBLIC, DIMENSION(:) ::   str_adv, str_ldf   !: Salt TRansports (adv, diff, overturn.) 
    5644    
    57    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   btmsk                  ! T-point basin interior masks 
    58    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   btm30                  ! mask out Southern Ocean (=0 south of 30°S) 
    59    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   htr  , str             ! adv heat and salt transports (approx) 
    60    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   tn_jk, sn_jk , v_msf   ! i-mean T and S, j-Stream-Function 
    61    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   sjk  , r1_sjk          ! i-mean i-k-surface and its inverse         
    62    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   htr_eiv, str_eiv       ! bolus adv heat ans salt transports ('key_diaeiv') 
    63    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   v_msf_eiv              ! bolus j-streamfuction              ('key_diaeiv') 
    64  
    65  
    66    INTEGER ::   niter       ! 
    67    INTEGER ::   nidom_ptr   ! 
    68    INTEGER ::   numptr      ! logical unit for Poleward TRansports 
    69    INTEGER ::   nptr        ! = 1 (ln_subbas=F) or = 5 (glo, atl, pac, ind, ipc) (ln_subbas=T)  
     45 
     46   LOGICAL, PUBLIC ::   ln_diaptr   !  Poleward transport flag (T) or not (F) 
     47   LOGICAL, PUBLIC ::   ln_subbas   !  Atlantic/Pacific/Indian basins calculation 
     48   INTEGER         ::   nptr        ! = 1 (l_subbas=F) or = 5 (glo, atl, pac, ind, ipc) (l_subbas=T)  
    7049 
    7150   REAL(wp) ::   rc_sv    = 1.e-6_wp   ! conversion from m3/s to Sverdrup 
    … …  
    7352   REAL(wp) ::   rc_ggram = 1.e-6_wp   ! conversion from g    to Pg 
    7453 
    75    REAL(wp), TARGET, DIMENSION(:),   ALLOCATABLE, SAVE :: p_fval1d 
    76    REAL(wp), TARGET, DIMENSION(:,:), ALLOCATABLE, SAVE :: p_fval2d 
    77  
    78    !! Integer, 1D workspace arrays. Not common enough to be implemented in  
    79    !! wrk_nemo module. 
    80    INTEGER, ALLOCATABLE, SAVE, DIMENSION(:) :: ndex  , ndex_atl     , ndex_pac     , ndex_ind     , ndex_ipc 
    81    INTEGER, ALLOCATABLE, SAVE, DIMENSION(:) ::         ndex_atl_30  , ndex_pac_30  , ndex_ind_30  , ndex_ipc_30 
    82    INTEGER, ALLOCATABLE, SAVE, DIMENSION(:) :: ndex_h, ndex_h_atl_30, ndex_h_pac_30, ndex_h_ind_30, ndex_h_ipc_30 
     54   CHARACTER(len=3), ALLOCATABLE, SAVE, DIMENSION(:)     :: clsubb 
     55   REAL(wp),         ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: btmsk   ! T-point basin interior masks 
     56   REAL(wp),         ALLOCATABLE, SAVE, DIMENSION(:,:)   :: btm30   ! mask out Southern Ocean (=0 south of 30°S) 
     57 
     58   REAL(wp), TARGET, ALLOCATABLE, SAVE, DIMENSION(:)     :: p_fval1d 
     59   REAL(wp), TARGET, ALLOCATABLE, SAVE, DIMENSION(:,:)   :: p_fval2d 
     60 
    8361 
    8462   !! * Substitutions 
    … …  
    9270CONTAINS 
    9371 
     72   SUBROUTINE dia_ptr( pvtr ) 
     73      !!---------------------------------------------------------------------- 
     74      !!                  ***  ROUTINE dia_ptr  *** 
     75      !!---------------------------------------------------------------------- 
     76      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in), OPTIONAL ::   pvtr   ! j-effective transport 
     77      ! 
     78      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices 
     79      REAL(wp) ::   zv, zsfc               ! local scalar 
     80      REAL(wp), DIMENSION(jpi,jpj)     ::  z2d   ! 2D workspace 
     81      REAL(wp), DIMENSION(jpi,jpj,jpk) ::  z3d   ! 3D workspace 
     82      REAL(wp), DIMENSION(jpi,jpj,jpk) ::  zmask   ! 3D workspace 
     83      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts) ::  zts   ! 3D workspace 
     84      CHARACTER( len = 10 )  :: cl1 
     85      !!---------------------------------------------------------------------- 
     86      ! 
     87      IF( nn_timing == 1 )   CALL timing_start('dia_ptr') 
     88 
     89      ! 
     90      IF( PRESENT( pvtr ) ) THEN 
     91         IF( iom_use("zomsfglo") ) THEN    ! effective MSF 
     92            z3d(1,:,:) = ptr_sjk( pvtr(:,:,:) )  ! zonal cumulative effective transport 
     93            DO jk = 2, jpkm1  
     94              z3d(1,:,jk) = z3d(1,:,jk-1) + z3d(1,:,jk)   ! effective j-Stream-Function (MSF) 
     95            END DO 
     96            DO ji = 1, jpi 
     97               z3d(ji,:,:) = z3d(1,:,:) 
     98            ENDDO 
     99            cl1 = TRIM('zomsf'//clsubb(1) ) 
     100            CALL iom_put( cl1, z3d * rc_sv ) 
     101            DO jn = 2, nptr                                    ! by sub-basins 
     102               z3d(1,:,:) =  ptr_sjk( pvtr(:,:,:), btmsk(:,:,jn)*btm30(:,:) )  
     103               DO jk = 2, jpkm1  
     104                  z3d(1,:,jk) = z3d(1,:,jk-1) + z3d(1,:,jk)    ! effective j-Stream-Function (MSF) 
     105               END DO 
     106               DO ji = 1, jpi 
     107                  z3d(ji,:,:) = z3d(1,:,:) 
     108               ENDDO 
     109               cl1 = TRIM('zomsf'//clsubb(jn) ) 
     110               CALL iom_put( cl1, z3d * rc_sv ) 
     111            END DO 
     112         ENDIF 
     113         ! 
     114      ELSE 
     115         ! 
     116         IF( iom_use("zotemglo") ) THEN    ! i-mean i-k-surface  
     117            DO jk = 1, jpkm1 
     118               DO jj = 1, jpj 
     119                  DO ji = 1, jpi 
     120                     zsfc = e1t(ji,jj) * fse3t(ji,jj,jk) 
     121                     zmask(ji,jj,jk)      = tmask(ji,jj,jk)      * zsfc 
     122                     zts(ji,jj,jk,jp_tem) = tsn(ji,jj,jk,jp_tem) * zsfc 
     123                     zts(ji,jj,jk,jp_sal) = tsn(ji,jj,jk,jp_sal) * zsfc 
     124                  ENDDO 
     125               ENDDO 
     126            ENDDO 
     127            DO jn = 1, nptr 
     128               zmask(1,:,:) = ptr_sjk( zmask(:,:,:), btmsk(:,:,jn) ) 
     129               cl1 = TRIM('zosrf'//clsubb(jn) ) 
     130               CALL iom_put( cl1, zmask ) 
     131               ! 
     132               z3d(1,:,:) = ptr_sjk( zts(:,:,:,jp_tem), btmsk(:,:,jn) ) & 
     133                  &            / MAX( zmask(1,:,:), 10.e-15 ) 
     134               DO ji = 1, jpi 
     135                  z3d(ji,:,:) = z3d(1,:,:) 
     136               ENDDO 
     137               cl1 = TRIM('zotem'//clsubb(jn) ) 
     138               CALL iom_put( cl1, z3d ) 
     139               ! 
     140               z3d(1,:,:) = ptr_sjk( zts(:,:,:,jp_sal), btmsk(:,:,jn) ) & 
     141                  &            / MAX( zmask(1,:,:), 10.e-15 ) 
     142               DO ji = 1, jpi 
     143                  z3d(ji,:,:) = z3d(1,:,:) 
     144               ENDDO 
     145               cl1 = TRIM('zosal'//clsubb(jn) ) 
     146               CALL iom_put( cl1, z3d ) 
     147            END DO 
     148         ENDIF 
     149         ! 
     150         !                                ! Advective and diffusive heat and salt transport 
     151         IF( iom_use("sophtadv") .OR. iom_use("sopstadv") ) THEN    
     152            z2d(1,:) = htr_adv(:) * rc_pwatt        !  (conversion in PW) 
     153            DO ji = 1, jpi 
     154               z2d(ji,:) = z2d(1,:) 
     155            ENDDO 
     156            cl1 = 'sophtadv'                  
     157            CALL iom_put( TRIM(cl1), z2d ) 
     158            z2d(1,:) = str_adv(:) * rc_ggram        ! (conversion in Gg) 
     159            DO ji = 1, jpi 
     160               z2d(ji,:) = z2d(1,:) 
     161            ENDDO 
     162            cl1 = 'sopstadv' 
     163            CALL iom_put( TRIM(cl1), z2d ) 
     164         ENDIF 
     165         ! 
     166         IF( iom_use("sophtldf") .OR. iom_use("sopstldf") ) THEN    
     167            z2d(1,:) = htr_ldf(:) * rc_pwatt        !  (conversion in PW)  
     168            DO ji = 1, jpi 
     169               z2d(ji,:) = z2d(1,:) 
     170            ENDDO 
     171            cl1 = 'sophtldf' 
     172            CALL iom_put( TRIM(cl1), z2d ) 
     173            z2d(1,:) = str_ldf(:) * rc_ggram        !  (conversion in Gg) 
     174            DO ji = 1, jpi 
     175               z2d(ji,:) = z2d(1,:) 
     176            ENDDO 
     177            cl1 = 'sopstldf' 
     178            CALL iom_put( TRIM(cl1), z2d ) 
     179         ENDIF 
     180         ! 
     181      ENDIF 
     182      ! 
     183      IF( nn_timing == 1 )   CALL timing_stop('dia_ptr') 
     184      ! 
     185   END SUBROUTINE dia_ptr 
     186 
     187 
     188   SUBROUTINE dia_ptr_init 
     189      !!---------------------------------------------------------------------- 
     190      !!                  ***  ROUTINE dia_ptr_init  *** 
     191      !!                    
     192      !! ** Purpose :   Initialization, namelist read 
     193      !!---------------------------------------------------------------------- 
     194      INTEGER ::  jn           ! local integers 
     195      INTEGER ::  inum, ierr   ! local integers 
     196      INTEGER ::  ios          ! Local integer output status for namelist read 
     197      !! 
     198      NAMELIST/namptr/ ln_diaptr, ln_subbas 
     199      !!---------------------------------------------------------------------- 
     200 
     201      REWIND( numnam_ref )              ! Namelist namptr in reference namelist : Poleward transport 
     202      READ  ( numnam_ref, namptr, IOSTAT = ios, ERR = 901) 
     203901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namptr in reference namelist', lwp ) 
     204 
     205      REWIND( numnam_cfg )              ! Namelist namptr in configuration namelist : Poleward transport 
     206      READ  ( numnam_cfg, namptr, IOSTAT = ios, ERR = 902 ) 
     207902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namptr in configuration namelist', lwp ) 
     208      IF(lwm) WRITE ( numond, namptr ) 
     209 
     210      IF(lwp) THEN                     ! Control print 
     211         WRITE(numout,*) 
     212         WRITE(numout,*) 'dia_ptr_init : poleward transport and msf initialization' 
     213         WRITE(numout,*) '~~~~~~~~~~~~' 
     214         WRITE(numout,*) '   Namelist namptr : set ptr parameters' 
     215         WRITE(numout,*) '      Poleward heat & salt transport (T) or not (F)      ln_diaptr  = ', ln_diaptr 
     216         WRITE(numout,*) '      Global (F) or glo/Atl/Pac/Ind/Indo-Pac basins      ln_subbas  = ', ln_subbas 
     217      ENDIF 
     218 
     219      IF( ln_diaptr ) THEN   
     220         ! 
     221         IF( ln_subbas ) THEN  
     222            nptr = 5            ! Global, Atlantic, Pacific, Indian, Indo-Pacific 
     223            ALLOCATE( clsubb(nptr) ) 
     224            clsubb(1) = 'glo' ;  clsubb(2) = 'atl'  ;  clsubb(3) = 'pac'  ;  clsubb(4) = 'ind'  ;  clsubb(5) = 'ipc' 
     225         ELSE                
     226            nptr = 1       ! Global only 
     227            ALLOCATE( clsubb(nptr) ) 
     228            clsubb(1) = 'glo'  
     229         ENDIF 
     230 
     231         !                                      ! allocate dia_ptr arrays 
     232         IF( dia_ptr_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'dia_ptr_init : unable to allocate arrays' ) 
     233 
     234         rc_pwatt = rc_pwatt * rau0_rcp          ! conversion from K.s-1 to PetaWatt 
     235 
     236         IF( lk_mpp )   CALL mpp_ini_znl( numout )     ! Define MPI communicator for zonal sum 
     237 
     238         IF( ln_subbas ) THEN                ! load sub-basin mask 
     239            CALL iom_open( 'subbasins', inum,  ldstop = .FALSE.  ) 
     240            CALL iom_get( inum, jpdom_data, 'atlmsk', btmsk(:,:,2) )   ! Atlantic basin 
     241            CALL iom_get( inum, jpdom_data, 'pacmsk', btmsk(:,:,3) )   ! Pacific  basin 
     242            CALL iom_get( inum, jpdom_data, 'indmsk', btmsk(:,:,4) )   ! Indian   basin 
     243            CALL iom_close( inum ) 
     244            btmsk(:,:,5) = MAX ( btmsk(:,:,3), btmsk(:,:,4) )          ! Indo-Pacific basin 
     245            WHERE( gphit(:,:) < -30._wp)   ;   btm30(:,:) = 0._wp      ! mask out Southern Ocean 
     246            ELSE WHERE                     ;   btm30(:,:) = ssmask(:,:) 
     247            END WHERE 
     248         ENDIF 
     249    
     250         btmsk(:,:,1) = tmask_i(:,:)                                   ! global ocean 
     251       
     252         DO jn = 1, nptr 
     253            btmsk(:,:,jn) = btmsk(:,:,jn) * tmask_i(:,:)               ! interior domain only 
     254         END DO 
     255 
     256         ! Initialise arrays to zero because diatpr is called before they are first calculated 
     257         ! Note that this means diagnostics will not be exactly correct when model run is restarted. 
     258         htr_adv(:) = 0._wp  ;  str_adv(:) =  0._wp   
     259         htr_ldf(:) = 0._wp  ;  str_ldf(:) =  0._wp  
     260         ! 
     261      ENDIF  
     262      !  
     263   END SUBROUTINE dia_ptr_init 
     264 
     265 
    94266   FUNCTION dia_ptr_alloc() 
    95267      !!---------------------------------------------------------------------- 
    … …  
    97269      !!---------------------------------------------------------------------- 
    98270      INTEGER               ::   dia_ptr_alloc   ! return value 
    99       INTEGER, DIMENSION(6) ::   ierr 
     271      INTEGER, DIMENSION(3) ::   ierr 
    100272      !!---------------------------------------------------------------------- 
    101273      ierr(:) = 0 
    … …  
    103275      ALLOCATE( btmsk(jpi,jpj,nptr) ,           & 
    104276         &      htr_adv(jpj) , str_adv(jpj) ,   & 
    105          &      htr_ldf(jpj) , str_ldf(jpj) ,   & 
    106          &      htr_ove(jpj) , str_ove(jpj),    & 
    107          &      htr(jpj,nptr) , str(jpj,nptr) , & 
    108          &      tn_jk(jpj,jpk,nptr) , sn_jk (jpj,jpk,nptr) , v_msf(jpj,jpk,nptr) , & 
    109          &      sjk  (jpj,jpk,nptr) , r1_sjk(jpj,jpk,nptr) , STAT=ierr(1)  ) 
    110          ! 
    111 #if defined key_diaeiv 
    112       ALLOCATE( htr_eiv(jpj,nptr) , str_eiv(jpj,nptr) , & 
    113          &      v_msf_eiv(jpj,jpk,nptr) , STAT=ierr(2) ) 
    114 #endif 
    115       ALLOCATE( p_fval1d(jpj), p_fval2d(jpj,jpk), Stat=ierr(3)) 
    116       ! 
    117       ALLOCATE(ndex(jpj*jpk),        ndex_atl(jpj*jpk), ndex_pac(jpj*jpk), & 
    118          &     ndex_ind(jpj*jpk),    ndex_ipc(jpj*jpk),                    & 
    119          &     ndex_atl_30(jpj*jpk), ndex_pac_30(jpj*jpk), Stat=ierr(4)) 
    120  
    121       ALLOCATE(ndex_ind_30(jpj*jpk), ndex_ipc_30(jpj*jpk),                   & 
    122          &     ndex_h(jpj),          ndex_h_atl_30(jpj), ndex_h_pac_30(jpj), & 
    123          &     ndex_h_ind_30(jpj),   ndex_h_ipc_30(jpj), Stat=ierr(5) ) 
    124          ! 
    125      ALLOCATE( btm30(jpi,jpj) , STAT=ierr(6)  ) 
     277         &      htr_ldf(jpj) , str_ldf(jpj) , STAT=ierr(1)  ) 
     278         ! 
     279      ALLOCATE( p_fval1d(jpj), p_fval2d(jpj,jpk), Stat=ierr(2)) 
     280      ! 
     281      ALLOCATE( btm30(jpi,jpj), STAT=ierr(3)  ) 
     282 
    126283         ! 
    127284      dia_ptr_alloc = MAXVAL( ierr ) 
    … …  
    131288 
    132289 
    133    FUNCTION ptr_vj_3d( pva )   RESULT ( p_fval ) 
    134       !!---------------------------------------------------------------------- 
    135       !!                    ***  ROUTINE ptr_vj_3d  *** 
     290   FUNCTION ptr_sj_3d( pva, pmsk )   RESULT ( p_fval ) 
     291      !!---------------------------------------------------------------------- 
     292      !!                    ***  ROUTINE ptr_sj_3d  *** 
    136293      !! 
    137294      !! ** Purpose :   i-k sum computation of a j-flux array 
    … …  
    142299      !! ** Action  : - p_fval: i-k-mean poleward flux of pva 
    143300      !!---------------------------------------------------------------------- 
    144       REAL(wp) , INTENT(in), DIMENSION(jpi,jpj,jpk) ::   pva   ! mask flux array at V-point 
    145       !! 
     301      REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk)       ::   pva   ! mask flux array at V-point 
     302      REAL(wp), INTENT(in), DIMENSION(jpi,jpj), OPTIONAL ::   pmsk   ! Optional 2D basin mask 
     303      ! 
    146304      INTEGER                  ::   ji, jj, jk   ! dummy loop arguments 
    147305      INTEGER                  ::   ijpj         ! ??? 
    … …  
    153311      ijpj = jpj 
    154312      p_fval(:) = 0._wp 
    155       DO jk = 1, jpkm1 
    156          DO jj = 2, jpjm1 
    157             DO ji = fs_2, fs_jpim1   ! Vector opt. 
    158                p_fval(jj) = p_fval(jj) + pva(ji,jj,jk) * tmask_i(ji,jj)  
    159             END DO 
    160          END DO 
    161       END DO 
     313      IF( PRESENT( pmsk ) ) THEN  
     314         DO jk = 1, jpkm1 
     315            DO jj = 2, jpjm1 
     316               DO ji = fs_2, fs_jpim1   ! Vector opt. 
     317                  p_fval(jj) = p_fval(jj) + pva(ji,jj,jk) * tmask_i(ji,jj) * pmsk(ji,jj) 
     318               END DO 
     319            END DO 
     320         END DO 
     321      ELSE 
     322         DO jk = 1, jpkm1 
     323            DO jj = 2, jpjm1 
     324               DO ji = fs_2, fs_jpim1   ! Vector opt. 
     325                  p_fval(jj) = p_fval(jj) + pva(ji,jj,jk) * tmask_i(ji,jj)  
     326               END DO 
     327            END DO 
     328         END DO 
     329      ENDIF 
    162330#if defined key_mpp_mpi 
    163331      IF(lk_mpp)   CALL mpp_sum( p_fval, ijpj, ncomm_znl) 
    164332#endif 
    165333      ! 
    166    END FUNCTION ptr_vj_3d 
    167  
    168  
    169    FUNCTION ptr_vj_2d( pva )   RESULT ( p_fval ) 
    170       !!---------------------------------------------------------------------- 
    171       !!                    ***  ROUTINE ptr_vj_2d  *** 
     334   END FUNCTION ptr_sj_3d 
     335 
     336 
     337   FUNCTION ptr_sj_2d( pva, pmsk )   RESULT ( p_fval ) 
     338      !!---------------------------------------------------------------------- 
     339      !!                    ***  ROUTINE ptr_sj_2d  *** 
    172340      !! 
    173341      !! ** Purpose :   "zonal" and vertical sum computation of a i-flux array 
    … …  
    178346      !! ** Action  : - p_fval: i-k-mean poleward flux of pva 
    179347      !!---------------------------------------------------------------------- 
    180       IMPLICIT none 
    181       REAL(wp) , INTENT(in), DIMENSION(jpi,jpj) ::   pva   ! mask flux array at V-point 
    182       !! 
     348      REAL(wp) , INTENT(in), DIMENSION(jpi,jpj)           ::   pva   ! mask flux array at V-point 
     349      REAL(wp) , INTENT(in), DIMENSION(jpi,jpj), OPTIONAL ::   pmsk   ! Optional 2D basin mask 
     350      ! 
    183351      INTEGER                  ::   ji,jj       ! dummy loop arguments 
    184352      INTEGER                  ::   ijpj        ! ??? 
    … …  
    190358      ijpj = jpj 
    191359      p_fval(:) = 0._wp 
    192       DO jj = 2, jpjm1 
    193          DO ji = nldi, nlei   ! No vector optimisation here. Better use a mask ? 
    194             p_fval(jj) = p_fval(jj) + pva(ji,jj) * tmask_i(ji,jj) 
    195          END DO 
    196       END DO 
     360      IF( PRESENT( pmsk ) ) THEN  
     361         DO jj = 2, jpjm1 
     362            DO ji = nldi, nlei   ! No vector optimisation here. Better use a mask ? 
     363               p_fval(jj) = p_fval(jj) + pva(ji,jj) * tmask_i(ji,jj) * pmsk(ji,jj) 
     364            END DO 
     365         END DO 
     366      ELSE 
     367         DO jj = 2, jpjm1 
     368            DO ji = nldi, nlei   ! No vector optimisation here. Better use a mask ? 
     369               p_fval(jj) = p_fval(jj) + pva(ji,jj) * tmask_i(ji,jj) 
     370            END DO 
     371         END DO 
     372      ENDIF 
    197373#if defined key_mpp_mpi 
    198374      CALL mpp_sum( p_fval, ijpj, ncomm_znl ) 
    199375#endif 
    200376      !  
    201    END FUNCTION ptr_vj_2d 
    202  
    203  
    204    FUNCTION ptr_vjk( pva, pmsk )   RESULT ( p_fval ) 
    205       !!---------------------------------------------------------------------- 
    206       !!                    ***  ROUTINE ptr_vjk  *** 
    207       !! 
    208       !! ** Purpose :   i-sum computation of a j-velocity array 
     377   END FUNCTION ptr_sj_2d 
     378 
     379 
     380   FUNCTION ptr_sjk( pta, pmsk )   RESULT ( p_fval ) 
     381      !!---------------------------------------------------------------------- 
     382      !!                    ***  ROUTINE ptr_sjk  *** 
     383      !! 
     384      !! ** Purpose :   i-sum computation of an array 
    209385      !! 
    210386      !! ** Method  : - i-sum of pva using the interior 2D vmask (vmask_i). 
    211       !!              pva is supposed to be a masked flux (i.e. * vmask) 
    212387      !! 
    213388      !! ** Action  : - p_fval: i-mean poleward flux of pva 
    … …  
    215390      !! 
    216391      IMPLICIT none 
    217       REAL(wp) , INTENT(in), DIMENSION(jpi,jpj,jpk)           ::   pva    ! mask flux array at V-point 
     392      REAL(wp) , INTENT(in), DIMENSION(jpi,jpj,jpk)           ::   pta    ! mask flux array at V-point 
    218393      REAL(wp) , INTENT(in), DIMENSION(jpi,jpj)    , OPTIONAL ::   pmsk   ! Optional 2D basin mask 
    219394      !! 
    … …  
    224399      INTEGER, DIMENSION(2) ::   ish2 
    225400      INTEGER               ::   ijpjjpk 
    226 #endif 
    227 #if defined key_mpp_mpi 
    228       REAL(wp), POINTER, DIMENSION(:) ::   zwork    ! mask flux array at V-point 
     401      REAL(wp), DIMENSION(jpj*jpk) ::   zwork    ! mask flux array at V-point 
    229402#endif 
    230403      !!-------------------------------------------------------------------- 
    231404      ! 
    232 #if defined key_mpp_mpi 
    233       ijpjjpk = jpj*jpk 
    234       CALL wrk_alloc( jpj*jpk, zwork ) 
    235 #endif 
    236  
    237405      p_fval => p_fval2d 
    238406 
    … …  
    244412!!gm here, use of tmask_i  ==> no need of loop over nldi, nlei.... 
    245413               DO ji =  nldi, nlei   ! No vector optimisation here. Better use a mask ? 
    246                   p_fval(jj,jk) = p_fval(jj,jk) + pva(ji,jj,jk) * e1v(ji,jj) * fse3v(ji,jj,jk) * pmsk(ji,jj) 
     414                  p_fval(jj,jk) = p_fval(jj,jk) + pta(ji,jj,jk) * pmsk(ji,jj) 
    247415               END DO 
    248416            END DO 
    … …  
    252420            DO jj = 2, jpjm1 
    253421               DO ji =  nldi, nlei   ! No vector optimisation here. Better use a mask ? 
    254                   p_fval(jj,jk) = p_fval(jj,jk) + pva(ji,jj,jk) * e1v(ji,jj) * fse3v(ji,jj,jk) * tmask_i(ji,jj) 
     422                  p_fval(jj,jk) = p_fval(jj,jk) + pta(ji,jj,jk) * tmask_i(ji,jj) 
    255423               END DO 
    256424            END DO 
    … …  
    266434#endif 
    267435      ! 
    268 #if defined key_mpp_mpi 
    269       CALL wrk_dealloc( jpj*jpk, zwork ) 
    270 #endif 
    271       ! 
    272    END FUNCTION ptr_vjk 
    273  
    274  
    275    FUNCTION ptr_tjk( pta, pmsk )   RESULT ( p_fval ) 
    276       !!---------------------------------------------------------------------- 
    277       !!                    ***  ROUTINE ptr_tjk  *** 
    278       !! 
    279       !! ** Purpose :   i-sum computation of e1t*e3t * a tracer field 
    280       !! 
    281       !! ** Method  : - i-sum of mj(pta) using tmask 
    282       !! 
    283       !! ** Action  : - p_fval: i-sum of e1t*e3t*pta 
    284       !!---------------------------------------------------------------------- 
    285       !! 
    286       REAL(wp) , INTENT(in), DIMENSION(jpi,jpj,jpk) ::   pta    ! tracer flux array at T-point 
    287       REAL(wp) , INTENT(in), DIMENSION(jpi,jpj)     ::   pmsk   ! Optional 2D basin mask 
    288       !! 
    289       INTEGER                           :: ji, jj, jk   ! dummy loop arguments 
    290       REAL(wp), POINTER, DIMENSION(:,:) :: p_fval       ! return function value 
    291 #if defined key_mpp_mpi 
    292       INTEGER, DIMENSION(1) ::   ish 
    293       INTEGER, DIMENSION(2) ::   ish2 
    294       INTEGER               ::   ijpjjpk 
    295 #endif 
    296 #if defined key_mpp_mpi 
    297       REAL(wp), POINTER, DIMENSION(:) ::   zwork    ! mask flux array at V-point 
    298 #endif 
    299       !!--------------------------------------------------------------------  
    300       ! 
    301 #if defined key_mpp_mpi 
    302       ijpjjpk = jpj*jpk 
    303       CALL wrk_alloc( jpj*jpk, zwork ) 
    304 #endif 
    305  
    306       p_fval => p_fval2d 
    307  
    308       p_fval(:,:) = 0._wp 
    309       DO jk = 1, jpkm1 
    310          DO jj = 2, jpjm1 
    311             DO ji =  nldi, nlei   ! No vector optimisation here. Better use a mask ? 
    312                p_fval(jj,jk) = p_fval(jj,jk) + pta(ji,jj,jk) * e1t(ji,jj) * fse3t(ji,jj,jk) * pmsk(ji,jj) 
    313             END DO 
    314          END DO 
    315       END DO 
    316 #if defined key_mpp_mpi 
    317       ijpjjpk = jpj*jpk 
    318       ish(1) = jpj*jpk   ;   ish2(1) = jpj   ;   ish2(2) = jpk 
    319       zwork(1:ijpjjpk)= RESHAPE( p_fval, ish ) 
    320       CALL mpp_sum( zwork, ijpjjpk, ncomm_znl ) 
    321       p_fval(:,:)= RESHAPE( zwork, ish2 ) 
    322 #endif 
    323       ! 
    324 #if defined key_mpp_mpi 
    325       CALL wrk_dealloc( jpj*jpk, zwork ) 
    326 #endif 
    327       !     
    328    END FUNCTION ptr_tjk 
    329  
    330  
    331    SUBROUTINE dia_ptr( kt ) 
    332       !!---------------------------------------------------------------------- 
    333       !!                  ***  ROUTINE dia_ptr  *** 
    334       !!---------------------------------------------------------------------- 
    335       USE oce,     vt  =>   ua   ! use ua as workspace 
    336       USE oce,     vs  =>   va   ! use va as workspace 
    337       IMPLICIT none 
    338       !! 
    339       INTEGER, INTENT(in) ::   kt   ! ocean time step index 
    340       ! 
    341       INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices 
    342       REAL(wp) ::   zv               ! local scalar 
    343       !!---------------------------------------------------------------------- 
    344       ! 
    345       IF( nn_timing == 1 )   CALL timing_start('dia_ptr') 
    346       ! 
    347       IF( kt == nit000 .OR. MOD( kt, nn_fptr ) == 0 )   THEN 
    348          ! 
    349          IF( MOD( kt, nn_fptr ) == 0 ) THEN  
    350             ! 
    351             IF( ln_diaznl ) THEN               ! i-mean temperature and salinity 
    352                DO jn = 1, nptr 
    353                   tn_jk(:,:,jn) = ptr_tjk( tsn(:,:,:,jp_tem), btmsk(:,:,jn) ) * r1_sjk(:,:,jn) 
    354                   sn_jk(:,:,jn) = ptr_tjk( tsn(:,:,:,jp_sal), btmsk(:,:,jn) ) * r1_sjk(:,:,jn) 
    355                END DO 
    356             ENDIF 
    357             ! 
    358             !                          ! horizontal integral and vertical dz  
    359             !                                ! eulerian velocity 
    360             v_msf(:,:,1) = ptr_vjk( vn(:,:,:) )  
    361             DO jn = 2, nptr 
    362                v_msf(:,:,jn) = ptr_vjk( vn(:,:,:), btmsk(:,:,jn)*btm30(:,:) )  
    363             END DO 
    364 #if defined key_diaeiv 
    365             DO jn = 1, nptr                  ! bolus velocity 
    366                v_msf_eiv(:,:,jn) = ptr_vjk( v_eiv(:,:,:), btmsk(:,:,jn) )   ! here no btm30 for MSFeiv 
    367             END DO 
    368             !                                ! add bolus stream-function to the eulerian one 
    369             v_msf(:,:,:) = v_msf(:,:,:) + v_msf_eiv(:,:,:) 
    370 #endif 
    371             ! 
    372             !                          ! Transports 
    373             !                                ! local heat & salt transports at T-points  ( tsn*mj[vn+v_eiv] ) 
    374             vt(:,:,jpk) = 0._wp   ;   vs(:,:,jpk) = 0._wp 
    375             DO jk= 1, jpkm1 
    376                DO jj = 2, jpj 
    377                   DO ji = 1, jpi 
    378 #if defined key_diaeiv  
    379                      zv = ( vn(ji,jj,jk) + vn(ji,jj-1,jk) + v_eiv(ji,jj,jk) + v_eiv(ji,jj-1,jk) ) * 0.5_wp 
    380 #else 
    381                      zv = ( vn(ji,jj,jk) + vn(ji,jj-1,jk) ) * 0.5_wp 
    382 #endif  
    383                      vt(ji,jj,jk) = zv * tsn(ji,jj,jk,jp_tem) 
    384                      vs(ji,jj,jk) = zv * tsn(ji,jj,jk,jp_sal) 
    385                   END DO 
    386                END DO 
    387             END DO 
    388 !!gm useless as overlap areas are not used in ptr_vjk 
    389             CALL lbc_lnk( vs, 'V', -1. )   ;   CALL lbc_lnk( vt, 'V', -1. ) 
    390 !!gm 
    391             !                                ! heat & salt advective transports (approximation) 
    392             htr(:,1) = SUM( ptr_vjk( vt(:,:,:) ) , 2 ) * rc_pwatt   ! SUM over jk + conversion 
    393             str(:,1) = SUM( ptr_vjk( vs(:,:,:) ) , 2 ) * rc_ggram 
    394             DO jn = 2, nptr  
    395                htr(:,jn) = SUM( ptr_vjk( vt(:,:,:), btmsk(:,:,jn)*btm30(:,:) ) , 2 ) * rc_pwatt   ! mask Southern Ocean 
    396                str(:,jn) = SUM( ptr_vjk( vs(:,:,:), btmsk(:,:,jn)*btm30(:,:) ) , 2 ) * rc_ggram   ! mask Southern Ocean 
    397             END DO 
    398  
    399             IF( ln_ptrcomp ) THEN            ! overturning transport 
    400                htr_ove(:) = SUM( v_msf(:,:,1) * tn_jk(:,:,1), 2 ) * rc_pwatt   ! SUM over jk + conversion 
    401                str_ove(:) = SUM( v_msf(:,:,1) * sn_jk(:,:,1), 2 ) * rc_ggram 
    402             END IF 
    403             !                                ! Advective and diffusive transport 
    404             htr_adv(:) = htr_adv(:) * rc_pwatt        ! these are computed in tra_adv... and tra_ldf... routines  
    405             htr_ldf(:) = htr_ldf(:) * rc_pwatt        ! here just the conversion in PW and Gg 
    406             str_adv(:) = str_adv(:) * rc_ggram 
    407             str_ldf(:) = str_ldf(:) * rc_ggram 
    408  
    409 #if defined key_diaeiv 
    410             DO jn = 1, nptr                  ! Bolus component 
    411                htr_eiv(:,jn) = SUM( v_msf_eiv(:,:,jn) * tn_jk(:,:,jn), 2 ) * rc_pwatt   ! SUM over jk 
    412                str_eiv(:,jn) = SUM( v_msf_eiv(:,:,jn) * sn_jk(:,:,jn), 2 ) * rc_ggram   ! SUM over jk 
    413             END DO 
    414 #endif 
    415             !                                ! "Meridional" Stream-Function 
    416             DO jn = 1, nptr 
    417                DO jk = 2, jpk  
    418                   v_msf    (:,jk,jn) = v_msf    (:,jk-1,jn) + v_msf    (:,jk,jn)       ! Eulerian j-Stream-Function 
    419 #if defined key_diaeiv 
    420                   v_msf_eiv(:,jk,jn) = v_msf_eiv(:,jk-1,jn) + v_msf_eiv(:,jk,jn)       ! Bolus    j-Stream-Function 
    421  
    422 #endif 
    423                END DO 
    424             END DO 
    425             v_msf    (:,:,:) = v_msf    (:,:,:) * rc_sv       ! converte in Sverdrups 
    426 #if defined key_diaeiv 
    427             v_msf_eiv(:,:,:) = v_msf_eiv(:,:,:) * rc_sv 
    428 #endif 
    429          ENDIF 
    430          ! 
    431          CALL dia_ptr_wri( kt )                        ! outputs 
    432          ! 
    433       ENDIF 
    434       ! 
    435 #if defined key_mpp_mpi 
    436       IF( kt == nitend .AND. l_znl_root )   CALL histclo( numptr )      ! Close the file 
    437 #else 
    438       IF( kt == nitend )                    CALL histclo( numptr )      ! Close the file 
    439 #endif 
    440       ! 
    441       IF( nn_timing == 1 )   CALL timing_stop('dia_ptr') 
    442       ! 
    443    END SUBROUTINE dia_ptr 
    444  
    445  
    446    SUBROUTINE dia_ptr_init 
    447       !!---------------------------------------------------------------------- 
    448       !!                  ***  ROUTINE dia_ptr_init  *** 
    449       !!                    
    450       !! ** Purpose :   Initialization, namelist read 
    451       !!---------------------------------------------------------------------- 
    452       INTEGER ::   jn           ! dummy loop indices  
    453       INTEGER ::   inum, ierr   ! local integers 
    454       INTEGER ::   ios          ! Local integer output status for namelist read 
    455 #if defined key_mpp_mpi 
    456       INTEGER, DIMENSION(1) :: iglo, iloc, iabsf, iabsl, ihals, ihale, idid 
    457 #endif 
    458       !! 
    459       NAMELIST/namptr/ ln_diaptr, ln_diaznl, ln_subbas, ln_ptrcomp, nn_fptr, nn_fwri 
    460       !!---------------------------------------------------------------------- 
    461  
    462       REWIND( numnam_ref )              ! Namelist namptr in reference namelist : Poleward transport 
    463       READ  ( numnam_ref, namptr, IOSTAT = ios, ERR = 901) 
    464 901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namptr in reference namelist', lwp ) 
    465  
    466       REWIND( numnam_cfg )              ! Namelist namptr in configuration namelist : Poleward transport 
    467       READ  ( numnam_cfg, namptr, IOSTAT = ios, ERR = 902 ) 
    468 902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namptr in configuration namelist', lwp ) 
    469       IF(lwm) WRITE ( numond, namptr ) 
    470  
    471       IF(lwp) THEN                     ! Control print 
    472          WRITE(numout,*) 
    473          WRITE(numout,*) 'dia_ptr_init : poleward transport and msf initialization' 
    474          WRITE(numout,*) '~~~~~~~~~~~~' 
    475          WRITE(numout,*) '   Namelist namptr : set ptr parameters' 
    476          WRITE(numout,*) '      Poleward heat & salt transport (T) or not (F)      ln_diaptr  = ', ln_diaptr 
    477          WRITE(numout,*) '      Overturning heat & salt transport                  ln_ptrcomp = ', ln_ptrcomp 
    478          WRITE(numout,*) '      T & S zonal mean and meridional stream function    ln_diaznl  = ', ln_diaznl  
    479          WRITE(numout,*) '      Global (F) or glo/Atl/Pac/Ind/Indo-Pac basins      ln_subbas  = ', ln_subbas 
    480          WRITE(numout,*) '      Frequency of computation                           nn_fptr    = ', nn_fptr 
    481          WRITE(numout,*) '      Frequency of outputs                               nn_fwri    = ', nn_fwri 
    482       ENDIF 
    483        
    484       IF( ln_diaptr) THEN   
    485       
    486          IF( nn_timing == 1 )   CALL timing_start('dia_ptr_init') 
    487        
    488          IF( ln_subbas ) THEN   ;   nptr = 5       ! Global, Atlantic, Pacific, Indian, Indo-Pacific 
    489          ELSE                   ;   nptr = 1       ! Global only 
    490          ENDIF 
    491  
    492          !                                      ! allocate dia_ptr arrays 
    493          IF( dia_ptr_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'dia_ptr_init : unable to allocate arrays' ) 
    494  
    495          rc_pwatt = rc_pwatt * rau0 * rcp          ! conversion from K.s-1 to PetaWatt 
    496  
    497          IF( lk_mpp )   CALL mpp_ini_znl( numout )     ! Define MPI communicator for zonal sum 
    498  
    499          IF( ln_subbas ) THEN                ! load sub-basin mask 
    500             CALL iom_open( 'subbasins', inum ) 
    501             CALL iom_get( inum, jpdom_data, 'atlmsk', btmsk(:,:,2) )   ! Atlantic basin 
    502             CALL iom_get( inum, jpdom_data, 'pacmsk', btmsk(:,:,3) )   ! Pacific  basin 
    503             CALL iom_get( inum, jpdom_data, 'indmsk', btmsk(:,:,4) )   ! Indian   basin 
    504             CALL iom_close( inum ) 
    505             btmsk(:,:,5) = MAX ( btmsk(:,:,3), btmsk(:,:,4) )          ! Indo-Pacific basin 
    506             WHERE( gphit(:,:) < -30._wp)   ;   btm30(:,:) = 0._wp      ! mask out Southern Ocean 
    507             ELSE WHERE                     ;   btm30(:,:) = tmask(:,:,1) 
    508             END WHERE 
    509          ENDIF 
    510          btmsk(:,:,1) = tmask_i(:,:)                                   ! global ocean 
    511        
    512          DO jn = 1, nptr 
    513             btmsk(:,:,jn) = btmsk(:,:,jn) * tmask_i(:,:)               ! interior domain only 
    514          END DO 
    515        
    516          IF( lk_vvl )   CALL ctl_stop( 'diaptr: error in vvl case as constant i-mean surface is used' ) 
    517  
    518          !                                   ! i-sum of e1v*e3v surface and its inverse 
    519          DO jn = 1, nptr 
    520             sjk(:,:,jn) = ptr_tjk( tmask(:,:,:), btmsk(:,:,jn) ) 
    521             r1_sjk(:,:,jn) = 0._wp 
    522             WHERE( sjk(:,:,jn) /= 0._wp )   r1_sjk(:,:,jn) = 1._wp / sjk(:,:,jn) 
    523          END DO 
    524  
    525       ! Initialise arrays to zero because diatpr is called before they are first calculated 
    526       ! Note that this means diagnostics will not be exactly correct when model run is restarted. 
    527       htr_adv(:) = 0._wp ; str_adv(:) =  0._wp ;  htr_ldf(:) = 0._wp ; str_ldf(:) =  0._wp 
    528  
    529 #if defined key_mpp_mpi  
    530          iglo (1) = jpjglo                   ! MPP case using MPI  ('key_mpp_mpi') 
    531          iloc (1) = nlcj 
    532          iabsf(1) = njmppt(narea) 
    533          iabsl(:) = iabsf(:) + iloc(:) - 1 
    534          ihals(1) = nldj - 1 
    535          ihale(1) = nlcj - nlej 
    536          idid (1) = 2 
    537          CALL flio_dom_set( jpnj, nproc/jpni, idid, iglo, iloc, iabsf, iabsl, ihals, ihale, 'BOX', nidom_ptr ) 
    538 #else 
    539          nidom_ptr = FLIO_DOM_NONE 
    540 #endif 
    541       IF( nn_timing == 1 )   CALL timing_stop('dia_ptr_init') 
    542       ! 
    543       ENDIF  
    544       !  
    545    END SUBROUTINE dia_ptr_init 
    546  
    547  
    548    SUBROUTINE dia_ptr_wri( kt ) 
    549       !!--------------------------------------------------------------------- 
    550       !!                ***  ROUTINE dia_ptr_wri  *** 
    551       !! 
    552       !! ** Purpose :   output of poleward fluxes 
    553       !! 
    554       !! ** Method  :   NetCDF file 
    555       !!---------------------------------------------------------------------- 
    556       !! 
    557       INTEGER, INTENT(in) ::   kt   ! ocean time-step index 
    558       !! 
    559       INTEGER, SAVE ::   nhoridz, ndepidzt, ndepidzw 
    560       INTEGER, SAVE ::   ndim  , ndim_atl     , ndim_pac     , ndim_ind     , ndim_ipc 
    561       INTEGER, SAVE ::           ndim_atl_30  , ndim_pac_30  , ndim_ind_30  , ndim_ipc_30 
    562       INTEGER, SAVE ::   ndim_h, ndim_h_atl_30, ndim_h_pac_30, ndim_h_ind_30, ndim_h_ipc_30 
    563       !! 
    564       CHARACTER (len=40) ::   clhstnam, clop, clop_once, cl_comment   ! temporary names 
    565       INTEGER            ::   iline, it, itmod, ji, jj, jk            ! 
    566 #if defined key_iomput 
    567       INTEGER            ::   inum                                    ! temporary logical unit 
    568 #endif 
    569       REAL(wp)           ::   zsto, zout, zdt, zjulian                ! temporary scalars 
    570       !! 
    571       REAL(wp), POINTER, DIMENSION(:)   ::   zphi, zfoo    ! 1D workspace 
    572       REAL(wp), POINTER, DIMENSION(:,:) ::   z_1           ! 2D workspace 
    573       !!--------------------------------------------------------------------  
    574       ! 
    575       CALL wrk_alloc( jpj      , zphi , zfoo ) 
    576       CALL wrk_alloc( jpj , jpk, z_1 ) 
    577  
    578       ! define time axis 
    579       it    = kt / nn_fptr 
    580       itmod = kt - nit000 + 1 
    581        
    582       ! Initialization 
    583       ! -------------- 
    584       IF( kt == nit000 ) THEN 
    585          niter = ( nit000 - 1 ) / nn_fptr 
    586          zdt = rdt 
    587          IF( nacc == 1 )   zdt = rdtmin 
    588          ! 
    589          IF(lwp) THEN 
    590             WRITE(numout,*) 
    591             WRITE(numout,*) 'dia_ptr_wri : poleward transport and msf writing: initialization , niter = ', niter 
    592             WRITE(numout,*) '~~~~~~~~~~~~' 
    593          ENDIF 
    594  
    595          ! Reference latitude (used in plots) 
    596          ! ------------------ 
    597          !                                           ! ======================= 
    598          IF( cp_cfg == "orca" ) THEN                 !   ORCA configurations 
    599             !                                        ! ======================= 
    600             IF( jp_cfg == 05  )   iline = 192   ! i-line that passes near the North Pole 
    601             IF( jp_cfg == 025 )   iline = 384   ! i-line that passes near the North Pole 
    602             IF( jp_cfg == 1   )   iline =  96   ! i-line that passes near the North Pole 
    603             IF( jp_cfg == 2   )   iline =  48   ! i-line that passes near the North Pole 
    604             IF( jp_cfg == 4   )   iline =  24   ! i-line that passes near the North Pole 
    605             zphi(1:jpj) = 0._wp 
    606             DO ji = mi0(iline), mi1(iline)  
    607                zphi(1:jpj) = gphiv(ji,:)         ! if iline is in the local domain 
    608                ! Correct highest latitude for some configurations - will work if domain is parallelized in J ? 
    609                IF( jp_cfg == 05 ) THEN 
    610                   DO jj = mj0(jpjdta), mj1(jpjdta)  
    611                      zphi( jj ) = zphi(mj0(jpjdta-1)) + ( zphi(mj0(jpjdta-1))-zphi(mj0(jpjdta-2)) ) * 0.5_wp 
    612                      zphi( jj ) = MIN( zphi(jj), 90._wp ) 
    613                   END DO 
    614                END IF 
    615                IF( jp_cfg == 1 .OR. jp_cfg == 2 .OR. jp_cfg == 4 ) THEN 
    616                   DO jj = mj0(jpjdta-1), mj1(jpjdta-1)  
    617                      zphi( jj ) = 88.5_wp 
    618                   END DO 
    619                   DO jj = mj0(jpjdta  ), mj1(jpjdta  )  
    620                      zphi( jj ) = 89.5_wp 
    621                   END DO 
    622                END IF 
    623             END DO 
    624             ! provide the correct zphi to all local domains 
    625 #if defined key_mpp_mpi 
    626             CALL mpp_sum( zphi, jpj, ncomm_znl )         
    627 #endif 
    628             !                                        ! ======================= 
    629          ELSE                                        !   OTHER configurations  
    630             !                                        ! ======================= 
    631             zphi(1:jpj) = gphiv(1,:)             ! assume lat/lon coordinate, select the first i-line 
    632             ! 
    633          ENDIF 
    634          ! 
    635          ! Work only on westmost processor (will not work if mppini2 is used) 
    636 #if defined key_mpp_mpi 
    637          IF( l_znl_root ) THEN  
    638 #endif 
    639             ! 
    640             ! OPEN netcdf file  
    641             ! ---------------- 
    642             ! Define frequency of output and means 
    643             zsto = nn_fptr * zdt 
    644             IF( ln_mskland )   THEN    ! put 1.e+20 on land (very expensive!!) 
    645                clop      = "ave(only(x))" 
    646                clop_once = "once(only(x))" 
    647             ELSE                       ! no use of the mask value (require less cpu time) 
    648                clop      = "ave(x)"        
    649                clop_once = "once" 
    650             ENDIF 
    651  
    652             zout = nn_fwri * zdt 
    653             zfoo(1:jpj) = 0._wp 
    654  
    655             CALL ymds2ju( nyear, nmonth, nday, rdt, zjulian )  ! Compute julian date from starting date of the run 
    656             zjulian = zjulian - adatrj                         ! set calendar origin to the beginning of the experiment 
    657  
    658 #if defined key_iomput 
    659             ! Requested by IPSL people, use by their postpro... 
    660             IF(lwp) THEN 
    661                CALL dia_nam( clhstnam, nn_fwri,' ' ) 
    662                CALL ctl_opn( inum, 'date.file', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp, narea ) 
    663                WRITE(inum,*) clhstnam 
    664                CLOSE(inum) 
    665             ENDIF 
    666 #endif 
    667  
    668             CALL dia_nam( clhstnam, nn_fwri, 'diaptr' ) 
    669             IF(lwp)WRITE( numout,*)" Name of diaptr NETCDF file : ", clhstnam 
    670  
    671             ! Horizontal grid : zphi() 
    672             CALL histbeg(clhstnam, 1, zfoo, jpj, zphi,   & 
    673                1, 1, 1, jpj, niter, zjulian, zdt*nn_fptr, nhoridz, numptr, domain_id=nidom_ptr) 
    674             ! Vertical grids : gdept_1d, gdepw_1d 
    675             CALL histvert( numptr, "deptht", "Vertical T levels",   & 
    676                &                   "m", jpk, gdept_1d, ndepidzt, "down" ) 
    677             CALL histvert( numptr, "depthw", "Vertical W levels",   & 
    678                &                   "m", jpk, gdepw_1d, ndepidzw, "down" ) 
    679             ! 
    680             CALL wheneq ( jpj*jpk, MIN(sjk(:,:,1), 1._wp), 1, 1., ndex  , ndim  )      ! Lat-Depth 
    681             CALL wheneq ( jpj    , MIN(sjk(:,1,1), 1._wp), 1, 1., ndex_h, ndim_h )     ! Lat 
    682  
    683             IF( ln_subbas ) THEN 
    684                z_1(:,1) = 1._wp 
    685                WHERE ( gphit(jpi/2,:) < -30._wp )   z_1(:,1) = 0._wp 
    686                DO jk = 2, jpk 
    687                   z_1(:,jk) = z_1(:,1) 
    688                END DO 
    689                !                       ! Atlantic (jn=2) 
    690                CALL wheneq ( jpj*jpk, MIN(sjk(:,:,2)         , 1._wp), 1, 1., ndex_atl     , ndim_atl      ) ! Lat-Depth 
    691                CALL wheneq ( jpj*jpk, MIN(sjk(:,:,2)*z_1(:,:), 1._wp), 1, 1., ndex_atl_30  , ndim_atl_30   ) ! Lat-Depth 
    692                CALL wheneq ( jpj    , MIN(sjk(:,1,2)*z_1(:,1), 1._wp), 1, 1., ndex_h_atl_30, ndim_h_atl_30 ) ! Lat 
    693                !                       ! Pacific (jn=3) 
    694                CALL wheneq ( jpj*jpk, MIN(sjk(:,:,3)         , 1._wp), 1, 1., ndex_pac     , ndim_pac      ) ! Lat-Depth 
    695                CALL wheneq ( jpj*jpk, MIN(sjk(:,:,3)*z_1(:,:), 1._wp), 1, 1., ndex_pac_30  , ndim_pac_30   ) ! Lat-Depth 
    696                CALL wheneq ( jpj    , MIN(sjk(:,1,3)*z_1(:,1), 1._wp), 1, 1., ndex_h_pac_30, ndim_h_pac_30 ) ! Lat 
    697                !                       ! Indian (jn=4) 
    698                CALL wheneq ( jpj*jpk, MIN(sjk(:,:,4)         , 1._wp), 1, 1., ndex_ind     , ndim_ind      ) ! Lat-Depth 
    699                CALL wheneq ( jpj*jpk, MIN(sjk(:,:,4)*z_1(:,:), 1._wp), 1, 1., ndex_ind_30  , ndim_ind_30   ) ! Lat-Depth 
    700                CALL wheneq ( jpj    , MIN(sjk(:,1,4)*z_1(:,1), 1._wp), 1, 1., ndex_h_ind_30, ndim_h_ind_30 ) ! Lat 
    701                !                       ! Indo-Pacific (jn=5) 
    702                CALL wheneq ( jpj*jpk, MIN(sjk(:,:,5)         , 1._wp), 1, 1., ndex_ipc     , ndim_ipc      ) ! Lat-Depth 
    703                CALL wheneq ( jpj*jpk, MIN(sjk(:,:,5)*z_1(:,:), 1._wp), 1, 1., ndex_ipc_30  , ndim_ipc_30   ) ! Lat-Depth 
    704                CALL wheneq ( jpj    , MIN(sjk(:,1,5)*z_1(:,1), 1._wp), 1, 1., ndex_h_ipc_30, ndim_h_ipc_30 ) ! Lat 
    705             ENDIF 
    706             !  
    707 #if defined key_diaeiv 
    708             cl_comment = ' (Bolus part included)' 
    709 #else 
    710             cl_comment = '                      ' 
    711 #endif 
    712             IF( ln_diaznl ) THEN             !  Zonal mean T and S 
    713                CALL histdef( numptr, "zotemglo", "Zonal Mean Temperature","C" ,   & 
    714                   1, jpj, nhoridz, jpk, 1, jpk, ndepidzt, 32, clop, zsto, zout ) 
    715                CALL histdef( numptr, "zosalglo", "Zonal Mean Salinity","PSU"  ,   & 
    716                   1, jpj, nhoridz, jpk, 1, jpk, ndepidzt, 32, clop, zsto, zout ) 
    717  
    718                CALL histdef( numptr, "zosrfglo", "Zonal Mean Surface","m^2"   ,   & 
    719                   1, jpj, nhoridz, jpk, 1, jpk, ndepidzt, 32, clop_once, zsto, zout ) 
    720                ! 
    721                IF (ln_subbas) THEN  
    722                   CALL histdef( numptr, "zotematl", "Zonal Mean Temperature: Atlantic","C" ,   & 
    723                      1, jpj, nhoridz, jpk, 1, jpk, ndepidzt, 32, clop, zsto, zout ) 
    724                   CALL histdef( numptr, "zosalatl", "Zonal Mean Salinity: Atlantic","PSU"  ,   & 
    725                      1, jpj, nhoridz, jpk, 1, jpk, ndepidzt, 32, clop, zsto, zout ) 
    726                   CALL histdef( numptr, "zosrfatl", "Zonal Mean Surface: Atlantic","m^2"   ,   & 
    727                      1, jpj, nhoridz, jpk, 1, jpk, ndepidzt, 32, clop_once, zsto, zout ) 
    728  
    729                   CALL histdef( numptr, "zotempac", "Zonal Mean Temperature: Pacific","C"  ,   & 
    730                      1, jpj, nhoridz, jpk, 1, jpk, ndepidzt, 32, clop, zsto, zout ) 
    731                   CALL histdef( numptr, "zosalpac", "Zonal Mean Salinity: Pacific","PSU"   ,   & 
    732                      1, jpj, nhoridz, jpk, 1, jpk, ndepidzt, 32, clop, zsto, zout ) 
    733                   CALL histdef( numptr, "zosrfpac", "Zonal Mean Surface: Pacific","m^2"    ,   & 
    734                      1, jpj, nhoridz, jpk, 1, jpk, ndepidzt, 32, clop_once, zsto, zout ) 
    735  
    736                   CALL histdef( numptr, "zotemind", "Zonal Mean Temperature: Indian","C"   ,   & 
    737                      1, jpj, nhoridz, jpk, 1, jpk, ndepidzt, 32, clop, zsto, zout ) 
    738                   CALL histdef( numptr, "zosalind", "Zonal Mean Salinity: Indian","PSU"    ,   & 
    739                      1, jpj, nhoridz, jpk, 1, jpk, ndepidzt, 32, clop, zsto, zout ) 
    740                   CALL histdef( numptr, "zosrfind", "Zonal Mean Surface: Indian","m^2"     ,   & 
    741                      1, jpj, nhoridz, jpk, 1, jpk, ndepidzt, 32, clop_once, zsto, zout ) 
    742  
    743                   CALL histdef( numptr, "zotemipc", "Zonal Mean Temperature: Pacific+Indian","C" ,   & 
    744                      1, jpj, nhoridz, jpk, 1, jpk, ndepidzt, 32, clop, zsto, zout ) 
    745                   CALL histdef( numptr, "zosalipc", "Zonal Mean Salinity: Pacific+Indian","PSU"  ,   & 
    746                      1, jpj, nhoridz, jpk, 1, jpk, ndepidzt, 32, clop, zsto, zout ) 
    747                   CALL histdef( numptr, "zosrfipc", "Zonal Mean Surface: Pacific+Indian","m^2"   ,   & 
    748                      1, jpj, nhoridz, jpk, 1, jpk, ndepidzt, 32, clop_once, zsto, zout ) 
    749                ENDIF 
    750             ENDIF 
    751             ! 
    752             !  Meridional Stream-Function (Eulerian and Bolus) 
    753             CALL histdef( numptr, "zomsfglo", "Meridional Stream-Function: Global"//TRIM(cl_comment),"Sv" ,   & 
    754                1, jpj, nhoridz, jpk, 1, jpk, ndepidzw, 32, clop, zsto, zout ) 
    755             IF( ln_subbas .AND. ln_diaznl ) THEN 
    756                CALL histdef( numptr, "zomsfatl", "Meridional Stream-Function: Atlantic"//TRIM(cl_comment),"Sv" ,   & 
    757                   1, jpj, nhoridz, jpk, 1, jpk, ndepidzw, 32, clop, zsto, zout ) 
    758                CALL histdef( numptr, "zomsfpac", "Meridional Stream-Function: Pacific"//TRIM(cl_comment),"Sv"  ,   & 
    759                   1, jpj, nhoridz, jpk, 1, jpk, ndepidzw, 32, clop, zsto, zout ) 
    760                CALL histdef( numptr, "zomsfind", "Meridional Stream-Function: Indian"//TRIM(cl_comment),"Sv"   ,   & 
    761                   1, jpj, nhoridz, jpk, 1, jpk, ndepidzw, 32, clop, zsto, zout ) 
    762                CALL histdef( numptr, "zomsfipc", "Meridional Stream-Function: Indo-Pacific"//TRIM(cl_comment),"Sv" ,& 
    763                   1, jpj, nhoridz, jpk, 1, jpk, ndepidzw, 32, clop, zsto, zout ) 
    764             ENDIF 
    765             ! 
    766             !  Heat transport  
    767             CALL histdef( numptr, "sophtadv", "Advective Heat Transport"      ,   & 
    768                "PW", 1, jpj, nhoridz, 1, 1, 1, -99, 32, clop, zsto, zout ) 
    769             CALL histdef( numptr, "sophtldf", "Diffusive Heat Transport"      ,   & 
    770                "PW",1, jpj, nhoridz, 1, 1, 1, -99, 32, clop, zsto, zout ) 
    771             IF ( ln_ptrcomp ) THEN  
    772                CALL histdef( numptr, "sophtove", "Overturning Heat Transport"    ,   & 
    773                   "PW",1, jpj, nhoridz, 1, 1, 1, -99, 32, clop, zsto, zout ) 
    774             END IF 
    775             IF( ln_subbas ) THEN 
    776                CALL histdef( numptr, "sohtatl", "Heat Transport Atlantic"//TRIM(cl_comment),  & 
    777                   "PW", 1, jpj, nhoridz, 1, 1, 1, -99, 32, clop, zsto, zout ) 
    778                CALL histdef( numptr, "sohtpac", "Heat Transport Pacific"//TRIM(cl_comment) ,  & 
    779                   "PW", 1, jpj, nhoridz, 1, 1, 1, -99, 32, clop, zsto, zout ) 
    780                CALL histdef( numptr, "sohtind", "Heat Transport Indian"//TRIM(cl_comment)  ,  & 
    781                   "PW", 1, jpj, nhoridz, 1, 1, 1, -99, 32, clop, zsto, zout ) 
    782                CALL histdef( numptr, "sohtipc", "Heat Transport Pacific+Indian"//TRIM(cl_comment), & 
    783                   "PW", 1, jpj, nhoridz, 1, 1, 1, -99, 32, clop, zsto, zout ) 
    784             ENDIF 
    785             ! 
    786             !  Salt transport  
    787             CALL histdef( numptr, "sopstadv", "Advective Salt Transport"      ,   & 
    788                "Giga g/s", 1, jpj, nhoridz, 1, 1, 1, -99, 32, clop, zsto, zout ) 
    789             CALL histdef( numptr, "sopstldf", "Diffusive Salt Transport"      ,   & 
    790                "Giga g/s", 1, jpj, nhoridz, 1, 1, 1, -99, 32, clop, zsto, zout ) 
    791             IF ( ln_ptrcomp ) THEN  
    792                CALL histdef( numptr, "sopstove", "Overturning Salt Transport"    ,   & 
    793                   "Giga g/s", 1, jpj, nhoridz, 1, 1, 1, -99, 32, clop, zsto, zout ) 
    794             END IF 
    795 #if defined key_diaeiv 
    796             ! Eddy induced velocity 
    797             CALL histdef( numptr, "zomsfeiv", "Bolus Meridional Stream-Function: global",   & 
    798                "Sv"      , 1, jpj, nhoridz, jpk, 1, jpk, ndepidzw, 32, clop, zsto, zout ) 
    799             CALL histdef( numptr, "sophteiv", "Bolus Advective Heat Transport",   & 
    800                "PW"      , 1, jpj, nhoridz, 1, 1, 1, -99, 32, clop, zsto, zout ) 
    801             CALL histdef( numptr, "sopsteiv", "Bolus Advective Salt Transport",   & 
    802                "Giga g/s", 1, jpj, nhoridz, 1, 1, 1, -99, 32, clop, zsto, zout ) 
    803 #endif 
    804             IF( ln_subbas ) THEN 
    805                CALL histdef( numptr, "sostatl", "Salt Transport Atlantic"//TRIM(cl_comment)      ,  & 
    806                   "Giga g/s", 1, jpj, nhoridz, 1, 1, 1, -99, 32, clop, zsto, zout ) 
    807                CALL histdef( numptr, "sostpac", "Salt Transport Pacific"//TRIM(cl_comment)      ,   & 
    808                   "Giga g/s", 1, jpj, nhoridz, 1, 1, 1, -99, 32, clop, zsto, zout ) 
    809                CALL histdef( numptr, "sostind", "Salt Transport Indian"//TRIM(cl_comment)      ,    & 
    810                   "Giga g/s", 1, jpj, nhoridz, 1, 1, 1, -99, 32, clop, zsto, zout ) 
    811                CALL histdef( numptr, "sostipc", "Salt Transport Pacific+Indian"//TRIM(cl_comment),  & 
    812                   "Giga g/s", 1, jpj, nhoridz, 1, 1, 1, -99, 32, clop, zsto, zout ) 
    813             ENDIF 
    814             ! 
    815             CALL histend( numptr ) 
    816             ! 
    817          END IF 
    818 #if defined key_mpp_mpi 
    819       END IF 
    820 #endif 
    821  
    822 #if defined key_mpp_mpi 
    823       IF( MOD( itmod, nn_fptr ) == 0 .AND. l_znl_root ) THEN 
    824 #else 
    825       IF( MOD( itmod, nn_fptr ) == 0  ) THEN 
    826 #endif 
    827          niter = niter + 1 
    828  
    829          IF( ln_diaznl ) THEN  
    830             CALL histwrite( numptr, "zosrfglo", niter, sjk  (:,:,1) , ndim, ndex ) 
    831             CALL histwrite( numptr, "zotemglo", niter, tn_jk(:,:,1)  , ndim, ndex ) 
    832             CALL histwrite( numptr, "zosalglo", niter, sn_jk(:,:,1)  , ndim, ndex ) 
    833  
    834             IF (ln_subbas) THEN  
    835                CALL histwrite( numptr, "zosrfatl", niter, sjk(:,:,2), ndim_atl, ndex_atl ) 
    836                CALL histwrite( numptr, "zosrfpac", niter, sjk(:,:,3), ndim_pac, ndex_pac ) 
    837                CALL histwrite( numptr, "zosrfind", niter, sjk(:,:,4), ndim_ind, ndex_ind ) 
    838                CALL histwrite( numptr, "zosrfipc", niter, sjk(:,:,5), ndim_ipc, ndex_ipc ) 
    839  
    840                CALL histwrite( numptr, "zotematl", niter, tn_jk(:,:,2)  , ndim_atl, ndex_atl ) 
    841                CALL histwrite( numptr, "zosalatl", niter, sn_jk(:,:,2)  , ndim_atl, ndex_atl ) 
    842                CALL histwrite( numptr, "zotempac", niter, tn_jk(:,:,3)  , ndim_pac, ndex_pac ) 
    843                CALL histwrite( numptr, "zosalpac", niter, sn_jk(:,:,3)  , ndim_pac, ndex_pac ) 
    844                CALL histwrite( numptr, "zotemind", niter, tn_jk(:,:,4)  , ndim_ind, ndex_ind ) 
    845                CALL histwrite( numptr, "zosalind", niter, sn_jk(:,:,4)  , ndim_ind, ndex_ind ) 
    846                CALL histwrite( numptr, "zotemipc", niter, tn_jk(:,:,5)  , ndim_ipc, ndex_ipc ) 
    847                CALL histwrite( numptr, "zosalipc", niter, sn_jk(:,:,5)  , ndim_ipc, ndex_ipc ) 
    848             END IF 
    849          ENDIF 
    850  
    851          ! overturning outputs: 
    852          CALL histwrite( numptr, "zomsfglo", niter, v_msf(:,:,1), ndim, ndex ) 
    853          IF( ln_subbas .AND. ln_diaznl ) THEN 
    854             CALL histwrite( numptr, "zomsfatl", niter, v_msf(:,:,2) , ndim_atl_30, ndex_atl_30 ) 
    855             CALL histwrite( numptr, "zomsfpac", niter, v_msf(:,:,3) , ndim_pac_30, ndex_pac_30 ) 
    856             CALL histwrite( numptr, "zomsfind", niter, v_msf(:,:,4) , ndim_ind_30, ndex_ind_30 ) 
    857             CALL histwrite( numptr, "zomsfipc", niter, v_msf(:,:,5) , ndim_ipc_30, ndex_ipc_30 ) 
    858          ENDIF 
    859 #if defined key_diaeiv 
    860          CALL histwrite( numptr, "zomsfeiv", niter, v_msf_eiv(:,:,1), ndim  , ndex   ) 
    861 #endif 
    862  
    863          ! heat transport outputs: 
    864          IF( ln_subbas ) THEN 
    865             CALL histwrite( numptr, "sohtatl", niter, htr(:,2)  , ndim_h_atl_30, ndex_h_atl_30 ) 
    866             CALL histwrite( numptr, "sohtpac", niter, htr(:,3)  , ndim_h_pac_30, ndex_h_pac_30 ) 
    867             CALL histwrite( numptr, "sohtind", niter, htr(:,4)  , ndim_h_ind_30, ndex_h_ind_30 ) 
    868             CALL histwrite( numptr, "sohtipc", niter, htr(:,5)  , ndim_h_ipc_30, ndex_h_ipc_30 ) 
    869             CALL histwrite( numptr, "sostatl", niter, str(:,2)  , ndim_h_atl_30, ndex_h_atl_30 ) 
    870             CALL histwrite( numptr, "sostpac", niter, str(:,3)  , ndim_h_pac_30, ndex_h_pac_30 ) 
    871             CALL histwrite( numptr, "sostind", niter, str(:,4)  , ndim_h_ind_30, ndex_h_ind_30 ) 
    872             CALL histwrite( numptr, "sostipc", niter, str(:,5)  , ndim_h_ipc_30, ndex_h_ipc_30 ) 
    873          ENDIF 
    874  
    875          CALL histwrite( numptr, "sophtadv", niter, htr_adv     , ndim_h, ndex_h ) 
    876          CALL histwrite( numptr, "sophtldf", niter, htr_ldf     , ndim_h, ndex_h ) 
    877          CALL histwrite( numptr, "sopstadv", niter, str_adv     , ndim_h, ndex_h ) 
    878          CALL histwrite( numptr, "sopstldf", niter, str_ldf     , ndim_h, ndex_h ) 
    879          IF( ln_ptrcomp ) THEN  
    880             CALL histwrite( numptr, "sopstove", niter, str_ove(:) , ndim_h, ndex_h ) 
    881             CALL histwrite( numptr, "sophtove", niter, htr_ove(:) , ndim_h, ndex_h ) 
    882          ENDIF 
    883 #if defined key_diaeiv 
    884          CALL histwrite( numptr, "sophteiv", niter, htr_eiv(:,1)  , ndim_h, ndex_h ) 
    885          CALL histwrite( numptr, "sopsteiv", niter, str_eiv(:,1)  , ndim_h, ndex_h ) 
    886 #endif 
    887          ! 
    888       ENDIF 
    889       ! 
    890       CALL wrk_dealloc( jpj      , zphi , zfoo ) 
    891       CALL wrk_dealloc( jpj , jpk, z_1 ) 
    892       ! 
    893   END SUBROUTINE dia_ptr_wri 
     436   END FUNCTION ptr_sjk 
     437 
    894438 
    895439   !!====================================================================== 

  • branches/2014/dev_r4650_UKMO10_Tidally_Meaned_Diagnostics/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90

    r4756 r5260  
    4444   USE in_out_manager  ! I/O manager 
    4545   USE diadimg         ! dimg direct access file format output 
    46    USE diaar5, ONLY :   lk_diaar5 
    47    USE dynadv, ONLY :   ln_dynadv_vec 
    4846   USE diatmb          ! Top,middle,bottom output 
    4947   USE dia25h          ! 25h Mean output 
    … …  
    8280   !!---------------------------------------------------------------------- 
    8381   !! NEMO/OPA 3.3 , NEMO Consortium (2010) 
    84    !! $Id $ 
     82   !! $Id$ 
    8583   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
    8684   !!---------------------------------------------------------------------- 
    … …  
    9189      INTEGER, DIMENSION(2) :: ierr 
    9290      !!---------------------------------------------------------------------- 
    93       ! 
    9491      ierr = 0 
    95       ! 
    9692      ALLOCATE( ndex_hT(jpi*jpj) , ndex_T(jpi*jpj*jpk) ,     & 
    9793         &      ndex_hU(jpi*jpj) , ndex_U(jpi*jpj*jpk) ,     & 
    … …  
    133129      REAL(wp)                     ::   zztmp, zztmpx, zztmpy   !  
    134130      !! 
    135       REAL(wp), POINTER, DIMENSION(:,:)   :: z2d       ! 2D workspace 
     131      REAL(wp), POINTER, DIMENSION(:,:)   :: z2d      ! 2D workspace 
    136132      REAL(wp), POINTER, DIMENSION(:,:,:) :: z3d      ! 3D workspace 
    137133      !!---------------------------------------------------------------------- 
    … …  
    148144      ENDIF 
    149145 
    150       IF( lk_vvl ) THEN 
    151          z3d(:,:,:) = tsn(:,:,:,jp_tem) * fse3t_n(:,:,:) 
    152          CALL iom_put( "toce" , z3d                        )   ! heat content 
    153          CALL iom_put( "sst"  , z3d(:,:,1)                 )   ! sea surface heat content 
    154          z3d(:,:,1) = tsn(:,:,1,jp_tem) * z3d(:,:,1) 
    155          CALL iom_put( "sst2" , z3d(:,:,1)                 )   ! sea surface content of squared temperature 
    156          z3d(:,:,:) = tsn(:,:,:,jp_sal) * fse3t_n(:,:,:)             
    157          CALL iom_put( "soce" , z3d                        )   ! salinity content 
    158          CALL iom_put( "sss"  , z3d(:,:,1)                 )   ! sea surface salinity content 
    159          z3d(:,:,1) = tsn(:,:,1,jp_sal) * z3d(:,:,1) 
    160          CALL iom_put( "sss2" , z3d(:,:,1)                 )   ! sea surface content of squared salinity 
    161       ELSE 
    162          CALL iom_put( "toce" , tsn(:,:,:,jp_tem)          )   ! temperature 
    163          CALL iom_put( "sst"  , tsn(:,:,1,jp_tem)          )   ! sea surface temperature 
    164          CALL iom_put( "sst2" , tsn(:,:,1,jp_tem) * tsn(:,:,1,jp_tem) ) ! square of sea surface temperature 
    165          CALL iom_put( "soce" , tsn(:,:,:,jp_sal)          )   ! salinity 
    166          CALL iom_put( "sss"  , tsn(:,:,1,jp_sal)          )   ! sea surface salinity 
    167          CALL iom_put( "sss2" , tsn(:,:,1,jp_sal) * tsn(:,:,1,jp_sal) ) ! square of sea surface salinity 
    168       END IF 
    169       IF( lk_vvl .AND. (.NOT. ln_dynadv_vec) ) THEN 
    170          CALL iom_put( "uoce" , un(:,:,:) * fse3u_n(:,:,:) )    ! i-transport 
    171          CALL iom_put( "voce" , vn(:,:,:) * fse3v_n(:,:,:) )    ! j-transport 
    172       ELSE 
    173          CALL iom_put( "uoce" , un                         )    ! i-current 
    174          CALL iom_put( "voce" , vn                         )    ! j-current 
    175       END IF 
    176       CALL iom_put(    "avt"  , avt                        )    ! T vert. eddy diff. coef. 
    177       CALL iom_put(    "avm"  , avmu                       )    ! T vert. eddy visc. coef. 
    178       IF( lk_zdfddm ) THEN 
    179          CALL iom_put( "avs" , fsavs(:,:,:)                          )    ! S vert. eddy diff. coef. 
    180       ENDIF 
    181  
    182       DO jj = 2, jpjm1                                    ! sst gradient 
    183          DO ji = fs_2, fs_jpim1   ! vector opt. 
    184             zztmp      = tsn(ji,jj,1,jp_tem) 
    185             zztmpx     = ( tsn(ji+1,jj  ,1,jp_tem) - zztmp ) / e1u(ji,jj) + ( zztmp - tsn(ji-1,jj  ,1,jp_tem) ) / e1u(ji-1,jj  ) 
    186             zztmpy     = ( tsn(ji  ,jj+1,1,jp_tem) - zztmp ) / e2v(ji,jj) + ( zztmp - tsn(ji  ,jj-1,1,jp_tem) ) / e2v(ji  ,jj-1) 
    187             z2d(ji,jj) = 0.25 * ( zztmpx * zztmpx + zztmpy * zztmpy )   & 
    188                &              * umask(ji,jj,1) * umask(ji-1,jj,1) * vmask(ji,jj,1) * umask(ji,jj-1,1) 
    189          END DO 
    190       END DO 
    191       CALL lbc_lnk( z2d, 'T', 1. ) 
    192       CALL iom_put( "sstgrad2",  z2d               )    ! square of module of sst gradient 
    193 !CDIR NOVERRCHK 
    194       z2d(:,:) = SQRT( z2d(:,:) ) 
    195       CALL iom_put( "sstgrad" ,  z2d               )    ! module of sst gradient 
    196  
    197       IF( lk_diaar5 ) THEN 
     146      IF( .NOT.lk_vvl ) THEN 
     147         CALL iom_put( "e3t" , fse3t_n(:,:,:) ) 
     148         CALL iom_put( "e3u" , fse3u_n(:,:,:) ) 
     149         CALL iom_put( "e3v" , fse3v_n(:,:,:) ) 
     150         CALL iom_put( "e3w" , fse3w_n(:,:,:) ) 
     151      ENDIF 
     152       
     153      CALL iom_put( "toce", tsn(:,:,:,jp_tem) )    ! 3D temperature 
     154      CALL iom_put(  "sst", tsn(:,:,1,jp_tem) )    ! surface temperature 
     155      IF ( iom_use("sbt") ) THEN 
     156         DO jj = 1, jpj 
     157            DO ji = 1, jpi 
     158               z2d(ji,jj) = tsn(ji,jj,MAX(mbathy(ji,jj),1),jp_tem) 
     159            END DO 
     160         END DO 
     161         CALL iom_put( "sbt", z2d )                ! bottom temperature 
     162      ENDIF 
     163       
     164      CALL iom_put( "soce", tsn(:,:,:,jp_sal) )    ! 3D salinity 
     165      CALL iom_put(  "sss", tsn(:,:,1,jp_sal) )    ! surface salinity 
     166      IF ( iom_use("sbs") ) THEN 
     167         DO jj = 1, jpj 
     168            DO ji = 1, jpi 
     169               z2d(ji,jj) = tsn(ji,jj,MAX(mbathy(ji,jj),1),jp_sal) 
     170            END DO 
     171         END DO 
     172         CALL iom_put( "sbs", z2d )                ! bottom salinity 
     173      ENDIF 
     174          
     175      CALL iom_put( "uoce", un(:,:,:)         )    ! 3D i-current 
     176      CALL iom_put(  "ssu", un(:,:,1)         )    ! surface i-current 
     177      IF ( iom_use("sbu") ) THEN 
     178         DO jj = 1, jpj 
     179            DO ji = 1, jpi 
     180               z2d(ji,jj) = un(ji,jj,MAX(mbathy(ji,jj),1)) 
     181            END DO 
     182         END DO 
     183         CALL iom_put( "sbu", z2d )                ! bottom i-current 
     184      ENDIF 
     185       
     186      CALL iom_put( "voce", vn(:,:,:)         )    ! 3D j-current 
     187      CALL iom_put(  "ssv", vn(:,:,1)         )    ! surface j-current 
     188      IF ( iom_use("sbv") ) THEN 
     189         DO jj = 1, jpj 
     190            DO ji = 1, jpi 
     191               z2d(ji,jj) = vn(ji,jj,MAX(mbathy(ji,jj),1)) 
     192            END DO 
     193         END DO 
     194         CALL iom_put( "sbv", z2d )                ! bottom j-current 
     195      ENDIF 
     196 
     197      CALL iom_put( "avt" , avt                        )    ! T vert. eddy diff. coef. 
     198      CALL iom_put( "avm" , avmu                       )    ! T vert. eddy visc. coef. 
     199      CALL iom_put( "avs" , fsavs(:,:,:)               )    ! S vert. eddy diff. coef. (useful only with key_zdfddm) 
     200 
     201      IF ( iom_use("sstgrad") .OR. iom_use("sstgrad2") ) THEN 
     202         DO jj = 2, jpjm1                                    ! sst gradient 
     203            DO ji = fs_2, fs_jpim1   ! vector opt. 
     204               zztmp      = tsn(ji,jj,1,jp_tem) 
     205               zztmpx     = ( tsn(ji+1,jj  ,1,jp_tem) - zztmp ) / e1u(ji,jj) + ( zztmp - tsn(ji-1,jj  ,1,jp_tem) ) / e1u(ji-1,jj  ) 
     206               zztmpy     = ( tsn(ji  ,jj+1,1,jp_tem) - zztmp ) / e2v(ji,jj) + ( zztmp - tsn(ji  ,jj-1,1,jp_tem) ) / e2v(ji  ,jj-1) 
     207               z2d(ji,jj) = 0.25 * ( zztmpx * zztmpx + zztmpy * zztmpy )   & 
     208                  &              * umask(ji,jj,1) * umask(ji-1,jj,1) * vmask(ji,jj,1) * umask(ji,jj-1,1) 
     209            END DO 
     210         END DO 
     211         CALL lbc_lnk( z2d, 'T', 1. ) 
     212         CALL iom_put( "sstgrad2",  z2d               )    ! square of module of sst gradient 
     213         z2d(:,:) = SQRT( z2d(:,:) ) 
     214         CALL iom_put( "sstgrad" ,  z2d               )    ! module of sst gradient 
     215      ENDIF 
     216          
     217      ! clem: heat and salt content 
     218      IF( iom_use("heatc") ) THEN 
     219         z2d(:,:)  = 0._wp  
     220         DO jk = 1, jpkm1 
     221            DO jj = 1, jpj 
     222               DO ji = 1, jpi 
     223                  z2d(ji,jj) = z2d(ji,jj) + fse3t(ji,jj,jk) * tsn(ji,jj,jk,jp_tem) * tmask(ji,jj,jk) 
     224               END DO 
     225            END DO 
     226         END DO 
     227         CALL iom_put( "heatc", (rau0 * rcp) * z2d )    ! vertically integrated heat content (J/m2) 
     228      ENDIF 
     229 
     230      IF( iom_use("saltc") ) THEN 
     231         z2d(:,:)  = 0._wp  
     232         DO jk = 1, jpkm1 
     233            DO jj = 1, jpj 
     234               DO ji = 1, jpi 
     235                  z2d(ji,jj) = z2d(ji,jj) + fse3t(ji,jj,jk) * tsn(ji,jj,jk,jp_sal) * tmask(ji,jj,jk) 
     236               END DO 
     237            END DO 
     238         END DO 
     239         CALL iom_put( "saltc", rau0 * z2d )   ! vertically integrated salt content (PSU*kg/m2) 
     240      ENDIF 
     241      ! 
     242      IF ( iom_use("eken") ) THEN 
     243         rke(:,:,jk) = 0._wp                               !      kinetic energy  
     244         DO jk = 1, jpkm1 
     245            DO jj = 2, jpjm1 
     246               DO ji = fs_2, fs_jpim1   ! vector opt. 
     247                  zztmp   = 1._wp / ( e1e2t(ji,jj) * fse3t(ji,jj,jk) ) 
     248                  zztmpx  = 0.5 * (  un(ji-1,jj,jk) * un(ji-1,jj,jk) * e2u(ji-1,jj) * fse3u(ji-1,jj,jk)    & 
     249                     &             + un(ji  ,jj,jk) * un(ji  ,jj,jk) * e2u(ji  ,jj) * fse3u(ji  ,jj,jk) )  & 
     250                     &          *  zztmp  
     251                  ! 
     252                  zztmpy  = 0.5 * (  vn(ji,jj-1,jk) * vn(ji,jj-1,jk) * e1v(ji,jj-1) * fse3v(ji,jj-1,jk)    & 
     253                     &             + vn(ji,jj  ,jk) * vn(ji,jj  ,jk) * e1v(ji,jj  ) * fse3v(ji,jj  ,jk) )  & 
     254                     &          *  zztmp  
     255                  ! 
     256                  rke(ji,jj,jk) = 0.5_wp * ( zztmpx + zztmpy ) 
     257                  ! 
     258               ENDDO 
     259            ENDDO 
     260         ENDDO 
     261         CALL lbc_lnk( rke, 'T', 1. ) 
     262         CALL iom_put( "eken", rke )            
     263      ENDIF 
     264          
     265      IF( iom_use("u_masstr") .OR. iom_use("u_heattr") .OR. iom_use("u_salttr") ) THEN 
    198266         z3d(:,:,jpk) = 0.e0 
    199267         DO jk = 1, jpkm1 
    200             z3d(:,:,jk) = rau0 * un(:,:,jk) * e2u(:,:) * fse3u(:,:,jk) 
     268            z3d(:,:,jk) = rau0 * un(:,:,jk) * e2u(:,:) * fse3u(:,:,jk) * umask(:,:,jk) 
    201269         END DO 
    202270         CALL iom_put( "u_masstr", z3d )                  ! mass transport in i-direction 
    203          zztmp = 0.5 * rcp 
     271      ENDIF 
     272       
     273      IF( iom_use("u_heattr") ) THEN 
    204274         z2d(:,:) = 0.e0  
    205275         DO jk = 1, jpkm1 
    206276            DO jj = 2, jpjm1 
    207277               DO ji = fs_2, fs_jpim1   ! vector opt. 
    208                   z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * zztmp * ( tsn(ji,jj,jk,jp_tem) + tsn(ji+1,jj,jk,jp_tem) ) 
     278                  z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( tsn(ji,jj,jk,jp_tem) + tsn(ji+1,jj,jk,jp_tem) ) 
    209279               END DO 
    210280            END DO 
    211281         END DO 
    212282         CALL lbc_lnk( z2d, 'U', -1. ) 
    213          CALL iom_put( "u_heattr", z2d )                  ! heat transport in i-direction 
    214          DO jk = 1, jpkm1 
    215             z3d(:,:,jk) = rau0 * vn(:,:,jk) * e1v(:,:) * fse3v(:,:,jk) 
    216          END DO 
    217          CALL iom_put( "v_masstr", z3d )                  ! mass transport in j-direction 
     283         CALL iom_put( "u_heattr", (0.5 * rcp) * z2d )    ! heat transport in i-direction 
     284      ENDIF 
     285 
     286      IF( iom_use("u_salttr") ) THEN 
    218287         z2d(:,:) = 0.e0  
    219288         DO jk = 1, jpkm1 
    220289            DO jj = 2, jpjm1 
    221290               DO ji = fs_2, fs_jpim1   ! vector opt. 
    222                   z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * zztmp * ( tsn(ji,jj,jk,jp_tem) + tsn(ji,jj+1,jk,jp_tem) ) 
     291                  z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( tsn(ji,jj,jk,jp_sal) + tsn(ji+1,jj,jk,jp_sal) ) 
    223292               END DO 
    224293            END DO 
    225294         END DO 
     295         CALL lbc_lnk( z2d, 'U', -1. ) 
     296         CALL iom_put( "u_salttr", 0.5 * z2d )            ! heat transport in i-direction 
     297      ENDIF 
     298 
     299       
     300      IF( iom_use("v_masstr") .OR. iom_use("v_heattr") .OR. iom_use("v_salttr") ) THEN 
     301         z3d(:,:,jpk) = 0.e0 
     302         DO jk = 1, jpkm1 
     303            z3d(:,:,jk) = rau0 * vn(:,:,jk) * e1v(:,:) * fse3v(:,:,jk) * vmask(:,:,jk) 
     304         END DO 
     305         CALL iom_put( "v_masstr", z3d )                  ! mass transport in j-direction 
     306      ENDIF 
     307       
     308      IF( iom_use("v_heattr") ) THEN 
     309         z2d(:,:) = 0.e0  
     310         DO jk = 1, jpkm1 
     311            DO jj = 2, jpjm1 
     312               DO ji = fs_2, fs_jpim1   ! vector opt. 
     313                  z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( tsn(ji,jj,jk,jp_tem) + tsn(ji,jj+1,jk,jp_tem) ) 
     314               END DO 
     315            END DO 
     316         END DO 
    226317         CALL lbc_lnk( z2d, 'V', -1. ) 
    227          CALL iom_put( "v_heattr", z2d )                  !  heat transport in i-direction 
     318         CALL iom_put( "v_heattr", (0.5 * rcp) * z2d )    !  heat transport in j-direction 
     319      ENDIF 
     320 
     321      IF( iom_use("v_salttr") ) THEN 
     322         z2d(:,:) = 0.e0  
     323         DO jk = 1, jpkm1 
     324            DO jj = 2, jpjm1 
     325               DO ji = fs_2, fs_jpim1   ! vector opt. 
     326                  z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( tsn(ji,jj,jk,jp_sal) + tsn(ji,jj+1,jk,jp_sal) ) 
     327               END DO 
     328            END DO 
     329         END DO 
     330         CALL lbc_lnk( z2d, 'V', -1. ) 
     331         CALL iom_put( "v_salttr", 0.5 * z2d )            !  heat transport in j-direction 
    228332      ENDIF 
    229333      ! 
    … …  
    500604         ENDIF 
    501605 
    502 #if ! defined key_coupled  
    503          CALL histdef( nid_T, "sohefldp", "Surface Heat Flux: Damping"         , "W/m2"   ,   &  ! qrp 
    504             &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout ) 
    505          CALL histdef( nid_T, "sowafldp", "Surface Water Flux: Damping"        , "Kg/m2/s",   &  ! erp 
    506             &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout ) 
    507          CALL histdef( nid_T, "sosafldp", "Surface salt flux: damping"         , "Kg/m2/s",   &  ! erp * sn 
    508             &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout ) 
    509 #endif 
    510  
    511  
    512  
    513 #if ( defined key_coupled && ! defined key_lim3 && ! defined key_lim2 )  
    514          CALL histdef( nid_T, "sohefldp", "Surface Heat Flux: Damping"         , "W/m2"   ,   &  ! qrp 
    515             &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout ) 
    516          CALL histdef( nid_T, "sowafldp", "Surface Water Flux: Damping"        , "Kg/m2/s",   &  ! erp 
    517             &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout ) 
    518          CALL histdef( nid_T, "sosafldp", "Surface salt flux: Damping"         , "Kg/m2/s",   &  ! erp * sn 
    519             &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout ) 
    520 #endif 
     606         IF( .NOT. lk_cpl ) THEN 
     607            CALL histdef( nid_T, "sohefldp", "Surface Heat Flux: Damping"         , "W/m2"   ,   &  ! qrp 
     608               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout ) 
     609            CALL histdef( nid_T, "sowafldp", "Surface Water Flux: Damping"        , "Kg/m2/s",   &  ! erp 
     610               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout ) 
     611            CALL histdef( nid_T, "sosafldp", "Surface salt flux: damping"         , "Kg/m2/s",   &  ! erp * sn 
     612               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout ) 
     613         ENDIF 
     614 
     615         IF( lk_cpl .AND. nn_ice <= 1 ) THEN 
     616            CALL histdef( nid_T, "sohefldp", "Surface Heat Flux: Damping"         , "W/m2"   ,   &  ! qrp 
     617               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout ) 
     618            CALL histdef( nid_T, "sowafldp", "Surface Water Flux: Damping"        , "Kg/m2/s",   &  ! erp 
     619               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout ) 
     620            CALL histdef( nid_T, "sosafldp", "Surface salt flux: Damping"         , "Kg/m2/s",   &  ! erp * sn 
     621               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout ) 
     622         ENDIF 
     623          
    521624         clmx ="l_max(only(x))"    ! max index on a period 
    522625         CALL histdef( nid_T, "sobowlin", "Bowl Index"                         , "W-point",   &  ! bowl INDEX  
    … …  
    533636#endif 
    534637 
    535 #if defined key_coupled  
    536 # if defined key_lim3 
    537          Must be adapted to LIM3 
    538 # endif  
    539 # if defined key_lim2 
    540          CALL histdef( nid_T,"soicetem" , "Ice Surface Temperature"            , "K"      ,   &  ! tn_ice 
    541             &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout ) 
    542          CALL histdef( nid_T,"soicealb" , "Ice Albedo"                         , "[0,1]"  ,   &  ! alb_ice 
    543             &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout ) 
    544 # endif  
    545 #endif  
     638         IF( lk_cpl .AND. nn_ice == 2 ) THEN 
     639            CALL histdef( nid_T,"soicetem" , "Ice Surface Temperature"            , "K"      ,   &  ! tn_ice 
     640               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout ) 
     641            CALL histdef( nid_T,"soicealb" , "Ice Albedo"                         , "[0,1]"  ,   &  ! alb_ice 
     642               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout ) 
     643         ENDIF 
    546644 
    547645         CALL histend( nid_T, snc4chunks=snc4set ) 
    … …  
    634732      ENDIF 
    635733 
    636       ! Write fields on T grid 
    637734      IF( lk_vvl ) THEN 
    638735         CALL histwrite( nid_T, "votemper", it, tsn(:,:,:,jp_tem) * fse3t_n(:,:,:) , ndim_T , ndex_T  )   ! heat content 
    … …  
    645742         CALL histwrite( nid_T, "sosstsst", it, tsn(:,:,1,jp_tem) , ndim_hT, ndex_hT )   ! sea surface temperature 
    646743         CALL histwrite( nid_T, "sosaline", it, tsn(:,:,1,jp_sal) , ndim_hT, ndex_hT )   ! sea surface salinity 
    647  
    648744      ENDIF 
    649745      IF( lk_vvl ) THEN 
    … …  
    695791      ENDIF 
    696792 
    697 #if ! defined key_coupled 
    698       CALL histwrite( nid_T, "sohefldp", it, qrp           , ndim_hT, ndex_hT )   ! heat flux damping 
    699       CALL histwrite( nid_T, "sowafldp", it, erp           , ndim_hT, ndex_hT )   ! freshwater flux damping 
    700       IF( ln_ssr ) zw2d(:,:) = erp(:,:) * tsn(:,:,1,jp_sal) * tmask(:,:,1) 
    701       CALL histwrite( nid_T, "sosafldp", it, zw2d          , ndim_hT, ndex_hT )   ! salt flux damping 
    702 #endif 
    703 #if ( defined key_coupled && ! defined key_lim3 && ! defined key_lim2 )  
    704       CALL histwrite( nid_T, "sohefldp", it, qrp           , ndim_hT, ndex_hT )   ! heat flux damping 
    705       CALL histwrite( nid_T, "sowafldp", it, erp           , ndim_hT, ndex_hT )   ! freshwater flux damping 
     793      IF( .NOT. lk_cpl ) THEN 
     794         CALL histwrite( nid_T, "sohefldp", it, qrp           , ndim_hT, ndex_hT )   ! heat flux damping 
     795         CALL histwrite( nid_T, "sowafldp", it, erp           , ndim_hT, ndex_hT )   ! freshwater flux damping 
    706796         IF( ln_ssr ) zw2d(:,:) = erp(:,:) * tsn(:,:,1,jp_sal) * tmask(:,:,1) 
    707       CALL histwrite( nid_T, "sosafldp", it, zw2d          , ndim_hT, ndex_hT )   ! salt flux damping 
    708 #endif 
    709       zw2d(:,:) = FLOAT( nmln(:,:) ) * tmask(:,:,1) 
    710       CALL histwrite( nid_T, "sobowlin", it, zw2d          , ndim_hT, ndex_hT )   ! ??? 
     797         CALL histwrite( nid_T, "sosafldp", it, zw2d          , ndim_hT, ndex_hT )   ! salt flux damping 
     798      ENDIF 
     799      IF( lk_cpl .AND. nn_ice <= 1 ) THEN 
     800         CALL histwrite( nid_T, "sohefldp", it, qrp           , ndim_hT, ndex_hT )   ! heat flux damping 
     801         CALL histwrite( nid_T, "sowafldp", it, erp           , ndim_hT, ndex_hT )   ! freshwater flux damping 
     802         IF( ln_ssr ) zw2d(:,:) = erp(:,:) * tsn(:,:,1,jp_sal) * tmask(:,:,1) 
     803         CALL histwrite( nid_T, "sosafldp", it, zw2d          , ndim_hT, ndex_hT )   ! salt flux damping 
     804      ENDIF 
     805!      zw2d(:,:) = FLOAT( nmln(:,:) ) * tmask(:,:,1) 
     806!      CALL histwrite( nid_T, "sobowlin", it, zw2d          , ndim_hT, ndex_hT )   ! ??? 
    711807 
    712808#if defined key_diahth 
    … …  
    717813#endif 
    718814 
    719 #if defined key_coupled  
    720 # if defined key_lim3 
    721       Must be adapted for LIM3 
    722       CALL histwrite( nid_T, "soicetem", it, tn_ice        , ndim_hT, ndex_hT )   ! surf. ice temperature 
    723       CALL histwrite( nid_T, "soicealb", it, alb_ice       , ndim_hT, ndex_hT )   ! ice albedo 
    724 # endif 
    725 # if defined key_lim2 
    726       CALL histwrite( nid_T, "soicetem", it, tn_ice(:,:,1) , ndim_hT, ndex_hT )   ! surf. ice temperature 
    727       CALL histwrite( nid_T, "soicealb", it, alb_ice(:,:,1), ndim_hT, ndex_hT )   ! ice albedo 
    728 # endif 
    729 #endif 
    730          ! Write fields on U grid 
     815      IF( lk_cpl .AND. nn_ice == 2 ) THEN 
     816         CALL histwrite( nid_T, "soicetem", it, tn_ice(:,:,1) , ndim_hT, ndex_hT )   ! surf. ice temperature 
     817         CALL histwrite( nid_T, "soicealb", it, alb_ice(:,:,1), ndim_hT, ndex_hT )   ! ice albedo 
     818      ENDIF 
     819 
    731820      CALL histwrite( nid_U, "vozocrtx", it, un            , ndim_U , ndex_U )    ! i-current 
    732821      IF( ln_traldf_gdia ) THEN 
    … …  
    750839      CALL histwrite( nid_U, "sozotaux", it, utau          , ndim_hU, ndex_hU )   ! i-wind stress 
    751840 
    752          ! Write fields on V grid 
    753841      CALL histwrite( nid_V, "vomecrty", it, vn            , ndim_V , ndex_V  )   ! j-current 
    754842      IF( ln_traldf_gdia ) THEN 
    … …  
    765853      CALL histwrite( nid_V, "sometauy", it, vtau          , ndim_hV, ndex_hV )   ! j-wind stress 
    766854 
    767          ! Write fields on W grid 
    768855      CALL histwrite( nid_W, "vovecrtz", it, wn             , ndim_T, ndex_T )    ! vert. current 
    769856      IF( ln_traldf_gdia ) THEN 
Note: See TracChangeset for help on using the changeset viewer.