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 6225 for branches/2014/dev_r4704_NOC5_MPP_BDY_UPDATE/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90 – NEMO

Ignore:
Timestamp:
2016-01-08T10:35:19+01:00 (8 years ago)
Author:
jamesharle
Message:

Update MPP_BDY_UPDATE branch to be consistent with head of trunk

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/2014/dev_r4704_NOC5_MPP_BDY_UPDATE/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90

    r4688 r6225  
    2323   !!   lim_sbc_tau   : update i- and j-stresses, and its modulus at the ocean surface 
    2424   !!---------------------------------------------------------------------- 
    25    USE par_oce          ! ocean parameters 
    26    USE phycst           ! physical constants 
    27    USE par_ice          ! ice parameters 
    28    USE dom_oce          ! ocean domain 
    29    USE dom_ice,    ONLY : tms, area 
    30    USE ice              ! LIM sea-ice variables 
    31    USE sbc_ice          ! Surface boundary condition: sea-ice fields 
    32    USE sbc_oce          ! Surface boundary condition: ocean fields 
    33    USE sbccpl 
    34    USE cpl_oasis3, ONLY : lk_cpl 
    35    USE oce       , ONLY : iatte, oatte, sshn, sshb, snwice_mass, snwice_mass_b, snwice_fmass 
    36    USE albedo           ! albedo parameters 
    37    USE lbclnk           ! ocean lateral boundary condition - MPP exchanges 
    38    USE lib_mpp          ! MPP library 
    39    USE wrk_nemo         ! work arrays 
    40    USE in_out_manager   ! I/O manager 
    41    USE prtctl           ! Print control 
    42    USE lib_fortran      ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)   
    43    USE traqsr           ! clem: add penetration of solar flux into the calculation of heat budget 
    44    USE iom 
    45    USE domvvl           ! Variable volume 
     25   USE par_oce        ! ocean parameters 
     26   USE oce     , ONLY : sshn, sshb, snwice_mass, snwice_mass_b, snwice_fmass 
     27   USE phycst         ! physical constants 
     28   USE dom_oce        ! ocean domain 
     29   USE ice            ! LIM sea-ice variables 
     30   USE sbc_ice        ! Surface boundary condition: sea-ice fields 
     31   USE sbc_oce        ! Surface boundary condition: ocean fields 
     32   USE sbccpl         ! Surface boundary condition: coupled interface 
     33   USE albedo         ! albedo parameters 
     34   USE traqsr         ! add penetration of solar flux in the calculation of heat budget 
     35   USE domvvl         ! Variable volume 
     36   USE limctl         !  
     37   USE limcons        !  
     38   ! 
     39   USE in_out_manager ! I/O manager 
     40   USE iom            ! xIO server 
     41   USE lbclnk         ! ocean lateral boundary condition - MPP exchanges 
     42   USE lib_mpp        ! MPP library 
     43   USE wrk_nemo       ! work arrays 
     44   USE prtctl         ! Print control 
     45   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)   
    4646 
    4747   IMPLICIT NONE 
    4848   PRIVATE 
    4949 
    50    PUBLIC   lim_sbc_init   ! called by ice_init 
     50   PUBLIC   lim_sbc_init   ! called by sbcice_lim 
    5151   PUBLIC   lim_sbc_flx    ! called by sbc_ice_lim 
    5252   PUBLIC   lim_sbc_tau    ! called by sbc_ice_lim 
    53  
    54    REAL(wp)  ::   epsi10 = 1.e-10   ! 
    55    REAL(wp)  ::   epsi20 = 1.e-20   ! 
    5653 
    5754   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   utau_oce, vtau_oce   ! air-ocean surface i- & j-stress     [N/m2] 
     
    6158   !! * Substitutions 
    6259#  include "vectopt_loop_substitute.h90" 
    63 #  include "domzgr_substitute.h90" 
    6460   !!---------------------------------------------------------------------- 
    6561   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011) 
     
    9894      !!              - fr_i    : ice fraction 
    9995      !!              - tn_ice  : sea-ice surface temperature 
    100       !!              - alb_ice : sea-ice alberdo (lk_cpl=T) 
     96      !!              - alb_ice : sea-ice albedo (only useful in coupled mode) 
    10197      !! 
    10298      !! References : Goosse, H. et al. 1996, Bul. Soc. Roy. Sc. Liege, 65, 87-90. 
    10399      !!              Tartinville et al. 2001 Ocean Modelling, 3, 95-108. 
     100      !!              These refs are now obsolete since everything has been revised 
     101      !!              The ref should be Rousset et al., 2015 
    104102      !!--------------------------------------------------------------------- 
    105       INTEGER, INTENT(in) ::   kt    ! number of iteration 
    106       ! 
    107       INTEGER  ::   ji, jj, jl, jk           ! dummy loop indices 
    108       REAL(wp) ::   zinda, zemp      ! local scalars 
    109       REAL(wp) ::   zf_mass         ! Heat flux associated with mass exchange ice->ocean (W.m-2) 
    110       REAL(wp) ::   zfcm1           ! New solar flux received by the ocean 
    111       REAL(wp), POINTER, DIMENSION(:,:,:) ::   zalb, zalbp     ! 2D/3D workspace 
     103      INTEGER, INTENT(in) ::   kt   ! number of iteration 
     104      ! 
     105      INTEGER  ::   ji, jj, jl, jk   ! dummy loop indices 
     106      REAL(wp) ::   zqmass           ! Heat flux associated with mass exchange ice->ocean (W.m-2) 
     107      REAL(wp) ::   zqsr             ! New solar flux received by the ocean 
     108      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zalb_cs, zalb_os     ! 3D workspace 
    112109      !!--------------------------------------------------------------------- 
    113        
    114       IF( lk_cpl )   CALL wrk_alloc( jpi, jpj, jpl, zalb, zalbp ) 
    115  
     110      ! 
    116111      ! make calls for heat fluxes before it is modified 
    117       CALL iom_put( "qsr_oce" , qsr(:,:) * pfrld(:,:) )   !     solar flux at ocean surface 
    118       CALL iom_put( "qns_oce" , qns(:,:) * pfrld(:,:) )   ! non-solar flux at ocean surface 
    119       CALL iom_put( "qsr_ice" , SUM( qsr_ice(:,:,:) * old_a_i(:,:,:), dim=3 ) )  !     solar flux at ice surface 
    120       CALL iom_put( "qns_ice" , SUM( qns_ice(:,:,:) * old_a_i(:,:,:), dim=3 ) )  ! non-solar flux at ice surface 
    121       CALL iom_put( "qtr_ice" , SUM( ftr_ice(:,:,:) * old_a_i(:,:,:), dim=3 ) )  !     solar flux transmitted thru ice 
    122       CALL iom_put( "qt_oce"  , ( qsr(:,:) + qns(:,:) ) * pfrld(:,:) )   
    123       CALL iom_put( "qt_ice"  , SUM( ( qns_ice(:,:,:) + qsr_ice(:,:,:) ) * old_a_i(:,:,:), dim=3 ) ) 
     112      IF( iom_use('qsr_oce') )   CALL iom_put( "qsr_oce" , qsr_oce(:,:) * pfrld(:,:) )                                   !     solar flux at ocean surface 
     113      IF( iom_use('qns_oce') )   CALL iom_put( "qns_oce" , qns_oce(:,:) * pfrld(:,:) + qemp_oce(:,:) )                   ! non-solar flux at ocean surface 
     114      IF( iom_use('qsr_ice') )   CALL iom_put( "qsr_ice" , SUM( qsr_ice(:,:,:) * a_i_b(:,:,:), dim=3 ) )                 !     solar flux at ice surface 
     115      IF( iom_use('qns_ice') )   CALL iom_put( "qns_ice" , SUM( qns_ice(:,:,:) * a_i_b(:,:,:), dim=3 ) + qemp_ice(:,:) ) ! non-solar flux at ice surface 
     116      IF( iom_use('qtr_ice') )   CALL iom_put( "qtr_ice" , SUM( ftr_ice(:,:,:) * a_i_b(:,:,:), dim=3 ) )                 !     solar flux transmitted thru ice 
     117      IF( iom_use('qt_oce' ) )   CALL iom_put( "qt_oce"  , ( qsr_oce(:,:) + qns_oce(:,:) ) * pfrld(:,:) + qemp_oce(:,:) )   
     118      IF( iom_use('qt_ice' ) )   CALL iom_put( "qt_ice"  , SUM( ( qns_ice(:,:,:) + qsr_ice(:,:,:) )   & 
     119         &                                                      * a_i_b(:,:,:), dim=3 ) + qemp_ice(:,:) ) 
     120      IF( iom_use('qemp_oce' ) )   CALL iom_put( "qemp_oce"  , qemp_oce(:,:) )   
     121      IF( iom_use('qemp_ice' ) )   CALL iom_put( "qemp_ice"  , qemp_ice(:,:) )   
    124122 
    125123      ! pfrld is the lead fraction at the previous time step (actually between TRP and THD) 
     
    130128            !      heat flux at the ocean surface      ! 
    131129            !------------------------------------------! 
    132             zinda   = 1._wp - MAX( 0._wp , SIGN( 1._wp , - ( 1._wp - pfrld(ji,jj) ) ) ) ! 1 if ice 
    133  
    134             ! Solar heat flux reaching the ocean = zfcm1 (W.m-2)  
     130            ! Solar heat flux reaching the ocean = zqsr (W.m-2)  
    135131            !--------------------------------------------------- 
    136             IF( lk_cpl ) THEN ! be carfeful: not been tested yet 
    137                ! original line 
    138                zfcm1 = qsr_tot(ji,jj) 
    139                !!!zfcm1 = qsr_tot(ji,jj) + ftr_ice(ji,jj) * ( 1._wp - pfrld(ji,jj) ) / ( 1._wp - zinda + zinda * iatte(ji,jj) ) 
    140                DO jl = 1, jpl 
    141                   zfcm1 = zfcm1 - ( qsr_ice(ji,jj,jl) - ftr_ice(ji,jj,jl) ) * old_a_i(ji,jj,jl) 
    142                END DO 
    143             ELSE 
    144                !!!zfcm1   = pfrld(ji,jj) * qsr(ji,jj)  + & 
    145                !!!     &    ( 1._wp - pfrld(ji,jj) ) * ftr_ice(ji,jj) / ( 1._wp - zinda + zinda * iatte(ji,jj) ) 
    146                zfcm1   = pfrld(ji,jj) * qsr(ji,jj) 
    147                DO jl = 1, jpl 
    148                   zfcm1   = zfcm1 + old_a_i(ji,jj,jl) * ftr_ice(ji,jj,jl) 
    149                END DO 
    150             ENDIF 
     132            zqsr = qsr_tot(ji,jj) 
     133            DO jl = 1, jpl 
     134               zqsr = zqsr - a_i_b(ji,jj,jl) * (  qsr_ice(ji,jj,jl) - ftr_ice(ji,jj,jl) )  
     135            END DO 
    151136 
    152137            ! Total heat flux reaching the ocean = hfx_out (W.m-2)  
    153138            !--------------------------------------------------- 
    154             zf_mass        = hfx_thd(ji,jj) + hfx_dyn(ji,jj) + hfx_res(ji,jj) ! heat flux from snow is 0 (T=0 degC) 
    155             hfx_out(ji,jj) = hfx_out(ji,jj) + zf_mass + zfcm1 
     139            zqmass         = hfx_thd(ji,jj) + hfx_dyn(ji,jj) + hfx_res(ji,jj) ! heat flux from snow is 0 (T=0 degC) 
     140            hfx_out(ji,jj) = hfx_out(ji,jj) + zqmass + zqsr 
     141 
     142            ! Add the residual from heat diffusion equation (W.m-2) 
     143            !------------------------------------------------------- 
     144            hfx_out(ji,jj) = hfx_out(ji,jj) + hfx_err_dif(ji,jj) 
    156145 
    157146            ! New qsr and qns used to compute the oceanic heat flux at the next time step 
    158147            !--------------------------------------------------- 
    159             qsr(ji,jj) = zfcm1                                       
    160             qns(ji,jj) = hfx_out(ji,jj) - zfcm1               
     148            qsr(ji,jj) = zqsr                                       
     149            qns(ji,jj) = hfx_out(ji,jj) - zqsr               
    161150 
    162151            !------------------------------------------! 
     
    171160            !                     Even if i see Ice melting as a FW and SALT flux 
    172161            !         
    173             !  computing freshwater exchanges at the ice/ocean interface 
    174             IF( lk_cpl ) THEN  
    175                zemp = - emp_tot(ji,jj) + emp_ice(ji,jj) * ( 1. - pfrld(ji,jj) )    &   ! 
    176                   &   + wfx_snw(ji,jj) 
    177             ELSE 
    178                zemp =   emp(ji,jj)     *           pfrld(ji,jj)            &   ! evaporation over oceanic fraction 
    179                   &   - tprecip(ji,jj) * ( 1._wp - pfrld(ji,jj) )          &   ! all precipitation reach the ocean 
    180                   &   + sprecip(ji,jj) * ( 1._wp - pfrld(ji,jj)**betas )       ! except solid precip intercepted by sea-ice 
    181             ENDIF 
    182  
    183162            ! mass flux from ice/ocean 
    184             wfx_ice(ji,jj) = wfx_bog(ji,jj) + wfx_bom(ji,jj) + wfx_sum(ji,jj) + wfx_sni(ji,jj) + wfx_opw(ji,jj) + wfx_dyn(ji,jj) + wfx_res(ji,jj) 
     163            wfx_ice(ji,jj) = wfx_bog(ji,jj) + wfx_bom(ji,jj) + wfx_sum(ji,jj) + wfx_sni(ji,jj)   & 
     164                           + wfx_opw(ji,jj) + wfx_dyn(ji,jj) + wfx_res(ji,jj) 
    185165 
    186166            ! mass flux at the ocean/ice interface 
    187             fmmflx(ji,jj) = - wfx_ice(ji,jj) * rdt_ice                   ! F/M mass flux save at least for biogeochemical model 
    188             emp(ji,jj)    = zemp - wfx_ice(ji,jj) - wfx_snw(ji,jj) - wfx_sub(ji,jj)   ! mass flux + F/M mass flux (always ice/ocean mass exchange) 
     167            fmmflx(ji,jj) = - ( wfx_ice(ji,jj) + wfx_snw(ji,jj) ) * r1_rdtice  ! F/M mass flux save at least for biogeochemical model 
     168            emp(ji,jj)    = emp_oce(ji,jj) - wfx_ice(ji,jj) - wfx_snw(ji,jj)   ! mass flux + F/M mass flux (always ice/ocean mass exchange) 
    189169             
    190170         END DO 
     
    194174      !      salt flux at the ocean surface      ! 
    195175      !------------------------------------------! 
    196       sfx(:,:) = sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) + sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:) + sfx_bri(:,:) 
     176      sfx(:,:) = sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) + sfx_opw(:,:)   & 
     177         &     + sfx_res(:,:) + sfx_dyn(:,:) + sfx_bri(:,:) 
    197178 
    198179      !-------------------------------------------------------------! 
     
    203184         snwice_mass_b(:,:) = snwice_mass(:,:)                   
    204185         ! new mass per unit area 
    205          snwice_mass  (:,:) = tms(:,:) * ( rhosn * vt_s(:,:) + rhoic * vt_i(:,:)  )  
     186         snwice_mass  (:,:) = tmask(:,:,1) * ( rhosn * vt_s(:,:) + rhoic * vt_i(:,:)  )  
    206187         ! time evolution of snow+ice mass 
    207188         snwice_fmass (:,:) = ( snwice_mass(:,:) - snwice_mass_b(:,:) ) * r1_rdtice 
     
    214195      tn_ice(:,:,:) = t_su(:,:,:)           ! Ice surface temperature                       
    215196 
    216       !------------------------------------------------! 
    217       !    Computation of snow/ice and ocean albedo    ! 
    218       !------------------------------------------------! 
    219       IF( lk_cpl ) THEN          ! coupled case 
    220          CALL albedo_ice( t_su, ht_i, ht_s, zalbp, zalb )                  ! snow/ice albedo 
    221          alb_ice(:,:,:) =  0.5_wp * zalbp(:,:,:) + 0.5_wp * zalb (:,:,:)   ! Ice albedo (mean clear and overcast skys) 
    222       ENDIF 
    223  
    224  
     197      !------------------------------------------------------------------------! 
     198      !    Snow/ice albedo (only if sent to coupler, useless in forced mode)   ! 
     199      !------------------------------------------------------------------------! 
     200      CALL wrk_alloc( jpi,jpj,jpl,   zalb_cs, zalb_os )     
     201      CALL albedo_ice( t_su, ht_i, ht_s, zalb_cs, zalb_os )  ! cloud-sky and overcast-sky ice albedos 
     202      alb_ice(:,:,:) = ( 1. - cldf_ice ) * zalb_cs(:,:,:) + cldf_ice * zalb_os(:,:,:) 
     203      CALL wrk_dealloc( jpi,jpj,jpl,   zalb_cs, zalb_os ) 
     204 
     205      ! conservation test 
     206      IF( ln_limdiahsb )   CALL lim_cons_final( 'limsbc' ) 
     207 
     208      ! control prints 
     209      IF( ln_icectl )   CALL lim_prt( kt, iiceprt, jiceprt, 3, ' - Final state lim_sbc - ' ) 
     210      ! 
    225211      IF(ln_ctl) THEN 
    226212         CALL prt_ctl( tab2d_1=qsr   , clinfo1=' lim_sbc: qsr    : ', tab2d_2=qns , clinfo2=' qns     : ' ) 
     
    230216      ENDIF 
    231217      ! 
    232       IF( lk_cpl )   CALL wrk_dealloc( jpi, jpj, jpl, zalb, zalbp ) 
    233       !  
    234218   END SUBROUTINE lim_sbc_flx 
    235219 
     
    262246      INTEGER ,                     INTENT(in) ::   kt               ! ocean time-step index 
    263247      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pu_oce, pv_oce   ! surface ocean currents 
    264       !! 
     248      ! 
    265249      INTEGER  ::   ji, jj   ! dummy loop indices 
    266250      REAL(wp) ::   zat_u, zutau_ice, zu_t, zmodt   ! local scalar 
     
    269253      ! 
    270254      IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN     !==  Ice time-step only  ==!   (i.e. surface module time-step) 
    271 !CDIR NOVERRCHK 
    272255         DO jj = 2, jpjm1                             !* update the modulus of stress at ocean surface (T-point) 
    273 !CDIR NOVERRCHK 
    274256            DO ji = fs_2, fs_jpim1 
    275257               !                                               ! 2*(U_ice-U_oce) at T-point 
     
    321303      !! ** input   : Namelist namicedia 
    322304      !!------------------------------------------------------------------- 
    323       REAL(wp) :: zsum, zarea 
    324       ! 
    325       INTEGER  ::   ji, jj, jk                       ! dummy loop indices 
    326       REAL(wp) ::   zcoefu, zcoefv, zcoeff          ! local scalar 
     305      INTEGER  ::   ji, jj, jk               ! dummy loop indices 
     306      REAL(wp) ::   zcoefu, zcoefv, zcoeff   ! local scalar 
     307      !!------------------------------------------------------------------- 
     308      ! 
    327309      IF(lwp) WRITE(numout,*) 
    328310      IF(lwp) WRITE(numout,*) 'lim_sbc_init : LIM-3 sea-ice - surface boundary condition' 
     
    342324         END WHERE 
    343325      ENDIF 
    344       ! clem modif 
    345       IF( .NOT. ln_rstart ) THEN 
    346          iatte(:,:) = 1._wp 
    347          oatte(:,:) = 1._wp 
    348       ENDIF 
    349       ! 
    350       ! clem: snwice_mass in the restart file now 
     326      ! 
    351327      IF( .NOT. ln_rstart ) THEN 
    352328         !                                      ! embedded sea ice 
    353329         IF( nn_ice_embd /= 0 ) THEN            ! mass exchanges between ice and ocean (case 1 or 2) set the snow+ice mass 
    354             snwice_mass  (:,:) = tms(:,:) * ( rhosn * vt_s(:,:) + rhoic * vt_i(:,:)  ) 
     330            snwice_mass  (:,:) = tmask(:,:,1) * ( rhosn * vt_s(:,:) + rhoic * vt_i(:,:)  ) 
    355331            snwice_mass_b(:,:) = snwice_mass(:,:) 
    356332         ELSE 
     
    361337            sshn(:,:) = sshn(:,:) - snwice_mass(:,:) * r1_rau0 
    362338            sshb(:,:) = sshb(:,:) - snwice_mass(:,:) * r1_rau0 
    363 #if defined key_vvl             
    364            ! key_vvl necessary? clem: yes for compilation purpose 
    365             DO jk = 1,jpkm1                     ! adjust initial vertical scale factors 
    366                fse3t_n(:,:,jk) = e3t_0(:,:,jk)*( 1._wp + sshn(:,:)*tmask(:,:,1)/(ht_0(:,:) + 1.0 - tmask(:,:,1)) ) 
    367                fse3t_b(:,:,jk) = e3t_0(:,:,jk)*( 1._wp + sshb(:,:)*tmask(:,:,1)/(ht_0(:,:) + 1.0 - tmask(:,:,1)) ) 
    368             ENDDO 
    369             fse3t_a(:,:,:) = fse3t_b(:,:,:) 
    370             ! Reconstruction of all vertical scale factors at now and before time 
    371             ! steps 
    372             ! ============================================================================= 
    373             ! Horizontal scale factor interpolations 
    374             ! -------------------------------------- 
    375             CALL dom_vvl_interpol( fse3t_b(:,:,:), fse3u_b(:,:,:), 'U' ) 
    376             CALL dom_vvl_interpol( fse3t_b(:,:,:), fse3v_b(:,:,:), 'V' ) 
    377             CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3u_n(:,:,:), 'U' ) 
    378             CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3v_n(:,:,:), 'V' ) 
    379             CALL dom_vvl_interpol( fse3u_n(:,:,:), fse3f_n(:,:,:), 'F' ) 
    380             ! Vertical scale factor interpolations 
    381             ! ------------------------------------ 
    382             CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3w_n (:,:,:), 'W'  ) 
    383             CALL dom_vvl_interpol( fse3u_n(:,:,:), fse3uw_n(:,:,:), 'UW' ) 
    384             CALL dom_vvl_interpol( fse3v_n(:,:,:), fse3vw_n(:,:,:), 'VW' ) 
    385             CALL dom_vvl_interpol( fse3u_b(:,:,:), fse3uw_b(:,:,:), 'UW' ) 
    386             CALL dom_vvl_interpol( fse3v_b(:,:,:), fse3vw_b(:,:,:), 'VW' ) 
    387             ! t- and w- points depth 
    388             ! ---------------------- 
    389             fsdept_n(:,:,1) = 0.5_wp * fse3w_n(:,:,1) 
    390             fsdepw_n(:,:,1) = 0.0_wp 
    391             fsde3w_n(:,:,1) = fsdept_n(:,:,1) - sshn(:,:) 
    392             DO jk = 2, jpk 
    393                fsdept_n(:,:,jk) = fsdept_n(:,:,jk-1) + fse3w_n(:,:,jk) 
    394                fsdepw_n(:,:,jk) = fsdepw_n(:,:,jk-1) + fse3t_n(:,:,jk-1) 
    395                fsde3w_n(:,:,jk) = fsdept_n(:,:,jk  ) - sshn   (:,:) 
    396             END DO 
    397 #endif 
     339 
     340!!gm I really don't like this staff here...  Find a way to put that elsewhere or differently 
     341!!gm 
     342            IF( .NOT.ln_linssh ) THEN 
     343               DO jk = 1,jpkm1                     ! adjust initial vertical scale factors 
     344                  e3t_n(:,:,jk) = e3t_0(:,:,jk)*( 1._wp + sshn(:,:)*tmask(:,:,1)/(ht_0(:,:) + 1.0 - tmask(:,:,1)) ) 
     345                  e3t_b(:,:,jk) = e3t_0(:,:,jk)*( 1._wp + sshb(:,:)*tmask(:,:,1)/(ht_0(:,:) + 1.0 - tmask(:,:,1)) ) 
     346               END DO 
     347               e3t_a(:,:,:) = e3t_b(:,:,:) 
     348               ! Reconstruction of all vertical scale factors at now and before time-steps 
     349               ! ========================================================================= 
     350               ! Horizontal scale factor interpolations 
     351               ! -------------------------------------- 
     352               CALL dom_vvl_interpol( e3t_b(:,:,:), e3u_b(:,:,:), 'U' ) 
     353               CALL dom_vvl_interpol( e3t_b(:,:,:), e3v_b(:,:,:), 'V' ) 
     354               CALL dom_vvl_interpol( e3t_n(:,:,:), e3u_n(:,:,:), 'U' ) 
     355               CALL dom_vvl_interpol( e3t_n(:,:,:), e3v_n(:,:,:), 'V' ) 
     356               CALL dom_vvl_interpol( e3u_n(:,:,:), e3f_n(:,:,:), 'F' ) 
     357               ! Vertical scale factor interpolations 
     358                 ! ------------------------------------ 
     359               CALL dom_vvl_interpol( e3t_n(:,:,:), e3w_n (:,:,:), 'W'  ) 
     360               CALL dom_vvl_interpol( e3u_n(:,:,:), e3uw_n(:,:,:), 'UW' ) 
     361               CALL dom_vvl_interpol( e3v_n(:,:,:), e3vw_n(:,:,:), 'VW' ) 
     362               CALL dom_vvl_interpol( e3u_b(:,:,:), e3uw_b(:,:,:), 'UW' ) 
     363               CALL dom_vvl_interpol( e3v_b(:,:,:), e3vw_b(:,:,:), 'VW' ) 
     364               ! t- and w- points depth 
     365               ! ---------------------- 
     366!!gm not sure of that.... 
     367               gdept_n(:,:,1) = 0.5_wp * e3w_n(:,:,1) 
     368               gdepw_n(:,:,1) = 0.0_wp 
     369               gde3w_n(:,:,1) = gdept_n(:,:,1) - sshn(:,:) 
     370               DO jk = 2, jpk 
     371                  gdept_n(:,:,jk) = gdept_n(:,:,jk-1) + e3w_n(:,:,jk) 
     372                  gdepw_n(:,:,jk) = gdepw_n(:,:,jk-1) + e3t_n(:,:,jk-1) 
     373                  gde3w_n(:,:,jk) = gdept_n(:,:,jk  ) - sshn   (:,:) 
     374               END DO 
     375            ENDIF 
    398376         ENDIF 
    399377      ENDIF ! .NOT. ln_rstart 
    400378      ! 
    401  
    402379   END SUBROUTINE lim_sbc_init 
    403380 
Note: See TracChangeset for help on using the changeset viewer.