Index: branches/UKMO/AMM15_v3_6_STABLE_package_reanalysis2/NEMOGCM/NEMO/OPA_SRC/DIA/diaar5.F90
===================================================================
--- branches/UKMO/AMM15_v3_6_STABLE_package_reanalysis2/NEMOGCM/NEMO/OPA_SRC/DIA/diaar5.F90	(revision 10941)
+++ branches/UKMO/AMM15_v3_6_STABLE_package_reanalysis2/NEMOGCM/NEMO/OPA_SRC/DIA/diaar5.F90	(revision 10962)
@@ -39,4 +39,9 @@
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:  ) ::   thick0       ! ocean thickness (interior domain)
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   sn0          ! initial salinity
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   tn0          ! initial temperature
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   sshthster_mat         ! ssh_thermosteric height
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   sshhlster_mat         ! ssh_halosteric height
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   sshsteric_mat         ! ssh_steric height
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   zbotpres_mat          ! bottom pressure
       
    !! * Substitutions
@@ -56,5 +61,7 @@
       !!----------------------------------------------------------------------
       !
-      ALLOCATE( area(jpi,jpj), thick0(jpi,jpj) , sn0(jpi,jpj,jpk) , STAT=dia_ar5_alloc )
+      ALLOCATE( area(jpi,jpj), thick0(jpi,jpj) , sn0(jpi,jpj,jpk), tn0(jpi,jpj,jpk) , &
+          & sshthster_mat(jpi,jpj),sshhlster_mat(jpi,jpj),sshsteric_mat(jpi,jpj), &
+          & zbotpres_mat(jpi,jpj),STAT=dia_ar5_alloc )
       !
       IF( lk_mpp             )   CALL mpp_sum ( dia_ar5_alloc )
@@ -85,4 +92,9 @@
       CALL wrk_alloc( jpi , jpj , jpk        , zrhd      , zrhop    )
       CALL wrk_alloc( jpi , jpj , jpk , jpts , ztsn                 )
+      
+      sshthster_mat(:,:) = 0._wp  
+      sshhlster_mat(:,:) = 0._wp  
+      sshsteric_mat(:,:) = 0._wp  
+      zbotpres_mat(:,:)  = 0._wp  
 
       zarea_ssh(:,:) = area(:,:) * sshn(:,:)
@@ -121,4 +133,36 @@
       zssh_steric = - zarho / area_tot
       CALL iom_put( 'sshthster', zssh_steric )
+      sshthster_mat(:,:) =  -zbotpres(:,:)
+      CALL iom_put( 'sshthster_mat', sshthster_mat )
+
+      !                     
+      ztsn(:,:,:,jp_tem) = tn0(:,:,:)                    ! thermohaline ssh
+      ztsn(:,:,:,jp_sal) = tsn(:,:,:,jp_sal) 
+      CALL eos( ztsn, zrhd, fsdept_n(:,:,:) )                       ! now in situ density using initial temperature
+      !
+      zbotpres(:,:) = 0._wp                        ! no atmospheric surface pressure, levitating sea-ice
+      DO jk = 1, jpkm1
+         zbotpres(:,:) = zbotpres(:,:) + fse3t(:,:,jk) * zrhd(:,:,jk)
+      END DO
+      IF( .NOT.lk_vvl ) THEN
+         IF ( ln_isfcav ) THEN
+            DO ji=1,jpi
+               DO jj=1,jpj
+                  zbotpres(ji,jj) = zbotpres(ji,jj) + sshn(ji,jj) * zrhd(ji,jj,mikt(ji,jj)) + riceload(ji,jj)
+               END DO
+            END DO
+         ELSE
+            zbotpres(:,:) = zbotpres(:,:) + sshn(:,:) * zrhd(:,:,1)
+         END IF
+      END IF
+      !                                         
+      zarho = SUM( area(:,:) * zbotpres(:,:) ) 
+      IF( lk_mpp )   CALL mpp_sum( zarho )
+      zssh_steric = - zarho / area_tot
+      CALL iom_put( 'sshhlster', zssh_steric )
+      sshhlster_mat(:,:) = -zbotpres(:,:)
+      CALL iom_put( 'sshhlster_mat', sshhlster_mat )
+      
+
       
       !                                         ! steric sea surface height
@@ -147,10 +191,13 @@
       zssh_steric = - zarho / area_tot
       CALL iom_put( 'sshsteric', zssh_steric )
+      sshsteric_mat(:,:) = -zbotpres(:,:) 
+      CALL iom_put( 'sshsteric_mat', sshsteric_mat )
       
       !                                         ! ocean bottom pressure
       zztmp = rau0 * grav * 1.e-4_wp               ! recover pressure from pressure anomaly and cover to dbar = 1.e4 Pa
       zbotpres(:,:) = zztmp * ( zbotpres(:,:) + sshn(:,:) + thick0(:,:) )
-      CALL iom_put( 'botpres', zbotpres )
-
+      zbotpres_mat(:,:) = zbotpres(:,:)
+      CALL iom_put( 'botpres', zbotpres_mat )
+      
       !                                         ! Mean density anomalie, temperature and salinity
       ztemp = 0._wp
@@ -211,4 +258,5 @@
       REAL(wp) ::   zztmp  
       REAL(wp), POINTER, DIMENSION(:,:,:,:) ::   zsaldta   ! Jan/Dec levitus salinity
+      REAL(wp), POINTER, DIMENSION(:,:,:,:) ::   ztemdta   ! Jan/Dec levitus salinity
       ! reading initial file
       LOGICAL  ::   ln_tsd_init      !: T & S data flag
@@ -234,4 +282,5 @@
       !
       CALL wrk_alloc( jpi , jpj , jpk, jpts, zsaldta )
+      CALL wrk_alloc( jpi , jpj , jpk, jpts, ztemdta )
       !                                      ! allocate dia_ar5 arrays
       IF( dia_ar5_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'dia_ar5_init : unable to allocate arrays' )
@@ -253,6 +302,13 @@
       CALL iom_get  ( inum, jpdom_data, TRIM(sn_sal%clvar), zsaldta(:,:,:,2), 12 )
       CALL iom_close( inum )
+
+      CALL iom_open ( TRIM( cn_dir )//TRIM(sn_tem%clname), inum )
+      CALL iom_get  ( inum, jpdom_data, TRIM(sn_tem%clvar), ztemdta(:,:,:,1), 1  )
+      CALL iom_get  ( inum, jpdom_data, TRIM(sn_tem%clvar), ztemdta(:,:,:,2), 12 )
+      CALL iom_close( inum )
       sn0(:,:,:) = 0.5_wp * ( zsaldta(:,:,:,1) + zsaldta(:,:,:,2) )        
-      sn0(:,:,:) = sn0(:,:,:) * tmask(:,:,:)
+      tn0(:,:,:) = 0.5_wp * ( ztemdta(:,:,:,1) + ztemdta(:,:,:,2) )        
+      sn0(:,:,:) = sn0(:,:,:) * tmask(:,:,:)  
+      tn0(:,:,:) = tn0(:,:,:) * tmask(:,:,:)
       IF( ln_zps ) THEN               ! z-coord. partial steps
          DO jj = 1, jpj               ! interpolation of salinity at the last ocean level (i.e. the partial step)
@@ -262,4 +318,5 @@
                   zztmp = ( gdept_1d(ik) - gdept_0(ji,jj,ik) ) / ( gdept_1d(ik) - gdept_1d(ik-1) )
                   sn0(ji,jj,ik) = ( 1._wp - zztmp ) * sn0(ji,jj,ik) + zztmp * sn0(ji,jj,ik-1)
+                  tn0(ji,jj,ik) = ( 1._wp - zztmp ) * tn0(ji,jj,ik) + zztmp * tn0(ji,jj,ik-1)
                ENDIF
             END DO
@@ -268,4 +325,5 @@
       !
       CALL wrk_dealloc( jpi , jpj , jpk, jpts, zsaldta )
+      CALL wrk_dealloc( jpi , jpj , jpk, jpts, ztemdta )
       !
       IF( nn_timing == 1 )   CALL timing_stop('dia_ar5_init')
