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 2222 – NEMO

Changeset 2222


Ignore:
Timestamp:
2010-10-12T15:53:09+02:00 (14 years ago)
Author:
cbricaud
Message:

commit changes from devmercator2010_1

Location:
branches/DEV_r2191_3partymerge2010/NEMO
Files:
1 added
33 edited
1 copied

Legend:

Unmodified
Added
Removed

  • branches/DEV_r2191_3partymerge2010/NEMO/LIM_SRC_2/dom_ice_2.F90

    r2208 r2222  
    44   !! LIM 2.0 Sea Ice :   Domain  variables 
    55   !!====================================================================== 
    6    !! History :   2.0  !  03-08  (C. Ethe)  Free form and module 
     6   !! History :   1.0  !  2003-08  (C. Ethe)  Free form and module 
     7   !!             3.3  !  2009-05 (G.Garric, C. Bricaud) addition of lim2_evp case 
    78   !!---------------------------------------------------------------------- 
    8 #if defined key_lim2 
     9#if defined   key_lim2 
    910   !!---------------------------------------------------------------------- 
    1011   !!   LIM 2.0, UCL-LOCEAN-IPSL (2005) 
    … …  
    2021 
    2122   INTEGER, PUBLIC ::   njeq , njeqm1          !: j-index of the equator if it is inside the domain 
    22       !                                        !  (otherwise = jpj+10 (SH) or -10 (SH) ) 
     23   !                                           !  (otherwise = jpj+10 (SH) or -10 (SH) ) 
    2324 
    24    REAL(wp), PUBLIC, DIMENSION(jpi,jpj)         ::   fs2cor , fcor,   &  !: coriolis factor and coeficient 
    25       &                                              covrai ,         &  !: sine of geographic latitude 
    26       &                                              area   ,         &  !: surface of grid cell  
    27       &                                              tms    , tmu        !: temperature and velocity points masks 
    28    REAL(wp), PUBLIC, DIMENSION(jpi,jpj,2,2)     ::   wght   ,         &  !: weight of the 4 neighbours to compute averages 
    29       &                                              akappa , bkappa     !: first and third group of metric coefficients 
    30    REAL(wp), PUBLIC, DIMENSION(jpi,jpj,2,2,2,2) ::   alambd   !: second group of metric coefficients 
     25   REAL(wp), PUBLIC, DIMENSION(jpi,jpj)         ::   fs2cor , fcor     !: coriolis factor and coeficient 
     26   REAL(wp), PUBLIC, DIMENSION(jpi,jpj)         ::   covrai            !: sine of geographic latitude 
     27   REAL(wp), PUBLIC, DIMENSION(jpi,jpj)         ::   area              !: surface of grid cell  
     28   REAL(wp), PUBLIC, DIMENSION(jpi,jpj)         ::   tms    , tmu      !: temperature and velocity points masks 
     29   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,2,2)     ::   wght              !: weight of the 4 neighbours to compute averages 
    3130 
     31# if defined key_lim2_vp 
     32   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,2,2)     ::   akappa , bkappa   !: first and third group of metric coefficients 
     33   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,2,2,2,2) ::   alambd            !: second group of metric coefficients 
     34# else 
     35   REAL(wp), PUBLIC, DIMENSION(jpi,jpj)         ::   tmv    , tmf      !: y-velocity and F-points masks 
     36   REAL(wp), PUBLIC, DIMENSION(jpi,jpj)         ::   tmi               !: ice mask: =1 if ice thick > 0 
     37# endif 
     38 
     39#else 
     40   !!---------------------------------------------------------------------- 
     41   !!   Default option          Empty module         NO LIM-2 sea-ice model 
     42   !!---------------------------------------------------------------------- 
     43#endif 
     44 
     45   !!---------------------------------------------------------------------- 
     46   !! NEMO/LIM2 3.3, UCL-LOCEAN-IPSL (2010) 
     47   !! $Id$ 
     48   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
    3249   !!====================================================================== 
    33 #endif 
    3450END MODULE dom_ice_2 

  • branches/DEV_r2191_3partymerge2010/NEMO/LIM_SRC_2/ice_2.F90

    r2208 r2222  
    44   !! Sea Ice physics:  diagnostics variables of ice defined in memory 
    55   !!===================================================================== 
    6    !! History :  2.0  !  03-08  (C. Ethe)  F90: Free form and module 
     6   !! History :  2.0  !  2003-08  (C. Ethe)  F90: Free form and module 
     7   !!            3.3  !  2009-05  (G.Garric) addition of the lim2_evp cas 
    78   !!---------------------------------------------------------------------- 
    89#if defined key_lim2 
    … …  
    2021   PRIVATE 
    2122 
    22    !!* Share parameters namelist (namicerun read in iceini) * 
     23   !                                                                     !!* namicerun namelist * 
    2324   CHARACTER(len=32)     , PUBLIC ::   cn_icerst_in  = "restart_ice_in"   !: suffix of ice restart name (input) 
    2425   CHARACTER(len=32)     , PUBLIC ::   cn_icerst_out = "restart_ice"      !: suffix of ice restart name (output) 
    2526   LOGICAL               , PUBLIC ::   ln_limdyn     = .TRUE.             !: flag for ice dynamics (T) or not (F) 
    2627   LOGICAL               , PUBLIC ::   ln_limdmp     = .FALSE.            !: Ice damping 
    27    REAL(wp)              , PUBLIC ::   hsndif        = 0.e0               !: computation of temp. in snow (0) or not (9999) 
    28    REAL(wp)              , PUBLIC ::   hicdif        = 0.e0               !: computation of temp. in ice (0) or not (9999) 
     28   LOGICAL               , PUBLIC ::   ln_nicep      = .TRUE.             !: flag grid points output (T) or not (F) 
     29   REAL(wp)              , PUBLIC ::   hsndif        = 0.e0               !: computation of snow temp. (0) or not (9999) 
     30   REAL(wp)              , PUBLIC ::   hicdif        = 0.e0               !: computation of ice  temp. (0) or not (9999) 
    2931   REAL(wp), DIMENSION(2), PUBLIC ::   acrit = (/ 1.e-06 , 1.e-06 /)      !: minimum fraction for leads in  
    3032   !                                                                      !: north and south hemisphere 
    31    !!* ice-dynamic namelist (namicedyn) * 
     33   !                                       !!* ice-dynamic namelist (namicedyn) * 
    3234   INTEGER , PUBLIC ::   nbiter = 1         !: number of sub-time steps for relaxation 
    3335   INTEGER , PUBLIC ::   nbitdr = 250       !: maximum number of iterations for relaxation 
    … …  
    4648   REAL(wp), PUBLIC ::   ecc    = 2.e0      !: eccentricity of the elliptical yield curve 
    4749   REAL(wp), PUBLIC ::   ahi0   = 350.e0    !: sea-ice hor. eddy diffusivity coeff. (m2/s) 
    48  
     50   INTEGER , PUBLIC ::   nevp   =   360     !: number of EVP subcycling iterations 
     51   INTEGER , PUBLIC ::   telast =  3600     !: timescale for EVP elastic waves 
     52   REAL(wp), PUBLIC ::   alphaevp = 1.e0    !: coefficient for the solution of EVP int. stresses 
    4953   REAL(wp), PUBLIC ::   usecc2             !:  = 1.0 / ( ecc * ecc ) 
    5054   REAL(wp), PUBLIC ::   rhoco              !: = rau0 * cw 
    … …  
    5458   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   ahiu , ahiv   !: hor. diffusivity coeff. at ocean U- and V-points (m2/s) 
    5559   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   pahu , pahv   !: ice hor. eddy diffusivity coef. at ocean U- and V-points 
    56    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   hsnm , hicm   !: mean snow and ice thicknesses 
    57    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   ust2s                 !: friction velocity 
     60   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   ust2s         !: friction velocity 
    5861 
    59    !!* diagnostic quantities 
    60    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   sist          !: Sea-Ice Surface Temperature (Kelvin) 
     62# if defined key_lim2_vp 
     63   !                                                        !!* VP rheology * 
     64   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   hsnm , hicm     !: mean snow and ice thicknesses 
     65   CHARACTER(len=1), PUBLIC :: cl_grid =   'B'                   !: type of grid used in ice dynamics, 'C' or 'B' 
     66# else 
     67   !                                                        !!* EVP rheology * 
     68   CHARACTER(len=1), PUBLIC             ::   cl_grid = 'C'   !: type of grid used in ice dynamics, 'C' or 'B' 
     69   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   stress1_i       !: first stress tensor element        
     70   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   stress2_i       !: second stress tensor element 
     71   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   stress12_i      !: diagonal stress tensor element 
     72   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   delta_i         !: rheology delta factor (see Flato and Hibler 95) [s-1] 
     73   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   divu_i          !: Divergence of the velocity field [s-1] -> limrhg.F90 
     74   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   shear_i         !: Shear of the velocity field [s-1] -> limrhg.F90 
     75   ! 
     76   REAL(wp), PUBLIC, DIMENSION(jpi,jpj)          :: at_i          !: 
     77   REAL(wp), PUBLIC, DIMENSION(:,:)    , POINTER :: vt_s ,vt_i    !: mean snow and ice thicknesses 
     78   REAL(wp), PUBLIC, DIMENSION(jpi,jpj), TARGET  :: hsnm , hicm   !: target vt_s ,vt_i pointers  
     79#endif 
     80 
     81   !                                                      !!* diagnostic quantities 
     82   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   rdvosif       !: ice volume change at ice surface (only used for outputs) 
     83   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   rdvobif       !: ice volume change at ice bottom  (only used for outputs) 
     84   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   fdvolif       !: Total   ice volume change (only used for outputs) 
     85   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   rdvonif       !: Lateral ice volume change (only used for outputs) 
     86   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   sist          !: Sea-Ice Surface Temperature [Kelvin] 
    6187   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   tfu           !: Freezing/Melting point temperature of sea water at SSS 
    6288   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   hicif         !: Ice thickness 
    … …  
    7197   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   rdmicif       !: Variation of ice mass 
    7298   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   qldif         !: heat balance of the lead (or of the open ocean) 
    73    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   qcmif         !: Energy needed to bring the ocean surface layer until its freezing  
     99   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   qcmif         !: Energy needed to bring  ocean surface layer to freezing  
    74100   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   fdtcn         !: net downward heat flux from the ice to the ocean 
    75101   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   qdtcn         !: energy from the ice to the ocean point (at a factor 2) 
    76102   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   thcm          !: part of the solar energy used in the lead heat budget 
    77103   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   fstric        !: Solar flux transmitted trough the ice 
    78    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   ffltbif       !: Array linked with the max heat contained in brine pockets (?) 
     104   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   ffltbif       !: Array linked with the max brine pockets heat content 
    79105   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   fscmbq        !: Linked with the solar flux below the ice (?) 
    80106   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   fsbbq         !: Also linked with the solar flux below the ice (?) 
    81    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   qfvbq         !: Array used to store energy in case of toral lateral ablation (?) 
     107   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   qfvbq         !: Array used to store energy (total lateral ablation case) 
    82108   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   dmgwi         !: Variation of the mass of snow ice 
    83109 
    84    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   u_ice, v_ice   !: two components of the ice   velocity at I-point (m/s) 
    85    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   u_oce, v_oce   !: two components of the ocean velocity at I-point (m/s) 
     110   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   u_ice, v_ice   !: two components of the ice   velocity at I-point [m/s] 
     111   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   u_oce, v_oce   !: two components of the ocean velocity at I-point [m/s] 
    86112 
    87113   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jplayersp1) ::   tbif  !: Temperature inside the ice/snow layer 
    88114 
    89    !!* moment used in the advection scheme 
     115   !                                                                               !!* moment used in the advection 
    90116   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   sxice, syice, sxxice, syyice, sxyice   !: for ice  volume 
    91117   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   sxsn,  sysn,  sxxsn,  syysn,  sxysn    !: for snow volume 
    … …  
    98124#else 
    99125   !!---------------------------------------------------------------------- 
    100    !!   Default option         Empty module        NO LIM 2.0 sea-ice model 
     126   !!   Default option          Empty module         NO LIM-2 sea-ice model 
    101127   !!---------------------------------------------------------------------- 
    102128#endif 
    103129 
     130   !!---------------------------------------------------------------------- 
     131   !! NEMO/LIM2 3.3, UCL-LOCEAN-IPSL (2010) 
     132   !! $Id$ 
     133   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
    104134   !!====================================================================== 
    105135END MODULE ice_2 

  • branches/DEV_r2191_3partymerge2010/NEMO/LIM_SRC_2/iceini_2.F90

    r1581 r2222  
    66   !! History :   1.0  !  02-08  (G. Madec)  F90: Free form and modules 
    77   !!             2.0  !  03-08  (C. Ethe)  add ice_run 
     8   !!             3.3  !  09-05  (G.Garric, C. Bricaud) addition of the lim2_evp case 
    89   !!---------------------------------------------------------------------- 
    910#if defined key_lim2 
    … …  
    1112   !!   'key_lim2' :                                  LIM 2.0 sea-ice model 
    1213   !!---------------------------------------------------------------------- 
    13    !!---------------------------------------------------------------------- 
    1414   !!   ice_init_2       : sea-ice model initialization 
    1515   !!   ice_run_2        : Definition some run parameter for ice model 
    1616   !!---------------------------------------------------------------------- 
    17    USE dom_oce 
    18    USE dom_ice_2 
    19    USE sbc_oce         ! surface boundary condition: ocean 
    20    USE sbc_ice         ! surface boundary condition: ice 
    21    USE phycst          ! Define parameters for the routines 
    22    USE ice_2 
    23    USE limmsh_2 
    24    USE limistate_2 
    25    USE limrst_2    
    26    USE in_out_manager 
     17   USE dom_oce          ! ocean domain 
     18   USE dom_ice_2        ! LIM2: ice domain 
     19   USE sbc_oce          ! surface boundary condition: ocean 
     20   USE sbc_ice          ! surface boundary condition: ice 
     21   USE phycst           ! Define parameters for the routines 
     22   USE ice_2            ! LIM2: ice variable 
     23   USE limmsh_2         ! LIM2: mesh 
     24   USE limistate_2      ! LIM2: initial state 
     25   USE limrst_2         ! LIM2: restart 
     26   USE in_out_manager   ! I/O manager 
    2727       
    2828   IMPLICIT NONE 
    2929   PRIVATE 
    3030 
    31    PUBLIC   ice_init_2               ! called by sbcice_lim_2.F90 
     31   PUBLIC   ice_init_2   ! called by sbcice_lim_2.F90 
     32 
     33   INTEGER, PUBLIC ::   numit   !: iteration number 
    3234 
    3335   !!---------------------------------------------------------------------- 
    34    !!   LIM 2.0,  UCL-LOCEAN-IPSL (2005)  
    35    !! $Id$  
    36    !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt  
     36   !! NEMO/LIM2 3.3, UCL-LOCEAN-IPSL (2010) 
     37   !! $Id$ 
     38   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
    3739   !!---------------------------------------------------------------------- 
    38  
    3940CONTAINS 
    4041 
    … …  
    4344      !!                  ***  ROUTINE ice_init_2  *** 
    4445      !! 
    45       !! ** purpose :    
     46      !! ** purpose :   initialisation of LIM-2 domain and variables  
    4647      !!---------------------------------------------------------------------- 
    4748      ! 
    48       ! Open the namelist file  
     49      !                               ! Open the namelist file  
    4950      CALL ctl_opn( numnam_ice, 'namelist_ice', 'OLD', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp )       
    50       CALL ice_run_2                    !  read in namelist some run parameters 
     51      CALL ice_run_2                  !  read in namelist some run parameters 
    5152                  
    52       ! Louvain la Neuve Ice model 
    53       rdt_ice = nn_fsbc * rdttra(1) 
     53      rdt_ice = nn_fsbc * rdttra(1)   ! sea-ice time step 
     54      numit   = nit000 - 1               
    5455 
    5556      CALL lim_msh_2                  ! ice mesh initialization 
    5657      
    57       ! Initial sea-ice state 
    58       IF( .NOT.ln_rstart ) THEN 
    59          CALL lim_istate_2            ! start from rest: sea-ice deduced from sst 
    60       ELSE 
    61          CALL lim_rst_read_2          ! start from a restart file 
     58      !                               ! ice initialisation (start from rest or from a restart) 
     59      IF( .NOT.ln_rstart ) THEN   ;   CALL lim_istate_2      
     60      ELSE                        ;   CALL lim_rst_read_2 
    6261      ENDIF 
    6362       
    64       tn_ice(:,:,1) = sist(:,:)         ! initialisation of ice temperature    
     63      tn_ice(:,:,1) = sist(:,:)       ! initialisation of ice temperature    
    6564      fr_i  (:,:) = 1.0 - frld(:,:)   ! initialisation of sea-ice fraction     
    6665      ! 
    … …  
    7574      !! 
    7675      !! ** Method  :   Read the namicerun namelist and check the parameter  
    77       !!       values called at the first timestep (nit000) 
     76      !!              values called at the first timestep (nit000) 
    7877      !! 
    7978      !! ** input   :   Namelist namicerun 
    … …  
    8281      !!------------------------------------------------------------------- 
    8382      !                     
    84       REWIND ( numnam_ice )                       ! Read Namelist namicerun  
    85       READ   ( numnam_ice , namicerun ) 
    86  
    87       IF(lwp) THEN 
     83      REWIND( numnam_ice )                   ! Read namicerun namelist 
     84      READ  ( numnam_ice , namicerun ) 
     85      ! 
     86      IF(lwp) THEN                           ! control print 
    8887         WRITE(numout,*) 
    8988         WRITE(numout,*) 'ice_run : ice share parameters for dynamics/advection/thermo of sea-ice' 

  • branches/DEV_r2191_3partymerge2010/NEMO/LIM_SRC_2/limdyn_2.F90

    r2219 r2222  
    88   !!             2.0  !  03-08  (C. Ethe) add lim_dyn_init 
    99   !!             2.0  !  06-07  (G. Madec)  Surface module 
     10   !!             3.3  !  09-05  (G.Garric, C. Bricaud) addition of the lim2_evp case 
    1011   !!--------------------------------------------------------------------- 
    1112#if defined key_lim2 
    … …  
    1617   !!    lim_dyn_init_2 : initialization and namelist read 
    1718   !!---------------------------------------------------------------------- 
    18    USE dom_oce        ! ocean space and time domain 
    19    USE sbc_oce        ! 
    20    USE phycst         ! 
    21    USE ice_2          ! 
    22    USE dom_ice_2      ! 
    23    USE limistate_2    ! 
    24    USE limrhg_2       ! ice rheology 
    25  
    26    USE lbclnk         ! 
    27    USE lib_mpp        ! 
    28    USE in_out_manager ! I/O manager 
    29    USE prtctl         ! Print control 
     19   USE dom_oce          ! ocean domain 
     20   USE sbc_oce          ! surface boundary condition variables 
     21   USE phycst           ! physical constant 
     22   USE ice_2            ! LIM2: ice variables 
     23   USE dom_ice_2        ! LIM2: ice domain 
     24   USE limistate_2      ! LIM2: ice initial state 
     25#if defined key_lim2_vp 
     26   USE limrhg_2         ! LIM2:  VP ice rheology 
     27#else 
     28   USE limrhg           ! LIM : EVP ice rheology 
     29#endif 
     30   USE lbclnk           ! lateral boundary condition - MPP exchanges 
     31   USE lib_mpp          ! MPP library 
     32   USE in_out_manager   ! I/O manager 
     33   USE prtctl           ! Print control 
    3034 
    3135   IMPLICIT NONE 
    3236   PRIVATE 
    3337 
    34    PUBLIC   lim_dyn_2 ! routine called by sbc_ice_lim 
    35  
    36    !! * Module variables 
    37    REAL(wp)  ::  rone    = 1.e0   ! constant value 
     38   PUBLIC   lim_dyn_2   ! routine called by sbc_ice_lim module 
     39 
     40   REAL(wp) ::  rone = 1.e0   ! constant value 
    3841 
    3942#  include "vectopt_loop_substitute.h90" 
    4043   !!---------------------------------------------------------------------- 
    41    !!   LIM 2.0,  UCL-LOCEAN-IPSL (2006)  
     44   !! NEMO/LIM2 3.3, UCL-LOCEAN-IPSL (2010) 
    4245   !! $Id$ 
    4346   !! Software governed by the CeCILL licence  (NEMOGCM/License_CeCILL.txt) 
    4447   !!---------------------------------------------------------------------- 
    45  
    4648CONTAINS 
    4749 
    … …  
    7173      IF( kt == nit000 )   CALL lim_dyn_init_2   ! Initialization (first time-step only) 
    7274       
    73       IF( ln_limdyn ) THEN 
    74          ! 
    75          ! Mean ice and snow thicknesses.           
    76          hsnm(:,:)  = ( 1.0 - frld(:,:) ) * hsnif(:,:) 
     75      IF( ln_limdyn ) THEN                     ! Rheology (ice dynamics) 
     76         !                                     ! ======== 
     77         !          
     78         hsnm(:,:)  = ( 1.0 - frld(:,:) ) * hsnif(:,:)   ! Mean ice and snow thicknesses 
    7779         hicm(:,:)  = ( 1.0 - frld(:,:) ) * hicif(:,:) 
    7880         ! 
    79          !                                     ! Rheology (ice dynamics) 
    80          !                                     ! ======== 
     81         !                                     ! Define the j-limits where ice rheology is computed 
    8182          
    8283         !  Define the j-limits where ice rheology is computed 
    … …  
    8788            i_jpj = jpj 
    8889            IF(ln_ctl)   CALL prt_ctl_info( 'lim_dyn  :    i_j1 = ', ivar1=i_j1, clinfo2=' ij_jpj = ', ivar2=i_jpj ) 
    89             CALL lim_rhg_2( i_j1, i_jpj ) 
     90#if defined key_lim2_vp 
     91            CALL lim_rhg_2( i_j1, i_jpj )             !  VP rheology 
     92#else 
     93            CALL lim_rhg  ( i_j1, i_jpj )             ! EVP rheology 
     94#endif 
     95         ELSE                                      !==  optimization of the computational area  ==! 
     96            DO jj = 1, jpj 
     97               zind(jj) = SUM( frld (:,jj  ) )        ! = FLOAT(jpj) if ocean everywhere on a j-line 
     98               zmsk(jj) = SUM( tmask(:,jj,1) )        ! = 0          if land  everywhere on a j-line 
     99            END DO 
    90100            ! 
    91          ELSE                                 ! optimization of the computational area 
    92             ! 
    93             DO jj = 1, jpj 
    94                zind(jj) = SUM( frld (:,jj  ) )   ! = FLOAT(jpj) if ocean everywhere on a j-line 
    95                zmsk(jj) = SUM( tmask(:,jj,1) )   ! = 0          if land  everywhere on a j-line 
    96             END DO 
    97             ! 
    98             IF( l_jeq ) THEN                     ! local domain include both hemisphere 
    99                !                                 ! Rheology is computed in each hemisphere 
    100                !                                 ! only over the ice cover latitude strip 
    101                ! Northern hemisphere 
    102                i_j1  = njeq 
     101            IF( l_jeq ) THEN                          ! local domain include both hemisphere: rheology is computed 
     102               !                                      ! in each hemisphere only over the ice cover latitude strip 
     103               i_j1  = njeq                  ! Northern hemisphere 
    103104               i_jpj = jpj 
    104105               DO WHILE ( i_j1 <= jpj .AND. zind(i_j1) == FLOAT(jpi) .AND. zmsk(i_j1) /=0 ) 
    105106                  i_j1 = i_j1 + 1 
    106107               END DO 
     108#if defined key_lim2_vp 
    107109               i_j1 = MAX( 1, i_j1-1 ) 
    108110               IF(ln_ctl)   WRITE(numout,*) 'lim_dyn : NH i_j1 = ', i_j1, ' ij_jpj = ', i_jpj 
    109111               !  
    110112               CALL lim_rhg_2( i_j1, i_jpj ) 
    111                !  
    112                ! Southern hemisphere 
    113                i_j1  =  1  
     113#else 
     114               i_j1 = MAX( 1, i_j1-2 ) 
     115               IF(ln_ctl)   WRITE(numout,*) 'lim_dyn : NH i_j1 = ', i_j1, ' ij_jpj = ', i_jpj 
     116               CALL lim_rhg( i_j1, i_jpj ) 
     117#endif 
     118               i_j1  =  1                    ! Southern hemisphere 
    114119               i_jpj = njeq 
    115120               DO WHILE ( i_jpj >= 1 .AND. zind(i_jpj) == FLOAT(jpi) .AND. zmsk(i_jpj) /=0 ) 
    116121                  i_jpj = i_jpj - 1 
    117122               END DO 
     123#if defined key_lim2_vp 
    118124               i_jpj = MIN( jpj, i_jpj+2 ) 
    119125               IF(ln_ctl)   WRITE(numout,*) 'lim_dyn : SH i_j1 = ', i_j1, ' ij_jpj = ', i_jpj 
    120126               !  
    121127               CALL lim_rhg_2( i_j1, i_jpj ) 
     128#else 
     129               i_jpj = MIN( jpj, i_jpj+1 ) 
     130               IF(ln_ctl)   WRITE(numout,*) 'lim_dyn : SH i_j1 = ', i_j1, ' ij_jpj = ', i_jpj 
     131               CALL lim_rhg( i_j1, i_jpj ) !!!!cbr CALL lim_rhg( i_j1, i_jpj, kt ) 
     132#endif 
    122133               !  
    123             ELSE                                 ! local domain extends over one hemisphere only 
    124                !                                 ! Rheology is computed only over the ice cover 
    125                !                                 ! latitude strip 
     134            ELSE                                      ! local domain extends over one hemisphere only: rheology is 
     135               !                                      ! computed only over the ice cover latitude strip 
    126136               i_j1  = 1 
    127137               DO WHILE ( i_j1 <= jpj .AND. zind(i_j1) == FLOAT(jpi) .AND. zmsk(i_j1) /=0 ) 
    … …  
    129139               END DO 
    130140               i_j1 = MAX( 1, i_j1-1 ) 
    131      
    132141               i_jpj  = jpj 
    133142               DO WHILE ( i_jpj >= 1  .AND. zind(i_jpj) == FLOAT(jpi) .AND. zmsk(i_jpj) /=0 ) 
    134143                  i_jpj = i_jpj - 1 
    135144               END DO 
     145               IF(ln_ctl)   WRITE(numout,*) 'lim_dyn : one hemisphere: i_j1 = ', i_j1, ' ij_jpj = ', i_jpj 
     146#if defined key_lim2_vp 
    136147               i_jpj = MIN( jpj, i_jpj+2) 
    137      
    138                IF(ln_ctl)   WRITE(numout,*) 'lim_dyn : one hemisphere: i_j1 = ', i_j1, ' ij_jpj = ', i_jpj 
    139                !  
    140                CALL lim_rhg_2( i_j1, i_jpj ) 
     148               CALL lim_rhg_2( i_j1, i_jpj )                !  VP rheology 
     149#else 
     150               i_j1 = MAX( 1, i_j1-2 ) 
     151               i_jpj = MIN( jpj, i_jpj+1) 
     152               CALL lim_rhg  ( i_j1, i_jpj )                ! EVP rheology 
     153#endif 
    141154               ! 
    142155            ENDIF 
    143156            ! 
    144157         ENDIF 
    145  
     158         ! 
    146159         IF(ln_ctl)   CALL prt_ctl(tab2d_1=u_ice , clinfo1=' lim_dyn  : u_ice :', tab2d_2=v_ice , clinfo2=' v_ice :') 
     160         !  
    147161          
    148          ! computation of friction velocity 
    149          ! -------------------------------- 
    150          ! ice-ocean velocity at U & V-points (u_ice v_ice at I-point ; ssu_m, ssv_m at U- & V-points) 
    151           
    152          DO jj = 1, jpjm1 
    153             DO ji = 1, jpim1   ! NO vector opt. 
    154                zu_io(ji,jj) = 0.5 * ( u_ice(ji+1,jj+1) + u_ice(ji+1,jj  ) ) - ssu_m(ji,jj) 
    155                zv_io(ji,jj) = 0.5 * ( v_ice(ji+1,jj+1) + v_ice(ji  ,jj+1) ) - ssv_m(ji,jj) 
    156             END DO 
    157          END DO 
    158          ! frictional velocity at T-point 
    159          DO jj = 2, jpjm1 
     162         !                                     ! friction velocity 
     163         !                                     ! ================= 
     164         SELECT CASE( cl_grid ) 
     165         CASE( 'C' )                          ! C-grid ice dynamics (EVP) 
     166            zu_io(:,:) = u_ice(:,:) - ssu_m(:,:)      ! ice-ocean & ice velocity at ocean velocity points  
     167            zv_io(:,:) = v_ice(:,:) - ssv_m(:,:) 
     168            ! 
     169         CASE( 'B' )                          ! B-grid ice dynamics (VP) 
     170            DO jj = 1, jpjm1                          ! ice velocity at I-point, ice-ocean velocity at ocean points 
     171               DO ji = 1, jpim1   ! NO vector opt. 
     172                  zu_io(ji,jj) = 0.5 * ( u_ice(ji+1,jj+1) + u_ice(ji+1,jj  ) ) - ssu_m(ji,jj) 
     173                  zv_io(ji,jj) = 0.5 * ( v_ice(ji+1,jj+1) + v_ice(ji  ,jj+1) ) - ssv_m(ji,jj) 
     174               END DO 
     175            END DO 
     176         END SELECT 
     177         ! 
     178         DO jj = 2, jpjm1                     ! frictional velocity at T-point 
    160179            DO ji = 2, jpim1   ! NO vector opt. because of zu_io 
    161180               ust2s(ji,jj) = 0.5 * cw                                                          & 
    … …  
    165184         END DO 
    166185         ! 
    167       ELSE      ! no ice dynamics : transmit directly the atmospheric stress to the ocean 
    168          ! 
     186      ELSE                                     ! no ice dynamics : transmit directly the atmospheric stress to the ocean 
     187         !                                     ! =============== 
    169188         zcoef = SQRT( 0.5 ) / rau0 
    170189         DO jj = 2, jpjm1 
    … …  
    177196      ENDIF 
    178197      ! 
    179       CALL lbc_lnk( ust2s, 'T',  1. )   ! T-point 
     198      CALL lbc_lnk( ust2s, 'T',  1. )   ! lateral boundary condition 
    180199      ! 
    181200      IF(ln_ctl)   CALL prt_ctl(tab2d_1=ust2s , clinfo1=' lim_dyn  : ust2s :') 
    182  
     201      ! 
    183202   END SUBROUTINE lim_dyn_2 
    184203 
    … …  
    188207      !!                  ***  ROUTINE lim_dyn_init_2  *** 
    189208      !! 
    190       !! ** Purpose :   Physical constants and parameters linked to the ice 
    191       !!              dynamics 
    192       !! 
    193       !! ** Method  :   Read the namicedyn namelist and check the ice-dynamic 
    194       !!              parameter values 
     209      !! ** Purpose :   initialisation of  the ice dynamics variables 
     210      !!               
     211      !! ** Method  :   Read the namicedyn namelist and check their values 
    195212      !! 
    196213      !! ** input   :   Namelist namicedyn 
    197214      !!------------------------------------------------------------------- 
    198       NAMELIST/namicedyn/ epsd, alpha,     & 
    199          &                dm, nbiter, nbitdr, om, resl, cw, angvg, pstar,   & 
    200          &                c_rhg, etamn, creepl, ecc, ahi0 
    201       !!------------------------------------------------------------------- 
    202  
     215      NAMELIST/namicedyn/ epsd, alpha, dm, nbiter, nbitdr, om, resl, cw, angvg, pstar,   & 
     216         &                c_rhg, etamn, creepl, ecc, ahi0, nevp, telast, alphaevp 
     217      !!------------------------------------------------------------------- 
     218      ! 
    203219      REWIND ( numnam_ice )                       ! Read Namelist namicedyn 
    204220      READ   ( numnam_ice  , namicedyn ) 
    205  
     221      ! 
    206222      IF(lwp) THEN                                ! Control print 
    207223         WRITE(numout,*) 
    … …  
    223239         WRITE(numout,*) '       eccentricity of the elliptical yield curve       ecc    = ', ecc 
    224240         WRITE(numout,*) '       horizontal diffusivity coeff. for sea-ice        ahi0   = ', ahi0 
     241         WRITE(numout,*) '       number of iterations for subcycling              nevp   = ', nevp 
     242         WRITE(numout,*) '       timescale for elastic waves                      telast = ', telast 
     243         WRITE(numout,*) '       coefficient for the solution of int. stresses  alphaevp = ', alphaevp 
    225244      ENDIF 
    226  
     245      ! 
    227246      !  Initialization 
    228247      usecc2 = 1.0 / ( ecc * ecc ) 
    … …  
    240259#else 
    241260   !!---------------------------------------------------------------------- 
    242    !!   Default option          Empty module       NO LIM 2.0 sea-ice model 
     261   !!   Default option          Dummy module       NO LIM 2.0 sea-ice model 
    243262   !!---------------------------------------------------------------------- 
    244263CONTAINS 
    245    SUBROUTINE lim_dyn_2         ! Empty routine 
     264   SUBROUTINE lim_dyn_2         ! Dummy routine 
    246265   END SUBROUTINE lim_dyn_2 
    247266#endif  

  • branches/DEV_r2191_3partymerge2010/NEMO/LIM_SRC_2/limmsh_2.F90

    r1923 r2222  
    44   !! LIM 2.0 ice model :   definition of the ice mesh parameters 
    55   !!====================================================================== 
    6 #if defined key_lim2 
     6   !! History :  LIM  ! 2001-04  (Louvain-la-Neuve) Original code 
     7   !!            1.0  ! 2002-08  (C. Ethe, G. Madec) 
     8   !!            3.3  ! 2009-05  (G. Garric, C. Bricaud) addition of the lim2_evp case 
     9   !!---------------------------------------------------------------------- 
     10#if defined   key_lim2 
    711   !!---------------------------------------------------------------------- 
    812   !!   'key_lim2'                                     LIM 2.0sea-ice model 
    … …  
    1014   !!   lim_msh_2   : definition of the ice mesh 
    1115   !!---------------------------------------------------------------------- 
    12    !! * Modules used 
    13    USE phycst 
    14    USE dom_oce 
    15    USE dom_ice_2 
    16    USE lbclnk 
    17    USE in_out_manager 
     16   USE phycst           ! physical constants 
     17   USE dom_oce          ! ocean domain 
     18   USE dom_ice_2        ! LIM2: ice domain 
     19   USE lbclnk           ! lateral boundary condition - MPP exchanges 
     20   USE in_out_manager   ! I/O manager 
    1821 
    1922   IMPLICIT NONE 
    2023   PRIVATE 
    2124 
    22    !! * Accessibility 
    23    PUBLIC lim_msh_2      ! routine called by ice_ini_2.F90 
    24  
    25    !!---------------------------------------------------------------------- 
    26    !!   LIM 2.0,  UCL-LOCEAN-IPSL (2005)  
     25   PUBLIC   lim_msh_2   ! routine called by ice_ini_2.F90 
     26 
     27   !!---------------------------------------------------------------------- 
     28   !! NEMO/LIM2 3.3, UCL-LOCEAN-IPSL (2010) 
    2729   !! $Id$ 
    28    !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt  
    29    !!---------------------------------------------------------------------- 
    30  
     30   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
     31   !!---------------------------------------------------------------------- 
    3132CONTAINS 
    3233 
    … …  
    4243      !!              - Initialization of the ice masks (tmsk, umsk) 
    4344      !!  
    44       !! ** Refer.  : Deleersnijder et al. Ocean Modelling 100, 7-10  
    45       !! 
    46       !! ** History : 
    47       !!         original    : 01-04 (LIM) 
    48       !!         addition    : 02-08 (C. Ethe, G. Madec) 
     45      !! References : Deleersnijder et al. Ocean Modelling 100, 7-10  
    4946      !!---------------------------------------------------------------------  
    50       !! * Local variables 
    5147      INTEGER :: ji, jj      ! dummy loop indices 
    52  
    53       REAL(wp), DIMENSION(jpi,jpj) ::  & 
    54          zd2d1 , zd1d2       ! Derivative of zh2 (resp. zh1) in the x direction 
    55          !                   ! (resp. y direction) (defined at the center) 
    56       REAL(wp) ::         & 
    57          zh1p  , zh2p   , &  ! Idem zh1, zh2 for the bottom left corner of the grid 
    58          zd2d1p, zd1d2p , &  ! Idem zd2d1, zd1d2 for the bottom left corner of the grid 
    59          zusden, zusden2     ! temporary scalars 
     48      REAL(wp) ::   zusden   ! local scalars 
     49#if defined key_lim2_vp 
     50      REAL(wp) ::   zusden2           ! local scalars 
     51      REAL(wp) ::   zh1p  , zh2p      !   -      - 
     52      REAL(wp) ::   zd2d1p, zd1d2p    !   -      - 
     53      REAL(wp), DIMENSION(jpi,jpj) ::   zd2d1 , zd1d2   ! 2D workspace 
     54#endif 
    6055      !!--------------------------------------------------------------------- 
    61  
     56      ! 
    6257      IF(lwp) THEN 
    6358         WRITE(numout,*) 
    … …  
    7671      njeqm1 = njeq - 1  
    7772 
    78       fcor(:,:) = 2. * omega * SIN( gphit(:,:) * rad )   !  coriolis factor 
     73      fcor(:,:) = 2. * omega * SIN( gphit(:,:) * rad )   !  coriolis factor at T-point 
    7974  
    80 !i    DO jj = 1, jpj 
    81 !i       zmsk(jj) = SUM( tmask(:,jj,:) )   ! = 0          if land  everywhere on a j-line 
    82 !!ii     write(numout,*) jj, zind(jj) 
    83 !i    END DO 
    84  
    8575      IF( fcor(1,1) * fcor(1,nlcj) < 0.e0 ) THEN   ! local domain include both hemisphere 
    8676         l_jeq = .TRUE. 
    … …  
    115105      !------------------- 
    116106!!ibug ??? 
    117       akappa(:,:,:,:) = 0.e0 
    118107      wght(:,:,:,:) = 0.e0 
     108      tmu(:,:)      = 0.e0 
     109#if defined key_lim2_vp  
     110      akappa(:,:,:,:)     = 0.e0 
    119111      alambd(:,:,:,:,:,:) = 0.e0 
    120       tmu(:,:) = 0.e0 
     112#else 
     113      tmv(:,:) = 0.e0 
     114      tmf(:,:) = 0.e0 
     115#endif 
    121116!!i 
    122117       
    123        
    124       ! metric coefficients for sea ice dynamic 
     118#if defined key_lim2_vp 
     119      ! metric coefficients for sea ice dynamic (VP rheology) 
    125120      !---------------------------------------- 
    126121      !                                                       ! akappa 
    … …  
    128123         zd1d2(:,jj) = e1v(:,jj) - e1v(:,jj-1) 
    129124      END DO 
    130       CALL lbc_lnk( zd1d2, 'T', -1. ) 
    131  
    132125      DO ji = 2, jpi 
    133126         zd2d1(ji,:) = e2u(ji,:) - e2u(ji-1,:) 
    134127      END DO 
    135       CALL lbc_lnk( zd2d1, 'T', -1. ) 
    136  
     128      CALL lbc_lnk( zd1d2, 'T', -1. )   ;   CALL lbc_lnk( zd2d1, 'T', -1. )      ! lateral boundary condition 
     129      ! 
    137130      akappa(:,:,1,1) =        1.0 / ( 2.0 * e1t(:,:) ) 
    138131      akappa(:,:,1,2) = zd1d2(:,:) / ( 4.0 * e1t(:,:) * e2t(:,:) ) 
    … …  
    140133      akappa(:,:,2,2) =        1.0 / ( 2.0 * e2t(:,:) ) 
    141134       
    142       !                                                      ! weights (wght) 
    143       DO jj = 2, jpj 
     135      !              
     136      DO jj = 2, jpj                                         ! weights (wght) at I-points 
    144137         DO ji = 2, jpi 
    145138            zusden = 1. / (  ( e1t(ji,jj) + e1t(ji-1,jj  ) )   & 
    … …  
    155148      CALL lbc_lnk( wght(:,:,2,1), 'I', 1. )      ! but it is never used 
    156149      CALL lbc_lnk( wght(:,:,2,2), 'I', 1. ) 
     150#else 
     151      ! metric coefficients for sea ice dynamic (EVP rheology) 
     152      !---------------------------------------- 
     153      DO jj = 1, jpjm1                                       ! weights (wght) at F-points 
     154         DO ji = 1, jpim1 
     155            zusden = 1. / (  ( e1t(ji+1,jj  ) + e1t(ji,jj) )   & 
     156               &           * ( e2t(ji  ,jj+1) + e2t(ji,jj) ) )  
     157            wght(ji,jj,1,1) = zusden * e1t(ji+1,jj) * e2t(ji,jj+1) 
     158            wght(ji,jj,1,2) = zusden * e1t(ji+1,jj) * e2t(ji,jj  ) 
     159            wght(ji,jj,2,1) = zusden * e1t(ji  ,jj) * e2t(ji,jj+1) 
     160            wght(ji,jj,2,2) = zusden * e1t(ji  ,jj) * e2t(ji,jj  ) 
     161         END DO 
     162      END DO 
     163      CALL lbc_lnk( wght(:,:,1,1), 'F', 1. )   ;   CALL lbc_lnk( wght(:,:,1,2), 'F', 1. )       ! lateral boundary cond.    
     164      CALL lbc_lnk( wght(:,:,2,1), 'F', 1. )   ;   CALL lbc_lnk( wght(:,:,2,2), 'F', 1. ) 
     165#endif 
    157166     
    158167      ! Coefficients for divergence of the stress tensor 
    159168      !------------------------------------------------- 
    160  
    161       DO jj = 2, jpj 
     169#if defined key_lim2_vp 
     170      DO jj = 2, jpj                                     ! VP rheology 
    162171         DO ji = 2, jpi   ! NO vector opt. 
    163             zh1p  =  e1t(ji  ,jj  ) * wght(ji,jj,2,2)   & 
    164                &   + e1t(ji-1,jj  ) * wght(ji,jj,1,2)   & 
    165                &   + e1t(ji  ,jj-1) * wght(ji,jj,2,1)   & 
    166                &   + e1t(ji-1,jj-1) * wght(ji,jj,1,1) 
    167  
    168             zh2p  =  e2t(ji  ,jj  ) * wght(ji,jj,2,2)   & 
    169                &   + e2t(ji-1,jj  ) * wght(ji,jj,1,2)   & 
    170                &   + e2t(ji  ,jj-1) * wght(ji,jj,2,1)   & 
    171                &   + e2t(ji-1,jj-1) * wght(ji,jj,1,1) 
    172  
    173 ! better written but change the last digit and thus solver in less than 100 timestep 
    174 !           zh1p  = e1t(ji-1,jj  ) * wght(ji,jj,1,2) + e1t(ji,jj  ) * wght(ji,jj,2,2)   & 
    175 !              &  + e1t(ji-1,jj-1) * wght(ji,jj,1,1) + e1t(ji,jj-1) * wght(ji,jj,2,1)  
    176  
    177 !           zh2p  = e2t(ji-1,jj  ) * wght(ji,jj,1,2) + e2t(ji,jj  ) * wght(ji,jj,2,2)   & 
    178 !              &  + e2t(ji-1,jj-1) * wght(ji,jj,1,1) + e2t(ji,jj-1) * wght(ji,jj,2,1) 
    179  
    180 !!ibug =0   zusden = 1.0 / ( zh1p * zh2p * 4.e0 ) 
    181             zusden = 1.0 / MAX( zh1p * zh2p * 4.e0 , 1.e-20 ) 
     172            zh1p  =  e1t(ji  ,jj  ) * wght(ji,jj,2,2) + e1t(ji-1,jj  ) * wght(ji,jj,1,2)   & 
     173               &   + e1t(ji  ,jj-1) * wght(ji,jj,2,1) + e1t(ji-1,jj-1) * wght(ji,jj,1,1) 
     174            zh2p  =  e2t(ji  ,jj  ) * wght(ji,jj,2,2) + e2t(ji-1,jj  ) * wght(ji,jj,1,2)   & 
     175               &   + e2t(ji  ,jj-1) * wght(ji,jj,2,1) + e2t(ji-1,jj-1) * wght(ji,jj,1,1) 
     176            ! 
     177            zusden  = 1.e0 / MAX( zh1p * zh2p * 4.e0 , 1.e-20 ) 
    182178            zusden2 = zusden * 2.0  
    183  
    184             zd1d2p = zusden * 0.5 * ( -e1t(ji-1,jj-1) + e1t(ji-1,jj  ) - e1t(ji,jj-1) + e1t(ji  ,jj)   ) 
    185             zd2d1p = zusden * 0.5 * (  e2t(ji  ,jj-1) - e2t(ji-1,jj-1) + e2t(ji,jj  ) - e2t(ji-1,jj)   ) 
    186  
     179            zd1d2p  = zusden * 0.5 * ( -e1t(ji-1,jj-1) + e1t(ji-1,jj  ) - e1t(ji,jj-1) + e1t(ji  ,jj)   ) 
     180            zd2d1p  = zusden * 0.5 * (  e2t(ji  ,jj-1) - e2t(ji-1,jj-1) + e2t(ji,jj  ) - e2t(ji-1,jj)   ) 
     181            ! 
    187182            alambd(ji,jj,2,2,2,1) = zusden2 * e2t(ji  ,jj-1) 
    188183            alambd(ji,jj,2,2,2,2) = zusden2 * e2t(ji  ,jj  ) 
    189184            alambd(ji,jj,2,2,1,1) = zusden2 * e2t(ji-1,jj-1) 
    190185            alambd(ji,jj,2,2,1,2) = zusden2 * e2t(ji-1,jj  ) 
    191  
     186            ! 
    192187            alambd(ji,jj,1,1,2,1) = zusden2 * e1t(ji  ,jj-1) 
    193188            alambd(ji,jj,1,1,2,2) = zusden2 * e1t(ji  ,jj  ) 
    194189            alambd(ji,jj,1,1,1,1) = zusden2 * e1t(ji-1,jj-1) 
    195190            alambd(ji,jj,1,1,1,2) = zusden2 * e1t(ji-1,jj  ) 
    196  
     191            ! 
    197192            alambd(ji,jj,1,2,2,1) = zd1d2p 
    198193            alambd(ji,jj,1,2,2,2) = zd1d2p 
    199194            alambd(ji,jj,1,2,1,1) = zd1d2p 
    200195            alambd(ji,jj,1,2,1,2) = zd1d2p 
    201  
     196            ! 
    202197            alambd(ji,jj,2,1,2,1) = zd2d1p 
    203198            alambd(ji,jj,2,1,2,2) = zd2d1p 
    … …  
    206201         END DO 
    207202      END DO 
    208  
    209       CALL lbc_lnk( alambd(:,:,2,2,2,1), 'I', 1. )      ! CAUTION: even with the lbc_lnk at ice U-V point 
    210       CALL lbc_lnk( alambd(:,:,2,2,2,2), 'I', 1. )      ! the value of wght at jpj is wrong 
    211       CALL lbc_lnk( alambd(:,:,2,2,1,1), 'I', 1. )      ! but it is never used 
    212       CALL lbc_lnk( alambd(:,:,2,2,1,2), 'I', 1. )      !  
    213  
    214       CALL lbc_lnk( alambd(:,:,1,1,2,1), 'I', 1. )      ! CAUTION: idem 
    215       CALL lbc_lnk( alambd(:,:,1,1,2,2), 'I', 1. )      !  
    216       CALL lbc_lnk( alambd(:,:,1,1,1,1), 'I', 1. )      ! 
    217       CALL lbc_lnk( alambd(:,:,1,1,1,2), 'I', 1. )      ! 
    218  
    219       CALL lbc_lnk( alambd(:,:,1,2,2,1), 'I', 1. )      ! CAUTION: idem 
    220       CALL lbc_lnk( alambd(:,:,1,2,2,2), 'I', 1. )      ! 
    221       CALL lbc_lnk( alambd(:,:,1,2,1,1), 'I', 1. )      ! 
    222       CALL lbc_lnk( alambd(:,:,1,2,1,2), 'I', 1. )      ! 
    223  
    224       CALL lbc_lnk( alambd(:,:,2,1,2,1), 'I', 1. )      ! CAUTION: idem 
    225       CALL lbc_lnk( alambd(:,:,2,1,2,2), 'I', 1. )      ! 
    226       CALL lbc_lnk( alambd(:,:,2,1,1,1), 'I', 1. )      ! 
    227       CALL lbc_lnk( alambd(:,:,2,1,1,2), 'I', 1. )      ! 
     203      ! 
     204      ! lateral boundary conditions 
     205      ! CAUTION: even with the lbc_lnk at ice U-V point, the value of wght at jpj is wrong but it is never used 
     206      CALL lbc_lnk( alambd(:,:,2,2,1,1), 'I', 1. )   ;   CALL lbc_lnk( alambd(:,:,2,2,2,1), 'I', 1. )  
     207      CALL lbc_lnk( alambd(:,:,2,2,1,2), 'I', 1. )   ;   CALL lbc_lnk( alambd(:,:,2,2,2,2), 'I', 1. )  
     208      ! 
     209      CALL lbc_lnk( alambd(:,:,1,1,2,2), 'I', 1. )   ;   CALL lbc_lnk( alambd(:,:,1,1,2,1), 'I', 1. ) 
     210      CALL lbc_lnk( alambd(:,:,1,1,1,1), 'I', 1. )   ;   CALL lbc_lnk( alambd(:,:,1,1,1,2), 'I', 1. ) 
     211      ! 
     212      CALL lbc_lnk( alambd(:,:,1,2,1,1), 'I', 1. )   ;   CALL lbc_lnk( alambd(:,:,1,2,2,1), 'I', 1. ) 
     213      CALL lbc_lnk( alambd(:,:,1,2,1,2), 'I', 1. )   ;   CALL lbc_lnk( alambd(:,:,1,2,2,2), 'I', 1. ) 
     214      ! 
     215      CALL lbc_lnk( alambd(:,:,2,1,1,1), 'I', 1. )   ;   CALL lbc_lnk( alambd(:,:,2,1,2,2), 'I', 1. ) 
     216      CALL lbc_lnk( alambd(:,:,2,1,1,2), 'I', 1. )   ;   CALL lbc_lnk( alambd(:,:,2,1,2,1), 'I', 1. ) 
     217#endif 
    228218             
    229219 
    230220      ! Initialization of ice masks 
    231221      !---------------------------- 
    232        
     222      !       
    233223      tms(:,:) = tmask(:,:,1)      ! ice T-point  : use surface tmask 
    234  
    235 !i here we can use umask with a i and j shift of -1,-1 
    236       tmu(:,1) = 0.e0 
     224      ! 
     225#if defined key_lim2_vp 
     226      ! VP rheology : ice velocity point is I-point 
     227      tmu(:,1) = 0.e0              !  
    237228      tmu(1,:) = 0.e0 
    238       DO jj = 2, jpj               ! ice U.V-point: computed from ice T-point mask 
     229      DO jj = 2, jpj 
    239230         DO ji = 2, jpim1   ! NO vector opt. 
    240231            tmu(ji,jj) =  tms(ji,jj) * tms(ji-1,jj) * tms(ji,jj-1) * tms(ji-1,jj-1)             
    241232         END DO 
    242233      END DO 
    243        
    244       !--lateral boundary conditions     
    245       CALL lbc_lnk( tmu(:,:), 'I', 1. ) 
    246        
    247       ! unmasked and masked area of T-grid cell 
    248       area(:,:) = e1t(:,:) * e2t(:,:) 
    249        
     234      CALL lbc_lnk( tmu(:,:), 'I', 1. )      ! lateral boundary conditions     
     235#else 
     236      ! EVP rheology : ice velocity point are U- & V-points ; ice vorticity point is F-point 
     237      tmu(:,:) = umask(:,:,1) 
     238      tmv(:,:) = vmask(:,:,1) 
     239      tmf(:,:) = 0.e0                        ! used of fmask except its special value along the coast (rn_shlat) 
     240      WHERE( fmask(:,:,1) == 1.e0 )   tmf(:,:) = 1.e0 
     241#endif 
     242      ! 
     243      area(:,:) = e1t(:,:) * e2t(:,:)        ! unmasked and masked area of T-grid cell 
     244      ! 
    250245   END SUBROUTINE lim_msh_2 
    251246 

  • branches/DEV_r2191_3partymerge2010/NEMO/LIM_SRC_2/limrhg_2.F90

    r2208 r2222  
    44   !!   Ice rheology :  performs sea ice rheology 
    55   !!====================================================================== 
    6    !! History :  0.0  !  93-12  (M.A. Morales Maqueda.)  Original code 
    7    !!            1.0  !  94-12  (H. Goosse)  
    8    !!            2.0  !  03-12  (C. Ethe, G. Madec)  F90, mpp 
    9    !!            " "  !  06-08  (G. Madec)  surface module, ice-stress at I-point 
    10    !!            " "  !  09-09  (G. Madec)  Huge verctor optimisation 
    11    !!---------------------------------------------------------------------- 
    12 #if defined key_lim2 
    13    !!---------------------------------------------------------------------- 
    14    !!   'key_lim2'                                    LIM 2.0 sea-ice model 
    15    !!---------------------------------------------------------------------- 
     6   !! History :  LIM  !  1993-12  (M.A. Morales Maqueda.)  Original code 
     7   !!            1.0  !  1994-12  (H. Goosse)  
     8   !!            2.0  !  2003-12  (C. Ethe, G. Madec)  F90, mpp 
     9   !!             -   !  2006-08  (G. Madec)  surface module, ice-stress at I-point 
     10   !!             -   !  2009-09  (G. Madec)  Huge verctor optimisation 
     11   !!            3.3  !  2009-05  (G.Garric, C. Bricaud) addition of the lim2_evp case 
     12   !!---------------------------------------------------------------------- 
     13#if defined   key_lim2   &&   defined key_lim2_vp 
     14   !!---------------------------------------------------------------------- 
     15   !!   'key_lim2'                and                 LIM 2.0 sea-ice model 
     16   !!   'key_lim2_vp'                                       VP ice rheology 
    1617   !!---------------------------------------------------------------------- 
    1718   !!   lim_rhg_2   : computes ice velocities 
    … …  
    2122   USE sbc_oce        ! surface boundary condition: ocean variables 
    2223   USE sbc_ice        ! surface boundary condition: ice variables 
    23    USE dom_ice_2      ! domaine: ice variables 
    24    USE phycst         ! physical constant 
    25    USE ice_2          ! ice variables 
    26    USE lbclnk         ! lateral boundary condition 
     24   USE dom_ice_2      ! LIM2: ice domain 
     25   USE phycst         ! physical constants 
     26   USE ice_2          ! LIM2: ice variables 
     27   USE lbclnk         ! lateral boundary condition - MPP exchanges 
    2728   USE lib_mpp        ! MPP library 
    2829   USE in_out_manager ! I/O manager 
    … …  
    3233   PRIVATE 
    3334 
    34    PUBLIC   lim_rhg_2 ! routine called by lim_dyn 
     35   PUBLIC   lim_rhg_2   ! routine called by lim_dyn 
    3536 
    3637   REAL(wp) ::   rzero   = 0.e0   ! constant value: zero 
    … …  
    4041#  include "vectopt_loop_substitute.h90" 
    4142   !!---------------------------------------------------------------------- 
    42    !!   LIM 2.0,  UCL-LOCEAN-IPSL (2006)  
     43   !! NEMO/LIM2 3.3, UCL-LOCEAN-IPSL (2010) 
    4344   !! $Id$ 
    4445   !! Software governed by the CeCILL licence  (NEMOGCM/License_CeCILL.txt) 
    … …  
    8990      REAL(wp), DIMENSION(jpi,0:jpj+1) ::   zzfrld, zztms 
    9091      REAL(wp), DIMENSION(jpi,0:jpj+1) ::   zi1, zi2, zmasst, zpresh 
    91  
    9292      !!------------------------------------------------------------------- 
    93  
    94 !!bug 
    95 !!    u_oce(:,:) = 0.e0 
    96 !!    v_oce(:,:) = 0.e0 
    97 !!    write(*,*) 'rhg min, max u & v', maxval(u_oce), minval(u_oce), maxval(v_oce), minval(v_oce) 
    98 !!bug 
    9993       
    10094      !  Store initial velocities 
    … …  
    106100      zztms(:,1:jpj) = tms(:,:)   ;    zzfrld(:,1:jpj) = frld(:,:) 
    107101      zu0(:,1:jpj) = u_ice(:,:)   ;    zv0(:,1:jpj) = v_ice(:,:) 
    108  
     102      ! 
    109103      zu_a(:,:)    = zu0(:,:)     ;   zv_a(:,:) = zv0(:,:) 
    110104      zu_n(:,:)    = zu0(:,:)     ;   zv_n(:,:) = zv0(:,:) 
    … …  
    126120      zviszeta(:,:) = 0.e0 
    127121      zviseta (:,:) = 0.e0 
    128  
    129 !i    zviszeta(:,0    ) = 0.e0    ;    zviseta(:,0    ) = 0.e0 
    130 !i    zviszeta(:,jpj  ) = 0.e0    ;    zviseta(:,jpj  ) = 0.e0 
    131 !i    zviszeta(:,jpj+1) = 0.e0    ;    zviseta(:,jpj+1) = 0.e0 
    132122 
    133123 
    … …  
    372362            END DO 
    373363         END DO 
    374  
    375          ! GAUSS-SEIDEL method 
     364         !  
    376365         !                                                      ! ================ ! 
    377 iflag:   DO jter = 1 , nbitdr                                   !    Relaxation    ! 
     366iflag:   DO jter = 1 , nbitdr                                   !    Relaxation    ! GAUSS-SEIDEL method 
    378367            !                                                   ! ================ ! 
    379368!CDIR NOVERRCHK 

  • branches/DEV_r2191_3partymerge2010/NEMO/LIM_SRC_2/limrst_2.F90

    r2208 r2222  
    66   !! History :  2.0  !  01-04  (C. Ethe, G. Madec)  Original code 
    77   !!                 !  06-07  (S. Masson)  use IOM for restart read/write 
     8   !!            3.3  !  09-05  (G.Garric) addition of the lim2_evp case 
    89   !!---------------------------------------------------------------------- 
    910#if defined key_lim2 
    … …  
    108109      CALL iom_rstput( iter, nitrst, numriw, 'kt_ice', REAL( iter, wp) )  
    109110       
    110       CALL iom_rstput( iter, nitrst, numriw, 'hicif' , hicif (:,:)   )      ! prognostic variables  
    111       CALL iom_rstput( iter, nitrst, numriw, 'hsnif' , hsnif (:,:)   ) 
    112       CALL iom_rstput( iter, nitrst, numriw, 'frld'  , frld  (:,:)   ) 
    113       CALL iom_rstput( iter, nitrst, numriw, 'sist'  , sist  (:,:)   ) 
    114       CALL iom_rstput( iter, nitrst, numriw, 'tbif1' , tbif  (:,:,1) ) 
    115       CALL iom_rstput( iter, nitrst, numriw, 'tbif2' , tbif  (:,:,2) ) 
    116       CALL iom_rstput( iter, nitrst, numriw, 'tbif3' , tbif  (:,:,3) ) 
    117       CALL iom_rstput( iter, nitrst, numriw, 'u_ice' , u_ice (:,:)   ) 
    118       CALL iom_rstput( iter, nitrst, numriw, 'v_ice' , v_ice (:,:)   ) 
    119       CALL iom_rstput( iter, nitrst, numriw, 'qstoif', qstoif(:,:)   ) 
    120       CALL iom_rstput( iter, nitrst, numriw, 'fsbbq' , fsbbq (:,:)   ) 
    121       CALL iom_rstput( iter, nitrst, numriw, 'sxice' , sxice (:,:)   ) 
    122       CALL iom_rstput( iter, nitrst, numriw, 'syice' , syice (:,:)   ) 
    123       CALL iom_rstput( iter, nitrst, numriw, 'sxxice', sxxice(:,:)   ) 
    124       CALL iom_rstput( iter, nitrst, numriw, 'syyice', syyice(:,:)   ) 
    125       CALL iom_rstput( iter, nitrst, numriw, 'sxyice', sxyice(:,:)   ) 
    126       CALL iom_rstput( iter, nitrst, numriw, 'sxsn'  , sxsn  (:,:)   ) 
    127       CALL iom_rstput( iter, nitrst, numriw, 'sysn'  , sysn  (:,:)   ) 
    128       CALL iom_rstput( iter, nitrst, numriw, 'sxxsn' , sxxsn (:,:)   ) 
    129       CALL iom_rstput( iter, nitrst, numriw, 'syysn' , syysn (:,:)   ) 
    130       CALL iom_rstput( iter, nitrst, numriw, 'sxysn' , sxysn (:,:)   ) 
    131       CALL iom_rstput( iter, nitrst, numriw, 'sxa'   , sxa   (:,:)   ) 
    132       CALL iom_rstput( iter, nitrst, numriw, 'sya'   , sya   (:,:)   ) 
    133       CALL iom_rstput( iter, nitrst, numriw, 'sxxa'  , sxxa  (:,:)   ) 
    134       CALL iom_rstput( iter, nitrst, numriw, 'syya'  , syya  (:,:)   ) 
    135       CALL iom_rstput( iter, nitrst, numriw, 'sxya'  , sxya  (:,:)   ) 
    136       CALL iom_rstput( iter, nitrst, numriw, 'sxc0'  , sxc0  (:,:)   ) 
    137       CALL iom_rstput( iter, nitrst, numriw, 'syc0'  , syc0  (:,:)   ) 
    138       CALL iom_rstput( iter, nitrst, numriw, 'sxxc0' , sxxc0 (:,:)   ) 
    139       CALL iom_rstput( iter, nitrst, numriw, 'syyc0' , syyc0 (:,:)   ) 
    140       CALL iom_rstput( iter, nitrst, numriw, 'sxyc0' , sxyc0 (:,:)   ) 
    141       CALL iom_rstput( iter, nitrst, numriw, 'sxc1'  , sxc1  (:,:)   ) 
    142       CALL iom_rstput( iter, nitrst, numriw, 'syc1'  , syc1  (:,:)   ) 
    143       CALL iom_rstput( iter, nitrst, numriw, 'sxxc1' , sxxc1 (:,:)   ) 
    144       CALL iom_rstput( iter, nitrst, numriw, 'syyc1' , syyc1 (:,:)   ) 
    145       CALL iom_rstput( iter, nitrst, numriw, 'sxyc1' , sxyc1 (:,:)   ) 
    146       CALL iom_rstput( iter, nitrst, numriw, 'sxc2'  , sxc2  (:,:)   ) 
    147       CALL iom_rstput( iter, nitrst, numriw, 'syc2'  , syc2  (:,:)   ) 
    148       CALL iom_rstput( iter, nitrst, numriw, 'sxxc2' , sxxc2 (:,:)   ) 
    149       CALL iom_rstput( iter, nitrst, numriw, 'syyc2' , syyc2 (:,:)   ) 
    150       CALL iom_rstput( iter, nitrst, numriw, 'sxyc2' , sxyc2 (:,:)   ) 
    151       CALL iom_rstput( iter, nitrst, numriw, 'sxst'  , sxst  (:,:)   ) 
    152       CALL iom_rstput( iter, nitrst, numriw, 'syst'  , syst  (:,:)   ) 
    153       CALL iom_rstput( iter, nitrst, numriw, 'sxxst' , sxxst (:,:)   ) 
    154       CALL iom_rstput( iter, nitrst, numriw, 'syyst' , syyst (:,:)   ) 
    155       CALL iom_rstput( iter, nitrst, numriw, 'sxyst' , sxyst (:,:)   ) 
     111      CALL iom_rstput( iter, nitrst, numriw, 'hicif'      , hicif     (:,:)   )      ! prognostic variables  
     112      CALL iom_rstput( iter, nitrst, numriw, 'hsnif'      , hsnif     (:,:)   ) 
     113      CALL iom_rstput( iter, nitrst, numriw, 'frld'       , frld      (:,:)   ) 
     114      CALL iom_rstput( iter, nitrst, numriw, 'sist'       , sist      (:,:)   ) 
     115      CALL iom_rstput( iter, nitrst, numriw, 'tbif1'      , tbif      (:,:,1) ) 
     116      CALL iom_rstput( iter, nitrst, numriw, 'tbif2'      , tbif      (:,:,2) ) 
     117      CALL iom_rstput( iter, nitrst, numriw, 'tbif3'      , tbif      (:,:,3) ) 
     118      CALL iom_rstput( iter, nitrst, numriw, 'u_ice'      , u_ice     (:,:)   ) 
     119      CALL iom_rstput( iter, nitrst, numriw, 'v_ice'      , v_ice     (:,:)   ) 
     120      CALL iom_rstput( iter, nitrst, numriw, 'qstoif'     , qstoif    (:,:)   ) 
     121      CALL iom_rstput( iter, nitrst, numriw, 'fsbbq'      , fsbbq     (:,:)   ) 
     122#if ! defined key_lim2_vp 
     123      CALL iom_rstput( iter, nitrst, numriw, 'stress1_i'  , stress1_i (:,:)   ) 
     124      CALL iom_rstput( iter, nitrst, numriw, 'stress2_i'  , stress2_i (:,:)   ) 
     125      CALL iom_rstput( iter, nitrst, numriw, 'stress12_i' , stress12_i(:,:)   ) 
     126#endif 
     127      CALL iom_rstput( iter, nitrst, numriw, 'sxice' ,      sxice     (:,:)   ) 
     128      CALL iom_rstput( iter, nitrst, numriw, 'syice' ,      syice     (:,:)   ) 
     129      CALL iom_rstput( iter, nitrst, numriw, 'sxxice',      sxxice    (:,:)   ) 
     130      CALL iom_rstput( iter, nitrst, numriw, 'syyice',      syyice    (:,:)   ) 
     131      CALL iom_rstput( iter, nitrst, numriw, 'sxyice',      sxyice    (:,:)   ) 
     132      CALL iom_rstput( iter, nitrst, numriw, 'sxsn'  ,      sxsn      (:,:)   ) 
     133      CALL iom_rstput( iter, nitrst, numriw, 'sysn'  ,      sysn      (:,:)   ) 
     134      CALL iom_rstput( iter, nitrst, numriw, 'sxxsn' ,      sxxsn     (:,:)   ) 
     135      CALL iom_rstput( iter, nitrst, numriw, 'syysn' ,      syysn     (:,:)   ) 
     136      CALL iom_rstput( iter, nitrst, numriw, 'sxysn' ,      sxysn     (:,:)   ) 
     137      CALL iom_rstput( iter, nitrst, numriw, 'sxa'   ,      sxa       (:,:)   ) 
     138      CALL iom_rstput( iter, nitrst, numriw, 'sya'   ,      sya       (:,:)   ) 
     139      CALL iom_rstput( iter, nitrst, numriw, 'sxxa'  ,      sxxa      (:,:)   ) 
     140      CALL iom_rstput( iter, nitrst, numriw, 'syya'  ,      syya      (:,:)   ) 
     141      CALL iom_rstput( iter, nitrst, numriw, 'sxya'  ,      sxya      (:,:)   ) 
     142      CALL iom_rstput( iter, nitrst, numriw, 'sxc0'  ,      sxc0      (:,:)   ) 
     143      CALL iom_rstput( iter, nitrst, numriw, 'syc0'  ,      syc0      (:,:)   ) 
     144      CALL iom_rstput( iter, nitrst, numriw, 'sxxc0' ,      sxxc0     (:,:)   ) 
     145      CALL iom_rstput( iter, nitrst, numriw, 'syyc0' ,      syyc0     (:,:)   ) 
     146      CALL iom_rstput( iter, nitrst, numriw, 'sxyc0' ,      sxyc0     (:,:)   ) 
     147      CALL iom_rstput( iter, nitrst, numriw, 'sxc1'  ,      sxc1      (:,:)   ) 
     148      CALL iom_rstput( iter, nitrst, numriw, 'syc1'  ,      syc1      (:,:)   ) 
     149      CALL iom_rstput( iter, nitrst, numriw, 'sxxc1' ,      sxxc1     (:,:)   ) 
     150      CALL iom_rstput( iter, nitrst, numriw, 'syyc1' ,      syyc1     (:,:)   ) 
     151      CALL iom_rstput( iter, nitrst, numriw, 'sxyc1' ,      sxyc1     (:,:)   ) 
     152      CALL iom_rstput( iter, nitrst, numriw, 'sxc2'  ,      sxc2      (:,:)   ) 
     153      CALL iom_rstput( iter, nitrst, numriw, 'syc2'  ,      syc2      (:,:)   ) 
     154      CALL iom_rstput( iter, nitrst, numriw, 'sxxc2' ,      sxxc2     (:,:)   ) 
     155      CALL iom_rstput( iter, nitrst, numriw, 'syyc2' ,      syyc2     (:,:)   ) 
     156      CALL iom_rstput( iter, nitrst, numriw, 'sxyc2' ,      sxyc2     (:,:)   ) 
     157      CALL iom_rstput( iter, nitrst, numriw, 'sxst'  ,      sxst      (:,:)   ) 
     158      CALL iom_rstput( iter, nitrst, numriw, 'syst'  ,      syst      (:,:)   ) 
     159      CALL iom_rstput( iter, nitrst, numriw, 'sxxst' ,      sxxst     (:,:)   ) 
     160      CALL iom_rstput( iter, nitrst, numriw, 'syyst' ,      syyst     (:,:)   ) 
     161      CALL iom_rstput( iter, nitrst, numriw, 'sxyst' ,      sxyst     (:,:)   ) 
    156162       
    157163      IF( iter == nitrst ) THEN 
    … …  
    218224      ENDIF 
    219225 
    220       CALL iom_get( numrir, jpdom_autoglo, 'qstoif', qstoif )     
    221       CALL iom_get( numrir, jpdom_autoglo, 'fsbbq' , fsbbq  )     
    222       CALL iom_get( numrir, jpdom_autoglo, 'sxice' , sxice  ) 
    223       CALL iom_get( numrir, jpdom_autoglo, 'syice' , syice  ) 
    224       CALL iom_get( numrir, jpdom_autoglo, 'sxxice', sxxice ) 
    225       CALL iom_get( numrir, jpdom_autoglo, 'syyice', syyice ) 
    226       CALL iom_get( numrir, jpdom_autoglo, 'sxyice', sxyice ) 
    227       CALL iom_get( numrir, jpdom_autoglo, 'sxsn'  , sxsn   ) 
    228       CALL iom_get( numrir, jpdom_autoglo, 'sysn'  , sysn   ) 
    229       CALL iom_get( numrir, jpdom_autoglo, 'sxxsn' , sxxsn  ) 
    230       CALL iom_get( numrir, jpdom_autoglo, 'syysn' , syysn  ) 
    231       CALL iom_get( numrir, jpdom_autoglo, 'sxysn' , sxysn  ) 
    232       CALL iom_get( numrir, jpdom_autoglo, 'sxa'   , sxa    ) 
    233       CALL iom_get( numrir, jpdom_autoglo, 'sya'   , sya    ) 
    234       CALL iom_get( numrir, jpdom_autoglo, 'sxxa'  , sxxa   ) 
    235       CALL iom_get( numrir, jpdom_autoglo, 'syya'  , syya   ) 
    236       CALL iom_get( numrir, jpdom_autoglo, 'sxya'  , sxya   ) 
    237       CALL iom_get( numrir, jpdom_autoglo, 'sxc0'  , sxc0   ) 
    238       CALL iom_get( numrir, jpdom_autoglo, 'syc0'  , syc0   ) 
    239       CALL iom_get( numrir, jpdom_autoglo, 'sxxc0' , sxxc0  ) 
    240       CALL iom_get( numrir, jpdom_autoglo, 'syyc0' , syyc0  ) 
    241       CALL iom_get( numrir, jpdom_autoglo, 'sxyc0' , sxyc0  ) 
    242       CALL iom_get( numrir, jpdom_autoglo, 'sxc1'  , sxc1   ) 
    243       CALL iom_get( numrir, jpdom_autoglo, 'syc1'  , syc1   ) 
    244       CALL iom_get( numrir, jpdom_autoglo, 'sxxc1' , sxxc1  ) 
    245       CALL iom_get( numrir, jpdom_autoglo, 'syyc1' , syyc1  ) 
    246       CALL iom_get( numrir, jpdom_autoglo, 'sxyc1' , sxyc1  ) 
    247       CALL iom_get( numrir, jpdom_autoglo, 'sxc2'  , sxc2   ) 
    248       CALL iom_get( numrir, jpdom_autoglo, 'syc2'  , syc2   ) 
    249       CALL iom_get( numrir, jpdom_autoglo, 'sxxc2' , sxxc2  ) 
    250       CALL iom_get( numrir, jpdom_autoglo, 'syyc2' , syyc2  ) 
    251       CALL iom_get( numrir, jpdom_autoglo, 'sxyc2' , sxyc2  ) 
    252       CALL iom_get( numrir, jpdom_autoglo, 'sxst'  , sxst   ) 
    253       CALL iom_get( numrir, jpdom_autoglo, 'syst'  , syst   ) 
    254       CALL iom_get( numrir, jpdom_autoglo, 'sxxst' , sxxst  ) 
    255       CALL iom_get( numrir, jpdom_autoglo, 'syyst' , syyst  ) 
    256       CALL iom_get( numrir, jpdom_autoglo, 'sxyst' , sxyst  ) 
     226      CALL iom_get( numrir, jpdom_autoglo, 'qstoif'     , qstoif     )     
     227      CALL iom_get( numrir, jpdom_autoglo, 'fsbbq'      , fsbbq      )     
     228#if ! defined key_lim2_vp 
     229      CALL iom_get( numrir, jpdom_autoglo, 'stress1_i'  , stress1_i  ) 
     230      CALL iom_get( numrir, jpdom_autoglo, 'stress2_i'  , stress2_i  ) 
     231      CALL iom_get( numrir, jpdom_autoglo, 'stress12_i' , stress12_i ) 
     232#endif 
     233      CALL iom_get( numrir, jpdom_autoglo, 'sxice'      , sxice      ) 
     234      CALL iom_get( numrir, jpdom_autoglo, 'syice'      , syice      ) 
     235      CALL iom_get( numrir, jpdom_autoglo, 'sxxice'     , sxxice     ) 
     236      CALL iom_get( numrir, jpdom_autoglo, 'syyice'     , syyice     ) 
     237      CALL iom_get( numrir, jpdom_autoglo, 'sxyice'     , sxyice     ) 
     238      CALL iom_get( numrir, jpdom_autoglo, 'sxsn'       , sxsn       ) 
     239      CALL iom_get( numrir, jpdom_autoglo, 'sysn'       , sysn       ) 
     240      CALL iom_get( numrir, jpdom_autoglo, 'sxxsn'      , sxxsn      ) 
     241      CALL iom_get( numrir, jpdom_autoglo, 'syysn'      , syysn      ) 
     242      CALL iom_get( numrir, jpdom_autoglo, 'sxysn'      , sxysn      ) 
     243      CALL iom_get( numrir, jpdom_autoglo, 'sxa'        , sxa        ) 
     244      CALL iom_get( numrir, jpdom_autoglo, 'sya'        , sya        ) 
     245      CALL iom_get( numrir, jpdom_autoglo, 'sxxa'       , sxxa       ) 
     246      CALL iom_get( numrir, jpdom_autoglo, 'syya'       , syya       ) 
     247      CALL iom_get( numrir, jpdom_autoglo, 'sxya'       , sxya       ) 
     248      CALL iom_get( numrir, jpdom_autoglo, 'sxc0'       , sxc0       ) 
     249      CALL iom_get( numrir, jpdom_autoglo, 'syc0'       , syc0       ) 
     250      CALL iom_get( numrir, jpdom_autoglo, 'sxxc0'      , sxxc0      ) 
     251      CALL iom_get( numrir, jpdom_autoglo, 'syyc0'      , syyc0      ) 
     252      CALL iom_get( numrir, jpdom_autoglo, 'sxyc0'      , sxyc0      ) 
     253      CALL iom_get( numrir, jpdom_autoglo, 'sxc1'       , sxc1       ) 
     254      CALL iom_get( numrir, jpdom_autoglo, 'syc1'       , syc1       ) 
     255      CALL iom_get( numrir, jpdom_autoglo, 'sxxc1'      , sxxc1      ) 
     256      CALL iom_get( numrir, jpdom_autoglo, 'syyc1'      , syyc1      ) 
     257      CALL iom_get( numrir, jpdom_autoglo, 'sxyc1'      , sxyc1      ) 
     258      CALL iom_get( numrir, jpdom_autoglo, 'sxc2'       , sxc2       ) 
     259      CALL iom_get( numrir, jpdom_autoglo, 'syc2'       , syc2       ) 
     260      CALL iom_get( numrir, jpdom_autoglo, 'sxxc2'      , sxxc2      ) 
     261      CALL iom_get( numrir, jpdom_autoglo, 'syyc2'      , syyc2      ) 
     262      CALL iom_get( numrir, jpdom_autoglo, 'sxyc2'      , sxyc2      ) 
     263      CALL iom_get( numrir, jpdom_autoglo, 'sxst'       , sxst       ) 
     264      CALL iom_get( numrir, jpdom_autoglo, 'syst'       , syst       ) 
     265      CALL iom_get( numrir, jpdom_autoglo, 'sxxst'      , sxxst      ) 
     266      CALL iom_get( numrir, jpdom_autoglo, 'syyst'      , syyst      ) 
     267      CALL iom_get( numrir, jpdom_autoglo, 'sxyst'      , sxyst      ) 
    257268       
    258269      CALL iom_close( numrir ) 

  • branches/DEV_r2191_3partymerge2010/NEMO/LIM_SRC_2/limsbc_2.F90

    r2208 r2222  
    44   !!           computation of the flux at the sea ice/ocean interface 
    55   !!====================================================================== 
    6    !! History : 00-01 (H. Goosse) Original code 
    7    !!           02-07 (C. Ethe, G. Madec) re-writing F90 
    8    !!           06-07 (G. Madec) surface module 
     6   !! History :  LIM  ! 2000-01 (H. Goosse) Original code 
     7   !!            1.0  ! 2002-07 (C. Ethe, G. Madec) re-writing F90 
     8   !!            3.0  ! 2006-07 (G. Madec) surface module 
     9   !!            3.3  ! 2009-05 (G.Garric, C. Bricaud) addition of the lim2_evp case 
    910   !!---------------------------------------------------------------------- 
    1011#if defined key_lim2 
    1112   !!---------------------------------------------------------------------- 
    1213   !!   'key_lim2'                                    LIM 2.0 sea-ice model 
    13    !!---------------------------------------------------------------------- 
    1414   !!---------------------------------------------------------------------- 
    1515   !!   lim_sbc_2  : flux at the ice / ocean interface 
    … …  
    1717   USE par_oce          ! ocean parameters 
    1818   USE dom_oce          ! ocean domain 
    19    USE sbc_ice          ! surface boundary condition 
    20    USE sbc_oce          ! surface boundary condition 
     19   USE sbc_ice          ! surface boundary condition: ice 
     20   USE sbc_oce          ! surface boundary condition: ocean 
    2121   USE phycst           ! physical constants 
    22    USE ice_2            ! LIM sea-ice variables 
    23  
    24    USE lbclnk           ! ocean lateral boundary condition 
     22   USE ice_2            ! LIM2: ice variables 
     23 
     24   USE lbclnk           ! ocean lateral boundary condition - MPP exchanges 
    2525   USE in_out_manager   ! I/O manager 
    2626   USE diaar5, ONLY :   lk_diaar5 
    27    USE iom              !  
     27   USE iom              ! IOM library 
    2828   USE albedo           ! albedo parameters 
    2929   USE prtctl           ! Print control 
    … …  
    3333   PRIVATE 
    3434 
    35    PUBLIC lim_sbc_2     ! called by sbc_ice_lim_2 
    36  
    37    REAL(wp)  ::   epsi16 = 1.e-16  ! constant values 
    38    REAL(wp)  ::   rzero  = 0.e0     
    39    REAL(wp)  ::   rone   = 1.e0 
    40    REAL(wp), DIMENSION(jpi,jpj)  ::   soce_r 
    41    REAL(wp), DIMENSION(jpi,jpj)  ::   sice_r 
     35   PUBLIC   lim_sbc_2   ! called by sbc_ice_lim_2 
     36 
     37   REAL(wp)  ::   r1_rdtice                    ! constant values 
     38   REAL(wp)  ::   epsi16 = 1.e-16              !     -      - 
     39   REAL(wp)  ::   rzero  = 0.e0                !     -      - 
     40   REAL(wp)  ::   rone   = 1.e0                !     -      - 
     41   ! 
     42   REAL(wp), DIMENSION(jpi,jpj) ::   soce_r, sice_r   ! constant SSS and ice salinity used in levitating sea-ice case 
    4243 
    4344   !! * Substitutions 
    4445#  include "vectopt_loop_substitute.h90" 
    4546   !!---------------------------------------------------------------------- 
    46    !!   LIM 2.0,  UCL-LOCEAN-IPSL (2006)  
     47   !! NEMO/LIM2 3.3, UCL-LOCEAN-IPSL (2010) 
    4748   !! $Id$ 
    4849   !! Software governed by the CeCILL licence  (NEMOGCM/License_CeCILL.txt) 
    4950   !!---------------------------------------------------------------------- 
    50  
    5151CONTAINS 
    5252 
    … …  
    9898      !!--------------------------------------------------------------------- 
    9999      
    100       zrdtir = 1. / rdt_ice 
    101100       
    102101      IF( kt == nit000 ) THEN 
    … …  
    104103         IF(lwp) WRITE(numout,*) 'lim_sbc_2 : LIM 2.0 sea-ice - surface boundary condition' 
    105104         IF(lwp) WRITE(numout,*) '~~~~~~~~~   ' 
    106  
     105         ! 
     106         r1_rdtice = 1. / rdt_ice 
     107         ! 
    107108         soce_r(:,:) = soce 
    108109         sice_r(:,:) = sice 
    … …  
    177178!!$ 
    178179 
    179             !   computation the solar flux at ocean surface 
     180            ! solar flux at ocean surface 
    180181#if defined key_coupled  
    181182            zqsr = qsr_tot(ji,jj) + ( fstric(ji,jj) - qsr_ice(ji,jj,1) ) * ( 1.0 - pfrld(ji,jj) ) 
    … …  
    183184            zqsr = pfrld(ji,jj) * qsr(ji,jj)  + ( 1.  - pfrld(ji,jj) ) * fstric(ji,jj) 
    184185#endif             
    185             !  computation the non solar heat flux at ocean surface 
     186            ! non solar heat flux at ocean surface 
    186187            zqns    =  - ( 1. - thcm(ji,jj) ) * zqsr   &   ! part of the solar energy used in leads 
    187                &       + iflt    * ( fscmbq(ji,jj) + ffltbif(ji,jj) )                            & 
    188                &       + ifral   * ( ial * qcmif(ji,jj) + (1 - ial) * qldif(ji,jj) ) * zrdtir    & 
    189                &       + ifrdv   * ( qfvbq(ji,jj) + qdtcn(ji,jj) ) * zrdtir 
    190  
     188               &       + iflt    * ( fscmbq(ji,jj) + ffltbif(ji,jj) )                              & 
     189               &       + ifral   * ( ial * qcmif(ji,jj) + (1 - ial) * qldif(ji,jj) ) * r1_rdtice   & 
     190               &       + ifrdv   * ( qfvbq(ji,jj) + qdtcn(ji,jj) )                   * r1_rdtice 
     191            ! 
    191192            fsbbq(ji,jj) = ( 1.0 - ( ifvt + iflt ) ) * fscmbq(ji,jj)     ! ??? 
    192              
     193            ! 
    193194            qsr  (ji,jj) = zqsr                                          ! solar heat flux  
    194195            qns  (ji,jj) = zqns - fdtcn(ji,jj)                           ! non solar heat flux 
    … …  
    196197      END DO 
    197198 
    198       CALL iom_put( 'hflx_ice_cea', - fdtcn(:,:) )       
    199       CALL iom_put( 'qns_io_cea', qns(:,:) - zqnsoce(:,:) * pfrld(:,:) )       
    200       CALL iom_put( 'qsr_io_cea', fstric(:,:) * (1. - pfrld(:,:)) ) 
     199      CALL iom_put( 'hflx_ice_cea', - fdtcn(:,:)                         )       
     200      CALL iom_put( 'qns_io_cea'  , qns(:,:) - zqnsoce(:,:) * pfrld(:,:) )       
     201      CALL iom_put( 'qsr_io_cea'  , fstric(:,:) * ( 1.e0 - pfrld(:,:) ) ) 
    201202 
    202203      !------------------------------------------! 
    203204      !      mass flux at the ocean surface      ! 
    204205      !------------------------------------------! 
    205  
    206 !!gm 
    207 !!gm CAUTION    
    208 !!gm re-verifies the emp & emps expression, especially the absence of 1-frld on zfm 
    209 !!gm 
    210206      DO jj = 1, jpj 
    211207         DO ji = 1, jpi 
    212              
    213208#if defined key_coupled 
    214           zemp = emp_tot(ji,jj) - emp_ice(ji,jj) * ( 1. - pfrld(ji,jj) )    &   !  
    215              &   + rdmsnif(ji,jj) * zrdtir                                      !  freshwaterflux due to snow melting  
     209            ! freshwater exchanges at the ice-atmosphere / ocean interface (coupled mode) 
     210            zemp = emp_tot(ji,jj) - emp_ice(ji,jj) * ( 1. - pfrld(ji,jj) )   &   ! atmosphere-ocean freshwater flux 
     211               &                  + rdmsnif(ji,jj) * r1_rdtice                   ! freshwater flux due to snow melting  
    216212#else 
    217 !!$            !  computing freshwater exchanges at the ice/ocean interface 
    218 !!$            zpme = - evap(ji,jj)    *   frld(ji,jj)           &   !  evaporation over oceanic fraction 
    219 !!$               &   + tprecip(ji,jj)                           &   !  total precipitation 
    220 !!$               &   - sprecip(ji,jj) * ( 1. - pfrld(ji,jj) )   &   !  remov. snow precip over ice 
    221 !!$               &   - rdmsnif(ji,jj) / rdt_ice                     !  freshwaterflux due to snow melting  
    222             !  computing freshwater exchanges at the ice/ocean interface 
    223             zemp = + emp(ji,jj)     *         frld(ji,jj)      &   !  e-p budget over open ocean fraction  
    224                &   - tprecip(ji,jj) * ( 1. -  frld(ji,jj) )    &   !  liquid precipitation reaches directly the ocean 
    225                &   + sprecip(ji,jj) * ( 1. - pfrld(ji,jj) )    &   !  taking into account change in ice cover within the time step 
    226                &   + rdmsnif(ji,jj) * zrdtir                       !  freshwaterflux due to snow melting  
    227                !                                                   !  ice-covered fraction: 
     213            ! freshwater exchanges at the ice-atmosphere / ocean interface (forced mode) 
     214            zemp = + emp(ji,jj)     *         frld(ji,jj)     &   ! e-p budget over open ocean fraction  
     215               &   - tprecip(ji,jj) * ( 1. -  frld(ji,jj) )   &   ! liquid precipitation reaches directly the ocean 
     216               &   + sprecip(ji,jj) * ( 1. - pfrld(ji,jj) )   &   ! (account for change in ice cover within the timestep 
     217               &   + rdmsnif(ji,jj) * r1_rdtice                   ! freshwaterflux due to snow melting  
    228218#endif             
    229  
    230             !  computing salt exchanges at the ice/ocean interface 
    231             zfons =  ( soce_r(ji,jj) - sice_r(ji,jj) ) * ( rdmicif(ji,jj) * zrdtir )  
    232              
    233             !  converting the salt flux from ice to a freshwater flux from ocean 
     219            ! salt exchanges at the ice/ocean interface 
     220            zfons =  ( soce_r(ji,jj) - sice_r(ji,jj) ) * ( rdmicif(ji,jj) * r1_rdtice )  
     221            ! 
     222            ! convert the salt flux from ice into a freshwater flux from ocean 
    234223            zfm  = zfons / ( sss_m(ji,jj) + epsi16 ) 
    235              
     224            ! 
    236225            emps(ji,jj) = zemp + zfm      ! surface ocean concentration/dilution effect (use on SSS evolution) 
    237226            emp (ji,jj) = zemp            ! surface ocean volume flux (use on sea-surface height evolution) 
    238  
    239227         END DO 
    240228      END DO 
    241  
     229      ! 
    242230      IF( lk_diaar5 ) THEN 
    243          CALL iom_put( 'isnwmlt_cea'  ,                 rdmsnif(:,:) * zrdtir ) 
    244          CALL iom_put( 'fsal_virt_cea',   soce_r(:,:) * rdmicif(:,:) * zrdtir ) 
    245          CALL iom_put( 'fsal_real_cea', - sice_r(:,:) * rdmicif(:,:) * zrdtir ) 
     231         CALL iom_put( 'isnwmlt_cea'  ,                 rdmsnif(:,:) * r1_rdtice ) 
     232         CALL iom_put( 'fsal_virt_cea',   soce_r(:,:) * rdmicif(:,:) * r1_rdtice ) 
     233         CALL iom_put( 'fsal_real_cea', - sice_r(:,:) * rdmicif(:,:) * r1_rdtice ) 
    246234      ENDIF 
    247235 
    … …  
    277265            DO ji = 2, jpim1   ! NO vector opt. 
    278266               ! ... components of ice-ocean stress at U and V-points  (from I-point values) 
    279                zutau  = 0.5 * ( ztio_u(ji+1,jj) + ztio_u(ji+1,jj+1) ) 
     267#if defined key_lim2_vp 
     268               zutau  = 0.5 * ( ztio_u(ji+1,jj) + ztio_u(ji+1,jj+1) )      ! VP rheology 
    280269               zvtau  = 0.5 * ( ztio_v(ji,jj+1) + ztio_v(ji+1,jj+1) ) 
     270#else 
     271               zutau  = ztio_u(ji,jj)                                      ! EVP rheology 
     272               zvtau  = ztio_v(ji,jj) 
     273#endif 
    281274               ! ... open-ocean (lead) fraction at U- & V-points (from T-point values) 
    282275               zfrldu = 0.5 * ( frld(ji,jj) + frld(ji+1,jj  ) ) 
    … …  
    292285            END DO 
    293286         END DO 
    294  
    295          ! boundary condition on the stress (utau,vtau,taum) 
    296          CALL lbc_lnk( utau, 'U', -1. ) 
    297          CALL lbc_lnk( vtau, 'V', -1. ) 
     287         CALL lbc_lnk( utau, 'U', -1. )   ;   CALL lbc_lnk( vtau, 'V', -1. )     ! lateral boundary condition  
    298288         CALL lbc_lnk( taum, 'T',  1. ) 
    299289 
    300290      ENDIF 
    301291 
     292      IF( lk_cpl ) THEN            
    302293      !-----------------------------------------------! 
    303294      !   Coupling variables                          ! 
    304295      !-----------------------------------------------! 
    305  
    306       IF ( lk_cpl ) THEN            
    307          ! Ice surface temperature  
    308          tn_ice(:,:,1) = sist(:,:)          ! sea-ice surface temperature        
    309          ! Computation of snow/ice and ocean albedo 
     296         tn_ice(:,:,1) = sist(:,:)           ! sea-ice surface temperature        
     297         !                                   ! snow/ice and ocean albedo 
    310298         CALL albedo_ice( tn_ice, reshape( hicif, (/jpi,jpj,1/) ), reshape( hsnif, (/jpi,jpj,1/) ), zalbp, zalb ) 
    311299         alb_ice(:,:,1) =  0.5 * ( zalbp(:,:,1) + zalb (:,:,1) )   ! Ice albedo (mean clear and overcast skys) 
    … …  
    320308         CALL prt_ctl(tab2d_1=fr_i  , clinfo1=' lim_sbc: fr_i   : ', tab2d_2=tn_ice(:,:,1), clinfo2=' tn_ice  : ') 
    321309      ENDIF  
    322     
    323     END SUBROUTINE lim_sbc_2 
     310      ! 
     311   END SUBROUTINE lim_sbc_2 
    324312 
    325313#else 

  • branches/DEV_r2191_3partymerge2010/NEMO/LIM_SRC_2/limtrp_2.F90

    r2208 r2222  
    44   !! LIM 2.0 transport ice model : sea-ice advection/diffusion 
    55   !!====================================================================== 
    6    !! History :  LIM  !  2000-01 (UCL)  Original code 
    7    !!            2.0  !  2001-05 (G. Madec, R. Hordoir) opa norm 
    8    !!             -   !  2004-01 (G. Madec, C. Ethe)  F90, mpp 
     6   !! History :  LIM  !  2000-01  (LIM)  Original code 
     7   !!            2.0  !  2001-05  (G. Madec, R. Hordoir) opa norm 
     8   !!             -   !  2004-01  (G. Madec, C. Ethe)  F90, mpp 
     9   !!            3.3  !  2009-05  (G.Garric, C. Bricaud) addition of the lim2_evp case 
    910   !!---------------------------------------------------------------------- 
    1011#if defined key_lim2 
    … …  
    3435   REAL(wp), PUBLIC  ::   bound  = 0.e0   !: boundary condit. (0.0 no-slip, 1.0 free-slip) 
    3536 
    36    REAL(wp)  ::           &  ! constant values 
    37       epsi06 = 1.e-06  ,  & 
    38       epsi03 = 1.e-03  ,  & 
    39       epsi16 = 1.e-16  ,  & 
    40       rzero  = 0.e0    ,  & 
    41       rone   = 1.e0 
     37   REAL(wp)  ::   epsi06 = 1.e-06   ! constant values 
     38   REAL(wp)  ::   epsi03 = 1.e-03   
     39   REAL(wp)  ::   epsi16 = 1.e-16   
     40   REAL(wp)  ::   rzero  = 0.e0    
     41   REAL(wp)  ::   rone   = 1.e0 
    4242 
    4343   !! * Substitution 
    4444#  include "vectopt_loop_substitute.h90" 
    4545   !!---------------------------------------------------------------------- 
    46    !! NEMO/LIM 3.2,  UCL-LOCEAN-IPSL (2010)  
     46   !! NEMO/LIM2 3.3,  UCL-LOCEAN-IPSL (2010)  
    4747   !! $Id$ 
    4848   !! Software governed by the CeCILL licence  (NEMOGCM/License_CeCILL.txt) 
    4949   !!---------------------------------------------------------------------- 
    50  
    5150CONTAINS 
    5251 
    … …  
    6766      INTEGER  ::   ji, jj, jk   ! dummy loop indices 
    6867      INTEGER  ::   initad       ! number of sub-timestep for the advection 
    69       REAL(wp) ::   zindb  , zindsn , zindic, zacrith   ! local scalars 
    70       REAL(wp) ::   zusvosn, zusvoic, zignm , zindhe    !   -      - 
    71       REAL(wp) ::   zvbord , zcfl   , zusnit            !   -      - 
    72       REAL(wp) ::   zrtt   , ztsn   , ztic1 , ztic2     !   -      - 
    73       REAL(wp), DIMENSION(jpi,jpj)  ::   zui_u , zvi_v , zsm             ! 2D workspace 
    74       REAL(wp), DIMENSION(jpi,jpj)  ::   zs0ice, zs0sn , zs0a            !  -      - 
    75       REAL(wp), DIMENSION(jpi,jpj)  ::   zs0c0 , zs0c1 , zs0c2 , zs0st   !  -      - 
     68 
     69      REAL(wp) ::   zindb  , zacrith, zindsn , zindic , zusvosn   ! local scalars 
     70      REAL(wp) ::   zusvoic, zignm  , zindhe , zvbord , zcfl      !   -      - 
     71      REAL(wp) ::   zusnit , zrtt   , ztsn   , ztic1  , ztic2     !   -      - 
     72 
     73      REAL(wp), DIMENSION(jpi,jpj) ::   zui_u , zvi_v , zsm             ! 2D workspace 
     74      REAL(wp), DIMENSION(jpi,jpj) ::   zs0ice, zs0sn , zs0a            !  -      - 
     75      REAL(wp), DIMENSION(jpi,jpj) ::   zs0c0 , zs0c1 , zs0c2 , zs0st   !  -      - 
    7676      !--------------------------------------------------------------------- 
    7777 
    … …  
    8888         ! --------------------------------------- 
    8989         zvbord = 1.0 + ( 1.0 - bound )      ! zvbord=2 no-slip, =0 free slip boundary conditions         
    90          DO jj = 1, jpjm1 
     90#if defined key_lim2_vp 
     91         DO jj = 1, jpjm1     ! VP rheology: ice (u,v) at I-point 
    9192            DO ji = 1, jpim1   ! NO vector opt. 
    92                zui_u(ji,jj) = ( u_ice(ji+1,jj  ) + u_ice(ji+1,jj+1) ) / ( MAX( tmu(ji+1,jj  ) + tmu(ji+1,jj+1), zvbord ) ) 
    93                zvi_v(ji,jj) = ( v_ice(ji  ,jj+1) + v_ice(ji+1,jj+1) ) / ( MAX( tmu(ji  ,jj+1) + tmu(ji+1,jj+1), zvbord ) ) 
     93               zui_u(ji,jj) = ( u_ice(ji+1,jj) + u_ice(ji+1,jj+1) ) / ( MAX( tmu(ji+1,jj  ) + tmu(ji+1,jj+1), zvbord ) ) 
     94               zvi_v(ji,jj) = ( v_ice(ji,jj+1) + v_ice(ji+1,jj+1) ) / ( MAX( tmu(ji  ,jj+1) + tmu(ji+1,jj+1), zvbord ) ) 
    9495            END DO 
    9596         END DO 
    96          CALL lbc_lnk( zui_u, 'U', -1. )   ;   CALL lbc_lnk( zvi_v, 'V', -1. )         ! Lateral boundary conditions 
     97         CALL lbc_lnk( zui_u, 'U', -1. )   ;   CALL lbc_lnk( zvi_v, 'V', -1. )      ! Lateral boundary conditions  
     98#else 
     99        zui_u(:,:) = u_ice(:,:)      ! EVP rheology: ice (u,v) at u- and v-points 
     100        zvi_v(:,:) = v_ice(:,:) 
     101#endif 
    97102 
    98103 
    … …  
    107112         IF( zcfl > 0.5 .AND. lwp )   WRITE(numout,*) 'lim_trp_2 : violation of cfl criterion the ',nday,'th day, cfl = ', zcfl 
    108113 
     114         IF(lk_mpp ) CALL mpp_max(zcfl) 
     115 
     116         IF( zcfl > 0.5 .AND. lwp )   WRITE(numout,*) 'lim_trp_2 : cfl criterion violation. day ',nday,' cfl = ',zcfl 
     117 
    109118         ! content of properties 
    110119         ! --------------------- 
    111          zs0sn (:,:) =  hsnm(:,:) * area(:,:)                 ! Snow volume. 
    112          zs0ice(:,:) =  hicm(:,:) * area(:,:)                 ! Ice volume. 
    113          zs0a  (:,:) =  ( 1.0 - frld(:,:) )    * area  (:,:)  ! Surface covered by ice. 
    114          zs0c0 (:,:) =  tbif(:,:,1) / rt0_snow * zs0sn (:,:)  ! Heat content of the snow layer. 
    115          zs0c1 (:,:) =  tbif(:,:,2) / rt0_ice  * zs0ice(:,:)  ! Heat content of the first ice layer. 
    116          zs0c2 (:,:) =  tbif(:,:,3) / rt0_ice  * zs0ice(:,:)  ! Heat content of the second ice layer. 
    117          zs0st (:,:) =  qstoif(:,:) / xlic     * zs0a  (:,:)  ! Heat reservoir for brine pockets. 
    118           
     120         zs0sn (:,:) =  hsnm(:,:)              * area(:,:)     ! Snow volume. 
     121         zs0ice(:,:) =  hicm (:,:)             * area(:,:)     ! Ice volume. 
     122         zs0a  (:,:) =  ( 1.0 - frld(:,:) )    * area(:,:)     ! Surface covered by ice. 
     123         zs0c0 (:,:) =  tbif(:,:,1) / rt0_snow * zs0sn(:,:)    ! Heat content of the snow layer. 
     124         zs0c1 (:,:) =  tbif(:,:,2) / rt0_ice  * zs0ice(:,:)   ! Heat content of the first ice layer. 
     125         zs0c2 (:,:) =  tbif(:,:,3) / rt0_ice  * zs0ice(:,:)   ! Heat content of the second ice layer. 
     126         zs0st (:,:) =  qstoif(:,:) / xlic     * zs0a  (:,:)   ! Heat reservoir for brine pockets. 
    119127  
    120128         ! Advection (Prather scheme) 

  • branches/DEV_r2191_3partymerge2010/NEMO/LIM_SRC_3/limrhg.F90

    r2208 r2222  
    44   !!   Ice rheology : sea ice rheology 
    55   !!====================================================================== 
    6    !! History :   -   !  2007-03  (M.A. Morales Maqueda, S. Bouillon) Original code 
    7    !!            3.0  !  2008-03  (M. Vancoppenolle) LIM3 
     6   !! History :  LIM  !  2007-03  (M.A. Morales Maqueda, S. Bouillon) Original code 
     7   !!            3.0  !  2008-03  (M. Vancoppenolle) LIM 3 
    88   !!             -   !  2008-11  (M. Vancoppenolle, S. Bouillon, Y. Aksenov) add surface tilt in ice rheolohy  
     9   !!            3.3  !  2009-05  (G.Garric) addition of the lim2_evp cas 
    910   !!---------------------------------------------------------------------- 
    10 #if defined key_lim3 
     11#if defined key_lim3 || (  defined key_lim2 && ! defined key_lim2_vp ) 
    1112   !!---------------------------------------------------------------------- 
    1213   !!   'key_lim3'                                      LIM3 sea-ice model 
    … …  
    1415   !!   lim_rhg   : computes ice velocities 
    1516   !!---------------------------------------------------------------------- 
    16    !! * Modules used 
    17    USE phycst 
    18    USE par_oce 
    19    USE dom_oce 
    20    USE dom_ice 
    21    USE sbc_oce         ! Surface boundary condition: ocean fields 
    22    USE sbc_ice         ! Surface boundary condition: ice fields 
    23    USE ice 
    24    USE iceini 
    25    USE lbclnk 
    26    USE lib_mpp 
    27    USE in_out_manager  ! I/O manager 
    28    USE limitd_me 
    29    USE prtctl          ! Print control 
    30  
     17   USE phycst           ! physical constants 
     18   USE par_oce          ! ocean parameters 
     19   USE dom_oce          ! ocean domain 
     20   USE sbc_oce          ! Surface boundary condition: ocean fields 
     21   USE sbc_ice          ! Surface boundary condition: ice fields 
     22   USE lbclnk           ! lateral boundary condition - MPP exchanges 
     23   USE lib_mpp          ! MPP library 
     24   USE in_out_manager   ! I/O manager 
     25   USE limitd_me        ! LIM3:  
     26   USE prtctl           ! control print 
     27#if defined key_lim3 
     28   USE ice              ! LIM3: ice variables 
     29   USE dom_ice          ! LIM3: ice domain 
     30   USE iceini           ! LIM3: ice initialisation 
     31#endif 
     32#if defined key_lim2 && ! defined key_lim2_vp 
     33   USE ice_2            ! LIM2: ice variables 
     34   USE dom_ice_2        ! LIM2: ice domain 
     35   USE iceini_2         ! LIM2: ice initialisation 
     36#endif 
    3137 
    3238   IMPLICIT NONE 
    3339   PRIVATE 
    3440 
    35    !! * Routine accessibility 
    36    PUBLIC lim_rhg  ! routine called by lim_dyn 
    37  
    38    !! * Substitutions 
    39 #  include "vectopt_loop_substitute.h90" 
     41   PUBLIC   lim_rhg   ! routine called by lim_dyn module 
    4042 
    4143   !! * Module variables 
    … …  
    4345      rzero   = 0.e0   ,  & 
    4446      rone    = 1.e0 
     47 
     48   !! * Substitutions 
     49#  include "vectopt_loop_substitute.h90" 
    4550   !!---------------------------------------------------------------------- 
    46    !!   LIM 3.0,  UCL-LOCEAN-IPSL (2008)  
     51   !! NEMO/LIM-3 3.3,  UCL-LOCEAN-IPSL (2010)  
    4752   !! $Id$ 
    4853   !! Software governed by the CeCILL licence  (NEMOGCM/License_CeCILL.txt) 
    4954   !!---------------------------------------------------------------------- 
    50  
    5155CONTAINS 
    5256 
    5357   SUBROUTINE lim_rhg( k_j1, k_jpj ) 
    54  
    5558      !!------------------------------------------------------------------- 
    5659      !!                 ***  SUBROUTINE lim_rhg  *** 
    … …  
    100103      !! 
    101104      !! ** References : Hunke and Dukowicz, JPO97 
    102       !!                 Bouillon et al., 08, in prep (update this when 
    103       !!                 published) 
    104       !!                 Vancoppenolle et al., OM08 
     105      !!                 Bouillon et al., 2009, Ocean. Modelling, 27, 174-184. 
     106      !!                 Vancoppenolle et al. 2009, Ocean Modelling, 27, 33-53. 
     107      !!------------------------------------------------------------------- 
     108      INTEGER, INTENT(in) ::   k_j1    ! southern j-index for ice computation 
     109      INTEGER, INTENT(in) ::   k_jpj   ! northern j-index for ice computation 
    105110      !! 
    106       !!------------------------------------------------------------------- 
    107       ! * Arguments 
    108       ! 
    109       INTEGER, INTENT(in) :: & 
    110          k_j1 ,                      & !: southern j-index for ice computation 
    111          k_jpj                         !: northern j-index for ice computation 
    112  
    113       ! * Local variables 
    114       INTEGER ::   ji, jj              !: dummy loop indices 
    115  
    116       INTEGER  :: & 
    117          jter                          !: temporary integers 
    118  
    119       CHARACTER (len=50) ::   charout 
    120  
    121       REAL(wp) :: & 
    122          zt11, zt12, zt21, zt22,     & !: temporary scalars 
    123          ztagnx, ztagny,             & !: wind stress on U/V points                        
    124          delta                         ! 
    125  
    126       REAL(wp) :: & 
    127          za,                         & !: 
    128          zstms,                      & !: temporary scalar for ice strength 
    129          zsang,                      & !: temporary scalar for coriolis term 
    130          zmask                         !: mask for the computation of ice mass 
    131  
    132       REAL(wp),DIMENSION(jpi,jpj) :: & 
    133          zpresh        ,             & !: temporary array for ice strength 
    134          zpreshc       ,             & !: Ice strength on grid cell corners (zpreshc) 
    135          zfrld1, zfrld2,             & !: lead fraction on U/V points                                     
    136          zmass1, zmass2,             & !: ice/snow mass on U/V points                                     
    137          zcorl1, zcorl2,             & !: coriolis parameter on U/V points 
    138          za1ct, za2ct  ,             & !: temporary arrays 
    139          zc1           ,             & !: ice mass 
    140          zusw          ,             & !: temporary weight for the computation 
    141                                 !: of ice strength 
    142          u_oce1, v_oce1,             & !: ocean u/v component on U points                            
    143          u_oce2, v_oce2,             & !: ocean u/v component on V points 
    144          u_ice2,                     & !: ice u component on V point 
    145          v_ice1                        !: ice v component on U point 
     111      INTEGER  ::   ji, jj   ! dummy loop indices 
     112      INTEGER  ::   jter     ! local integers 
     113      CHARACTER (len=50) ::   charout   ! local character 
     114      REAL(wp) ::   zt11, zt12, zt21, zt22   ! local scalars 
     115      REAL(wp) ::   ztagnx, ztagny, delta    !   -      - 
     116      REAL(wp) ::   za, zstms, zsang, zmask  !   -      - 
     117      REAL(wp) ::   zresm, zindb, zdummy     !   -    - 
     118 
     119      REAL(wp),DIMENSION(jpi,jpj) ::   zpresh , zfrld1, zmass1, zcorl1, za1ct    ! 2D workspace 
     120      REAL(wp),DIMENSION(jpi,jpj) ::   zpreshc, zfrld2, zmass2, zcorl2, za2ct    !  -      - 
     121      REAL(wp),DIMENSION(jpi,jpj) ::   u_oce1, v_oce1, u_ice2, zc1               !  -      -                           
     122      REAL(wp),DIMENSION(jpi,jpj) ::   u_oce2, v_oce2, v_ice1, zusw              !  -      - 
     123      REAL(wp),DIMENSION(jpi,jpj) ::   zf1, zf2                                  !  -      - 
    146124 
    147125      REAL(wp) :: & 
    … …  
    159137         sigma1, sigma2                !: internal ice stress 
    160138 
    161       REAL(wp),DIMENSION(jpi,jpj) :: & 
    162          zf1, zf2                      !: arrays for internal stresses 
    163139 
    164140      REAL(wp),DIMENSION(jpi,jpj) :: & 
    … …  
    170146         zs12                          !: Non-diagonal stress tensor component zs12 
    171147 
    172       REAL(wp) :: & 
    173          zresm            ,          & !: Maximal error on ice velocity 
    174          zindb            ,          & !: ice (1) or not (0)       
    175          zdummy                        !: dummy argument 
    176  
    177       REAL(wp),DIMENSION(jpi,jpj) :: & 
    178          zu_ice           ,          & !: Ice velocity on previous time step 
    179          zv_ice           ,          & 
    180          zresr                         !: Local error on velocity 
    181  
     148 
     149      REAL(wp),DIMENSION(jpi,jpj) ::   zu_ice, zv_ice, zresr   ! 
     150      !!------------------------------------------------------------------- 
     151 
     152#if  defined key_lim2 && ! defined key_lim2_vp 
     153     vt_s => hsnm 
     154     vt_i => hicm 
     155     at_i(:,:) = 1. - frld(:,:) 
     156#endif 
    182157      ! 
    183158      !------------------------------------------------------------------------------! 
    … …  
    190165      u_ice2(:,:)  = 0.0 ; v_ice1(:,:)  = 0.0 
    191166      zdd(:,:)     = 0.0 ; zdt(:,:)     = 0.0 ; zds(:,:)     = 0.0 
    192  
     167#if defined key_lim3 
    193168      ! Ice strength on T-points 
    194169      CALL lim_itd_me_icestrength(ridge_scheme_swi) 
     170#endif 
    195171 
    196172      ! Ice mass and temp variables 
    … …  
    200176         DO ji = 1 , jpi 
    201177            zc1(ji,jj)    = tms(ji,jj) * ( rhosn * vt_s(ji,jj) + rhoic * vt_i(ji,jj) ) 
     178#if defined key_lim3 
    202179            zpresh(ji,jj) = tms(ji,jj) *  strength(ji,jj) / 2. 
     180#else 
     181            zpresh(ji,jj) = tms(ji,jj) *  2. * pstar * hicm(ji,jj) * EXP( -c_rhg * frld(ji,jj) ) 
     182#endif 
    203183            ! tmi = 1 where there is ice or on land 
    204184            tmi(ji,jj)    = 1.0 - ( 1.0 - MAX( 0.0 , SIGN ( 1.0 , vt_i(ji,jj) - & 
    … …  
    269249               / ( e2t(ji,jj+1) + e2t(ji,jj) + epsd ) 
    270250            ! 
     251            ! Mass, coriolis coeff. and currents 
    271252            u_oce1(ji,jj)  = u_oce(ji,jj) 
    272253            v_oce2(ji,jj)  = v_oce(ji,jj) 

  • branches/DEV_r2191_3partymerge2010/NEMO/OPA_SRC/DIA/diawri.F90

    r2207 r2222  
    3030   USE limwri_2  
    3131#endif 
     32   USE dtatem 
     33   USE dtasal 
     34 
    3235   IMPLICIT NONE 
    3336   PRIVATE 

  • branches/DEV_r2191_3partymerge2010/NEMO/OPA_SRC/DOM/daymod.F90

    r2208 r2222  
    6767      !!              - nmonth_len, nyear_len, nmonth_half, nmonth_end through day_mth 
    6868      !!---------------------------------------------------------------------- 
     69      INTEGER :: inbday, irest 
     70      REAL(wp) :: zjul 
     71      !!---------------------------------------------------------------------- 
    6972 
    7073      ! all calendar staff is based on the fact that MOD( rday, rdttra(1) ) == 0 
    … …  
    105108      ! day since january 1st 
    106109      nday_year = nday + SUM( nmonth_len(1:nmonth - 1) ) 
    107        
     110 
     111      !compute number of days between last monday and today       
     112      IF( nn_leapy==1 )THEN 
     113         CALL ymds2ju( 1900, 01, 01, 0.0, zjul )  ! compute julian day value of 01.01.1900 (monday) 
     114         inbday = INT(fjulday) - NINT(zjul)       ! compute nb day between  01.01.1900 and current day fjulday  
     115         irest = MOD(inbday,7)                    ! compute nb day between last monday and current day fjulday  
     116         IF(irest==0 )irest = 7  
     117      ENDIF 
     118 
    108119      ! number of seconds since the beginning of current year/month at the middle of the time-step 
    109120      nsec_year  = nday_year * nsecd - ndt05   ! 1 time step before the middle of the first time step 
    110121      nsec_month = nday      * nsecd - ndt05   ! because day will be called at the beginning of step 
    111122      nsec_day   =             nsecd - ndt05 
     123      nsec_week  = 0 
     124      IF( nn_leapy==1 ) nsec_week  = irest     * nsecd - ndt05 
    112125 
    113126      ! control print 
    114127      IF(lwp) WRITE(numout,'(a,i6,a,i2,a,i2,a,i6)')' ==============>> 1/2 time step before the start of the run DATE Y/M/D = ',   & 
    115            &                   nyear, '/', nmonth, '/', nday, '  nsec_day:', nsec_day 
     128           &                   nyear, '/', nmonth, '/', nday, '  nsec_day:', nsec_day, '  nsec_week:', nsec_week 
    116129 
    117130      ! Up to now, calendar parameters are related to the end of previous run (nit000-1) 
    … …  
    200213      nsec_year  = nsec_year  + ndt  
    201214      nsec_month = nsec_month + ndt                  
     215      IF( nn_leapy==1 ) nsec_week  = nsec_week  + ndt 
    202216      nsec_day   = nsec_day   + ndt                 
    203217      adatrj  = adatrj  + rdttra(1) / rday 
    … …  
    228242         ndastp = nyear * 10000 + nmonth * 100 + nday   ! NEW date 
    229243         ! 
     244         !compute first day of the year in julian days 
     245         CALL ymds2ju( nyear, 01, 01, 0.0, fjulstartyear ) 
     246         ! 
    230247         IF(lwp) WRITE(numout,'(a,i8,a,i4.4,a,i2.2,a,i2.2,a,i3.3)') '======>> time-step =', kt,   & 
    231248              &   '      New day, DATE Y/M/D = ', nyear, '/', nmonth, '/', nday, '      nday_year = ', nday_year 
    232249         IF(lwp) WRITE(numout,'(a,i8,a,i7,a,i5)') '         nsec_year = ', nsec_year,   & 
    233               &   '   nsec_month = ', nsec_month, '   nsec_day = ', nsec_day 
    234       ENDIF 
     250              &   '   nsec_month = ', nsec_month, '   nsec_day = ', nsec_day, '   nsec_week = ', nsec_week 
     251      ENDIF 
     252 
     253      IF( nsec_week .GT. 7*86400 ) nsec_week = ndt05 
    235254       
    236255      IF(ln_ctl) THEN 

  • branches/DEV_r2191_3partymerge2010/NEMO/OPA_SRC/DOM/dom_oce.F90

    r2219 r2222  
    195195   !! calendar variables 
    196196   !! --------------------------------------------------------------------- 
    197    INTEGER , PUBLIC ::   nyear       !: current year 
    198    INTEGER , PUBLIC ::   nmonth      !: current month 
    199    INTEGER , PUBLIC ::   nday        !: current day of the month 
    200    INTEGER , PUBLIC ::   ndastp      !: time step date in yyyymmdd format 
    201    INTEGER , PUBLIC ::   nday_year   !: current day counted from jan 1st of the current year 
    202    INTEGER , PUBLIC ::   nsec_year   !: current time step counted in second since 00h jan 1st of the current year 
    203    INTEGER , PUBLIC ::   nsec_month  !: current time step counted in second since 00h 1st day of the current month 
    204    INTEGER , PUBLIC ::   nsec_day    !: current time step counted in second since 00h of the current day 
    205    REAL(wp), PUBLIC ::   fjulday     !: julian day  
    206    REAL(wp), PUBLIC ::   adatrj      !: number of elapsed days since the begining of the whole simulation 
    207    !                                 !: (cumulative duration of previous runs that may have used different time-step size) 
     197   INTEGER , PUBLIC ::   nyear         !: current year 
     198   INTEGER , PUBLIC ::   nmonth        !: current month 
     199   INTEGER , PUBLIC ::   nday          !: current day of the month 
     200   INTEGER , PUBLIC ::   ndastp        !: time step date in yyyymmdd format 
     201   INTEGER , PUBLIC ::   nday_year     !: current day counted from jan 1st of the current year 
     202   INTEGER , PUBLIC ::   nsec_year     !: current time step counted in second since 00h jan 1st of the current year 
     203   INTEGER , PUBLIC ::   nsec_month    !: current time step counted in second since 00h 1st day of the current month 
     204   INTEGER , PUBLIC ::   nsec_week     !: current time step counted in second since 00h of last monday 
     205   INTEGER , PUBLIC ::   nsec_day      !: current time step counted in second since 00h of the current day 
     206   REAL(wp), PUBLIC ::   fjulday       !: current julian day  
     207   REAL(wp), PUBLIC ::   fjulstartyear !: first day of the current year in julian days 
     208   REAL(wp), PUBLIC ::   adatrj        !: number of elapsed days since the begining of the whole simulation 
     209   !                                   !: (cumulative duration of previous runs that may have used different time-step size) 
    208210   INTEGER , PUBLIC, DIMENSION(0: 1) ::   nyear_len     !: length in days of the previous/current year 
    209211   INTEGER , PUBLIC, DIMENSION(0:13) ::   nmonth_len    !: length in days of the months of the current year 

  • branches/DEV_r2191_3partymerge2010/NEMO/OPA_SRC/DTA/dtasal.F90

    r1715 r2222  
    1313   USE oce             ! ocean dynamics and tracers 
    1414   USE dom_oce         ! ocean space and time domain 
     15   USE fldread         ! read input fields 
    1516   USE in_out_manager  ! I/O manager 
    1617   USE phycst          ! physical constants 
    … …  
    2728   !! * Shared module variables 
    2829   LOGICAL , PUBLIC, PARAMETER ::   lk_dtasal = .TRUE.    !: salinity data flag 
    29    REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   &  !: 
    30       s_dta       !: salinity data at given time-step 
     30   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   s_dta    !: salinity data at given time-step 
    3131 
    3232   !! * Module variables 
    33    INTEGER ::   & 
    34       numsdt,           &  !: logical unit for data salinity 
    35       nsal1, nsal2         ! first and second record used 
    36    REAL(wp), DIMENSION(jpi,jpj,jpk,2) ::   & 
    37       saldta    ! salinity data at two consecutive times 
     33   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_sal       ! structure of input SST (file informations, fields read) 
    3834 
    3935   !! * Substitutions 
    … …  
    5248 
    5349   SUBROUTINE dta_sal( kt ) 
    54      !!---------------------------------------------------------------------- 
    55      !!                   ***  ROUTINE dta_sal  *** 
    56      !!         
    57      !! ** Purpose :   Reads monthly salinity data 
    58      !!              
    59      !! ** Method  : - Read on unit numsdt the monthly salinity data interpo- 
    60      !!     lated onto the model grid. 
    61      !!              - At each time step, a linear interpolation is applied 
    62      !!     between two monthly values. 
    63      !! 
    64      !! History : 
    65      !!        !  91-03  ()  Original code 
    66      !!        !  92-07  (M. Imbard) 
    67      !!   9.0  !  02-06  (G. Madec)  F90: Free form and module  
    68      !!---------------------------------------------------------------------- 
    69      !! * Modules used 
    70      USE iom 
    71       
    72      !! * Arguments 
    73      INTEGER, INTENT(in) ::   kt             ! ocean time step 
    74       
    75      !! * Local declarations 
    76       
    77      INTEGER ::   ji, jj, jk, jl, jkk   ! dummy loop indicies 
    78      INTEGER ::   & 
    79           imois, iman, i15, ik           ! temporary integers 
    80 #  if defined key_tradmp 
    81      INTEGER ::   & 
     50      !!---------------------------------------------------------------------- 
     51      !!                   ***  ROUTINE dta_sal  *** 
     52      !!         
     53      !! ** Purpose :   Reads monthly salinity data 
     54      !!              
     55      !! ** Method  : - Read on unit numsdt the monthly salinity data interpo- 
     56      !!     lated onto the model grid. 
     57      !!              - At each time step, a linear interpolation is applied 
     58      !!     between two monthly values. 
     59      !! 
     60      !! History : 
     61      !!        !  91-03  ()  Original code 
     62      !!        !  92-07  (M. Imbard) 
     63      !!   9.0  !  02-06  (G. Madec)  F90: Free form and module  
     64      !!---------------------------------------------------------------------- 
     65      
     66      !! * Arguments 
     67      INTEGER, INTENT(in) ::   kt             ! ocean time step 
     68       
     69      !! * Local declarations 
     70      
     71      INTEGER ::   ji, jj, jk, jl, jkk   ! dummy loop indicies 
     72      INTEGER ::   & 
     73           imois, iman, i15, ik           ! temporary integers 
     74      INTEGER            :: ierror 
     75#if defined key_tradmp 
     76      INTEGER ::   & 
    8277          il0, il1, ii0, ii1, ij0, ij1   ! temporary integers          
    83 # endif 
    84      REAL(wp) ::   zxy, zl 
    85 #if defined key_orca_lev10 
    86      REAL(wp), DIMENSION(jpi,jpj,jpkdta,2) :: zsal 
    87      INTEGER   :: ikr, ikw, ikt, jjk 
    88      REAL(wp)  :: zfac 
    89 #endif 
    90      REAL(wp), DIMENSION(jpk,2) ::   & 
     78#endif 
     79      REAL(wp) ::   zxy, zl 
     80#if defined key_orca_lev10 
     81      INTEGER   :: ikr, ikw, ikt, jjk 
     82      REAL(wp)  :: zfac 
     83#endif 
     84      REAL(wp), DIMENSION(jpk) ::   & 
    9185          zsaldta            ! auxiliary array for interpolation 
    92      !!---------------------------------------------------------------------- 
    93       
    94      ! 0. Initialization 
    95      ! ----------------- 
    96       
    97      iman  = INT( raamo ) 
    98 !!! better but change the results     i15 = INT( 2*FLOAT( nday ) / ( FLOAT( nobis(nmonth) ) + 0.5 ) ) 
    99      i15   = nday / 16 
    100      imois = nmonth + i15 - 1 
    101      IF( imois == 0 ) imois = iman 
    102       
    103      ! 1. First call kt=nit000 
    104      ! ----------------------- 
    105       
    106      IF( kt == nit000 ) THEN 
    107          
    108         nsal1 = 0   ! initializations 
    109         IF(lwp) WRITE(numout,*) ' dta_sal : monthly salinity data in NetCDF file' 
    110         CALL iom_open ( 'data_1m_salinity_nomask', numsdt )  
    111          
    112      ENDIF 
    113       
    114       
    115      ! 2. Read monthly file 
    116      ! ------------------- 
    117       
    118      IF( kt == nit000 .OR. imois /= nsal1 ) THEN 
    119          
    120         ! 2.1 Calendar computation 
    121          
    122         nsal1 = imois        ! first file record used  
    123         nsal2 = nsal1 + 1    ! last  file record used 
    124         nsal1 = MOD( nsal1, iman ) 
    125         IF( nsal1 == 0 ) nsal1 = iman 
    126         nsal2 = MOD( nsal2, iman ) 
    127         IF( nsal2 == 0 ) nsal2 = iman 
    128         IF(lwp) WRITE(numout,*) 'first record file used nsal1 ', nsal1 
    129         IF(lwp) WRITE(numout,*) 'last  record file used nsal2 ', nsal2 
    130          
    131         ! 2.3 Read monthly salinity data Levitus  
    132          
    133 #if defined key_orca_lev10 
    134         if (ln_zps) stop 
    135         zsal(:,:,:,:) = 0. 
    136         CALL iom_get (numsdt,jpdom_data,'vosaline',zsal(:,:,:,1),nsal1) 
    137         CALL iom_get (numsdt,jpdom_data,'vosaline',zsal(:,:,:,2),nsal2) 
     86      CHARACTER(len=100) :: cn_dir          ! Root directory for location of ssr files 
     87      TYPE(FLD_N)        :: sn_sal 
     88      LOGICAL , SAVE     :: linit_sal = .FALSE. 
     89      !!---------------------------------------------------------------------- 
     90      NAMELIST/namdta_sal/cn_dir,sn_sal 
     91      
     92      ! 1. Initialization 
     93      ! ----------------------- 
     94      
     95      IF( kt == nit000 .AND. ( .NOT. linit_sal ) ) THEN 
     96         
     97         !                         ! set file information 
     98         cn_dir = './'             ! directory in which the model is executed 
     99         ! ... default values (NB: frequency positive => hours, negative => months) 
     100         !            !   file    ! frequency !  variable  ! time intep !  clim   ! 'yearly' or ! weights  ! rotation   ! 
     101         !            !   name    !  (hours)  !   name     !   (T/F)    !  (T/F)  !  'monthly'  ! filename ! pairs      ! 
     102         sn_sal = FLD_N( 'salinity',  -1.  ,  'vosaline',  .false.   , .true.  ,  'monthly'  , ''       , ''         ) 
     103 
     104         REWIND ( numnam )         ! ... read in namlist namdta_sal  
     105         READ( numnam, namdta_sal )  
     106 
     107         IF(lwp) THEN              ! control print 
     108            WRITE(numout,*) 
     109            WRITE(numout,*) 'dta_sal : Salinity Climatology ' 
     110            WRITE(numout,*) '~~~~~~~ ' 
     111         ENDIF 
     112         ALLOCATE( sf_sal(1), STAT=ierror ) 
     113         IF( ierror > 0 ) THEN 
     114             CALL ctl_stop( 'dta_sal: unable to allocate sf_sal structure' )   ;   RETURN 
     115         ENDIF 
     116 
     117#if defined key_orca_lev10 
     118         ALLOCATE( sf_sal(1)%fnow(jpi,jpj,jpkdta)   ) 
     119         IF( sn_sal%ln_tint ) ALLOCATE( sf_sal(1)%fdta(jpi,jpj,jpkdta,2) ) 
    138120#else 
    139         CALL iom_get (numsdt,jpdom_data,'vosaline',saldta(:,:,:,1),nsal1) 
    140         CALL iom_get (numsdt,jpdom_data,'vosaline',saldta(:,:,:,2),nsal2) 
    141 #endif 
    142          
    143         IF(lwp) THEN 
    144            WRITE(numout,*) 
    145            WRITE(numout,*) ' read Levitus salinity ok' 
    146            WRITE(numout,*) 
    147         ENDIF 
    148          
     121         ALLOCATE( sf_sal(1)%fnow(jpi,jpj,jpk)   ) 
     122         IF( sn_sal%ln_tint ) ALLOCATE( sf_sal(1)%fdta(jpi,jpj,jpk,2) ) 
     123#endif 
     124         ! fill sf_sal with sn_sal and control print 
     125         CALL fld_fill( sf_sal, (/ sn_sal /), cn_dir, 'dta_sal', 'Salinity data', 'namdta_sal' ) 
     126         linit_sal = .TRUE.         
     127      ENDIF 
     128      
     129      
     130      ! 2. Read monthly file 
     131      ! ------------------- 
     132      
     133      CALL fld_read( kt, 1, sf_sal ) 
     134 
     135      IF( lwp .AND. kt==nn_it000 ) THEN 
     136         WRITE(numout,*) 
     137         WRITE(numout,*) ' read Levitus salinity ok' 
     138         WRITE(numout,*) 
     139      ENDIF 
     140 
    149141#if defined key_tradmp 
    150         IF( cp_cfg == "orca"  .AND. jp_cfg == 2 ) THEN 
     142      IF( cp_cfg == "orca"  .AND. jp_cfg == 2 ) THEN 
     143    
     144         !                                        ! ======================= 
     145         !                                        !  ORCA_R2 configuration 
     146         !                                        ! ======================= 
     147         ij0 = 101   ;   ij1 = 109 
     148         ii0 = 141   ;   ii1 = 155    
     149         DO jj = mj0(ij0), mj1(ij1)                  ! Reduced salinity in the Alboran Sea 
     150            DO ji = mi0(ii0), mi1(ii1) 
     151               sf_sal(1)%fnow(ji,jj,13:13) = sf_sal(1)%fnow(ji,jj,13:13) - 0.15 
     152               sf_sal(1)%fnow(ji,jj,14:15) = sf_sal(1)%fnow(ji,jj,14:15) - 0.25 
     153               sf_sal(1)%fnow(ji,jj,16:17) = sf_sal(1)%fnow(ji,jj,16:17) - 0.30 
     154               sf_sal(1)%fnow(ji,jj,18:25) = sf_sal(1)%fnow(ji,jj,18:25) - 0.35 
     155            END DO 
     156         END DO 
     157 
     158         IF( n_cla == 1 ) THEN  
     159            !                                         ! New salinity profile at Gibraltar 
     160            il0 = 138   ;   il1 = 138    
     161            ij0 = 101   ;   ij1 = 102 
     162            ii0 = 139   ;   ii1 = 139    
     163            DO jl = mi0(il0), mi1(il1) 
     164               DO jj = mj0(ij0), mj1(ij1) 
     165                  DO ji = mi0(ii0), mi1(ii1) 
     166                        sf_sal(1)%fnow(ji,jj,:) = sf_sal(1)%fnow(jl,jj,:) 
     167                  END DO 
     168               END DO 
     169            END DO 
     170            !                                         ! New salinity profile at Bab el Mandeb 
     171            il0 = 164   ;   il1 = 164    
     172            ij0 =  87   ;   ij1 =  88 
     173            ii0 = 161   ;   ii1 = 163    
     174            DO jl = mi0(il0), mi1(il1) 
     175               DO jj = mj0(ij0), mj1(ij1) 
     176                  DO ji = mi0(ii0), mi1(ii1) 
     177                     sf_sal(1)%fnow(ji,jj,:) = sf_sal(1)%fnow(jl,jj,:) 
     178                  END DO 
     179               END DO 
     180            END DO 
     181            ! 
     182         ENDIF 
     183            ! 
     184      ENDIF 
     185#endif    
     186         
     187#if defined key_orca_lev10 
     188      DO jjk = 1, 5 
     189         s_dta(:,:,jjk) = sf_sal(1)%fnow(:,:,1) 
     190      ENDDO 
     191      DO jk = 1, jpk-20,10 
     192         ikr =  INT(jk/10) + 1 
     193         ikw =  (ikr-1) *10 + 1 
     194         ikt =  ikw + 5 
     195         DO jjk=ikt,ikt+9 
     196            zfac = ( gdept_0(jjk   ) - gdepw_0(ikt) ) / ( gdepw_0(ikt+10) - gdepw_0(ikt) ) 
     197            s_dta(:,:,jjk) = sf_sal(1)%fnow(:,:,ikr) + ( sf_sal(1)%fnow(:,:,ikr+1) - sf_sal(1)%fnow(:,:,ikr) ) * zfac 
     198         END DO 
     199      END DO 
     200      DO jjk = jpk-5, jpk 
     201         s_dta(:,:,jjk) = sf_sal(1)%fnow(:,:,jpkdta-1) 
     202      END DO 
     203      ! fill the overlap areas 
     204      CALL lbc_lnk (s_dta(:,:,:),'Z',-999.,'no0')         
     205#else 
     206      s_dta(:,:,:)=sf_sal(1)%fnow(:,:,:) 
     207#endif 
     208         
     209      IF( ln_sco ) THEN 
     210         DO jj = 1, jpj                  ! interpolation of salinites 
     211            DO ji = 1, jpi 
     212               DO jk = 1, jpk 
     213                  zl=fsdept_0(ji,jj,jk) 
     214                  IF(zl < gdept_0(1)  ) zsaldta(jk) =  s_dta(ji,jj,1    )  
     215                  IF(zl > gdept_0(jpk)) zsaldta(jk) =  s_dta(ji,jj,jpkm1)  
     216                  DO jkk = 1, jpkm1 
     217                     IF((zl-gdept_0(jkk))*(zl-gdept_0(jkk+1)).le.0.0) THEN 
     218                          zsaldta(jk) = s_dta(ji,jj,jkk)                                 & 
     219                                     &           + (zl-gdept_0(jkk))/(gdept_0(jkk+1)-gdept_0(jkk))      & 
     220                                     &                              *(s_dta(ji,jj,jkk+1) - s_dta(ji,jj,jkk)) 
     221                     ENDIF 
     222                  END DO 
     223               END DO 
     224               DO jk = 1, jpkm1 
     225                  s_dta(ji,jj,jk) = zsaldta(jk)  
     226               END DO 
     227               s_dta(ji,jj,jpk) = 0.0  
     228            END DO 
     229         END DO 
    151230            
    152            !                                        ! ======================= 
    153            !                                        !  ORCA_R2 configuration 
    154            !                                        ! ======================= 
    155            ij0 = 101   ;   ij1 = 109 
    156            ii0 = 141   ;   ii1 = 155    
    157            DO jj = mj0(ij0), mj1(ij1)                  ! Reduced salinity in the Alboran Sea 
    158               DO ji = mi0(ii0), mi1(ii1) 
    159 #if defined key_orca_lev10 
    160                  zsal  (ji,jj,13:13,:) = zsal  (ji,jj,13:13,:) - 0.15 
    161                  zsal  (ji,jj,14:15,:) = zsal  (ji,jj,14:15,:) - 0.25 
    162                  zsal  (ji,jj,16:17,:) = zsal  (ji,jj,16:17,:) - 0.30 
    163                  zsal  (ji,jj,18:25,:) = zsal  (ji,jj,18:25,:) - 0.35 
    164 #else 
    165                  saldta(ji,jj,13:13,:) = saldta(ji,jj,13:13,:) - 0.15 
    166                  saldta(ji,jj,14:15,:) = saldta(ji,jj,14:15,:) - 0.25 
    167                  saldta(ji,jj,16:17,:) = saldta(ji,jj,16:17,:) - 0.30 
    168                  saldta(ji,jj,18:25,:) = saldta(ji,jj,18:25,:) - 0.35 
    169 #endif 
    170               END DO 
    171            END DO 
    172  
    173            IF( n_cla == 1 ) THEN  
    174               !                                         ! New salinity profile at Gibraltar 
    175               il0 = 138   ;   il1 = 138    
    176               ij0 = 101   ;   ij1 = 102 
    177               ii0 = 139   ;   ii1 = 139    
    178               DO jl = mi0(il0), mi1(il1) 
    179                  DO jj = mj0(ij0), mj1(ij1) 
    180                     DO ji = mi0(ii0), mi1(ii1) 
    181 #if defined key_orca_lev10 
    182                        zsal  (ji,jj,:,:) = zsal  (jl,jj,:,:) 
    183 #else 
    184                        saldta(ji,jj,:,:) = saldta(jl,jj,:,:) 
    185 #endif 
    186                     END DO 
    187                  END DO 
    188               END DO 
    189               !                                         ! New salinity profile at Bab el Mandeb 
    190               il0 = 164   ;   il1 = 164    
    191               ij0 =  87   ;   ij1 =  88 
    192               ii0 = 161   ;   ii1 = 163    
    193               DO jl = mi0(il0), mi1(il1) 
    194                  DO jj = mj0(ij0), mj1(ij1) 
    195                     DO ji = mi0(ii0), mi1(ii1) 
    196 #if defined key_orca_lev10 
    197                        zsal  (ji,jj,:,:) = zsal  (jl,jj,:,:) 
    198 #else 
    199                        saldta(ji,jj,:,:) = saldta(jl,jj,:,:) 
    200 #endif 
    201                     END DO 
    202                  END DO 
    203               END DO 
    204               ! 
    205            ENDIF 
    206            ! 
    207         ENDIF 
    208 #endif    
    209          
    210 #if defined key_orca_lev10 
    211         !  interpolate from 31 to 301 level the zsal field result in saldta 
    212         DO jl = 1, 2 
    213            DO jjk = 1, 5 
    214               saldta(:,:,jjk,jl) = zsal(:,:,1,jl) 
    215            ENDDO 
    216            DO jk = 1, jpk - 20, 10 
    217               ikr = INT( jk / 10 ) + 1 
    218               ikw = (ikr-1) * 10 + 1 
    219               ikt = ikw + 5 
    220               DO jjk = ikt , ikt + 9 
    221                  zfac = ( gdept_0(jjk) - gdepw_0(ikt) ) / ( gdepw_0(ikt+10) - gdepw_0(ikt) ) 
    222                  saldta(:,:,jjk,jl) = zsal(:,:,ikr,jl) + ( zsal(:,:,ikr+1,jl) - zsal(:,:,ikr,jl) ) * zfac 
    223               END DO 
    224            END DO 
    225            DO jjk = jpk-5, jpk 
    226               saldta(:,:,jjk,jl) = zsal(:,:,jpkdta-1,jl) 
    227            END DO 
    228            ! fill the overlap areas 
    229            CALL lbc_lnk (saldta(:,:,:,jl),'Z',-999.,'no0') 
    230         END DO 
    231          
    232 #endif 
    233          
    234         IF( ln_sco ) THEN 
    235            DO jl = 1, 2 
    236               DO jj = 1, jpj                  ! interpolation of salinites 
    237                  DO ji = 1, jpi 
    238                     DO jk = 1, jpk 
    239                        zl=fsdept_0(ji,jj,jk) 
    240                        IF(zl <  gdept_0(1)) zsaldta(jk,jl) =  saldta(ji,jj,1,jl) 
    241                        IF(zl >  gdept_0(jpk)) zsaldta(jk,jl) =  saldta(ji,jj,jpkm1,jl) 
    242                        DO jkk = 1, jpkm1 
    243                           IF((zl-gdept_0(jkk))*(zl-gdept_0(jkk+1)).le.0.0) THEN 
    244                              zsaldta(jk,jl) = saldta(ji,jj,jkk,jl)                                  & 
    245                                   &           + (zl-gdept_0(jkk))/(gdept_0(jkk+1)-gdept_0(jkk))       & 
    246                                   &                              *(saldta(ji,jj,jkk+1,jl) - saldta(ji,jj,jkk,jl)) 
    247                           ENDIF 
    248                        END DO 
    249                     END DO 
    250                     DO jk = 1, jpkm1 
    251                        saldta(ji,jj,jk,jl) = zsaldta(jk,jl) 
    252                     END DO 
    253                     saldta(ji,jj,jpk,jl) = 0.0 
    254                  END DO 
    255               END DO 
    256            END DO 
    257             
    258            IF(lwp) WRITE(numout,*) 
    259            IF(lwp) WRITE(numout,*) ' Levitus salinity data interpolated to s-coordinate' 
    260            IF(lwp) WRITE(numout,*) 
    261             
    262         ELSE 
    263            !                                  ! Mask 
    264            DO jl = 1, 2 
    265               saldta(:,:,:,jl) = saldta(:,:,:,jl)*tmask(:,:,:) 
    266               saldta(:,:,jpk,jl) = 0. 
    267               IF( ln_zps ) THEN               ! z-coord. partial steps 
    268                  DO jj = 1, jpj               ! interpolation of salinity at the last ocean level (i.e. the partial step) 
    269                     DO ji = 1, jpi 
    270                        ik = mbathy(ji,jj) - 1 
    271                        IF( ik > 2 ) THEN 
    272                           zl = ( gdept_0(ik) - fsdept_0(ji,jj,ik) ) / ( gdept_0(ik) - gdept_0(ik-1) ) 
    273                           saldta(ji,jj,ik,jl) = (1.-zl) * saldta(ji,jj,ik,jl) +zl * saldta(ji,jj,ik-1,jl) 
    274                        ENDIF 
    275                     END DO 
    276                  END DO 
    277               ENDIF 
    278            END DO 
    279         ENDIF 
    280          
    281          
    282         IF(lwp) THEN 
    283            WRITE(numout,*)' salinity Levitus month ',nsal1,nsal2 
    284            WRITE(numout,*) 
    285            WRITE(numout,*) ' Levitus month = ',nsal1,'  level = 1' 
    286            CALL prihre(saldta(:,:,1,1),jpi,jpj,1,jpi,20,1,jpj,20,1.,numout) 
    287            WRITE(numout,*) ' Levitus month = ',nsal1,'  level = ',jpk/2 
    288            CALL prihre(saldta(:,:,jpk/2,1),jpi,jpj,1,jpi,20,1,jpj,20,1.,numout) 
    289            WRITE(numout,*) ' Levitus month = ',nsal1,'  level = ',jpkm1 
    290            CALL prihre(saldta(:,:,jpkm1,1),jpi,jpj,1,jpi,20,1,jpj,20,1.,numout) 
    291         ENDIF 
    292      ENDIF 
    293       
    294       
    295      ! 3. At every time step compute salinity data 
    296      ! ------------------------------------------- 
    297       
    298      zxy = FLOAT(nday + 15 - 30*i15)/30. 
    299      s_dta(:,:,:) = ( 1.- zxy ) * saldta(:,:,:,1) + zxy * saldta(:,:,:,2) 
    300       
    301      ! Close the file 
    302      ! -------------- 
    303       
    304      IF( kt == nitend )   CALL iom_close (numsdt) 
     231         IF( lwp .AND. kt==nn_it000 ) THEN 
     232            WRITE(numout,*) 
     233            WRITE(numout,*) ' Levitus salinity data interpolated to s-coordinate' 
     234            WRITE(numout,*) 
     235         ENDIF 
     236 
     237      ELSE 
     238         !                                  ! Mask 
     239         s_dta(:,:,:) = s_dta(:,:,:) * tmask(:,:,:) 
     240         s_dta(:,:,jpk) = 0.  
     241         IF( ln_zps ) THEN               ! z-coord. partial steps 
     242            DO jj = 1, jpj               ! interpolation of salinity at the last ocean level (i.e. the partial step) 
     243               DO ji = 1, jpi 
     244                  ik = mbathy(ji,jj) - 1 
     245                  IF( ik > 2 ) THEN 
     246                     zl = ( gdept_0(ik) - fsdept_0(ji,jj,ik) ) / ( gdept_0(ik) - gdept_0(ik-1) ) 
     247                     s_dta(ji,jj,ik) = (1.-zl) * s_dta(ji,jj,ik) + zl * s_dta(ji,jj,ik-1) 
     248                  ENDIF 
     249               END DO 
     250            END DO 
     251         ENDIF 
     252      ENDIF 
     253         
     254      IF( lwp .AND. kt==nn_it000 ) THEN 
     255         WRITE(numout,*)' salinity Levitus ' 
     256         WRITE(numout,*) 
     257         WRITE(numout,*)'  level = 1' 
     258         CALL prihre(s_dta(:,:,1),    jpi,jpj,1,jpi,20,1,jpj,20,1.,numout) 
     259         WRITE(numout,*)'  level = ',jpk/2 
     260         CALL prihre(s_dta(:,:,jpk/2),jpi,jpj,1,jpi,20,1,jpj,20,1.,numout)            
     261         WRITE(numout,*) '  level = ',jpkm1 
     262         CALL prihre(s_dta(:,:,jpkm1),jpi,jpj,1,jpi,20,1,jpj,20,1.,numout) 
     263      ENDIF 
    305264 
    306265   END SUBROUTINE dta_sal 

  • branches/DEV_r2191_3partymerge2010/NEMO/OPA_SRC/DTA/dtatem.F90

    r1715 r2222  
    1313   USE oce             ! ocean dynamics and tracers 
    1414   USE dom_oce         ! ocean space and time domain 
     15   USE fldread         ! read input fields 
    1516   USE in_out_manager  ! I/O manager 
    1617   USE phycst          ! physical constants 
    … …  
    2627   !! * Shared module variables 
    2728   LOGICAL , PUBLIC, PARAMETER ::   lk_dtatem = .TRUE.   !: temperature data flag 
    28    REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   &  !: 
    29       t_dta             !: temperature data at given time-step 
     29   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::  t_dta    !: temperature data at given time-step 
    3030 
    3131   !! * Module variables 
    32    INTEGER ::   & 
    33       numtdt,        &  !: logical unit for data temperature 
    34       ntem1, ntem2  ! first and second record used 
    35    REAL(wp), DIMENSION(jpi,jpj,jpk,2) ::   & 
    36       temdta            ! temperature data at two consecutive times 
     32   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_tem      ! structure of input SST (file informations, fields read) 
    3733 
    3834   !! * Substitutions 
    … …  
    7369      !!   8.5  !  02-09  (G. Madec)  F90: Free form and module 
    7470      !!---------------------------------------------------------------------- 
    75       !! * Modules used 
    76       USE iom 
    77  
    7871      !! * Arguments 
    7972      INTEGER, INTENT( in ) ::   kt     ! ocean time-step 
    8073 
    8174      !! * Local declarations 
    82       INTEGER ::   ji, jj, jl, jk, jkk       ! dummy loop indicies 
     75      INTEGER ::   ji, jj, jk, jl, jkk       ! dummy loop indicies 
    8376      INTEGER ::   & 
    84          imois, iman, i15 , ik      ! temporary integers 
    85 #  if defined key_tradmp 
     77        imois, iman, i15 , ik      ! temporary integers 
     78      INTEGER            :: ierror 
     79#if defined key_tradmp 
    8680      INTEGER ::   & 
    8781         il0, il1, ii0, ii1, ij0, ij1   ! temporary integers 
    88 # endif 
     82#endif 
    8983      REAL(wp) ::   zxy, zl 
    9084#if defined key_orca_lev10 
    91       REAL(wp), DIMENSION(jpi,jpj,jpkdta,2) :: ztem 
     85      !!!REAL(wp), DIMENSION(jpi,jpj,jpkdta,2) :: ztem 
    9286      INTEGER   :: ikr, ikw, ikt, jjk  
    9387      REAL(wp)  :: zfac 
    9488#endif 
    95       REAL(wp), DIMENSION(jpk,2) ::   & 
     89      REAL(wp), DIMENSION(jpk) ::   & 
    9690         ztemdta            ! auxiliary array for interpolation 
     91      CHARACTER(len=100) :: cn_dir          ! Root directory for location of ssr files 
     92      TYPE(FLD_N)        :: sn_tem 
     93      LOGICAL , SAVE     :: linit_tem = .FALSE. 
    9794      !!---------------------------------------------------------------------- 
    98        
    99       ! 0. Initialization 
    100       ! ----------------- 
    101        
    102       iman  = INT( raamo ) 
    103 !!! better but change the results     i15 = INT( 2*FLOAT( nday ) / ( FLOAT( nobis(nmonth) ) + 0.5 ) ) 
    104       i15   = nday / 16 
    105       imois = nmonth + i15 - 1 
    106       IF( imois == 0 ) imois = iman 
    107        
    108       ! 1. First call kt=nit000 
     95      NAMELIST/namdta_tem/cn_dir,sn_tem 
     96  
     97      ! 1. Initialization  
    10998      ! ----------------------- 
    11099       
    111       IF( kt == nit000 ) THEN 
    112           
    113          ntem1= 0   ! initializations 
    114          IF(lwp) WRITE(numout,*) ' dta_tem : Levitus monthly fields' 
    115          CALL iom_open ( 'data_1m_potential_temperature_nomask', numtdt )  
    116           
    117       ENDIF 
    118        
     100      IF( kt == nit000 .AND. (.NOT. linit_tem ) ) THEN 
     101 
     102         !                   ! set file information 
     103         cn_dir = './'       ! directory in which the model is executed 
     104         ! ... default values (NB: frequency positive => hours, negative => months) 
     105         !            !   file    ! frequency !  variable  ! time intep !  clim   ! 'yearly' or ! weights  ! rotation   ! 
     106         !            !   name    !  (hours)  !   name     !   (T/F)    !  (T/F)  !  'monthly'  ! filename ! pairs      ! 
     107         sn_tem = FLD_N( 'temperature',  -1.  ,  'votemper',  .false.   , .true.  ,  'yearly'  , ''       , ''         ) 
     108 
     109         REWIND( numnam )         ! ... read in namlist namdta_tem  
     110         READ( numnam, namdta_tem )  
     111 
     112         IF(lwp) THEN              ! control print 
     113            WRITE(numout,*) 
     114            WRITE(numout,*) 'dta_tem : Temperature Climatology ' 
     115            WRITE(numout,*) '~~~~~~~ ' 
     116         ENDIF 
     117         ALLOCATE( sf_tem(1), STAT=ierror ) 
     118         IF( ierror > 0 ) THEN 
     119             CALL ctl_stop( 'dta_tem: unable to allocate sf_tem structure' )   ;   RETURN 
     120         ENDIF 
     121 
     122#if defined key_orca_lev10 
     123         ALLOCATE( sf_tem(1)%fnow(jpi,jpj,jpkdta)   ) 
     124         IF( sn_tem%ln_tint ) ALLOCATE( sf_tem(1)%fdta(jpi,jpj,jpkdta,2) ) 
     125#else 
     126         ALLOCATE( sf_tem(1)%fnow(jpi,jpj,jpk)   ) 
     127         IF( sn_tem%ln_tint ) ALLOCATE( sf_tem(1)%fdta(jpi,jpj,jpk,2) ) 
     128#endif 
     129         ! fill sf_tem with sn_tem and control print 
     130         CALL fld_fill( sf_tem, (/ sn_tem /), cn_dir, 'dta_tem', 'Temperature data', 'namdta_tem' ) 
     131         linit_tem = .TRUE. 
     132 
     133      ENDIF 
    119134       
    120135      ! 2. Read monthly file 
    121136      ! ------------------- 
    122        
    123       IF( kt == nit000 .OR. imois /= ntem1 ) THEN 
    124           
    125          ! Calendar computation 
    126           
    127          ntem1 = imois        ! first file record used  
    128          ntem2 = ntem1 + 1    ! last  file record used 
    129          ntem1 = MOD( ntem1, iman ) 
    130          IF( ntem1 == 0 )   ntem1 = iman 
    131          ntem2 = MOD( ntem2, iman ) 
    132          IF( ntem2 == 0 )   ntem2 = iman 
    133          IF(lwp) WRITE(numout,*) 'first record file used ntem1 ', ntem1 
    134          IF(lwp) WRITE(numout,*) 'last  record file used ntem2 ', ntem2 
    135           
    136          ! Read monthly temperature data Levitus  
    137           
    138 #if defined key_orca_lev10 
    139          if (ln_zps) stop 
    140          ztem(:,:,:,:) = 0. 
    141          CALL iom_get (numtdt,jpdom_data,'votemper',ztem(:,:,:,1),ntem1) 
    142          CALL iom_get (numtdt,jpdom_data,'votemper',ztem(:,:,:,2),ntem2) 
    143 #else          
    144          CALL iom_get (numtdt,jpdom_data,'votemper',temdta(:,:,:,1),ntem1) 
    145          CALL iom_get (numtdt,jpdom_data,'votemper',temdta(:,:,:,2),ntem2) 
    146 #endif 
    147           
    148          IF(lwp) WRITE(numout,*) 
    149          IF(lwp) WRITE(numout,*) ' read Levitus temperature ok' 
    150          IF(lwp) WRITE(numout,*) 
     137          
     138      CALL fld_read( kt, 1, sf_tem ) 
     139        
     140      IF( lwp .AND. kt==nn_it000 )THEN  
     141         WRITE(numout,*) 
     142         WRITE(numout,*) ' read Levitus temperature ok' 
     143         WRITE(numout,*) 
     144      ENDIF 
    151145          
    152146#if defined key_tradmp 
    153          IF( cp_cfg == "orca"  .AND. jp_cfg == 2 ) THEN 
    154              
    155             !                                        ! ======================= 
    156             !                                        !  ORCA_R2 configuration 
    157             !                                        ! =======================  
    158             ij0 = 101   ;   ij1 = 109 
    159             ii0 = 141   ;   ii1 = 155 
    160             DO jj = mj0(ij0), mj1(ij1)                      ! Reduced temperature in the Alboran Sea 
    161                DO ji = mi0(ii0), mi1(ii1) 
    162 #if defined key_orca_lev10 
    163                   ztem(  ji,jj, 13:13 ,:) = ztem  (ji,jj, 13:13 ,:) - 0.20 
    164                   ztem  (ji,jj, 14:15 ,:) = ztem  (ji,jj, 14:15 ,:) - 0.35 
    165                   ztem  (ji,jj, 16:25 ,:) = ztem  (ji,jj, 16:25 ,:) - 0.40 
     147      IF( cp_cfg == "orca"  .AND. jp_cfg == 2 ) THEN 
     148             
     149         !                                        ! ======================= 
     150         !                                        !  ORCA_R2 configuration 
     151         !                                        ! =======================  
     152         ij0 = 101   ;   ij1 = 109 
     153         ii0 = 141   ;   ii1 = 155 
     154         DO jj = mj0(ij0), mj1(ij1)                      ! Reduced temperature in the Alboran Sea 
     155            DO ji = mi0(ii0), mi1(ii1) 
     156               sf_tem(1)%fnow(ji,jj, 13:13 ) = sf_tem(1)%fnow(ji,jj, 13:13 ) - 0.20 
     157               sf_tem(1)%fnow(ji,jj, 14:15 ) = sf_tem(1)%fnow(ji,jj, 14:15 ) - 0.35   
     158               sf_tem(1)%fnow(ji,jj, 16:25 ) = sf_tem(1)%fnow(ji,jj, 16:25 ) - 0.40 
     159            END DO 
     160         END DO 
     161             
     162         IF( n_cla == 1 ) THEN  
     163            !                                         ! New temperature profile at Gibraltar 
     164            il0 = 138   ;   il1 = 138 
     165            ij0 = 101   ;   ij1 = 102 
     166            ii0 = 139   ;   ii1 = 139 
     167            DO jl = mi0(il0), mi1(il1) 
     168               DO jj = mj0(ij0), mj1(ij1) 
     169                  DO ji = mi0(ii0), mi1(ii1) 
     170                     sf_tem(1)%fnow(ji,jj,:) = sf_tem(1)%fnow(jl,jj,:) 
     171                  END DO 
     172               END DO 
     173            END DO 
     174            !                                         ! New temperature profile at Bab el Mandeb 
     175            il0 = 164   ;   il1 = 164 
     176            ij0 =  87   ;   ij1 =  88 
     177            ii0 = 161   ;   ii1 = 163 
     178            DO jl = mi0(il0), mi1(il1) 
     179               DO jj = mj0(ij0), mj1(ij1) 
     180                  DO ji = mi0(ii0), mi1(ii1) 
     181                     sf_tem(1)%fnow(ji,jj,:) = sf_tem(1)%fnow(jl,jj,:) 
     182                  END DO 
     183               END DO 
     184            END DO 
     185            ! 
     186         ELSE 
     187            !                                         ! Reduced temperature at Red Sea 
     188            ij0 =  87   ;   ij1 =  96 
     189            ii0 = 148   ;   ii1 = 160 
     190            sf_tem(1)%fnow( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ,  4:10 ) = 7.0 
     191            sf_tem(1)%fnow( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 11:13 ) = 6.5 
     192            sf_tem(1)%fnow( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 14:20 ) = 6.0 
     193         ENDIF 
     194            ! 
     195      ENDIF 
     196#endif 
     197          
     198#if defined key_orca_lev10 
     199      DO jjk = 1, 5 
     200         t_dta(:,:,jjk) = sf_tem(1)%fnow(:,:,1) 
     201      END DO 
     202      DO jk = 1, jpk-20,10 
     203         ik = jk+5 
     204         ikr =  INT(jk/10) + 1 
     205         ikw =  (ikr-1) *10 + 1 
     206         ikt =  ikw + 5 
     207         DO jjk=ikt,ikt+9 
     208            zfac = ( gdept_0(jjk   ) - gdepw_0(ikt) ) / ( gdepw_0(ikt+10) - gdepw_0(ikt) ) 
     209            t_dta(:,:,jjk) = sf_tem(1)%fnow(:,:,ikr) + ( sf_tem(1)%fnow(:,:,ikr+1) - sf_tem(1)%fnow(:,:,ikr) ) * zfac 
     210         END DO 
     211      END DO 
     212      DO jjk = jpk-5, jpk 
     213         t_dta(:,:,jjk) = sf_tem(1)%fnow(:,:,jpkdta-1) 
     214      END DO 
     215      ! fill the overlap areas 
     216      CALL lbc_lnk (t_dta(:,:,:),'Z',-999.,'no0') 
    166217#else 
    167                   temdta(ji,jj, 13:13 ,:) = temdta(ji,jj, 13:13 ,:) - 0.20 
    168                   temdta(ji,jj, 14:15 ,:) = temdta(ji,jj, 14:15 ,:) - 0.35 
    169                   temdta(ji,jj, 16:25 ,:) = temdta(ji,jj, 16:25 ,:) - 0.40 
    170 #endif 
    171                END DO 
    172             END DO 
    173              
    174             IF( n_cla == 1 ) THEN  
    175                !                                         ! New temperature profile at Gibraltar 
    176                il0 = 138   ;   il1 = 138 
    177                ij0 = 101   ;   ij1 = 102 
    178                ii0 = 139   ;   ii1 = 139 
    179                DO jl = mi0(il0), mi1(il1) 
    180                   DO jj = mj0(ij0), mj1(ij1) 
    181                      DO ji = mi0(ii0), mi1(ii1) 
    182 #if defined key_orca_lev10 
    183                         ztem  (ji,jj,:,:) = ztem  (jl,jj,:,:) 
    184 #else 
    185                         temdta(ji,jj,:,:) = temdta(jl,jj,:,:) 
    186 #endif 
    187                      END DO 
     218      t_dta(:,:,:) = sf_tem(1)%fnow(:,:,:)  
     219#endif 
     220          
     221      IF( ln_sco ) THEN 
     222         DO jj = 1, jpj                  ! interpolation of temperatures 
     223            DO ji = 1, jpi 
     224               DO jk = 1, jpk 
     225                  zl=fsdept_0(ji,jj,jk) 
     226                  IF(zl < gdept_0(1))   ztemdta(jk) =  t_dta(ji,jj,1) 
     227                  IF(zl > gdept_0(jpk)) ztemdta(jk) =  t_dta(ji,jj,jpkm1)  
     228                  DO jkk = 1, jpkm1 
     229                     IF((zl-gdept_0(jkk))*(zl-gdept_0(jkk+1)).le.0.0) THEN 
     230                        ztemdta(jk) = t_dta(ji,jj,jkk)                                 & 
     231                                  &    + (zl-gdept_0(jkk))/(gdept_0(jkk+1)-gdept_0(jkk))  & 
     232                                  &    * (t_dta(ji,jj,jkk+1) - t_dta(ji,jj,jkk)) 
     233                     ENDIF 
    188234                  END DO 
    189235               END DO 
    190                !                                         ! New temperature profile at Bab el Mandeb 
    191                il0 = 164   ;   il1 = 164 
    192                ij0 =  87   ;   ij1 =  88 
    193                ii0 = 161   ;   ii1 = 163 
    194                DO jl = mi0(il0), mi1(il1) 
    195                   DO jj = mj0(ij0), mj1(ij1) 
    196                      DO ji = mi0(ii0), mi1(ii1) 
    197 #if defined key_orca_lev10 
    198                         ztem  (ji,jj,:,:) = ztem  (jl,jj,:,:) 
    199 #else 
    200                         temdta(ji,jj,:,:) = temdta(jl,jj,:,:) 
    201 #endif 
    202                      END DO 
    203                   END DO 
    204                END DO 
    205                ! 
    206             ELSE 
    207                !                                         ! Reduced temperature at Red Sea 
    208                ij0 =  87   ;   ij1 =  96 
    209                ii0 = 148   ;   ii1 = 160 
    210 #if defined key_orca_lev10 
    211                ztem  ( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ,  4:10 , : ) = 7.0  
    212                ztem  ( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 11:13 , : ) = 6.5  
    213                ztem  ( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 14:20 , : ) = 6.0 
    214 #else 
    215                temdta( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ,  4:10 , : ) = 7.0  
    216                temdta( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 11:13 , : ) = 6.5  
    217                temdta( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 14:20 , : ) = 6.0 
    218 #endif 
    219             ENDIF 
    220             ! 
    221          ENDIF 
    222 #endif 
    223           
    224 #if defined key_orca_lev10 
    225          ! interpolate from 31 to 301 level the ztem field result in temdta 
    226          DO jl = 1, 2 
    227             DO jjk = 1, 5 
    228                temdta(:,:,jjk,jl) = ztem(:,:,1,jl) 
    229             END DO 
    230             DO jk = 1, jpk-20,10 
    231                ik = jk+5 
    232                ikr =  INT(jk/10) + 1 
    233                ikw =  (ikr-1) *10 + 1 
    234                ikt =  ikw + 5 
    235                DO jjk=ikt,ikt+9 
    236                   zfac = ( gdept_0(jjk   ) - gdepw_0(ikt) ) / ( gdepw_0(ikt+10) - gdepw_0(ikt) ) 
    237                   temdta(:,:,jjk,jl) = ztem(:,:,ikr,jl) + ( ztem(:,:,ikr+1,jl) - ztem(:,:,ikr,jl) ) * zfac 
    238                END DO 
    239             END DO 
    240             DO jjk = jpk-5, jpk 
    241                temdta(:,:,jjk,jl) = ztem(:,:,jpkdta-1,jl) 
    242             END DO 
    243             ! fill the overlap areas 
    244             CALL lbc_lnk (temdta(:,:,:,jl),'Z',-999.,'no0') 
    245          END DO 
    246 #endif 
    247           
    248          IF( ln_sco ) THEN 
    249             DO jl = 1, 2 
    250                DO jj = 1, jpj                  ! interpolation of temperatures 
    251                   DO ji = 1, jpi 
    252                      DO jk = 1, jpk 
    253                         zl=fsdept_0(ji,jj,jk) 
    254                         IF(zl < gdept_0(1)) ztemdta(jk,jl) =  temdta(ji,jj,1,jl) 
    255                         IF(zl > gdept_0(jpk)) ztemdta(jk,jl) =  temdta(ji,jj,jpkm1,jl) 
    256                         DO jkk = 1, jpkm1 
    257                            IF((zl-gdept_0(jkk))*(zl-gdept_0(jkk+1)).le.0.0) THEN 
    258                               ztemdta(jk,jl) = temdta(ji,jj,jkk,jl)                                 & 
    259                                    &           + (zl-gdept_0(jkk))/(gdept_0(jkk+1)-gdept_0(jkk))      & 
    260                                    &                              *(temdta(ji,jj,jkk+1,jl) - temdta(ji,jj,jkk,jl)) 
    261                            ENDIF 
    262                         END DO 
    263                      END DO 
    264                      DO jk = 1, jpkm1 
    265                         temdta(ji,jj,jk,jl) = ztemdta(jk,jl) 
    266                      END DO 
    267                      temdta(ji,jj,jpk,jl) = 0.0 
    268                   END DO 
    269                END DO 
    270             END DO 
    271              
    272             IF(lwp) WRITE(numout,*) 
    273             IF(lwp) WRITE(numout,*) ' Levitus temperature data interpolated to s-coordinate' 
    274             IF(lwp) WRITE(numout,*) 
    275              
    276          ELSE 
    277              
    278             !                                  ! Mask 
    279             DO jl = 1, 2 
    280                temdta(:,:,:,jl) = temdta(:,:,:,jl) * tmask(:,:,:) 
    281                temdta(:,:,jpk,jl) = 0. 
    282                IF( ln_zps ) THEN                ! z-coord. with partial steps 
    283                   DO jj = 1, jpj                  ! interpolation of temperature at the last level 
    284                      DO ji = 1, jpi 
    285                         ik = mbathy(ji,jj) - 1 
    286                         IF( ik > 2 ) THEN 
    287                            zl = ( gdept_0(ik) - fsdept_0(ji,jj,ik) ) / ( gdept_0(ik) - gdept_0(ik-1) ) 
    288                            temdta(ji,jj,ik,jl) = (1.-zl) * temdta(ji,jj,ik,jl) + zl * temdta(ji,jj,ik-1,jl) 
    289                         ENDIF 
    290                      END DO 
    291                   END DO 
    292                ENDIF 
    293             END DO 
    294              
    295          ENDIF 
    296           
    297          IF(lwp) THEN 
    298             WRITE(numout,*) ' temperature Levitus month ', ntem1, ntem2 
     236               DO jk = 1, jpkm1 
     237                  t_dta(ji,jj,jk) = ztemdta(jk) 
     238               END DO 
     239               t_dta(ji,jj,jpk) = 0.0 
     240            END DO 
     241         END DO 
     242             
     243         IF( lwp .AND. kt==nn_it000 )THEN 
    299244            WRITE(numout,*) 
    300             WRITE(numout,*) ' Levitus month = ', ntem1, '  level = 1' 
    301             CALL prihre( temdta(:,:,1,1), jpi, jpj, 1, jpi, 20, 1, jpj, 20, 1., numout ) 
    302             WRITE(numout,*) ' Levitus month = ', ntem1, '  level = ', jpk/2 
    303             CALL prihre( temdta(:,:,jpk/2,1), jpi, jpj, 1, jpi, 20, 1, jpj, 20, 1., numout ) 
    304             WRITE(numout,*) ' Levitus month = ',ntem1,'  level = ', jpkm1 
    305             CALL prihre( temdta(:,:,jpkm1,1), jpi, jpj, 1, jpi, 20, 1, jpj, 20, 1., numout ) 
    306          ENDIF 
    307       ENDIF 
    308        
    309        
    310       ! 2. At every time step compute temperature data 
    311       ! ---------------------------------------------- 
    312        
    313       zxy = FLOAT( nday + 15 - 30 * i15 ) / 30. 
    314       t_dta(:,:,:) = (1.-zxy) * temdta(:,:,:,1) + zxy * temdta(:,:,:,2) 
    315        
    316       ! Close the file 
    317       ! -------------- 
    318        
    319       IF( kt == nitend )   CALL iom_close (numtdt) 
    320        
    321     END SUBROUTINE dta_tem 
     245            WRITE(numout,*) ' Levitus temperature data interpolated to s-coordinate' 
     246            WRITE(numout,*) 
     247         ENDIF 
     248             
     249      ELSE 
     250         !                                  ! Mask 
     251         t_dta(:,:,:  ) = t_dta(:,:,:) * tmask(:,:,:) 
     252         t_dta(:,:,jpk) = 0. 
     253         IF( ln_zps ) THEN                ! z-coord. with partial steps 
     254            DO jj = 1, jpj                ! interpolation of temperature at the last level 
     255               DO ji = 1, jpi 
     256                  ik = mbathy(ji,jj) - 1 
     257                  IF( ik > 2 ) THEN 
     258                     zl = ( gdept_0(ik) - fsdept_0(ji,jj,ik) ) / ( gdept_0(ik) - gdept_0(ik-1) ) 
     259                     t_dta(ji,jj,ik) = (1.-zl) * t_dta(ji,jj,ik) + zl * t_dta(ji,jj,ik-1) 
     260                  ENDIF 
     261            END DO 
     262         END DO 
     263      ENDIF 
     264 
     265   ENDIF 
     266          
     267   IF( lwp .AND. kt==nn_it000 ) THEN 
     268      WRITE(numout,*) ' temperature Levitus ' 
     269      WRITE(numout,*) 
     270      WRITE(numout,*)'  level = 1' 
     271      CALL prihre( t_dta(:,:,1    ), jpi, jpj, 1, jpi, 20, 1, jpj, 20, 1., numout ) 
     272      WRITE(numout,*)'  level = ', jpk/2 
     273      CALL prihre( t_dta(:,:,jpk/2), jpi, jpj, 1, jpi, 20, 1, jpj, 20, 1., numout ) 
     274      WRITE(numout,*)'  level = ', jpkm1 
     275      CALL prihre( t_dta(:,:,jpkm1), jpi, jpj, 1, jpi, 20, 1, jpj, 20, 1., numout ) 
     276   ENDIF 
     277 
     278   END SUBROUTINE dta_tem 
    322279 
    323280#else 

  • branches/DEV_r2191_3partymerge2010/NEMO/OPA_SRC/DYN/dynzdf.F90

    r2208 r2222  
    107107      USE zdftke_old 
    108108      USE zdftke 
     109      USE zdfgls 
    109110      USE zdfkpp 
    110111      !!---------------------------------------------------------------------- 
    … …  
    116117 
    117118      ! Force implicit schemes 
    118       IF( lk_zdftke_old .OR. lk_zdftke .OR. lk_zdfkpp )   nzdf = 1   ! TKE or KPP physics 
    119       IF( ln_dynldf_iso                               )   nzdf = 1   ! iso-neutral lateral physics 
    120       IF( ln_dynldf_hor .AND. ln_sco                  )   nzdf = 1   ! horizontal lateral physics in s-coordinate 
     119      IF( lk_zdftke_old .OR. lk_zdftke .OR. lk_zdfgls .OR. lk_zdfkpp )   nzdf = 1   ! TKE, GLS or KPP physics 
     120      IF( ln_dynldf_iso                                              )   nzdf = 1   ! iso-neutral lateral physics 
     121      IF( ln_dynldf_hor .AND. ln_sco                                 )   nzdf = 1   ! horizontal lateral physics in s-coordinate 
    121122 
    122123      IF( lk_esopa )    nzdf = -1                   ! Esopa key: All schemes used 

  • branches/DEV_r2191_3partymerge2010/NEMO/OPA_SRC/DYN/dynzdf_imp.F90

    r2208 r2222  
    2424   USE phycst          ! physical constants 
    2525   USE in_out_manager  ! I/O manager 
     26#if defined key_zdfgls 
     27   USE zdfbfr, ONLY : bfrua, bfrva, wbotu, wbotv ! bottom stresses 
     28   USE lbclnk         ! ocean lateral boundary conditions (or mpp link) 
     29#endif 
    2630 
    2731   IMPLICIT NONE 
    … …  
    7579      REAL(wp) ::   zzwi                               ! temporary scalars 
    7680      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwi        ! temporary workspace arrays 
     81#if defined key_zdfgls 
     82      INTEGER :: ikbu, ikbv, ikbum1, ikbvm1 
     83      REAL(wp) :: zcbcu, zcbcv 
     84#endif 
    7785      !!---------------------------------------------------------------------- 
    7886 
    … …  
    173181         END DO 
    174182      END DO 
     183 
     184#if defined key_zdfgls 
     185      ! Save bottom stress for next time step 
     186      DO jj = 2, jpjm1 
     187         DO ji = fs_2, fs_jpim1   ! vector opt. 
     188            ikbu   = MIN( mbathy(ji+1,jj  ), mbathy(ji,jj) ) 
     189            ikbum1 = MAX( ikbu-1, 1 ) 
     190            wbotu(ji,jj) = bfrua(ji,jj) * ua(ji,jj,ikbum1) * umask(ji,jj,ikbum1) 
     191         END DO 
     192      END DO 
     193      CALL lbc_lnk( wbotu(:,:), 'U', -1. ) 
     194#endif 
    175195 
    176196      ! Normalization to obtain the general momentum trend ua 
    … …  
    272292      END DO 
    273293 
     294#if defined key_zdfgls 
     295      ! Save bottom stress for next time step 
     296      DO jj = 2, jpjm1 
     297         DO ji = fs_2, fs_jpim1   ! vector opt. 
     298            ikbv   = MIN( mbathy(ji,jj+1), mbathy(ji,jj) ) 
     299            ikbvm1 = MAX( ikbv-1, 1 ) 
     300            wbotv(ji,jj) = bfrva(ji,jj) * va(ji,jj,ikbvm1) * vmask(ji,jj,ikbvm1) 
     301         END DO 
     302      END DO 
     303      CALL lbc_lnk( wbotv(:,:), 'V', -1. ) 
     304#endif 
     305 
    274306      ! Normalization to obtain the general momentum trend va 
    275307      DO jk = 1, jpkm1 

  • branches/DEV_r2191_3partymerge2010/NEMO/OPA_SRC/IOM/restart.F90

    r2208 r2222  
    2525   USE zdfmxl          ! mixed layer depth 
    2626   USE trdmld_oce      ! ocean active mixed layer tracers trends variables 
     27#if defined key_zdfgls 
     28   USE zdfbfr, ONLY : wbotu, wbotv ! bottom stresses 
     29   USE zdf_oce 
     30#endif 
    2731 
    2832   IMPLICIT NONE 
    … …  
    135139#endif 
    136140 
     141#if defined key_zdfgls 
     142      ! Save bottom stresses 
     143      CALL iom_rstput( kt, nitrst, numrow, 'wbotu' , wbotu ) 
     144      CALL iom_rstput( kt, nitrst, numrow, 'wbotv' , wbotv ) 
     145#endif 
     146 
    137147      IF( kt == nitrst ) THEN 
    138148         CALL iom_close( numrow )     ! close the restart file (only at last time step) 

  • branches/DEV_r2191_3partymerge2010/NEMO/OPA_SRC/OBC/obcini.F90

    r2208 r2222  
    427427               obcsurftot = obcsurftot+hu(ji,jj)*e2u(ji,jj)*uemsk(jj,1) * MAX(obctmsk(ji,jj),obctmsk(ji+1,jj) ) 
    428428            END DO 
    429 <<<<<<< .working 
    430429         END DO 
    431430      END IF 
    … …  
    435434            DO jj = 1, jpj  
    436435               obcsurftot = obcsurftot+hu(ji,jj)*e2u(ji,jj)*uwmsk(jj,1) * MAX(obctmsk(ji,jj),obctmsk(ji+1,jj) ) 
    437 ======= 
    438          END DO 
    439       END IF 
     436            END DO 
     437         END DO 
     438      END IF 
     439 
    440440      IF( lp_obc_west ) THEN ! ... West open boundary lateral surface 
    441441         DO ji = niw0, niw1 
    442442            DO jj = 1, jpj  
    443443               obcsurftot = obcsurftot+hu(ji,jj)*e2u(ji,jj)*uwmsk(jj,1) * MAX(obctmsk(ji,jj),obctmsk(ji+1,jj) ) 
    444 >>>>>>> .merge-right.r2207 
    445             END DO 
    446          END DO 
    447       END IF 
     444            END DO 
     445         END DO 
     446      END IF 
     447 
    448448      IF( lp_obc_north ) THEN ! ... North open boundary lateral surface 
    449449         DO jj = njn0, njn1 

  • branches/DEV_r2191_3partymerge2010/NEMO/OPA_SRC/SBC/fldread.F90

    r2208 r2222  
    1515   USE oce             ! ocean dynamics and tracers 
    1616   USE dom_oce         ! ocean space and time domain 
     17   USE ioipsl, ONLY :   ymds2ju, ju2ymds   ! for calendar 
    1718   USE phycst          ! ??? 
    1819   USE in_out_manager  ! I/O manager 
    … …  
    2930      LOGICAL              ::   ln_tint     ! time interpolation or not (T/F) 
    3031      LOGICAL              ::   ln_clim     ! climatology or not (T/F) 
    31       CHARACTER(len = 7)   ::   cltype      ! type of data file 'daily', 'monthly' or yearly' 
     32      CHARACTER(len = 8)   ::   cltype      ! type of data file 'daily', 'monthly' or yearly' 
    3233      CHARACTER(len = 34)  ::   wname       ! generic name of a NetCDF weights file to be used, blank if not 
    3334      CHARACTER(len = 34)  ::   vcomp       ! symbolic component name if a vector that needs rotation 
    … …  
    4344      LOGICAL                         ::   ln_tint      ! time interpolation or not (T/F) 
    4445      LOGICAL                         ::   ln_clim      ! climatology or not (T/F) 
    45       CHARACTER(len = 7)              ::   cltype       ! type of data file 'daily', 'monthly' or yearly' 
     46      CHARACTER(len = 8)              ::   cltype       ! type of data file 'daily', 'monthly' or yearly' 
    4647      INTEGER                         ::   num          ! iom id of the jpfld files to be read 
    4748      INTEGER                         ::   nswap_sec    ! swapping time in second since Jan. 1st 00h of nit000 year 
    4849      INTEGER , DIMENSION(2)          ::   nrec_b       ! before record (1: index, 2: second since Jan. 1st 00h of nit000 year) 
    4950      INTEGER , DIMENSION(2)          ::   nrec_a       ! after  record (1: index, 2: second since Jan. 1st 00h of nit000 year) 
    50       REAL(wp) , ALLOCATABLE, DIMENSION(:,:)   ::   fnow         ! input fields interpolated to now time step 
    51       REAL(wp) , ALLOCATABLE, DIMENSION(:,:,:) ::   fdta         ! 2 consecutive record of input fields 
     51      REAL(wp) , ALLOCATABLE, DIMENSION(:,:,:)   ::   fnow       ! input fields interpolated to now time step 
     52      REAL(wp) , ALLOCATABLE, DIMENSION(:,:,:,:) ::   fdta       ! 2 consecutive record of input fields 
    5253      CHARACTER(len = 256)            ::   wgtname      ! current name of the NetCDF weight file acting as a key 
    5354                                                        ! into the WGTLIST structure 
    … …  
    7879      INTEGER, DIMENSION(:,:,:), POINTER      ::   data_jpj     ! array of source integers 
    7980      REAL(wp), DIMENSION(:,:,:), POINTER     ::   data_wgt     ! array of weights on model grid 
    80       REAL(wp), DIMENSION(:,:), POINTER       ::   fly_dta      ! array of values on input grid 
    81       REAL(wp), DIMENSION(:,:), POINTER       ::   col2         ! temporary array for reading in columns 
     81      REAL(wp), DIMENSION(:,:,:), POINTER       ::   fly_dta      ! array of values on input grid 
     82      REAL(wp), DIMENSION(:,:,:), POINTER       ::   col2         ! temporary array for reading in columns 
    8283   END TYPE WGT 
    8384 
    … …  
    120121 
    121122      INTEGER  ::   jf         ! dummy indices 
     123      INTEGER  ::   jk         ! dummy indices 
     124      INTEGER  ::   ipk        ! number of vertical levels of sdjf%fdta ( 2D: ipk=1 ; 3D: ipk=jpk ) 
    122125      INTEGER  ::   kw         ! index into wgts array 
    123126      INTEGER  ::   ireclast   ! last record to be read in the current year file 
    … …  
    143146            IF( sd(jf)%ln_tint ) THEN         ! time interpolation: swap before record field 
    144147!CDIR COLLAPSE 
    145                sd(jf)%fdta(:,:,1) = sd(jf)%fdta(:,:,2) 
     148               sd(jf)%fdta(:,:,:,1) = sd(jf)%fdta(:,:,:,2) 
    146149               sd(jf)%rotn(1)     = sd(jf)%rotn(2) 
    147150            ENDIF 
    … …  
    157160 
    158161               ! last record to be read in the current file 
    159                IF( sd(jf)%nfreqh == -1 ) THEN                  ;   ireclast = 12 
     162               IF( sd(jf)%nfreqh == -1 ) THEN 
     163                  IF(     sd(jf)%cltype == 'monthly'   ) THEN  ;   ireclast = 1 
     164                  ELSE                                         ;   ireclast = 12 
     165                  ENDIF 
    160166               ELSE                              
    161                   IF(     sd(jf)%cltype == 'monthly'   ) THEN  ;   ireclast = 24 * nmonth_len(nmonth) / sd(jf)%nfreqh  
    162                   ELSEIF( sd(jf)%cltype == 'daily'     ) THEN  ;   ireclast = 24                      / sd(jf)%nfreqh 
    163                   ELSE                                         ;   ireclast = 24 * nyear_len(     1 ) / sd(jf)%nfreqh  
     167                  IF(     sd(jf)%cltype      == 'monthly' ) THEN  ;   ireclast = 24 * nmonth_len(nmonth) / sd(jf)%nfreqh  
     168                  ELSEIF( sd(jf)%cltype(1:4) == 'week'    ) THEN  ;   ireclast = 24.* 7                  / sd(jf)%nfreqh 
     169                  ELSEIF( sd(jf)%cltype      == 'daily'   ) THEN  ;   ireclast = 24                      / sd(jf)%nfreqh 
     170                  ELSE                                            ;   ireclast = 24 * nyear_len(     1 ) / sd(jf)%nfreqh  
    164171                  ENDIF 
    165172               ENDIF 
    … …  
    204211            IF( LEN(TRIM(sd(jf)%wgtname)) > 0 ) THEN 
    205212               CALL wgt_list( sd(jf), kw ) 
    206                CALL fld_interp( sd(jf)%num, sd(jf)%clvar, kw, sd(jf)%fdta(:,:,2), sd(jf)%nrec_a(1) ) 
     213               ipk = SIZE(sd(jf)%fnow,3) 
     214               IF( sd(jf)%ln_tint ) THEN 
     215                  CALL fld_interp( sd(jf)%num, sd(jf)%clvar , kw , ipk, sd(jf)%fdta(:,:,:,2) , sd(jf)%nrec_a(1) ) 
     216               ELSE 
     217                  CALL fld_interp( sd(jf)%num, sd(jf)%clvar , kw , ipk, sd(jf)%fnow(:,:,:)   , sd(jf)%nrec_a(1) ) 
     218               ENDIF 
    207219            ELSE 
    208                CALL iom_get( sd(jf)%num, jpdom_data, sd(jf)%clvar, sd(jf)%fdta(:,:,2), sd(jf)%nrec_a(1) ) 
     220               SELECT CASE( SIZE(sd(jf)%fnow,3) ) 
     221               CASE(1)    
     222                  IF( sd(jf)%ln_tint ) THEN 
     223                     CALL iom_get( sd(jf)%num, jpdom_data, sd(jf)%clvar, sd(jf)%fdta(:,:,1,2), sd(jf)%nrec_a(1) ) 
     224                  ELSE 
     225                     CALL iom_get( sd(jf)%num, jpdom_data, sd(jf)%clvar, sd(jf)%fnow(:,:,1)  , sd(jf)%nrec_a(1) ) 
     226                  ENDIF  
     227               CASE(jpk) 
     228                  IF( sd(jf)%ln_tint ) THEN 
     229                     CALL iom_get( sd(jf)%num, jpdom_data, sd(jf)%clvar, sd(jf)%fdta(:,:,:,2), sd(jf)%nrec_a(1) ) 
     230                  ELSE 
     231                     CALL iom_get( sd(jf)%num, jpdom_data, sd(jf)%clvar, sd(jf)%fnow(:,:,:)  , sd(jf)%nrec_a(1) ) 
     232                  ENDIF  
     233               END SELECT 
    209234            ENDIF 
    210235            sd(jf)%rotn(2) = .FALSE. 
    … …  
    240265                IF( kf > 0 ) THEN 
    241266                   !! fields jf,kf are two components which need to be rotated together 
    242                    DO nf = 1,2 
     267                   IF( sd(jf)%ln_tint )THEN 
     268                      DO nf = 1,2 
     269                         !! check each time level of this pair 
     270                         IF( .NOT. sd(jf)%rotn(nf) .AND. .NOT. sd(kf)%rotn(nf) ) THEN 
     271                            utmp(:,:) = 0.0 
     272                            vtmp(:,:) = 0.0 
     273                            ! 
     274                            ipk = SIZE( sd(kf)%fnow(:,:,:) ,3 ) 
     275                            DO jk = 1,ipk 
     276                               CALL rot_rep( sd(jf)%fdta(:,:,jk,nf),sd(kf)%fdta(:,:,jk,nf),'T', 'en->i', utmp(:,:) ) 
     277                               CALL rot_rep( sd(jf)%fdta(:,:,jk,nf),sd(kf)%fdta(:,:,jk,nf),'T', 'en->j', vtmp(:,:) ) 
     278                               sd(jf)%fdta(:,:,jk,nf) = utmp(:,:) 
     279                               sd(kf)%fdta(:,:,jk,nf) = vtmp(:,:) 
     280                            ENDDO 
     281                            ! 
     282                            sd(jf)%rotn(nf) = .TRUE. 
     283                            sd(kf)%rotn(nf) = .TRUE. 
     284                            IF( lwp .AND. kt == nit000 ) & 
     285                                      WRITE(numout,*) 'fld_read: vector pair (',  & 
     286                                                      TRIM(sd(jf)%clvar),',',TRIM(sd(kf)%clvar), & 
     287                                                      ') rotated on to model grid' 
     288                         ENDIF 
     289                      END DO 
     290                   ELSE  
    243291                      !! check each time level of this pair 
    244292                      IF( .NOT. sd(jf)%rotn(nf) .AND. .NOT. sd(kf)%rotn(nf) ) THEN 
    245293                         utmp(:,:) = 0.0 
    246294                         vtmp(:,:) = 0.0 
    247                          CALL rot_rep( sd(jf)%fdta(:,:,nf), sd(kf)%fdta(:,:,nf), 'T', 'en->i', utmp(:,:) ) 
    248                          CALL rot_rep( sd(jf)%fdta(:,:,nf), sd(kf)%fdta(:,:,nf), 'T', 'en->j', vtmp(:,:) ) 
    249                          sd(jf)%fdta(:,:,nf) = utmp(:,:) 
    250                          sd(kf)%fdta(:,:,nf) = vtmp(:,:) 
     295                         ! 
     296                         ipk = SIZE( sd(kf)%fnow(:,:,:) ,3 ) 
     297                         DO jk = 1,ipk 
     298                            CALL rot_rep( sd(jf)%fnow(:,:,jk),sd(kf)%fnow(:,:,jk),'T', 'en->i', utmp(:,:) ) 
     299                            CALL rot_rep( sd(jf)%fnow(:,:,jk),sd(kf)%fnow(:,:,jk),'T', 'en->j', vtmp(:,:) ) 
     300                            sd(jf)%fnow(:,:,jk) = utmp(:,:) 
     301                            sd(kf)%fnow(:,:,jk) = vtmp(:,:) 
     302                         ENDDO 
     303                         ! 
    251304                         sd(jf)%rotn(nf) = .TRUE. 
    252305                         sd(kf)%rotn(nf) = .TRUE. 
    … …  
    256309                                                   ') rotated on to model grid' 
    257310                      ENDIF 
    258                    END DO 
     311                   ENDIF 
    259312                ENDIF 
    260313             ENDIF 
    … …  
    280333               ztintb =  1. - ztinta 
    281334!CDIR COLLAPSE 
    282                sd(jf)%fnow(:,:) = ztintb * sd(jf)%fdta(:,:,1) + ztinta * sd(jf)%fdta(:,:,2) 
     335               sd(jf)%fnow(:,:,:) = ztintb * sd(jf)%fdta(:,:,:,1) + ztinta * sd(jf)%fdta(:,:,:,2) 
    283336            ELSE 
    284337               IF(lwp .AND. kt - nit000 <= 100 ) THEN 
    … …  
    288341               ENDIF 
    289342!CDIR COLLAPSE 
    290                sd(jf)%fnow(:,:) = sd(jf)%fdta(:,:,2)   ! piecewise constant field 
    291   
    292343            ENDIF 
    293344            ! 
    … …  
    313364      TYPE(FLD), INTENT(inout) ::   sdjf        ! input field related variables 
    314365      !! 
    315       LOGICAL :: llprevyr       ! are we reading previous year  file? 
    316       LOGICAL :: llprevmth      ! are we reading previous month file? 
    317       LOGICAL :: llprevday      ! are we reading previous day   file? 
    318       LOGICAL :: llprev         ! llprevyr .OR. llprevmth .OR. llprevday 
    319       INTEGER :: idvar          ! variable id  
    320       INTEGER :: inrec          ! number of record existing for this variable 
     366      LOGICAL :: llprevyr              ! are we reading previous year  file? 
     367      LOGICAL :: llprevmth             ! are we reading previous month file? 
     368      LOGICAL :: llprevweek            ! are we reading previous week file? 
     369      LOGICAL :: llprevday             ! are we reading previous day   file? 
     370      LOGICAL :: llprev                ! llprevyr .OR. llprevmth .OR. llprevday 
     371      INTEGER :: idvar                 ! variable id  
     372      INTEGER :: inrec                 ! number of record existing for this variable 
    321373      INTEGER :: kwgt 
     374      INTEGER :: jk             !vertical loop variable 
     375      INTEGER :: ipk            !number of vertical levels of sdjf%fdta ( 2D: ipk=1 ; 3D: ipk=jpk ) 
     376      INTEGER :: iyear, imonth, iday   ! first day of the current file in yyyy mm dd 
     377      INTEGER :: isec_week             ! number of seconds since start of the weekly file 
    322378      CHARACTER(LEN=1000) ::   clfmt   ! write format 
    323379      !!--------------------------------------------------------------------- 
    324  
     380       
    325381      ! some default definitions... 
    326382      sdjf%num = 0   ! default definition for non-opened file 
    327383      IF( sdjf%ln_clim )   sdjf%clname = TRIM( sdjf%clrootname )   ! file name defaut definition, never change in this case 
    328       llprevyr  = .FALSE. 
    329       llprevmth = .FALSE. 
    330       llprevday = .FALSE. 
     384      llprevyr   = .FALSE. 
     385      llprevmth  = .FALSE. 
     386      llprevweek = .FALSE. 
     387      llprevday  = .FALSE. 
     388      isec_week  = 0 
    331389             
    332390      ! define record informations 
    … …  
    339397               IF( sdjf%cltype == 'monthly' ) THEN   ! monthly file 
    340398                  sdjf%nrec_b(1) = 1                                                       ! force to read the unique record 
    341                   llprevmth = .NOT. sdjf%ln_clim                                           ! use previous month file? 
     399                  llprevmth = .TRUE.                                                       ! use previous month file? 
    342400                  llprevyr  = llprevmth .AND. nmonth == 1                                  ! use previous year  file? 
    343401               ELSE                                  ! yearly file 
    … …  
    350408                  llprevmth = .NOT. sdjf%ln_clim                                           ! use previous month file? 
    351409                  llprevyr  = llprevmth .AND. nmonth == 1                                  ! use previous year  file? 
     410               ELSE IF ( sdjf%cltype(1:4) == 'week' ) THEN !weekly file 
     411                  isec_week = 86400 * 7 
     412                  sdjf%nrec_b(1) = 24. / sdjf%nfreqh * 7                                   ! last record of previous weekly file 
    352413               ELSEIF( sdjf%cltype == 'daily' ) THEN ! daily file 
    353414                  sdjf%nrec_b(1) = 24 / sdjf%nfreqh                                        ! last record of previous day 
    … …  
    361422            ENDIF 
    362423         ENDIF 
    363          llprev = llprevyr .OR. llprevmth .OR. llprevday 
     424         llprev = llprevyr .OR. llprevmth .OR. llprevweek .OR. llprevday 
    364425 
    365426         CALL fld_clopn( sdjf, nyear  - COUNT((/llprevyr /))                                              ,               & 
    366427            &                  nmonth - COUNT((/llprevmth/)) + 12                   * COUNT((/llprevyr /)),               & 
    367428            &                  nday   - COUNT((/llprevday/)) + nmonth_len(nmonth-1) * COUNT((/llprevmth/)), .NOT. llprev ) 
    368           
     429 
     430         IF ( sdjf%cltype(1:4) == 'week' ) THEN 
     431            isec_week  = ksec_week( sdjf%cltype(6:8) ) 
     432            if(lwp)write(numout,*)'cbr test2 isec_week = ',isec_week 
     433            llprevmth  = ( isec_week .GT. nsec_month ) 
     434            llprevyr   = llprevmth  .AND. nmonth==1 
     435         ENDIF 
     436         ! 
     437         iyear  = nyear  - COUNT((/llprevyr /)) 
     438         imonth = nmonth - COUNT((/llprevmth/)) + 12* COUNT((/llprevyr /)) 
     439         iday   = nday   - COUNT((/llprevday/)) + nmonth_len(nmonth-1) * COUNT((/llprevmth/)) - INT( isec_week )/86400 
     440         ! 
     441         CALL fld_clopn( sdjf , iyear , imonth , iday , .NOT. llprev ) 
     442 
    369443         ! if previous year/month/day file does not exist, we switch to the current year/month/day 
    370444         IF( llprev .AND. sdjf%num <= 0 ) THEN 
    … …  
    384458 
    385459         ! read before data into sdjf%fdta(:,:,2) because we will swap data in the following part of fld_read 
     460          
    386461         IF( LEN(TRIM(sdjf%wgtname)) > 0 ) THEN 
    387462            CALL wgt_list( sdjf, kwgt ) 
    388             CALL fld_interp( sdjf%num, sdjf%clvar, kwgt, sdjf%fdta(:,:,2), sdjf%nrec_b(1) ) 
     463            ipk = SIZE(sdjf%fnow,3) 
     464            IF( sdjf%ln_tint ) THEN 
     465               CALL fld_interp( sdjf%num, sdjf%clvar, kwgt, ipk, sdjf%fdta(:,:,:,2), sdjf%nrec_a(1) ) 
     466            ELSE 
     467               CALL fld_interp( sdjf%num, sdjf%clvar, kwgt, ipk, sdjf%fnow(:,:,:)  , sdjf%nrec_a(1) ) 
     468            ENDIF 
    389469         ELSE 
    390             CALL iom_get( sdjf%num, jpdom_data, sdjf%clvar, sdjf%fdta(:,:,2), sdjf%nrec_b(1) ) 
     470            SELECT CASE( SIZE(sdjf%fnow,3) ) 
     471            CASE(1) 
     472               IF( sdjf%ln_tint ) THEN 
     473                  CALL iom_get( sdjf%num, jpdom_data, sdjf%clvar, sdjf%fdta(:,:,1,2), sdjf%nrec_b(1) ) 
     474               ELSE 
     475                  CALL iom_get( sdjf%num, jpdom_data, sdjf%clvar, sdjf%fnow(:,:,1)  , sdjf%nrec_b(1) ) 
     476               ENDIF 
     477            CASE(jpk) 
     478               IF( sdjf%ln_tint ) THEN 
     479                  CALL iom_get( sdjf%num, jpdom_data, sdjf%clvar, sdjf%fdta(:,:,:,2), sdjf%nrec_b(1) ) 
     480               ELSE 
     481                  CALL iom_get( sdjf%num, jpdom_data, sdjf%clvar, sdjf%fnow(:,:,:)  , sdjf%nrec_b(1) ) 
     482               ENDIF 
     483            END SELECT 
    391484         ENDIF 
    392485         sdjf%rotn(2) = .FALSE. 
    … …  
    400493 
    401494      IF( sdjf%num <= 0 )   CALL fld_clopn( sdjf, nyear, nmonth, nday )   ! make sure current year/month/day file is opened 
     495      ! make sure current year/month/day file is opened 
     496      IF( sdjf%num == 0 ) THEN 
     497         isec_week   = 0 
     498         llprevyr    = .FALSE. 
     499         llprevmth   = .FALSE. 
     500         llprevweek  = .FALSE. 
     501         ! 
     502         IF ( sdjf%cltype(1:4) == 'week' ) THEN 
     503            isec_week  = ksec_week( sdjf%cltype(6:8) ) 
     504            llprevmth  = ( isec_week .GT. nsec_month ) 
     505            llprevyr   = llprevmth  .AND. nmonth==1 
     506         ENDIF 
     507         ! 
     508         iyear  = nyear  - COUNT((/llprevyr /)) 
     509         imonth = nmonth - COUNT((/llprevmth/)) + 12* COUNT((/llprevyr /)) 
     510         iday   = nday   + nmonth_len(nmonth-1) * COUNT((/llprevmth/)) - isec_week/86400 
     511         ! 
     512         CALL fld_clopn( sdjf, iyear, imonth, iday ) 
     513      ENDIF  
    402514 
    403515      sdjf%nswap_sec = nsec_year + nsec1jan000 - 1   ! force read/update the after data in the following part of fld_read  
    404        
     516      
     517 
    405518   END SUBROUTINE fld_init 
    406519 
    … …  
    420533      REAL(wp) ::   ztmp        ! temporary variable 
    421534      INTEGER  ::   ifreq_sec   ! frequency mean (in seconds) 
     535      INTEGER  ::   isec_week   ! number of seconds since the start of the weekly file 
    422536      !!---------------------------------------------------------------------- 
    423537      ! 
    … …  
    436550            !       forcing record :  nmonth  
    437551            !                             
     552            ztmp  = 0.e0 
    438553            ztmp  = REAL( nday, wp ) / REAL( nmonth_len(nmonth), wp ) + 0.5 
    439554         ELSE 
    … …  
    446561         ENDIF 
    447562 
    448          sdjf%nrec_a(:) = (/ irec, nmonth_half(irec) + nsec1jan000 /)   ! define after  record number and time 
    449          irec = irec - 1                                                ! move back to previous record 
    450          sdjf%nrec_b(:) = (/ irec, nmonth_half(irec) + nsec1jan000 /)   ! define before record number and time 
     563         IF( sdjf%cltype == 'monthly' ) THEN 
     564 
     565            sdjf%nrec_b(:) = (/ 0, nmonth_half(irec - 1 ) + nsec1jan000 /) 
     566            sdjf%nrec_a(:) = (/ 1, nmonth_half(irec     ) + nsec1jan000 /) 
     567 
     568            IF( ztmp  == 1. ) THEN 
     569              sdjf%nrec_b(1) = 1 
     570              sdjf%nrec_a(1) = 2 
     571            ENDIF 
     572 
     573         ELSE 
     574 
     575            sdjf%nrec_a(:) = (/ irec, nmonth_half(irec) + nsec1jan000 /)   ! define after  record number and time 
     576            irec = irec - 1                                                ! move back to previous record 
     577            sdjf%nrec_b(:) = (/ irec, nmonth_half(irec) + nsec1jan000 /)   ! define before record number and time 
     578 
     579         ENDIF 
    451580         ! 
    452581      ELSE                              ! higher frequency mean (in hours) 
    453582         ! 
    454583         ifreq_sec = sdjf%nfreqh * 3600   ! frequency mean (in seconds) 
     584         IF( sdjf%cltype(1:4) == 'week'    ) isec_week = ksec_week( sdjf%cltype(6:8)) !since the first day of the current week 
    455585         ! number of second since the beginning of the file 
    456          IF(     sdjf%cltype == 'monthly' ) THEN   ;   ztmp = REAL(nsec_month,wp)   ! since 00h on the 1st day of the current month 
    457          ELSEIF( sdjf%cltype == 'daily'   ) THEN   ;   ztmp = REAL(nsec_day  ,wp)   ! since 00h of the current day 
    458          ELSE                                      ;   ztmp = REAL(nsec_year ,wp)   ! since 00h on Jan 1 of the current year 
     586         IF(     sdjf%cltype      == 'monthly' ) THEN   ;   ztmp = REAL(nsec_month ,wp)  ! since 00h on the 1st day of the current month 
     587         ELSEIF( sdjf%cltype(1:4) == 'week'    ) THEN   ;   ztmp = REAL(isec_week  ,wp)  ! since the first day of the current week 
     588         ELSEIF( sdjf%cltype      == 'daily'   ) THEN   ;   ztmp = REAL(nsec_day   ,wp)  ! since 00h of the current day 
     589         ELSE                                           ;   ztmp = REAL(nsec_year  ,wp)  ! since 00h on Jan 1 of the current year 
    459590         ENDIF 
    460591         IF( sdjf%ln_tint ) THEN                ! time interpolation, shift by 1/2 record 
    … …  
    492623         ! after record index and second since Jan. 1st 00h of nit000 year 
    493624         sdjf%nrec_a(:) = (/ irec, ifreq_sec * irec - ifreq_sec / 2 + nsec1jan000 /) 
    494          IF( sdjf%cltype == 'monthly' )   &   ! add the number of seconds between 00h Jan 1 and the end of previous month 
     625         IF( sdjf%cltype == 'monthly' )       &   ! add the number of seconds between 00h Jan 1 and the end of previous month 
    495626            sdjf%nrec_a(2) = sdjf%nrec_a(2) + isecd * SUM(nmonth_len(1:nmonth -1))   ! ok if nmonth=1 
    496          IF( sdjf%cltype == 'daily'   )   &   ! add the number of seconds between 00h Jan 1 and the end of previous day 
     627         IF( sdjf%cltype(1:4) == 'week'   )   &   ! add the number of seconds between 00h Jan 1 and the end of previous week  
     628            sdjf%nrec_a(2) = sdjf%nrec_a(2) + ( nsec_year - isec_week ) 
     629         IF( sdjf%cltype == 'daily'   )       &   ! add the number of seconds between 00h Jan 1 and the end of previous day 
    497630            sdjf%nrec_a(2) = sdjf%nrec_a(2) + isecd * ( nday_year - 1 ) 
    498631 
    … …  
    500633         irec = irec - 1.                           ! move back to previous record 
    501634         sdjf%nrec_b(:) = (/ irec, ifreq_sec * irec - ifreq_sec / 2 + nsec1jan000 /) 
    502          IF( sdjf%cltype == 'monthly' )   &   ! add the number of seconds between 00h Jan 1 and the end of previous month 
     635         IF( sdjf%cltype == 'monthly' )       &   ! add the number of seconds between 00h Jan 1 and the end of previous month 
    503636            sdjf%nrec_b(2) = sdjf%nrec_b(2) + isecd * SUM(nmonth_len(1:nmonth -1))   ! ok if nmonth=1 
    504          IF( sdjf%cltype == 'daily'   )   &   ! add the number of seconds between 00h Jan 1 and the end of previous day 
     637         IF( sdjf%cltype(1:4) == 'week'   )   &   ! add the number of seconds between 00h Jan 1 and the end of previous week 
     638            sdjf%nrec_b(2) = sdjf%nrec_b(2) + ( nsec_year - isec_week ) 
     639         IF( sdjf%cltype == 'daily'   )       &   ! add the number of seconds between 00h Jan 1 and the end of previous day 
    505640            sdjf%nrec_b(2) = sdjf%nrec_b(2) + isecd * ( nday_year - 1 ) 
    506641 
    … …  
    523658      !! ** Method  :    
    524659      !!---------------------------------------------------------------------- 
    525       TYPE(FLD), INTENT(inout)           ::   sdjf     ! input field related variables 
    526       INTEGER  , INTENT(in   )           ::   kyear    ! year value 
    527       INTEGER  , INTENT(in   )           ::   kmonth   ! month value 
    528       INTEGER  , INTENT(in   )           ::   kday     ! day value 
    529       LOGICAL  , INTENT(in   ), OPTIONAL ::   ldstop   ! stop if open to read a non-existing file (default = .TRUE.) 
     660      TYPE(FLD), INTENT(inout)           ::   sdjf                      ! input field related variables 
     661      INTEGER  , INTENT(in   )           ::   kyear                     ! year value 
     662      INTEGER  , INTENT(in   )           ::   kmonth                    ! month value 
     663      INTEGER  , INTENT(in   )           ::   kday                      ! day value 
     664      LOGICAL  , INTENT(in   ), OPTIONAL ::   ldstop                    ! stop if open to read a non-existing file (default = .TRUE.) 
     665      INTEGER                            ::   iyear, imonth, iday       ! firt day of the current week in yyyy mm dd 
     666      REAL(wp)                           ::   zsec, zjul                !temp variable 
    530667 
    531668      IF( sdjf%num /= 0 )   CALL iom_close( sdjf%num )   ! close file if already open 
    532669      ! build the new filename if not climatological data 
    533       IF( .NOT. sdjf%ln_clim ) THEN   ;   WRITE(sdjf%clname, '(a,"_y",i4.4)' ) TRIM( sdjf%clrootname ), kyear    ! add year 
    534          IF( sdjf%cltype /= 'yearly' )    WRITE(sdjf%clname, '(a,"m" ,i2.2)' ) TRIM( sdjf%clname     ), kmonth   ! add month 
    535          IF( sdjf%cltype == 'daily'  )    WRITE(sdjf%clname, '(a,"d" ,i2.2)' ) TRIM( sdjf%clname     ), kday     ! add day 
     670      sdjf%clname=TRIM(sdjf%clrootname) 
     671      ! 
     672      IF(  sdjf%cltype(1:4) == 'week' .AND. nn_leapy==0 )CALL ctl_stop( 'fld_clopn: weekly file and nn_leapy=0 are not compatible' ) 
     673      ! 
     674      IF( .NOT. sdjf%ln_clim ) THEN    
     675         WRITE(sdjf%clname, '(a,"_y",i4.4)' ) TRIM( sdjf%clrootname ), kyear    ! add year 
     676         IF( sdjf%cltype /= 'yearly'        )   &  
     677            &     WRITE(sdjf%clname, '(a,"m" ,i2.2)' ) TRIM( sdjf%clname ), kmonth   ! add month 
     678         IF( sdjf%cltype == 'daily'  .OR. sdjf%cltype(1:4) == 'week' ) & 
     679            &     WRITE(sdjf%clname, '(a,"d" ,i2.2)' ) TRIM( sdjf%clname ), kday     ! add day 
     680      ELSE 
     681         ! build the new filename if climatological data 
     682         IF( sdjf%cltype == 'monthly' )   WRITE(sdjf%clname, '(a,"_m" ,i2.2)' ) TRIM( sdjf%clrootname ), kmonth   ! add month 
    536683      ENDIF 
    537684      CALL iom_open( sdjf%clname, sdjf%num, ldstop = ldstop, ldiof =  LEN(TRIM(sdjf%wgtname)) > 0 ) 
    … …  
    564711         sdf(jf)%ln_tint    = sdf_n(jf)%ln_tint 
    565712         sdf(jf)%ln_clim    = sdf_n(jf)%ln_clim 
    566          IF( sdf(jf)%nfreqh == -1. ) THEN   ;   sdf(jf)%cltype = 'yearly' 
    567          ELSE                               ;   sdf(jf)%cltype = sdf_n(jf)%cltype 
    568          ENDIF 
     713         sdf(jf)%cltype     = sdf_n(jf)%cltype 
    569714         sdf(jf)%wgtname = " " 
    570715         IF( LEN( TRIM(sdf_n(jf)%wname) ) > 0 )   sdf(jf)%wgtname = TRIM( cdir )//TRIM( sdf_n(jf)%wname ) 
    … …  
    587732               &                          ' pairing    : '    , TRIM( sdf(jf)%vcomp      ),   & 
    588733               &                          ' data type: '      ,       sdf(jf)%cltype 
     734            call flush(numout) 
    589735         END DO 
    590736      ENDIF 
    … …  
    684830      INTEGER                                 ::   inum          ! temporary logical unit 
    685831      INTEGER                                 ::   id            ! temporary variable id 
     832      INTEGER                                 ::   ipk           ! temporary vertical dimension 
    686833      CHARACTER (len=5)                       ::   aname 
    687834      INTEGER , DIMENSION(3)                  ::   ddims 
    … …  
    846993         ! SA: +3 stencil is a patch to avoid out-of-bound computation in some configuration.  
    847994         ! a more robust solution will be given in next release 
    848          ALLOCATE( ref_wgts(nxt_wgt)%fly_dta(ref_wgts(nxt_wgt)%jpiwgt+3, ref_wgts(nxt_wgt)%jpjwgt+3) ) 
    849          IF( ref_wgts(nxt_wgt)%cyclic ) ALLOCATE( ref_wgts(nxt_wgt)%col2(2,ref_wgts(nxt_wgt)%jpjwgt+3) ) 
     995         ipk =  SIZE(sd%fnow,3) 
     996         ALLOCATE( ref_wgts(nxt_wgt)%fly_dta(ref_wgts(nxt_wgt)%jpiwgt+3, ref_wgts(nxt_wgt)%jpjwgt+3 ,ipk) ) 
     997         IF( ref_wgts(nxt_wgt)%cyclic ) ALLOCATE( ref_wgts(nxt_wgt)%col2(2,ref_wgts(nxt_wgt)%jpjwgt+3,ipk) ) 
    850998 
    851999         nxt_wgt = nxt_wgt + 1 
    … …  
    8571005   END SUBROUTINE fld_weight 
    8581006 
    859    SUBROUTINE fld_interp(num, clvar, kw, dta, nrec) 
     1007   SUBROUTINE fld_interp(num, clvar, kw, kk, dta, nrec) 
    8601008      !!--------------------------------------------------------------------- 
    8611009      !!                    ***  ROUTINE fld_interp  *** 
    … …  
    8691017      CHARACTER(LEN=*), INTENT(in)                        ::   clvar               ! variable name 
    8701018      INTEGER,          INTENT(in)                        ::   kw                  ! weights number 
    871       REAL(wp),         INTENT(inout), DIMENSION(jpi,jpj) ::   dta                 ! output field on model grid 
     1019      INTEGER,          INTENT(in)                        ::   kk                  ! vertical dimension of kk 
     1020      REAL(wp),         INTENT(inout), DIMENSION(jpi,jpj,kk) ::   dta              ! output field on model grid 
    8721021      INTEGER,          INTENT(in)                        ::   nrec                ! record number to read (ie time slice) 
    8731022      !!  
    874       INTEGER, DIMENSION(2)                               ::   rec1,recn           ! temporary arrays for start and length 
     1023      INTEGER, DIMENSION(3)                               ::   rec1,recn           ! temporary arrays for start and length 
    8751024      INTEGER                                             ::  jk, jn, jm           ! loop counters 
    8761025      INTEGER                                             ::  ni, nj               ! lengths 
    … …  
    8951044      rec1(1) = MAX( jpimin-1, 1 ) 
    8961045      rec1(2) = MAX( jpjmin-1, 1 ) 
     1046      rec1(3) = 1 
    8971047      recn(1) = MIN( jpiwid+2, ref_wgts(kw)%ddims(1)-rec1(1)+1 ) 
    8981048      recn(2) = MIN( jpjwid+2, ref_wgts(kw)%ddims(2)-rec1(2)+1 ) 
     1049      recn(3) = kk 
    8991050 
    9001051      !! where we need to read it to 
    … …  
    9041055      jpj2 = jpj1 + recn(2) - 1 
    9051056 
    906       ref_wgts(kw)%fly_dta(:,:) = 0.0 
    907       CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%fly_dta(jpi1:jpi2,jpj1:jpj2), nrec, rec1, recn) 
     1057      ref_wgts(kw)%fly_dta(:,:,:) = 0.0 
     1058      SELECT CASE( SIZE(ref_wgts(kw)%fly_dta(jpi1:jpi2,jpj1:jpj2,:),3) ) 
     1059      CASE(1) 
     1060           CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%fly_dta(jpi1:jpi2,jpj1:jpj2,1), nrec, rec1, recn) 
     1061      CASE(jpk)   
     1062           CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%fly_dta(jpi1:jpi2,jpj1:jpj2,:), nrec, rec1, recn) 
     1063      END SELECT  
    9081064 
    9091065      !! first four weights common to both bilinear and bicubic 
    9101066      !! note that we have to offset by 1 into fly_dta array because of halo 
    911       dta(:,:) = 0.0 
     1067      dta(:,:,:) = 0.0 
    9121068      DO jk = 1,4 
    913         DO jn = 1, jpj 
    914           DO jm = 1,jpi 
     1069        DO jn = 1, nlcj 
     1070          DO jm = 1,nlci 
    9151071            ni = ref_wgts(kw)%data_jpi(jm,jn,jk) 
    9161072            nj = ref_wgts(kw)%data_jpj(jm,jn,jk) 
    917             dta(jm,jn) = dta(jm,jn) + ref_wgts(kw)%data_wgt(jm,jn,jk) * ref_wgts(kw)%fly_dta(ni+1,nj+1) 
     1073            dta(jm,jn,:) = dta(jm,jn,:) + ref_wgts(kw)%data_wgt(jm,jn,jk) * ref_wgts(kw)%fly_dta(ni+1,nj+1,jk) 
    9181074          END DO 
    9191075        END DO 
    … …  
    9241080        !! fix up halo points that we couldnt read from file 
    9251081        IF( jpi1 == 2 ) THEN 
    926            ref_wgts(kw)%fly_dta(jpi1-1,:) = ref_wgts(kw)%fly_dta(jpi1,:) 
     1082           ref_wgts(kw)%fly_dta(jpi1-1,:,:) = ref_wgts(kw)%fly_dta(jpi1,:,:) 
    9271083        ENDIF 
    9281084        IF( jpi2 + jpimin - 1 == ref_wgts(kw)%ddims(1)+1 ) THEN 
    929            ref_wgts(kw)%fly_dta(jpi2+1,:) = ref_wgts(kw)%fly_dta(jpi2,:) 
     1085           ref_wgts(kw)%fly_dta(jpi2+1,:,:) = ref_wgts(kw)%fly_dta(jpi2,:,:) 
    9301086        ENDIF 
    9311087        IF( jpj1 == 2 ) THEN 
    932            ref_wgts(kw)%fly_dta(:,jpj1-1) = ref_wgts(kw)%fly_dta(:,jpj1) 
     1088           ref_wgts(kw)%fly_dta(:,jpj1-1,:) = ref_wgts(kw)%fly_dta(:,jpj1,:) 
    9331089        ENDIF 
    9341090        IF( jpj2 + jpjmin - 1 == ref_wgts(kw)%ddims(2)+1 .AND. jpj2 .lt. jpjwid+2 ) THEN 
    935            ref_wgts(kw)%fly_dta(:,jpj2+1) = 2.0*ref_wgts(kw)%fly_dta(:,jpj2) - ref_wgts(kw)%fly_dta(:,jpj2-1) 
     1091           ref_wgts(kw)%fly_dta(:,jpj2+1,:) = 2.0*ref_wgts(kw)%fly_dta(:,jpj2,:) - ref_wgts(kw)%fly_dta(:,jpj2-1,:) 
    9361092        ENDIF 
    9371093 
    … …  
    9461102           IF( jpi1 == 2 ) THEN 
    9471103              rec1(1) = ref_wgts(kw)%ddims(1) - 1 
    948               CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%col2(:,jpj1:jpj2), nrec, rec1, recn) 
    949               ref_wgts(kw)%fly_dta(jpi1-1,jpj1:jpj2) = ref_wgts(kw)%col2(ref_wgts(kw)%offset+1,jpj1:jpj2) 
     1104              SELECT CASE( SIZE( ref_wgts(kw)%col2(:,jpj1:jpj2,:),3) ) 
     1105              CASE(1) 
     1106                   CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%col2(:,jpj1:jpj2,1), nrec, rec1, recn) 
     1107              CASE(jpk)          
     1108                   CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%col2(:,jpj1:jpj2,:), nrec, rec1, recn) 
     1109              END SELECT       
     1110              ref_wgts(kw)%fly_dta(jpi1-1,jpj1:jpj2,:) = ref_wgts(kw)%col2(ref_wgts(kw)%offset+1,jpj1:jpj2,:) 
    9501111           ENDIF 
    9511112           IF( jpi2 + jpimin - 1 == ref_wgts(kw)%ddims(1)+1 ) THEN 
    9521113              rec1(1) = 1 
    953               CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%col2(:,jpj1:jpj2), nrec, rec1, recn) 
    954               ref_wgts(kw)%fly_dta(jpi2+1,jpj1:jpj2) = ref_wgts(kw)%col2(2-ref_wgts(kw)%offset,jpj1:jpj2) 
     1114              SELECT CASE( SIZE( ref_wgts(kw)%col2(:,jpj1:jpj2,:),3) ) 
     1115              CASE(1) 
     1116                   CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%col2(:,jpj1:jpj2,1), nrec, rec1, recn) 
     1117              CASE(jpk) 
     1118                   CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%col2(:,jpj1:jpj2,:), nrec, rec1, recn) 
     1119              END SELECT 
     1120              ref_wgts(kw)%fly_dta(jpi2+1,jpj1:jpj2,:) = ref_wgts(kw)%col2(2-ref_wgts(kw)%offset,jpj1:jpj2,:) 
    9551121           ENDIF 
    9561122        ENDIF 
    … …  
    9581124        ! gradient in the i direction 
    9591125        DO jk = 1,4 
    960           DO jn = 1, jpj 
    961             DO jm = 1,jpi 
     1126          DO jn = 1, nlcj 
     1127            DO jm = 1,nlci 
    9621128              ni = ref_wgts(kw)%data_jpi(jm,jn,jk) 
    9631129              nj = ref_wgts(kw)%data_jpj(jm,jn,jk) 
    964               dta(jm,jn) = dta(jm,jn) + ref_wgts(kw)%data_wgt(jm,jn,jk+4) * 0.5 *         & 
    965                                (ref_wgts(kw)%fly_dta(ni+2,nj+1) - ref_wgts(kw)%fly_dta(ni,nj+1)) 
     1130              dta(jm,jn,:) = dta(jm,jn,:) + ref_wgts(kw)%data_wgt(jm,jn,jk+4) * 0.5 *         & 
     1131                               (ref_wgts(kw)%fly_dta(ni+2,nj+1,:) - ref_wgts(kw)%fly_dta(ni,nj+1,:)) 
    9661132            END DO 
    9671133          END DO 
    … …  
    9701136        ! gradient in the j direction 
    9711137        DO jk = 1,4 
    972           DO jn = 1, jpj 
    973             DO jm = 1,jpi 
     1138          DO jn = 1, nlcj 
     1139            DO jm = 1,nlci 
    9741140              ni = ref_wgts(kw)%data_jpi(jm,jn,jk) 
    9751141              nj = ref_wgts(kw)%data_jpj(jm,jn,jk) 
    976               dta(jm,jn) = dta(jm,jn) + ref_wgts(kw)%data_wgt(jm,jn,jk+8) * 0.5 *         & 
    977                                (ref_wgts(kw)%fly_dta(ni+1,nj+2) - ref_wgts(kw)%fly_dta(ni+1,nj)) 
     1142              dta(jm,jn,:) = dta(jm,jn,:) + ref_wgts(kw)%data_wgt(jm,jn,jk+8) * 0.5 *         & 
     1143                               (ref_wgts(kw)%fly_dta(ni+1,nj+2,:) - ref_wgts(kw)%fly_dta(ni+1,nj,:)) 
    9781144            END DO 
    9791145          END DO 
    … …  
    9861152              ni = ref_wgts(kw)%data_jpi(jm,jn,jk) 
    9871153              nj = ref_wgts(kw)%data_jpj(jm,jn,jk) 
    988               dta(jm,jn) = dta(jm,jn) + ref_wgts(kw)%data_wgt(jm,jn,jk+12) * 0.25 * ( & 
    989                                (ref_wgts(kw)%fly_dta(ni+2,nj+2) - ref_wgts(kw)%fly_dta(ni  ,nj+2)) -   & 
    990                                (ref_wgts(kw)%fly_dta(ni+2,nj  ) - ref_wgts(kw)%fly_dta(ni  ,nj  ))) 
     1154              dta(jm,jn,:) = dta(jm,jn,:) + ref_wgts(kw)%data_wgt(jm,jn,jk+12) * 0.25 * ( & 
     1155                               (ref_wgts(kw)%fly_dta(ni+2,nj+2,:) - ref_wgts(kw)%fly_dta(ni  ,nj+2,:)) -   & 
     1156                               (ref_wgts(kw)%fly_dta(ni+2,nj  ,:) - ref_wgts(kw)%fly_dta(ni  ,nj  ,:))) 
    9911157            END DO 
    9921158          END DO 
    … …  
    9961162 
    9971163   END SUBROUTINE fld_interp 
    998    
     1164 
     1165   FUNCTION ksec_week( cdday ) 
     1166      !!--------------------------------------------------------------------- 
     1167      !!                    ***  FUNCTION kshift_week ***  
     1168      !! 
     1169      !! ** Purpose :   
     1170      !! 
     1171      !! ** Method  : 
     1172      !!--------------------------------------------------------------------- 
     1173      CHARACTER(len=*), INTENT(in)   ::   cdday   !3 first letters of the first day of the weekly file 
     1174      !! 
     1175      INTEGER                        ::   ksec_week  ! output variable 
     1176      INTEGER                        ::   ijul       !temp variable 
     1177      INTEGER                        ::   ishift     !temp variable 
     1178      CHARACTER(len=3),DIMENSION(7)  ::   cl_week  
     1179      !!---------------------------------------------------------------------- 
     1180      cl_week = (/"sun","sat","fri","thu","wed","tue","mon"/) 
     1181      DO ijul=1,7 
     1182         IF(  cl_week(ijul)==TRIM(cdday) ) EXIT 
     1183      ENDDO 
     1184      IF( ijul .GT. 7 ) CALL ctl_stop( 'ksec_week: wrong day for sdjf%cltype(6:8): ',TRIM(cdday) ) 
     1185      ! 
     1186      ishift = ( ijul  ) * 86400 
     1187      !  
     1188      ksec_week = nsec_week + ishift 
     1189      ksec_week = MOD( ksec_week , 86400*7 ) 
     1190      if(lwp)write(numout,*)'cbr ijul ksec_week ',ijul,ksec_week 
     1191      !  
     1192   END FUNCTION ksec_week 
     1193 
    9991194END MODULE fldread 

  • branches/DEV_r2191_3partymerge2010/NEMO/OPA_SRC/SBC/sbcblk_clio.F90

    r2208 r2222  
    160160            CALL ctl_stop( 'sbc_blk_clio: unable to allocate sf structure' )   ;   RETURN 
    161161         ENDIF 
    162  
    163162         DO ifpr= 1, jpfld 
    164             ALLOCATE( sf(ifpr)%fnow(jpi,jpj) ) 
    165             ALLOCATE( sf(ifpr)%fdta(jpi,jpj,2) ) 
    166          END DO 
    167  
    168  
     163            ALLOCATE( sf(ifpr)%fnow(jpi,jpj,1) ) 
     164            IF( slf_i(ifpr)%ln_tint ) ALLOCATE( sf(ifpr)%fdta(jpi,jpj,1,2) ) 
     165         END DO 
    169166         ! fill sf with slf_i and control print 
    170167         CALL fld_fill( sf, slf_i, cn_dir, 'sbc_blk_clio', 'flux formulation for ocean surface boundary condition', 'namsbc_clio' ) 
    … …  
    178175      ! 
    179176#if defined key_lim3       
    180       tatm_ice(:,:) = sf(jp_tair)%fnow(:,:)     !RB ugly patch 
     177      tatm_ice(:,:) = sf(jp_tair)%fnow(:,:,1)     !RB ugly patch 
    181178#endif 
    182179      ! 
    … …  
    272269      DO jj = 1 , jpj 
    273270         DO ji = 1, jpi 
    274             utau(ji,jj) = sf(jp_utau)%fnow(ji,jj) 
    275             vtau(ji,jj) = sf(jp_vtau)%fnow(ji,jj) 
     271            utau(ji,jj) = sf(jp_utau)%fnow(ji,jj,1) 
     272            vtau(ji,jj) = sf(jp_vtau)%fnow(ji,jj,1) 
    276273         END DO 
    277274      END DO 
    … …  
    297294      DO jj = 1 , jpj 
    298295         DO ji = 1, jpi 
    299             wndm(ji,jj) = sf(jp_wndm)%fnow(ji,jj) 
     296            wndm(ji,jj) = sf(jp_wndm)%fnow(ji,jj,1) 
    300297         END DO 
    301298      END DO 
    … …  
    317314            ! 
    318315            zsst  = pst(ji,jj)              + rt0           ! converte Celcius to Kelvin the SST 
    319             ztatm = sf(jp_tair)%fnow(ji,jj)                 ! and set minimum value far above 0 K (=rt0 over land) 
    320             zcco1 = 1.0 - sf(jp_ccov)%fnow(ji,jj)           ! fraction of clear sky ( 1 - cloud cover) 
     316            ztatm = sf(jp_tair)%fnow(ji,jj,1)               ! and set minimum value far above 0 K (=rt0 over land) 
     317            zcco1 = 1.0 - sf(jp_ccov)%fnow(ji,jj,1)         ! fraction of clear sky ( 1 - cloud cover) 
    321318            zrhoa = zpatm / ( 287.04 * ztatm )              ! air density (equation of state for dry air)  
    322319            ztamr = ztatm - rtt                             ! Saturation water vapour 
    … …  
    325322            zmt3  = SIGN( 28.200, -ztamr )                  !           \/ 
    326323            zes   = 611.0 * EXP(  ABS( ztamr ) * MIN ( zmt1, zmt2 ) / ( ztatm - 35.86  + MAX( 0.e0, zmt3 ) )  ) 
    327             zev    = sf(jp_humi)%fnow(ji,jj) * zes          ! vapour pressure   
     324            zev    = sf(jp_humi)%fnow(ji,jj,1) * zes        ! vapour pressure   
    328325            zevsqr = SQRT( zev * 0.01 )                     ! square-root of vapour pressure 
    329326            zqatm = 0.622 * zev / ( zpatm - 0.378 * zev )   ! specific humidity  
    … …  
    333330            !--------------------------------------! 
    334331            ztatm3  = ztatm * ztatm * ztatm 
    335             zcldeff = 1.0 - sbudyko(ji,jj) * sf(jp_ccov)%fnow(ji,jj) * sf(jp_ccov)%fnow(ji,jj)     
     332            zcldeff = 1.0 - sbudyko(ji,jj) * sf(jp_ccov)%fnow(ji,jj,1) * sf(jp_ccov)%fnow(ji,jj,1)     
    336333            ztaevbk = ztatm * ztatm3 * zcldeff * ( 0.39 - 0.05 * zevsqr )  
    337334            ! 
    … …  
    351348            zdeltaq = zqatm - zqsato 
    352349            ztvmoy  = ztatm * ( 1. + 2.2e-3 * ztatm * zqatm ) 
    353             zdenum  = MAX( sf(jp_wndm)%fnow(ji,jj) * sf(jp_wndm)%fnow(ji,jj) * ztvmoy, zeps ) 
     350            zdenum  = MAX( sf(jp_wndm)%fnow(ji,jj,1) * sf(jp_wndm)%fnow(ji,jj,1) * ztvmoy, zeps ) 
    354351            zdtetar = zdteta / zdenum 
    355352            ztvmoyr = ztvmoy * ztvmoy * zdeltaq / zdenum 
    … …  
    373370            zpsil   = zpsih 
    374371             
    375             zvatmg         = MAX( 0.032 * 1.5e-3 * sf(jp_wndm)%fnow(ji,jj) * sf(jp_wndm)%fnow(ji,jj) / grav, zeps ) 
     372            zvatmg         = MAX( 0.032 * 1.5e-3 * sf(jp_wndm)%fnow(ji,jj,1) * sf(jp_wndm)%fnow(ji,jj,1) / grav, zeps ) 
    376373            zcmn           = vkarmn / LOG ( 10. / zvatmg ) 
    377374            zchn           = 0.0327 * zcmn 
    … …  
    387384            zcleo          = zcln * zclcm  
    388385 
    389             zrhova         = zrhoa * sf(jp_wndm)%fnow(ji,jj) 
     386            zrhova         = zrhoa * sf(jp_wndm)%fnow(ji,jj,1) 
    390387 
    391388            ! sensible heat flux 
    … …  
    408405         DO ji = 1, jpi 
    409406            qns (ji,jj) = zqlw(ji,jj) - zqsb(ji,jj) - zqla(ji,jj)      ! Downward Non Solar flux 
    410             emp (ji,jj) = zqla(ji,jj) / cevap - sf(jp_prec)%fnow(ji,jj) / rday * tmask(ji,jj,1) 
     407            emp (ji,jj) = zqla(ji,jj) / cevap - sf(jp_prec)%fnow(ji,jj,1) / rday * tmask(ji,jj,1) 
    411408         END DO 
    412409      END DO 
    … …  
    530527!CDIR NOVERRCHK 
    531528         DO ji = 1, jpi 
    532             ztatm (ji,jj) = sf(jp_tair)%fnow(ji,jj)                  ! air temperature in Kelvins  
     529            ztatm (ji,jj) = sf(jp_tair)%fnow(ji,jj,1)                ! air temperature in Kelvins  
    533530       
    534531            zrhoa(ji,jj) = zpatm / ( 287.04 * ztatm(ji,jj) )         ! air density (equation of state for dry air)  
    … …  
    541538               &                / ( ztatm(ji,jj) - 35.86  + MAX( 0.e0, zmt3 ) )  ) 
    542539 
    543             zev = sf(jp_humi)%fnow(ji,jj) * zes                      ! vapour pressure   
     540            zev = sf(jp_humi)%fnow(ji,jj,1) * zes                      ! vapour pressure   
    544541            zevsqr(ji,jj) = SQRT( zev * 0.01 )                       ! square-root of vapour pressure 
    545542            zqatm(ji,jj) = 0.622 * zev / ( zpatm - 0.378 * zev )     ! specific humidity  
    … …  
    551548            zmt2  = ( 272.0 - ztatm(ji,jj) ) / 38.0   ;   zind2 = MAX( 0.e0, SIGN( 1.e0, zmt2 ) ) 
    552549            zmt3  = ( 281.0 - ztatm(ji,jj) ) / 18.0   ;   zind3 = MAX( 0.e0, SIGN( 1.e0, zmt3 ) ) 
    553             p_spr(ji,jj) = sf(jp_prec)%fnow(ji,jj) / rday   &        ! rday = converte mm/day to kg/m2/s 
     550            p_spr(ji,jj) = sf(jp_prec)%fnow(ji,jj,1) / rday   &      ! rday = converte mm/day to kg/m2/s 
    554551               &         * (          zind1      &                   ! solid  (snow) precipitation [kg/m2/s] 
    555552               &            + ( 1.0 - zind1 ) * (          zind2   * ( 0.5 + zmt2 )   & 
    … …  
    561558            ! fraction of qsr_ice which is NOT absorbed in the thin surface layer 
    562559            ! and thus which penetrates inside the ice cover ( Maykut and Untersteiner, 1971 ; Elbert anbd Curry, 1993 ) 
    563             p_fr1(ji,jj) = 0.18  * ( 1.e0 - sf(jp_ccov)%fnow(ji,jj) ) + 0.35 * sf(jp_ccov)%fnow(ji,jj)  
    564             p_fr2(ji,jj) = 0.82  * ( 1.e0 - sf(jp_ccov)%fnow(ji,jj) ) + 0.65 * sf(jp_ccov)%fnow(ji,jj) 
     560            p_fr1(ji,jj) = 0.18  * ( 1.e0 - sf(jp_ccov)%fnow(ji,jj,1) ) + 0.35 * sf(jp_ccov)%fnow(ji,jj,1)  
     561            p_fr2(ji,jj) = 0.82  * ( 1.e0 - sf(jp_ccov)%fnow(ji,jj,1) ) + 0.65 * sf(jp_ccov)%fnow(ji,jj,1) 
    565562         END DO 
    566563      END DO 
    … …  
    584581               !-------------------------------------------! 
    585582               ztatm3  = ztatm(ji,jj) * ztatm(ji,jj) * ztatm(ji,jj) 
    586                zcldeff = 1.0 - sbudyko(ji,jj) * sf(jp_ccov)%fnow(ji,jj) * sf(jp_ccov)%fnow(ji,jj)     
     583               zcldeff = 1.0 - sbudyko(ji,jj) * sf(jp_ccov)%fnow(ji,jj,1) * sf(jp_ccov)%fnow(ji,jj,1)     
    587584               ztaevbk = ztatm3 * ztatm(ji,jj) * zcldeff * ( 0.39 - 0.05 * zevsqr(ji,jj) )  
    588585               ! 
    … …  
    609606                
    610607               !  sensible and latent fluxes over ice 
    611                zrhova     = zrhoa(ji,jj) * sf(jp_wndm)%fnow(ji,jj)      ! computation of intermediate values 
     608               zrhova     = zrhoa(ji,jj) * sf(jp_wndm)%fnow(ji,jj,1)      ! computation of intermediate values 
    612609               zrhovaclei = zrhova * zcshi * 2.834e+06 
    613610               zrhovacshi = zrhova * zclei * 1004.0 
    … …  
    639636      p_qns(:,:,:) = z_qlw (:,:,:) - z_qsb (:,:,:) - p_qla (:,:,:)      ! Downward Non Solar flux 
    640637!CDIR COLLAPSE 
    641       p_tpr(:,:)   = sf(jp_prec)%fnow(:,:) / rday                       ! total precipitation [kg/m2/s] 
     638      p_tpr(:,:)   = sf(jp_prec)%fnow(:,:,1) / rday                       ! total precipitation [kg/m2/s] 
    642639      ! 
    643640!!gm : not necessary as all input data are lbc_lnk... 
    … …  
    735732!CDIR NOVERRCHK 
    736733         DO ji = 1, jpi 
    737             ztamr = sf(jp_tair)%fnow(ji,jj) - rtt 
     734            ztamr = sf(jp_tair)%fnow(ji,jj,1) - rtt 
    738735            zmt1  = SIGN( 17.269,  ztamr ) 
    739736            zmt2  = SIGN( 21.875,  ztamr ) 
    740737            zmt3  = SIGN( 28.200, -ztamr ) 
    741738            zes = 611.0 * EXP(  ABS( ztamr ) * MIN ( zmt1, zmt2 )   &              ! Saturation water vapour 
    742                &                     / ( sf(jp_tair)%fnow(ji,jj) - 35.86  + MAX( 0.e0, zmt3 ) )  ) 
    743             zev(ji,jj) = sf(jp_humi)%fnow(ji,jj) * zes * 1.0e-05                   ! vapour pressure   
     739               &                     / ( sf(jp_tair)%fnow(ji,jj,1) - 35.86  + MAX( 0.e0, zmt3 ) )  ) 
     740            zev(ji,jj) = sf(jp_humi)%fnow(ji,jj,1) * zes * 1.0e-05                   ! vapour pressure   
    744741         END DO 
    745742      END DO 
    … …  
    798795 
    799796               ! ocean albedo depending on the cloud cover (Payne, 1972) 
    800                za_oce     = ( 1.0 - sf(jp_ccov)%fnow(ji,jj) ) * 0.05 / ( 1.1 * zcmue**1.4 + 0.15 )   &   ! clear sky 
    801                   &       +         sf(jp_ccov)%fnow(ji,jj)   * 0.06                                     ! overcast 
     797               za_oce     = ( 1.0 - sf(jp_ccov)%fnow(ji,jj,1) ) * 0.05 / ( 1.1 * zcmue**1.4 + 0.15 )   &   ! clear sky 
     798                  &       +         sf(jp_ccov)%fnow(ji,jj,1)   * 0.06                                     ! overcast 
    802799 
    803800                  ! solar heat flux absorbed by the ocean (Zillman, 1972) 
    … …  
    814811         DO ji = 1, jpi 
    815812            zlmunoon = ASIN( zps(ji,jj) + zpc(ji,jj) ) / rad                         ! local noon solar altitude 
    816             zcldcor  = MIN(  1.e0, ( 1.e0 - 0.62 * sf(jp_ccov)%fnow(ji,jj)   &       ! cloud correction (Reed 1977) 
     813            zcldcor  = MIN(  1.e0, ( 1.e0 - 0.62 * sf(jp_ccov)%fnow(ji,jj,1)   &       ! cloud correction (Reed 1977) 
    817814               &                          + 0.0019 * zlmunoon )                 ) 
    818815            pqsr_oce(ji,jj) = zcoef1 * zcldcor * pqsr_oce(ji,jj) * tmask(ji,jj,1)   ! and zcoef1: ellipsity 
    … …  
    865862!CDIR NOVERRCHK 
    866863         DO ji = 1, jpi            
    867             ztamr = sf(jp_tair)%fnow(ji,jj) - rtt            
     864            ztamr = sf(jp_tair)%fnow(ji,jj,1) - rtt            
    868865            zmt1  = SIGN( 17.269,  ztamr ) 
    869866            zmt2  = SIGN( 21.875,  ztamr ) 
    870867            zmt3  = SIGN( 28.200, -ztamr ) 
    871868            zes = 611.0 * EXP(  ABS( ztamr ) * MIN ( zmt1, zmt2 )   &              ! Saturation water vapour 
    872                &                     / ( sf(jp_tair)%fnow(ji,jj) - 35.86  + MAX( 0.e0, zmt3 ) )  ) 
    873             zev(ji,jj) = sf(jp_humi)%fnow(ji,jj) * zes * 1.0e-05                   ! vapour pressure   
     869               &                     / ( sf(jp_tair)%fnow(ji,jj,1) - 35.86  + MAX( 0.e0, zmt3 ) )  ) 
     870            zev(ji,jj) = sf(jp_humi)%fnow(ji,jj,1) * zes * 1.0e-05                   ! vapour pressure   
    874871         END DO 
    875872      END DO 
    … …  
    938935                     &        / (  1.0 + 0.139  * stauc(ji,jj) * ( 1.0 - 0.9435 * pa_ice_os(ji,jj,jl) ) )        
    939936              
    940                   pqsr_ice(ji,jj,jl) = pqsr_ice(ji,jj,jl) + (  ( 1.0 - sf(jp_ccov)%fnow(ji,jj) ) * zqsr_ice_cs    & 
    941                      &                                       +         sf(jp_ccov)%fnow(ji,jj)   * zqsr_ice_os  ) 
     937                  pqsr_ice(ji,jj,jl) = pqsr_ice(ji,jj,jl) + (  ( 1.0 - sf(jp_ccov)%fnow(ji,jj,1) ) * zqsr_ice_cs    & 
     938                     &                                       +         sf(jp_ccov)%fnow(ji,jj,1)   * zqsr_ice_os  ) 
    942939               END DO 
    943940            END DO 

  • branches/DEV_r2191_3partymerge2010/NEMO/OPA_SRC/SBC/sbcblk_core.F90

    r2208 r2222  
    164164         ENDIF 
    165165         DO ifpr= 1, jfld 
    166             ALLOCATE( sf(ifpr)%fnow(jpi,jpj) ) 
    167             ALLOCATE( sf(ifpr)%fdta(jpi,jpj,2) ) 
     166            ALLOCATE( sf(ifpr)%fnow(jpi,jpj,1) ) 
     167            IF( slf_i(ifpr)%ln_tint ) ALLOCATE( sf(ifpr)%fdta(jpi,jpj,1,2) ) 
    168168         END DO 
    169169         ! 
    … …  
    176176 
    177177#if defined key_lim3 
    178       tatm_ice(:,:) = sf(jp_tair)%fnow(:,:) 
     178      tatm_ice(:,:) = sf(jp_tair)%fnow(:,:,1) 
    179179#endif 
    180180 
    … …  
    244244      DO jj = 2, jpjm1 
    245245         DO ji = fs_2, fs_jpim1   ! vect. opt. 
    246             zwnd_i(ji,jj) = (  sf(jp_wndi)%fnow(ji,jj) - 0.5 * ( pu(ji-1,jj  ) + pu(ji,jj) )  ) 
    247             zwnd_j(ji,jj) = (  sf(jp_wndj)%fnow(ji,jj) - 0.5 * ( pv(ji  ,jj-1) + pv(ji,jj) )  ) 
     246            zwnd_i(ji,jj) = (  sf(jp_wndi)%fnow(ji,jj,1) - 0.5 * ( pu(ji-1,jj  ) + pu(ji,jj) )  ) 
     247            zwnd_j(ji,jj) = (  sf(jp_wndj)%fnow(ji,jj,1) - 0.5 * ( pv(ji  ,jj-1) + pv(ji,jj) )  ) 
    248248         END DO 
    249249      END DO 
    … …  
    262262      ! ocean albedo assumed to be 0.066 
    263263!CDIR COLLAPSE 
    264       qsr (:,:) = ( 1. - 0.066 ) * sf(jp_qsr)%fnow(:,:) * tmask(:,:,1)                                 ! Short Wave 
    265 !CDIR COLLAPSE 
    266       zqlw(:,:) = (  sf(jp_qlw)%fnow(:,:) - Stef * zst(:,:)*zst(:,:)*zst(:,:)*zst(:,:)  ) * tmask(:,:,1)   ! Long  Wave 
     264      qsr (:,:) = ( 1. - 0.066 ) * sf(jp_qsr)%fnow(:,:,1) * tmask(:,:,1)                                 ! Short Wave 
     265!CDIR COLLAPSE 
     266      zqlw(:,:) = (  sf(jp_qlw)%fnow(:,:,1) - Stef * zst(:,:)*zst(:,:)*zst(:,:)*zst(:,:)  ) * tmask(:,:,1)   ! Long  Wave 
    267267                       
    268268      ! ----------------------------------------------------------------------------- ! 
    … …  
    307307      IF( lhftau ) THEN  
    308308!CDIR COLLAPSE 
    309          taum(:,:) = taum(:,:) + sf(jp_tdif)%fnow(:,:) 
     309         taum(:,:) = taum(:,:) + sf(jp_tdif)%fnow(:,:,1) 
    310310      ENDIF 
    311311      CALL iom_put( "taum_oce", taum )   ! output wind stress module 
    … …  
    330330      ELSE 
    331331!CDIR COLLAPSE 
    332          zevap(:,:) = MAX( 0.e0, rhoa    *Ce(:,:)*( zqsatw(:,:) - sf(jp_humi)%fnow(:,:) ) * wndm(:,:) )   ! Evaporation 
    333 !CDIR COLLAPSE 
    334          zqsb (:,:) =            rhoa*cpa*Ch(:,:)*( zst   (:,:) - sf(jp_tair)%fnow(:,:) ) * wndm(:,:)     ! Sensible Heat 
     332         zevap(:,:) = MAX( 0.e0, rhoa    *Ce(:,:)*( zqsatw(:,:) - sf(jp_humi)%fnow(:,:,1) ) * wndm(:,:) )   ! Evaporation 
     333!CDIR COLLAPSE 
     334         zqsb (:,:) =            rhoa*cpa*Ch(:,:)*( zst   (:,:) - sf(jp_tair)%fnow(:,:,1) ) * wndm(:,:)     ! Sensible Heat 
    335335      ENDIF 
    336336!CDIR COLLAPSE 
    … …  
    355355      qns(:,:) = zqlw(:,:) - zqsb(:,:) - zqla(:,:)      ! Downward Non Solar flux 
    356356!CDIR COLLAPSE 
    357       emp (:,:) = zevap(:,:) - sf(jp_prec)%fnow(:,:) * rn_pfac * tmask(:,:,1) 
     357      emp (:,:) = zevap(:,:) - sf(jp_prec)%fnow(:,:,1) * rn_pfac * tmask(:,:,1) 
    358358!CDIR COLLAPSE 
    359359      emps(:,:) = emp(:,:) 
    … …  
    453453            DO ji = 2, jpim1   ! B grid : no vector opt 
    454454               ! ... scalar wind at I-point (fld being at T-point) 
    455                zwndi_f = 0.25 * (  sf(jp_wndi)%fnow(ji-1,jj  ) + sf(jp_wndi)%fnow(ji  ,jj  )   & 
    456                   &              + sf(jp_wndi)%fnow(ji-1,jj-1) + sf(jp_wndi)%fnow(ji  ,jj-1)  ) - pui(ji,jj) 
    457                zwndj_f = 0.25 * (  sf(jp_wndj)%fnow(ji-1,jj  ) + sf(jp_wndj)%fnow(ji  ,jj  )   & 
    458                   &              + sf(jp_wndj)%fnow(ji-1,jj-1) + sf(jp_wndj)%fnow(ji  ,jj-1)  ) - pvi(ji,jj) 
     455               zwndi_f = 0.25 * (  sf(jp_wndi)%fnow(ji-1,jj  ,1) + sf(jp_wndi)%fnow(ji  ,jj  ,1)   & 
     456                  &              + sf(jp_wndi)%fnow(ji-1,jj-1,1) + sf(jp_wndi)%fnow(ji  ,jj-1,1)  ) - pui(ji,jj) 
     457               zwndj_f = 0.25 * (  sf(jp_wndj)%fnow(ji-1,jj  ,1) + sf(jp_wndj)%fnow(ji  ,jj  ,1)   & 
     458                  &              + sf(jp_wndj)%fnow(ji-1,jj-1,1) + sf(jp_wndj)%fnow(ji  ,jj-1,1)  ) - pvi(ji,jj) 
    459459               zwnorm_f = zcoef_wnorm * SQRT( zwndi_f * zwndi_f + zwndj_f * zwndj_f ) 
    460460               ! ... ice stress at I-point 
    … …  
    462462               p_tauj(ji,jj) = zwnorm_f * zwndj_f 
    463463               ! ... scalar wind at T-point (fld being at T-point) 
    464                zwndi_t = sf(jp_wndi)%fnow(ji,jj) - 0.25 * (  pui(ji,jj+1) + pui(ji+1,jj+1)   & 
     464               zwndi_t = sf(jp_wndi)%fnow(ji,jj,1) - 0.25 * (  pui(ji,jj+1) + pui(ji+1,jj+1)   & 
    465465                  &                                        + pui(ji,jj  ) + pui(ji+1,jj  )  ) 
    466                zwndj_t = sf(jp_wndj)%fnow(ji,jj) - 0.25 * (  pvi(ji,jj+1) + pvi(ji+1,jj+1)   & 
     466               zwndj_t = sf(jp_wndj)%fnow(ji,jj,1) - 0.25 * (  pvi(ji,jj+1) + pvi(ji+1,jj+1)   & 
    467467                  &                                        + pvi(ji,jj  ) + pvi(ji+1,jj  )  ) 
    468468               z_wnds_t(ji,jj)  = SQRT( zwndi_t * zwndi_t + zwndj_t * zwndj_t ) * tmask(ji,jj,1) 
    … …  
    479479         DO jj = 2, jpj 
    480480            DO ji = fs_2, jpi   ! vect. opt. 
    481                zwndi_t = (  sf(jp_wndi)%fnow(ji,jj) - 0.5 * ( pui(ji-1,jj  ) + pui(ji,jj) )  ) 
    482                zwndj_t = (  sf(jp_wndj)%fnow(ji,jj) - 0.5 * ( pvi(ji  ,jj-1) + pvi(ji,jj) )  ) 
     481               zwndi_t = (  sf(jp_wndi)%fnow(ji,jj,1) - 0.5 * ( pui(ji-1,jj  ) + pui(ji,jj) )  ) 
     482               zwndj_t = (  sf(jp_wndj)%fnow(ji,jj,1) - 0.5 * ( pvi(ji  ,jj-1) + pvi(ji,jj) )  ) 
    483483               z_wnds_t(ji,jj)  = SQRT( zwndi_t * zwndi_t + zwndj_t * zwndj_t ) * tmask(ji,jj,1) 
    484484            END DO 
    … …  
    490490            DO ji = fs_2, fs_jpim1   ! vect. opt. 
    491491               p_taui(ji,jj) = zcoef_wnorm2 * ( z_wnds_t(ji+1,jj) + z_wnds_t(ji,jj) )                          & 
    492                   &          * ( 0.5 * (sf(jp_wndi)%fnow(ji+1,jj) + sf(jp_wndi)%fnow(ji,jj) ) - pui(ji,jj) ) 
     492                  &          * ( 0.5 * (sf(jp_wndi)%fnow(ji+1,jj,1) + sf(jp_wndi)%fnow(ji,jj,1) ) - pui(ji,jj) ) 
    493493               p_tauj(ji,jj) = zcoef_wnorm2 * ( z_wnds_t(ji,jj+1) + z_wnds_t(ji,jj) )                          & 
    494                   &          * ( 0.5 * (sf(jp_wndj)%fnow(ji,jj+1) + sf(jp_wndj)%fnow(ji,jj) ) - pvi(ji,jj) ) 
     494                  &          * ( 0.5 * (sf(jp_wndj)%fnow(ji,jj+1,1) + sf(jp_wndj)%fnow(ji,jj,1) ) - pvi(ji,jj) ) 
    495495            END DO 
    496496         END DO 
    … …  
    515515               zst3 = pst(ji,jj,jl) * zst2 
    516516               ! Short Wave (sw) 
    517                p_qsr(ji,jj,jl) = ( 1. - palb(ji,jj,jl) ) * sf(jp_qsr)%fnow(ji,jj) * tmask(ji,jj,1) 
     517               p_qsr(ji,jj,jl) = ( 1. - palb(ji,jj,jl) ) * sf(jp_qsr)%fnow(ji,jj,1) * tmask(ji,jj,1) 
    518518               ! Long  Wave (lw) 
    519                z_qlw(ji,jj,jl) = 0.95 * (  sf(jp_qlw)%fnow(ji,jj)       &                          
     519               z_qlw(ji,jj,jl) = 0.95 * (  sf(jp_qlw)%fnow(ji,jj,1)       &                          
    520520                  &                   - Stef * pst(ji,jj,jl) * zst3  ) * tmask(ji,jj,1) 
    521521               ! lw sensitivity 
    … …  
    528528               ! ... turbulent heat fluxes 
    529529               ! Sensible Heat 
    530                z_qsb(ji,jj,jl) = rhoa * cpa * Cice * z_wnds_t(ji,jj) * ( pst(ji,jj,jl) - sf(jp_tair)%fnow(ji,jj) ) 
     530               z_qsb(ji,jj,jl) = rhoa * cpa * Cice * z_wnds_t(ji,jj) * ( pst(ji,jj,jl) - sf(jp_tair)%fnow(ji,jj,1) ) 
    531531               ! Latent Heat 
    532532               p_qla(ji,jj,jl) = MAX( 0.e0, rhoa * Ls  * Cice * z_wnds_t(ji,jj)   &                            
    533                   &                    * (  11637800. * EXP( -5897.8 / pst(ji,jj,jl) ) / rhoa - sf(jp_humi)%fnow(ji,jj)  ) ) 
     533                  &                    * (  11637800. * EXP( -5897.8 / pst(ji,jj,jl) ) / rhoa - sf(jp_humi)%fnow(ji,jj,1)  ) ) 
    534534               ! Latent heat sensitivity for ice (Dqla/Dt) 
    535535               p_dqla(ji,jj,jl) = zcoef_dqla * z_wnds_t(ji,jj) / ( zst2 ) * EXP( -5897.8 / pst(ji,jj,jl) ) 
    … …  
    561561        
    562562!CDIR COLLAPSE 
    563       p_tpr(:,:) = sf(jp_prec)%fnow(:,:) * rn_pfac      ! total precipitation [kg/m2/s] 
    564 !CDIR COLLAPSE 
    565       p_spr(:,:) = sf(jp_snow)%fnow(:,:) * rn_pfac      ! solid precipitation [kg/m2/s] 
     563      p_tpr(:,:) = sf(jp_prec)%fnow(:,:,1) * rn_pfac      ! total precipitation [kg/m2/s] 
     564!CDIR COLLAPSE 
     565      p_spr(:,:) = sf(jp_snow)%fnow(:,:,1) * rn_pfac      ! solid precipitation [kg/m2/s] 
    566566      CALL iom_put( 'snowpre', p_spr )                  ! Snow precipitation  
    567567      ! 

  • branches/DEV_r2191_3partymerge2010/NEMO/OPA_SRC/SBC/sbcflx.F90

    r2208 r2222  
    126126         ENDIF 
    127127         DO ji= 1, jpfld 
    128             ALLOCATE( sf(ji)%fnow(jpi,jpj) ) 
    129             ALLOCATE( sf(ji)%fdta(jpi,jpj,2) ) 
     128            ALLOCATE( sf(ji)%fnow(jpi,jpj,1) ) 
     129            IF( slf_i(ji)%ln_tint ) ALLOCATE( sf(ji)%fdta(jpi,jpj,1,2) ) 
    130130         END DO 
    131131 
    … …  
    145145         DO jj = 1, jpj 
    146146            DO ji = 1, jpi 
    147                utau(ji,jj) = sf(jp_utau)%fnow(ji,jj) 
    148                vtau(ji,jj) = sf(jp_vtau)%fnow(ji,jj) 
    149                qns (ji,jj) = sf(jp_qtot)%fnow(ji,jj) - sf(jp_qsr)%fnow(ji,jj) 
    150                qsr (ji,jj) = sf(jp_qsr )%fnow(ji,jj) 
    151                emp (ji,jj) = sf(jp_emp )%fnow(ji,jj) 
     147               utau(ji,jj) = sf(jp_utau)%fnow(ji,jj,1) 
     148               vtau(ji,jj) = sf(jp_vtau)%fnow(ji,jj,1) 
     149               qns (ji,jj) = sf(jp_qtot)%fnow(ji,jj,1) - sf(jp_qsr)%fnow(ji,jj,1) 
     150               qsr (ji,jj) = sf(jp_qsr )%fnow(ji,jj,1) 
     151               emp (ji,jj) = sf(jp_emp )%fnow(ji,jj,1) 
    152152            END DO 
    153153         END DO 

  • branches/DEV_r2191_3partymerge2010/NEMO/OPA_SRC/SBC/sbcfwb.F90

    r2219 r2222  
    168168            IF( lk_mpp )   CALL  mpp_sum( zsurf_neg ) 
    169169            IF( lk_mpp )   CALL  mpp_sum( zsurf_pos ) 
     170            IF( lk_mpp )   CALL  mpp_sum( zsurf_neg ) 
     171            IF( lk_mpp )   CALL  mpp_sum( zsurf_pos ) 
    170172             
    171173            IF( z_fwf < 0.e0 ) THEN 

  • branches/DEV_r2191_3partymerge2010/NEMO/OPA_SRC/SBC/sbcice_if.F90

    r2208 r2222  
    8181            CALL ctl_stop( 'sbc_ice_if: unable to allocate sf_ice structure' )   ;   RETURN 
    8282         ENDIF 
    83          ALLOCATE( sf_ice(1)%fnow(jpi,jpj) ) 
    84          ALLOCATE( sf_ice(1)%fdta(jpi,jpj,2) ) 
     83         ALLOCATE( sf_ice(1)%fnow(jpi,jpj,1) ) 
     84         IF( sn_ice%ln_tint ) ALLOCATE( sf_ice(1)%fdta(jpi,jpj,1,2) ) 
    8585 
    8686 
    … …  
    107107               ! 
    108108               zt_fzp  = fr_i(ji,jj)                        ! freezing point temperature 
    109                zfr_obs = sf_ice(1)%fnow(ji,jj)              ! observed ice cover 
     109               zfr_obs = sf_ice(1)%fnow(ji,jj,1)              ! observed ice cover 
    110110               !                                            ! ocean ice fraction (0/1) from the freezing point temperature 
    111111               IF( sst_m(ji,jj) <= zt_fzp ) THEN   ;   fr_i(ji,jj) = 1.e0 

  • branches/DEV_r2191_3partymerge2010/NEMO/OPA_SRC/SBC/sbcice_lim_2.F90

    r2208 r2222  
    88   !! History :  1.0   !  06-2006  (G. Madec)  from icestp_2.F90 
    99   !!            3.0   !  08-2008  (S. Masson, E. .... ) coupled interface 
     10   !!            3.3   !  05-2009  (G.Garric) addition of the lim2_evp case 
    1011   !!---------------------------------------------------------------------- 
    1112#if defined key_lim2 
    … …  
    5253   PUBLIC sbc_ice_lim_2 ! routine called by sbcmod.F90 
    5354    
    54    CHARACTER(len=1) ::   cl_grid = 'B'     ! type of grid used in ice dynamics 
    55  
    5655   !! * Substitutions 
    5756#  include "domzgr_substitute.h90" 
    … …  
    171170         !  Ice model step  ! 
    172171         ! ---------------- ! 
     172         numit = numit + nn_fsbc                                             ! Ice model time step 
    173173 
    174174                                        CALL lim_rst_opn_2  ( kt )      ! Open Ice restart file 

  • branches/DEV_r2191_3partymerge2010/NEMO/OPA_SRC/SBC/sbcrnf.F90

    r2208 r2222  
    113113            DO jj = 1, jpj 
    114114               DO ji = 1, jpi 
    115                   IF( gphit(ji,jj) > 40 .AND. gphit(ji,jj) < 65 )   sf_rnf(1)%fnow(ji,jj) = 0.85 * sf_rnf(1)%fnow(ji,jj) 
     115                  IF( gphit(ji,jj) > 40 .AND. gphit(ji,jj) < 65 )   sf_rnf(1)%fnow(ji,jj,1) = 0.85 * sf_rnf(1)%fnow(ji,jj,1) 
    116116               END DO 
    117117            END DO 
    … …  
    267267            CALL ctl_stop( 'sbc_rnf: unable to allocate sf_rnf structure' )   ;   RETURN 
    268268         ENDIF 
    269          ALLOCATE( sf_rnf(1)%fnow(jpi,jpj)   )   ;   ALLOCATE( sf_rnf(1)%fdta(jpi,jpj,2) ) 
     269         ALLOCATE( sf_rnf(1)%fnow(jpi,jpj,1)   ) 
     270         IF( sf_rnf(1)%ln_tint ) ALLOCATE( sf_rnf(1)%fdta(jpi,jpj,1,2) ) 
    270271         !                                          ! fill sf_rnf with the namelist (sn_rnf) and control print 
    271272         CALL fld_fill( sf_rnf, (/ sn_rnf /), cn_dir, 'sbc_rnf_init', 'read runoffs data', 'namsbc_rnf' ) 
    … …  
    278279               CALL ctl_stop( 'sbc_rnf_init: unable to allocate sf_t_rnf structure' )   ;   RETURN 
    279280            ENDIF 
    280             ALLOCATE( sf_t_rnf(1)%fnow(jpi,jpj)   )   ;   ALLOCATE( sf_t_rnf(1)%fdta(jpi,jpj,2) ) 
     281            ALLOCATE( sf_t_rnf(1)%fnow(jpi,jpj)   ) 
     282            IF( sf_t_rnf(1)%ln_tint ) ALLOCATE( sf_t_rnf(1)%fdta(jpi,jpj,1,2) ) 
    281283            CALL fld_fill (sf_t_rnf, (/ sn_t_rnf /), cn_dir, 'sbc_rnf_init', 'read runoff temperature data', 'namsbc_rnf' )   
    282284         ENDIF 
    … …  
    289291               CALL ctl_stop( 'sbc_rnf_init: unable to allocate sf_s_rnf structure' )   ;   RETURN 
    290292            ENDIF 
    291             ALLOCATE( sf_s_rnf(1)%fnow(jpi,jpj)   )   ;   ALLOCATE( sf_s_rnf(1)%fdta(jpi,jpj,2) ) 
     293            ALLOCATE( sf_s_rnf(1)%fnow(jpi,jpj,1)   ) 
     294            IF( sf_s_rnf(1)%ln_tint ) ALLOCATE( sf_s_rnf(1)%fdta(jpi,jpj,1,2) ) 
    292295            CALL fld_fill (sf_s_rnf, (/ sn_s_rnf /), cn_dir, 'sbc_rnf_init', 'read runoff salinity data', 'namsbc_rnf' )   
    293296         ENDIF 

  • branches/DEV_r2191_3partymerge2010/NEMO/OPA_SRC/SBC/sbcssr.F90

    r2208 r2222  
    115115               CALL ctl_stop( 'sbc_ssr: unable to allocate sf_sst structure' )   ;   RETURN 
    116116            ENDIF 
    117             ALLOCATE( sf_sst(1)%fnow(jpi,jpj) ) 
    118             ALLOCATE( sf_sst(1)%fdta(jpi,jpj,2) ) 
     117            ALLOCATE( sf_sst(1)%fnow(jpi,jpj,1) ) 
    119118            ! 
    120119            ! fill sf_sst with sn_sst and control print 
    121120            CALL fld_fill( sf_sst, (/ sn_sst /), cn_dir, 'sbc_ssr', 'SST restoring term toward SST data', 'namsbc_ssr' ) 
     121            IF( sf_sst(1)%ln_tint ) ALLOCATE( sf_sst(1)%fdta(jpi,jpj,1,2) ) 
    122122         ENDIF 
    123123         ! 
    … …  
    128128               CALL ctl_stop( 'sbc_ssr: unable to allocate sf_sss structure' )   ;   RETURN 
    129129            ENDIF 
    130             ALLOCATE( sf_sss(1)%fnow(jpi,jpj) ) 
    131             ALLOCATE( sf_sss(1)%fdta(jpi,jpj,2) ) 
     130            ALLOCATE( sf_sss(1)%fnow(jpi,jpj,1) ) 
    132131            ! 
    133132            ! fill sf_sss with sn_sss and control print 
    134133            CALL fld_fill( sf_sss, (/ sn_sss /), cn_dir, 'sbc_ssr', 'SSS restoring term toward SSS data', 'namsbc_ssr' ) 
     134            IF( sf_sss(1)%ln_tint )ALLOCATE( sf_sss(1)%fdta(jpi,jpj,1,2) ) 
    135135         ENDIF 
    136136         ! 
    … …  
    153153               DO jj = 1, jpj 
    154154                  DO ji = 1, jpi 
    155                      zqrp = rn_dqdt * ( sst_m(ji,jj) - sf_sst(1)%fnow(ji,jj) ) 
     155                     zqrp = rn_dqdt * ( sst_m(ji,jj) - sf_sst(1)%fnow(ji,jj,1) ) 
    156156                     qns(ji,jj) = qns(ji,jj) + zqrp 
    157157                     qrp(ji,jj) = zqrp 
    … …  
    167167                  DO ji = 1, jpi 
    168168                     zerp = zsrp * ( 1. - 2.*rnfmsk(ji,jj) )   &      ! No damping in vicinity of river mouths 
    169                         &        * ( sss_m(ji,jj) - sf_sss(1)%fnow(ji,jj) )   & 
     169                        &        * ( sss_m(ji,jj) - sf_sss(1)%fnow(ji,jj,1) )   & 
    170170                        &        / ( sss_m(ji,jj) + 1.e-20   ) 
    171171                     emps(ji,jj) = emps(ji,jj) + zerp 
    … …  
    182182                  DO ji = 1, jpi                             
    183183                     zerp = zsrp * ( 1. - 2.*rnfmsk(ji,jj) )   &      ! No damping in vicinity of river mouths 
    184                         &        * ( sss_m(ji,jj) - sf_sss(1)%fnow(ji,jj) )   & 
     184                        &        * ( sss_m(ji,jj) - sf_sss(1)%fnow(ji,jj,1) )   & 
    185185                        &        / ( sss_m(ji,jj) + 1.e-20   ) 
    186186                     IF( ln_sssr_bnd )   zerp = SIGN( 1., zerp ) * MIN( zerp_bnd, ABS(zerp) ) 

  • branches/DEV_r2191_3partymerge2010/NEMO/OPA_SRC/TRA/traqsr.F90

    r2208 r2222  
    144144!CDIR NOVERRCHK 
    145145                  DO ji = 1, jpi 
    146                      zchl = MIN( 10. , MAX( 0.03, sf_chl(1)%fnow(ji,jj) ) ) 
     146                     zchl = MIN( 10. , MAX( 0.03, sf_chl(1)%fnow(ji,jj,1) ) ) 
    147147                     irgb = NINT( 41 + 20.*LOG10(zchl) + 1.e-15 ) 
    148148                     zekb(ji,jj) = rkrgb(1,irgb) 
    … …  
    336336                  CALL ctl_stop( 'tra_qsr_init: unable to allocate sf_chl structure' )   ;   RETURN 
    337337               ENDIF 
    338                ALLOCATE( sf_chl(1)%fnow(jpi,jpj)   ) 
    339                ALLOCATE( sf_chl(1)%fdta(jpi,jpj,2) ) 
     338               ALLOCATE( sf_chl(1)%fnow(jpi,jpj,1)   ) 
     339               IF( sn_chl%ln_tint )ALLOCATE( sf_chl(1)%fdta(jpi,jpj,1,2) ) 
    340340               !                                        ! fill sf_chl with sn_chl and control print 
    341341               CALL fld_fill( sf_chl, (/ sn_chl /), cn_dir, 'tra_qsr_init',   & 

  • branches/DEV_r2191_3partymerge2010/NEMO/OPA_SRC/TRA/trazdf.F90

    r2208 r2222  
    125125      USE zdftke_old 
    126126      USE zdftke 
     127      USE zdfgls 
    127128      USE zdfkpp 
    128129      !!---------------------------------------------------------------------- 
    … …  
    139140 
    140141      ! Force implicit schemes 
    141       IF( lk_zdftke_old .OR. lk_zdftke .OR. lk_zdfkpp )   nzdf = 1      ! TKE or KPP physics 
    142       IF( ln_traldf_iso                               )   nzdf = 1      ! iso-neutral lateral physics 
    143       IF( ln_traldf_hor .AND. ln_sco                  )   nzdf = 1      ! horizontal lateral physics in s-coordinate 
     142      IF( lk_zdftke_old .OR. lk_zdftke .OR. lk_zdfgls .OR. lk_zdfkpp )   nzdf = 1      ! TKE, GLS or KPP physics 
     143      IF( ln_traldf_iso                                              )   nzdf = 1      ! iso-neutral lateral physics 
     144      IF( ln_traldf_hor .AND. ln_sco                                 )   nzdf = 1      ! horizontal lateral physics in s-coordinate 
    144145 
    145146      IF( ln_zdfexp .AND. nzdf == 1 )   THEN 

  • branches/DEV_r2191_3partymerge2010/NEMO/OPA_SRC/ZDF/zdfbfr.F90

    r2208 r2222  
    2929    
    3030   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   bfrua , bfrva   !: Bottom friction coefficients set in zdfbfr 
     31#if defined key_zdfgls 
     32   REAL(wp), PUBLIC                     ::   rn_hbro =  0.003_wp  ! Bottom roughness (m) 
     33   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   wbotu, wbotv    !  Bottom stresses 
     34#endif 
    3135 
    3236   !                                    !!* Namelist nambfr: bottom friction namelist * 

  • branches/DEV_r2191_3partymerge2010/NEMO/OPA_SRC/ZDF/zdfini.F90

    r2208 r2222  
    2020   USE zdftke_old      ! TKE vertical mixing  (old scheme) 
    2121   USE zdftke          ! TKE vertical mixing 
     22   USE zdfgls          ! GLS vertical mixing 
    2223   USE zdfkpp          ! KPP vertical mixing           
    2324   USE zdfddm          ! double diffusion mixing       
    … …  
    106107         ioptio = ioptio+1 
    107108      ENDIF 
     109      IF( lk_zdfgls ) THEN 
     110         IF(lwp) WRITE(numout,*) '      GLS dependent eddy coefficients' 
     111         ioptio = ioptio+1 
     112      ENDIF 
    108113      IF( lk_zdfkpp ) THEN 
    109114         IF(lwp) WRITE(numout,*) '      KPP dependent eddy coefficients' 
    … …  
    128133         IF(lwp) WRITE(numout,*) '      use the 1.5 turbulent closure' 
    129134      ENDIF 
     135      IF( lk_zdfgls ) THEN 
     136         IF(lwp) WRITE(numout,*) '      use the GLS closure scheme' 
     137      ENDIF 
    130138      IF( lk_zdfkpp ) THEN 
    131139         IF(lwp) WRITE(numout,*) '      use the KPP closure scheme' 
    … …  
    136144      ENDIF 
    137145      IF ( ioptio > 1 .AND. .NOT. lk_esopa )   CALL ctl_stop( ' chose between ln_zdfnpc and ln_zdfevd' ) 
    138       IF( ioptio == 0 .AND. .NOT.( lk_zdftke_old .OR. lk_zdftke .OR. lk_zdfkpp ) ) & 
    139          CALL ctl_stop( ' except for TKE sor KPP physics, a convection scheme is', & 
     146      IF( ioptio == 0 .AND. .NOT.( lk_zdftke_old .OR. lk_zdftke .OR. lk_zdfgls .OR. lk_zdfkpp ) ) & 
     147         CALL ctl_stop( ' except for TKE, GLS or KPP physics, a convection scheme is', & 
    140148         &              ' required: ln_zdfevd or ln_zdfnpc logicals' ) 
    141149 

  • branches/DEV_r2191_3partymerge2010/NEMO/OPA_SRC/step.F90

    r2208 r2222  
    9191   USE zdftke_old      ! old TKE vertical mixing      (zdf_tke_old routine) 
    9292   USE zdftke          ! TKE vertical mixing              (zdf_tke routine) 
     93   USE zdfgls          ! GLS vertical mixing              (zdf_gls routine)   
    9394   USE zdfkpp          ! KPP vertical mixing              (zdf_kpp routine) 
    9495   USE zdfddm          ! double diffusion mixing          (zdf_ddm routine) 
    … …  
    208209      IF( lk_zdftke_old )   CALL zdf_tke_old( kstp )  ! TKE closure scheme for Kz (old scheme) 
    209210      IF( lk_zdftke     )   CALL zdf_tke    ( kstp )  ! TKE closure scheme for Kz 
     211      IF( lk_zdfgls     )   CALL zdf_gls    ( kstp )  ! GLS closure scheme for Kz 
    210212      IF( lk_zdfkpp     )   CALL zdf_kpp    ( kstp )  ! KPP closure scheme for Kz 
    211213      IF( lk_zdfcst     )   THEN                      ! Constant Kz (reset avt, avm[uv] to the background value) 
    … …  
    227229                                                      ! write tke information in the restart file 
    228230      IF( lrst_oce .AND. lk_zdftke )   CALL tke_rst( kstp, 'WRITE' ) 
     231                                                      ! write gls information in the restart file 
     232      IF( lrst_oce .AND. lk_zdfgls )   CALL gls_rst( kstp, 'WRITE' ) 
    229233      ! 
    230234      !  LATERAL  PHYSICS  
Note: See TracChangeset for help on using the changeset viewer.