Changeset 4689


Ignore:
Timestamp:
2014-06-25T01:40:18+02:00 (6 years ago)
Author:
clem
Message:

new version of LIM3 with perfect conservation of heat, see ticket #1352

Location:
trunk/NEMOGCM/NEMO/OPA_SRC
Files:
12 edited

Legend:

Unmodified
Added
Removed
  • trunk/NEMOGCM/NEMO/OPA_SRC/BDY/bdy_oce.F90

    r4608 r4689  
    109109                                                              !: = 1 read it in a NetCDF file 
    110110#endif 
     111#if defined key_lim3 
     112   REAL,    DIMENSION(jp_bdy) ::   rn_ice_tem               !: choice of the temperature of incoming sea ice 
     113   REAL,    DIMENSION(jp_bdy) ::   rn_ice_sal               !: choice of the salinity    of incoming sea ice 
     114   REAL,    DIMENSION(jp_bdy) ::   rn_ice_age               !: choice of the age         of incoming sea ice 
     115#endif 
    111116   ! 
    112117    
  • trunk/NEMOGCM/NEMO/OPA_SRC/BDY/bdydta.F90

    r4673 r4689  
    676676               CALL iom_close ( inum ) 
    677677               !CALL fld_clopn ( bn_a_i, nyear, nmonth, nday, ldstop=.TRUE. ) 
    678                !CALL iom_open ( bn_a_i %clname, inum ) 
     678               !CALL iom_open ( bn_a_i%clname, inum ) 
    679679               !id1 = iom_varid ( bn_a_i%num, bn_a_i%clvar, kdimsz=zdimsz, kndims=zndims, ldstop = .FALSE. ) 
    680680                IF ( zndims == 4 ) THEN 
     
    907907   !!============================================================================== 
    908908END MODULE bdydta 
    909  
    910  
  • trunk/NEMOGCM/NEMO/OPA_SRC/BDY/bdydyn.F90

    r4370 r4689  
    3030   USE lbclnk          ! ocean lateral boundary conditions (or mpp link) 
    3131   USE in_out_manager  ! 
    32    USE domvvl 
     32   USE domvvl          ! variable volume 
    3333 
    3434   IMPLICIT NONE 
  • trunk/NEMOGCM/NEMO/OPA_SRC/BDY/bdyice_lim.F90

    r4333 r4689  
    9999      REAL(wp) ::   zinda, ztmelts, zdh 
    100100 
    101       REAL(wp), PARAMETER  ::   zsal = 6.3    ! arbitrary salinity    for incoming ice 
    102       REAL(wp), PARAMETER  ::   ztem = 270.0  ! arbitrary temperature for incoming ice 
    103       REAL(wp), PARAMETER  ::   zage = 30.0   ! arbitrary age         for incoming ice 
    104101      !!------------------------------------------------------------------------------ 
    105102      ! 
     
    233230 
    234231               ! Ice salinity, age, temperature 
    235                sm_i(ji,jj,jl)   = zinda * zsal  + ( 1.0 - zinda ) * s_i_min 
    236                o_i(ji,jj,jl)    = zinda * zage  + ( 1.0 - zinda ) 
    237                t_su(ji,jj,jl)   = zinda * ztem  + ( 1.0 - zinda ) * ztem 
     232               sm_i(ji,jj,jl)   = zinda * rn_ice_sal(ib_bdy)  + ( 1.0 - zinda ) * s_i_min 
     233               o_i(ji,jj,jl)    = zinda * rn_ice_age(ib_bdy)  + ( 1.0 - zinda ) 
     234               t_su(ji,jj,jl)   = zinda * rn_ice_tem(ib_bdy)  + ( 1.0 - zinda ) * rn_ice_tem(ib_bdy) 
    238235               DO jk = 1, nlay_s 
    239                   t_s(ji,jj,jk,jl) = zinda * ztem + ( 1.0 - zinda ) * rtt 
     236                  t_s(ji,jj,jk,jl) = zinda * rn_ice_tem(ib_bdy) + ( 1.0 - zinda ) * rtt 
    240237               END DO 
    241238               DO jk = 1, nlay_i 
    242                   t_i(ji,jj,jk,jl) = zinda * ztem + ( 1.0 - zinda ) * rtt  
    243                   s_i(ji,jj,jk,jl) = zinda * zsal + ( 1.0 - zinda ) * s_i_min 
     239                  t_i(ji,jj,jk,jl) = zinda * rn_ice_tem(ib_bdy) + ( 1.0 - zinda ) * rtt  
     240                  s_i(ji,jj,jk,jl) = zinda * rn_ice_sal(ib_bdy) + ( 1.0 - zinda ) * s_i_min 
    244241               END DO 
    245242                
     
    259256 
    260257            END SELECT 
     258 
     259            ! if salinity is constant, then overwrite rn_ice_sal 
     260            IF( num_sal == 1 ) THEN 
     261               sm_i(ji,jj,jl)   = bulk_sal 
     262               s_i (ji,jj,:,jl) = bulk_sal 
     263            ENDIF 
    261264 
    262265            ! contents 
     
    338341      DO ib_bdy=1, nb_bdy 
    339342         ! 
    340          SELECT CASE( nn_ice_lim(ib_bdy) ) 
     343         SELECT CASE( cn_ice_lim(ib_bdy) ) 
    341344 
    342345         CASE('none') 
     
    355358                  ji    = idx_bdy(ib_bdy)%nbi(jb,jgrd) 
    356359                  jj    = idx_bdy(ib_bdy)%nbj(jb,jgrd) 
    357                   zflag = idx_bdy(ib_bdy)%flagu(jb) 
     360                  zflag = idx_bdy(ib_bdy)%flagu(jb,jgrd) 
    358361                   
    359362                  IF ( ABS( zflag ) == 1. ) THEN  ! eastern and western boundaries 
     
    384387                  ji    = idx_bdy(ib_bdy)%nbi(jb,jgrd) 
    385388                  jj    = idx_bdy(ib_bdy)%nbj(jb,jgrd) 
    386                   zflag = idx_bdy(ib_bdy)%flagv(jb) 
     389                  zflag = idx_bdy(ib_bdy)%flagv(jb,jgrd) 
    387390                   
    388391                  IF ( ABS( zflag ) == 1. ) THEN  ! northern and southern boundaries 
  • trunk/NEMOGCM/NEMO/OPA_SRC/BDY/bdyini.F90

    r4624 r4689  
    102102#if ( defined key_lim2 || defined key_lim3 ) 
    103103         &             cn_ice_lim, nn_ice_lim_dta,                           & 
     104#endif 
     105#if defined key_lim3 
     106         &             rn_ice_tem, rn_ice_sal, rn_ice_age,                 & 
    104107#endif 
    105108         &             ln_vol, nn_volctl, nn_rimwidth 
     
    359362        ENDIF 
    360363        IF(lwp) WRITE(numout,*) 
     364        IF(lwp) WRITE(numout,*) '      tem of bdy sea-ice = ', rn_ice_tem(ib_bdy)          
     365        IF(lwp) WRITE(numout,*) '      sal of bdy sea-ice = ', rn_ice_sal(ib_bdy)          
     366        IF(lwp) WRITE(numout,*) '      age of bdy sea-ice = ', rn_ice_age(ib_bdy)          
    361367#endif 
    362368 
  • trunk/NEMOGCM/NEMO/OPA_SRC/DOM/phycst.F90

    r3625 r4689  
    5454   REAL(wp), PUBLIC ::   r1_rau0                     !: = 1. / rau0                   [m3/kg] 
    5555   REAL(wp), PUBLIC ::   rauw     = 1000._wp         !: volumic mass of pure water    [m3/kg] 
    56    REAL(wp), PUBLIC ::   rcp      =    4.e3_wp       !: ocean specific heat           [J/Kelvin] 
    57    REAL(wp), PUBLIC ::   r1_rcp                      !: = 1. / rcp                    [Kelvin/J] 
     56   REAL(wp), PUBLIC ::   rcp      =    4.e3_wp       !: ocean specific heat           [J/kg/K] 
     57   REAL(wp), PUBLIC ::   r1_rcp                      !: = 1. / rcp                    [kg.K/J] 
    5858   REAL(wp), PUBLIC ::   r1_rau0_rcp                 !: = 1. / ( rau0 * rcp ) 
    5959 
     
    6969#if defined key_lim3 || defined key_cice 
    7070   REAL(wp), PUBLIC ::   rhoic    =  917._wp         !: volumic mass of sea ice                               [kg/m3] 
    71    REAL(wp), PUBLIC ::   rcdic    =    2.034396_wp   !: thermal conductivity of fresh ice 
    72    REAL(wp), PUBLIC ::   rcdsn    =    0.31_wp       !: thermal conductivity of snow 
    73    REAL(wp), PUBLIC ::   cpic     = 2067.0_wp        !: specific heat for ice  
     71   REAL(wp), PUBLIC ::   rcdic    =    2.034396_wp   !: thermal conductivity of fresh ice                     [W/m/K] 
     72   REAL(wp), PUBLIC ::   rcdsn    =    0.31_wp       !: thermal conductivity of snow                          [W/m/K]  
     73   REAL(wp), PUBLIC ::   cpic     = 2067.0_wp        !: specific heat for ice                                 [J/kg/K] 
    7474   REAL(wp), PUBLIC ::   lsub     =    2.834e+6_wp   !: pure ice latent heat of sublimation                   [J/kg] 
    7575   REAL(wp), PUBLIC ::   lfus     =    0.334e+6_wp   !: latent heat of fusion of fresh ice                    [J/kg] 
    76    REAL(wp), PUBLIC ::   tmut     =    0.054_wp      !: decrease of seawater meltpoint with salinity 
     76   REAL(wp), PUBLIC ::   tmut     =    0.054_wp      !: decrease of seawater meltpoint with salinity          [degC/ppt] 
    7777   REAL(wp), PUBLIC ::   xlsn                        !: = lfus*rhosn (volumetric latent heat fusion of snow)  [J/m3] 
    7878#else 
  • trunk/NEMOGCM/NEMO/OPA_SRC/IOM/iom.F90

    r4650 r4689  
    3232   USE trc_oce, ONLY :   nn_dttrc        !  !: frequency of step on passive tracers 
    3333   USE icb_oce, ONLY :   nclasses, class_num       !  !: iceberg classes 
     34   USE par_ice 
    3435   USE domngb          ! ocean space and time domain 
    3536   USE phycst          ! physical constants 
     
    4950#endif 
    5051   PUBLIC iom_init, iom_swap, iom_open, iom_close, iom_setkt, iom_varid, iom_get, iom_gettime, iom_rstput, iom_put 
    51    PUBLIC iom_getatt, iom_context_finalize 
     52   PUBLIC iom_getatt, iom_use, iom_context_finalize 
    5253 
    5354   PRIVATE iom_rp0d, iom_rp1d, iom_rp2d, iom_rp3d 
     
    143144      CALL iom_set_axis_attr( "nfloat", (/ (REAL(ji,wp), ji=1,nfloat) /) ) 
    144145# endif 
     146#if defined key_lim3 
     147      CALL iom_set_axis_attr( "ncatice", (/ (REAL(ji,wp), ji=1,jpl) /) ) 
     148#endif 
    145149      CALL iom_set_axis_attr( "icbcla", class_num ) 
    146150       
     
    10151019      CHARACTER(LEN=*), INTENT(in) ::   cdname 
    10161020      REAL(wp)        , INTENT(in) ::   pfield0d 
     1021      REAL(wp)        , DIMENSION(jpi,jpj) ::   zz     ! masson 
    10171022#if defined key_iomput 
    1018       CALL xios_send_field(cdname, (/pfield0d/)) 
     1023      zz(:,:)=pfield0d 
     1024      CALL xios_send_field(cdname, zz) 
     1025      !CALL xios_send_field(cdname, (/pfield0d/))  
    10191026#else 
    10201027      IF( .FALSE. )   WRITE(numout,*) cdname, pfield0d   ! useless test to avoid compilation warnings 
     
    12071214      !! 
    12081215      !!---------------------------------------------------------------------- 
    1209       REAL(wp), DIMENSION(1,1) ::   zz = 1. 
     1216      REAL(wp), DIMENSION(1) ::   zz = 1. 
    12101217      !!---------------------------------------------------------------------- 
    12111218      CALL iom_set_domain_attr('scalarpoint', ni_glo=jpnij, nj_glo=1, ibegin=narea, jbegin=1, ni=1, nj=1) 
    1212       CALL iom_set_domain_attr('scalarpoint', data_dim=1) 
    1213       CALL iom_set_domain_attr('scalarpoint', lonvalue=(/ zz /), latvalue=(/ zz /)) 
     1219      CALL iom_set_domain_attr('scalarpoint', data_dim=2, data_ibegin = 1, data_ni = 1, data_jbegin = 1, data_nj = 1) 
     1220      zz=REAL(narea,wp) 
     1221      CALL iom_set_domain_attr('scalarpoint', lonvalue=zz, latvalue=zz) 
    12141222 
    12151223   END SUBROUTINE set_scalar 
     
    14991507 
    15001508#endif 
     1509 
     1510   LOGICAL FUNCTION iom_use( cdname ) 
     1511      CHARACTER(LEN=*), INTENT(in) ::   cdname 
     1512#if defined key_iomput 
     1513      iom_use = xios_field_is_active( cdname ) 
     1514#else 
     1515      iom_use = .FALSE. 
     1516#endif 
     1517   END FUNCTION iom_use 
    15011518    
    15021519   !!====================================================================== 
  • trunk/NEMOGCM/NEMO/OPA_SRC/IOM/iom_nf90.F90

    r4292 r4689  
    217217         CALL iom_nf90_check(NF90_Inquire_Variable(if90id, ivarid, dimids = idimid(1:i_nvd)), clinfo)   ! dimensions ids 
    218218         iom_file(kiomid)%luld(kiv) = .FALSE.   ! default value 
    219          iom_file(kiomid)%dimsz(:,kiv) = 0   ! reset dimsz in case previously used 
     219         iom_file(kiomid)%dimsz(:,kiv) = 0      ! reset dimsz in case previously used 
    220220         DO ji = 1, i_nvd                       ! dimensions size 
    221221            CALL iom_nf90_check(NF90_Inquire_Dimension(if90id, idimid(ji), len = iom_file(kiomid)%dimsz(ji,kiv)), clinfo)    
  • trunk/NEMOGCM/NEMO/OPA_SRC/IOM/restart.F90

    r4334 r4689  
    120120                     CALL iom_rstput( kt, nitrst, numrow, 'hdivb'  , hdivb     ) 
    121121                     CALL iom_rstput( kt, nitrst, numrow, 'sshb'   , sshb      ) 
     122      IF( lk_vvl )   CALL iom_rstput( kt, nitrst, numrow, 'fse3t_b', fse3t_b(:,:,:) ) 
    122123                     ! 
    123124                     CALL iom_rstput( kt, nitrst, numrow, 'un'     , un        )     ! now fields 
     
    210211         CALL iom_get( numror, jpdom_autoglo, 'hdivb'  , hdivb   ) 
    211212         CALL iom_get( numror, jpdom_autoglo, 'sshb'   , sshb    ) 
     213         IF( lk_vvl )   CALL iom_get( numror, jpdom_autoglo, 'fse3t_b', fse3t_b(:,:,:) ) 
    212214      ELSE 
    213215         neuler = 0 
     
    245247         hdivb(:,:,:)   = hdivn(:,:,:) 
    246248         sshb (:,:)     = sshn (:,:) 
    247       ENDIF 
    248       ! 
    249       IF( lk_lim3 ) THEN  
     249         IF( lk_vvl ) THEN 
     250            DO jk = 1, jpk 
     251               fse3t_b(:,:,jk) = fse3t_n(:,:,jk) 
     252            END DO 
     253         ENDIF 
     254      ENDIF 
     255      ! 
     256      IF( lk_lim3 ) THEN 
    250257         CALL iom_get( numror, jpdom_autoglo, 'iatte' , iatte ) ! clem modif 
    251258         CALL iom_get( numror, jpdom_autoglo, 'oatte' , oatte ) ! clem modif 
  • trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90

    r4682 r4689  
    563563      zcoef_dqsb   = rhoa * cpa * Cice 
    564564      zcoef_frca   = 1.0  - 0.3 
     565      ! MV 2014 the proper cloud fraction (mean summer months from the CLIO climato, NH+SH) is 0.19 
     566      zcoef_frca   = 1.0  - 0.19 
    565567 
    566568!!gm brutal.... 
     
    648650               p_qsr(ji,jj,jl) = zztmp * ( 1. - palb(ji,jj,jl) ) * qsr(ji,jj) 
    649651               ! Long  Wave (lw) 
    650                ! iovino 
    651                IF( ff(ji,jj) .GT. 0._wp ) THEN 
    652                   z_qlw(ji,jj,jl) = ( 0.95 * sf(jp_qlw)%fnow(ji,jj,1) - Stef * pst(ji,jj,jl) * zst3 ) * tmask(ji,jj,1) 
    653                ELSE 
    654                   z_qlw(ji,jj,jl) = 0.95 * ( sf(jp_qlw)%fnow(ji,jj,1) - Stef * pst(ji,jj,jl) * zst3 ) * tmask(ji,jj,1) 
    655                ENDIF 
     652               z_qlw(ji,jj,jl) = 0.95 * ( sf(jp_qlw)%fnow(ji,jj,1) - Stef * pst(ji,jj,jl) * zst3 ) * tmask(ji,jj,1) 
    656653               ! lw sensitivity 
    657654               z_dqlw(ji,jj,jl) = zcoef_dqlw * zst3                                                
     
    668665                  &                         * (  11637800. * EXP( -5897.8 / pst(ji,jj,jl) ) / rhoa - sf(jp_humi)%fnow(ji,jj,1)  ) ) 
    669666               ! Latent heat sensitivity for ice (Dqla/Dt) 
    670                p_dqla(ji,jj,jl) = rn_efac * zcoef_dqla * z_wnds_t(ji,jj) / ( zst2 ) * EXP( -5897.8 / pst(ji,jj,jl) ) 
     667               ! MV we also have to cap the sensitivity if the flux is zero 
     668               IF ( p_qla(ji,jj,jl) .GT. 0.0 ) THEN 
     669                  p_dqla(ji,jj,jl) = rn_efac * zcoef_dqla * z_wnds_t(ji,jj) / ( zst2 ) * EXP( -5897.8 / pst(ji,jj,jl) ) 
     670               ELSE 
     671                  p_dqla(ji,jj,jl) = 0.0 
     672               ENDIF 
     673                              
    671674               ! Sensible heat sensitivity (Dqsb_ice/Dtn_ice) 
    672675               z_dqsb(ji,jj,jl) = zcoef_dqsb * z_wnds_t(ji,jj) 
  • trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcfwb.F90

    r4347 r4689  
    5757      !!                =1 global mean of emp set to zero at each nn_fsbc time step 
    5858      !!                =2 annual global mean corrected from previous year 
     59      !!                =3 global mean of emp set to zero at each nn_fsbc time step 
     60      !!                   & spread out over erp area depending its sign 
    5961      !! Note: if sea ice is embedded it is taken into account when computing the budget  
    6062      !!---------------------------------------------------------------------- 
     
    8183            IF( kn_fwb == 1 )   WRITE(numout,*) '          instantaneously set to zero' 
    8284            IF( kn_fwb == 2 )   WRITE(numout,*) '          adjusted from previous year budget' 
    83          ENDIF 
     85            IF( kn_fwb == 3 )   WRITE(numout,*) '          fwf set to zero and spread out over erp area' 
     86         ENDIF 
     87         ! 
     88         IF( kn_fwb == 3 .AND. nn_sssr /= 2 )   CALL ctl_stop( 'sbc_fwb: nn_fwb = 3 requires nn_sssr = 2, we stop ' ) 
    8489         ! 
    8590         area = glob_sum( e1e2t(:,:) )           ! interior global domain surface 
     
    142147         ENDIF 
    143148         ! 
     149      CASE ( 3 )                             !==  global fwf set to zero and spread out over erp area  ==! 
     150         ! 
     151         IF( MOD( kt-1, kn_fsbc ) == 0 ) THEN 
     152            ztmsk_pos(:,:) = tmask_i(:,:)                      ! Select <0 and >0 area of erp 
     153            WHERE( erp < 0._wp )   ztmsk_pos = 0._wp 
     154            ztmsk_neg(:,:) = tmask_i(:,:) - ztmsk_pos(:,:) 
     155            ! 
     156            zsurf_neg = glob_sum( e1e2t(:,:)*ztmsk_neg(:,:) )  ! Area filled by <0 and >0 erp  
     157            zsurf_pos = glob_sum( e1e2t(:,:)*ztmsk_pos(:,:) ) 
     158            !                                                  ! fwf global mean (excluding ocean to ice/snow exchanges)  
     159            z_fwf     = glob_sum( e1e2t(:,:) * ( emp(:,:) - rnf(:,:) - snwice_fmass(:,:) ) ) / area 
     160            !             
     161            IF( z_fwf < 0._wp ) THEN         ! spread out over >0 erp area to increase evaporation 
     162                zsurf_tospread      = zsurf_pos 
     163                ztmsk_tospread(:,:) = ztmsk_pos(:,:) 
     164            ELSE                             ! spread out over <0 erp area to increase precipitation 
     165                zsurf_tospread      = zsurf_neg 
     166                ztmsk_tospread(:,:) = ztmsk_neg(:,:) 
     167            ENDIF 
     168            ! 
     169            zsum_fwf   = glob_sum( e1e2t(:,:) * z_fwf )         ! fwf global mean over <0 or >0 erp area 
     170!!gm :  zsum_fwf   = z_fwf * area   ???  it is right?  I think so.... 
     171            z_fwf_nsrf =  zsum_fwf / ( zsurf_tospread + rsmall ) 
     172            !                                                  ! weight to respect erp field 2D structure  
     173            zsum_erp   = glob_sum( ztmsk_tospread(:,:) * erp(:,:) * e1e2t(:,:) ) 
     174            z_wgt(:,:) = ztmsk_tospread(:,:) * erp(:,:) / ( zsum_erp + rsmall ) 
     175            !                                                  ! final correction term to apply 
     176            zerp_cor(:,:) = -1. * z_fwf_nsrf * zsurf_tospread * z_wgt(:,:) 
     177            ! 
     178!!gm   ===>>>>  lbc_lnk should be useless as all the computation is done over the whole domain ! 
     179            CALL lbc_lnk( zerp_cor, 'T', 1. ) 
     180            ! 
     181            emp(:,:) = emp(:,:) + zerp_cor(:,:) 
     182            qns(:,:) = qns(:,:) - zerp_cor(:,:) * rcp * sst_m(:,:)  ! account for change to the heat budget due to fw correction 
     183            erp(:,:) = erp(:,:) + zerp_cor(:,:) 
     184            ! 
     185            IF( nprint == 1 .AND. lwp ) THEN                   ! control print 
     186               IF( z_fwf < 0._wp ) THEN 
     187                  WRITE(numout,*)'   z_fwf < 0' 
     188                  WRITE(numout,*)'   SUM(erp+)     = ', SUM( ztmsk_tospread(:,:)*erp(:,:)*e1e2t(:,:) )*1.e-9,' Sv' 
     189               ELSE 
     190                  WRITE(numout,*)'   z_fwf >= 0' 
     191                  WRITE(numout,*)'   SUM(erp-)     = ', SUM( ztmsk_tospread(:,:)*erp(:,:)*e1e2t(:,:) )*1.e-9,' Sv' 
     192               ENDIF 
     193               WRITE(numout,*)'   SUM(empG)     = ', SUM( z_fwf*e1e2t(:,:) )*1.e-9,' Sv' 
     194               WRITE(numout,*)'   z_fwf         = ', z_fwf      ,' Kg/m2/s' 
     195               WRITE(numout,*)'   z_fwf_nsrf    = ', z_fwf_nsrf ,' Kg/m2/s' 
     196               WRITE(numout,*)'   MIN(zerp_cor) = ', MINVAL(zerp_cor)  
     197               WRITE(numout,*)'   MAX(zerp_cor) = ', MAXVAL(zerp_cor)  
     198            ENDIF 
     199         ENDIF 
     200         ! 
    144201      CASE DEFAULT                           !==  you should never be there  ==! 
    145          CALL ctl_stop( 'sbc_fwb : wrong nn_fwb value for the FreshWater Budget correction, choose either 1 or 2' ) 
     202         CALL ctl_stop( 'sbc_fwb : wrong nn_fwb value for the FreshWater Budget correction, choose either 1, 2 or 3' ) 
    146203         ! 
    147204      END SELECT 
  • trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90

    r4333 r4689  
    5959   USE prtctl          ! Print control 
    6060   USE lib_fortran     !  
     61   USE cpl_oasis3, ONLY : lk_cpl 
    6162 
    6263#if defined key_bdy  
     
    6869 
    6970   PUBLIC sbc_ice_lim  ! routine called by sbcmod.F90 
     71   PUBLIC lim_prt_state 
    7072    
    7173   !! * Substitutions 
     
    133135      INTEGER, INTENT(in) ::   kblk    ! type of bulk (=3 CLIO, =4 CORE) 
    134136      !! 
    135       INTEGER  ::   jl      ! dummy loop index 
     137      INTEGER  ::   ji, jj, jl, jk      ! dummy loop index 
    136138      REAL(wp) ::   zcoef   ! local scalar 
    137139      REAL(wp), POINTER, DIMENSION(:,:,:)   ::   zalb_ice_os, zalb_ice_cs  ! albedo of the ice under overcast/clear sky 
     
    146148      REAL(wp), POINTER, DIMENSION(:,:) :: z_dqns_ice_all  ! Mean d(qns)/dT over all categories 
    147149      REAL(wp), POINTER, DIMENSION(:,:) :: z_dqla_ice_all  ! Mean d(qla)/dT over all categories 
     150      REAL(wp) ::   ztmelts           ! clem 2014: for HC diags 
     151      REAL(wp) ::   epsi20 = 1.e-20   ! 
    148152      !!---------------------------------------------------------------------- 
    149153 
     
    152156      IF( nn_timing == 1 )  CALL timing_start('sbc_ice_lim') 
    153157 
    154       CALL wrk_alloc( jpi,jpj,jpl, zalb_ice_os, zalb_ice_cs ) 
    155  
    156 #if defined key_coupled 
    157       IF ( ln_cpl .OR. ln_iceflx_ave .OR. ln_iceflx_linear ) CALL wrk_alloc( jpi,jpj,jpl, zalb_ice) 
    158       IF ( ln_iceflx_ave .OR. ln_iceflx_linear ) & 
    159          &   CALL wrk_alloc( jpi,jpj, ztem_ice_all, zalb_ice_all, z_qsr_ice_all, z_qns_ice_all, z_qla_ice_all, z_dqns_ice_all, z_dqla_ice_all) 
    160 #endif 
     158      CALL wrk_alloc( jpi,jpj,jpl, zalb_ice_os, zalb_ice_cs, zalb_ice ) 
     159 
     160      IF( lk_cpl ) THEN 
     161         IF ( ln_iceflx_ave .OR. ln_iceflx_linear ) & 
     162            &   CALL wrk_alloc( jpi,jpj, ztem_ice_all, zalb_ice_all, z_qsr_ice_all, z_qns_ice_all, z_qla_ice_all, z_dqns_ice_all, z_dqla_ice_all) 
     163      ENDIF 
    161164 
    162165      IF( kt == nit000 ) THEN 
     
    168171         ! 
    169172         IF( ln_nicep ) THEN      ! control print at a given point 
    170             jiindx = 177   ;   jjindx = 112 
     173            jiindx = 15    ;   jjindx =  44 
    171174            IF(lwp) WRITE(numout,*) ' The debugging point is : jiindx : ',jiindx, ' jjindx : ',jjindx 
    172175         ENDIF 
     
    176179      IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN     !  Ice time-step only  ! 
    177180         !                                     !----------------------! 
    178          !                                           !  Bulk Formulea ! 
     181         !                                           !  Bulk Formulae ! 
    179182         !                                           !----------------! 
    180183         ! 
    181184         u_oce(:,:) = ssu_m(:,:)                     ! mean surface ocean current at ice velocity point 
    182185         v_oce(:,:) = ssv_m(:,:)                     ! (C-grid dynamics :  U- & V-points as the ocean) 
    183          ! 
    184          t_bo(:,:) = tfreez( sss_m ) +  rt0          ! masked sea surface freezing temperature [Kelvin] 
    185          !                                           ! (set to rt0 over land) 
     186 
     187         ! masked sea surface freezing temperature [Kelvin] 
     188         t_bo(:,:) = ( tfreez( sss_m ) +  rt0 ) * tmask(:,:,1) + rt0 * ( 1. - tmask(:,:,1) ) 
     189 
    186190         CALL albedo_ice( t_su, ht_i, ht_s, zalb_ice_cs, zalb_ice_os )  ! ... ice albedo 
    187191 
     
    192196         IF ( ln_cpl ) zalb_ice (:,:,:) = 0.5 * ( zalb_ice_cs (:,:,:) +  zalb_ice_os (:,:,:) ) 
    193197          
    194 #if defined key_coupled 
    195          IF ( ln_iceflx_ave .OR. ln_iceflx_linear ) THEN 
    196             ! 
    197             ! Compute mean albedo and temperature 
    198             zalb_ice_all (:,:) = fice_ice_ave ( zalb_ice (:,:,:) )  
    199             ztem_ice_all (:,:) = fice_ice_ave ( tn_ice   (:,:,:) )  
    200             ! 
     198         IF( lk_cpl ) THEN 
     199            IF ( ln_iceflx_ave .OR. ln_iceflx_linear ) THEN 
     200               ! 
     201               ! Compute mean albedo and temperature 
     202               zalb_ice_all (:,:) = fice_ice_ave ( zalb_ice (:,:,:) )  
     203               ztem_ice_all (:,:) = fice_ice_ave ( tn_ice   (:,:,:) )  
     204               ! 
     205            ENDIF 
    201206         ENDIF 
    202 #endif 
    203207                                               ! Bulk formulea - provides the following fields: 
    204208         ! utau_ice, vtau_ice : surface ice stress                     (U- & V-points)   [N/m2] 
     
    218222            !          
    219223         CASE( 4 )                                       ! CORE bulk formulation 
    220             CALL blk_ice_core( t_su , u_ice     , v_ice     , zalb_ice_cs,               & 
     224            ! MV 2014 
     225            ! We must account for cloud fraction in the computation of the albedo 
     226            ! The present ref just uses the clear sky value 
     227            ! The overcast sky value is 0.06 higher, and polar skies are mostly overcast 
     228            ! CORE has no cloud fraction, hence we must prescribe it 
     229            ! Mean summer cloud fraction computed from CLIO = 0.81 
     230            zalb_ice(:,:,:) = 0.19 * zalb_ice_cs(:,:,:) + 0.81 * zalb_ice_os(:,:,:) 
     231            ! Following line, we replace zalb_ice_cs by simply zalb_ice 
     232            CALL blk_ice_core( t_su , u_ice     , v_ice     , zalb_ice   ,               & 
    221233               &                      utau_ice  , vtau_ice  , qns_ice    , qsr_ice   ,   & 
    222234               &                      qla_ice   , dqns_ice  , dqla_ice   ,               & 
     
    239251 
    240252         ! Average over all categories 
    241 #if defined key_coupled 
     253         IF( lk_cpl ) THEN 
    242254         IF ( ln_iceflx_ave .OR. ln_iceflx_linear ) THEN 
    243255 
     
    269281            END IF 
    270282         END IF 
    271 #endif 
     283         ENDIF 
    272284         !                                           !----------------------! 
    273285         !                                           ! LIM-3  time-stepping ! 
     
    285297         old_smv_i(:,:,:)   = smv_i(:,:,:)     ! salt content 
    286298         old_oa_i (:,:,:)   = oa_i (:,:,:)     ! areal age content 
    287          ! 
    288          old_u_ice(:,:) = u_ice(:,:) 
    289          old_v_ice(:,:) = v_ice(:,:) 
    290          !                                           ! intialisation to zero    !!gm is it truly necessary ??? 
     299         old_u_ice(:,:)     = u_ice(:,:) 
     300         old_v_ice(:,:)     = v_ice(:,:) 
     301 
     302         ! trends    !!gm is it truly necessary ??? 
    291303         d_a_i_thd  (:,:,:)   = 0._wp   ;   d_a_i_trp  (:,:,:)   = 0._wp 
    292304         d_v_i_thd  (:,:,:)   = 0._wp   ;   d_v_i_trp  (:,:,:)   = 0._wp 
     
    296308         d_smv_i_thd(:,:,:)   = 0._wp   ;   d_smv_i_trp(:,:,:)   = 0._wp 
    297309         d_oa_i_thd (:,:,:)   = 0._wp   ;   d_oa_i_trp (:,:,:)   = 0._wp 
    298          ! 
    299          d_u_ice_dyn(:,:) = 0._wp 
    300          d_v_ice_dyn(:,:) = 0._wp 
    301          ! 
    302          sfx    (:,:) = 0._wp   ;   sfx_thd  (:,:) = 0._wp 
    303          sfx_bri(:,:) = 0._wp   ;   sfx_mec  (:,:) = 0._wp   ;   sfx_res  (:,:) = 0._wp 
    304          fhbri  (:,:) = 0._wp   ;   fheat_mec(:,:) = 0._wp   ;   fheat_res(:,:) = 0._wp 
    305          fhmec  (:,:) = 0._wp   ;    
    306          fmmec  (:,:) = 0._wp 
    307          fmmflx (:,:) = 0._wp      
    308          focea2D(:,:) = 0._wp 
    309          fsup2D (:,:) = 0._wp 
    310  
    311          ! used in limthd.F90 
    312          rdvosif(:,:) = 0._wp   ! variation of ice volume at surface 
    313          rdvobif(:,:) = 0._wp   ! variation of ice volume at bottom 
    314          fdvolif(:,:) = 0._wp   ! total variation of ice volume 
    315          rdvonif(:,:) = 0._wp   ! lateral variation of ice volume 
    316          fstric (:,:) = 0._wp   ! part of solar radiation transmitted through the ice 
    317          ffltbif(:,:) = 0._wp   ! linked with fstric 
    318          qfvbq  (:,:) = 0._wp   ! linked with fstric 
    319          rdm_snw(:,:) = 0._wp   ! variation of snow mass per unit area 
    320          rdm_ice(:,:) = 0._wp   ! variation of ice mass per unit area 
    321          hicifp (:,:) = 0._wp   ! daily thermodynamic ice production.  
    322          ! 
    323          diag_sni_gr(:,:) = 0._wp   ;   diag_lat_gr(:,:) = 0._wp 
    324          diag_bot_gr(:,:) = 0._wp   ;   diag_dyn_gr(:,:) = 0._wp 
    325          diag_bot_me(:,:) = 0._wp   ;   diag_sur_me(:,:) = 0._wp 
    326          diag_res_pr(:,:) = 0._wp   ;   diag_trp_vi(:,:) = 0._wp 
     310         d_u_ice_dyn(:,:)     = 0._wp   ;   d_v_ice_dyn(:,:)     = 0._wp 
     311 
     312         ! salt, heat and mass fluxes 
     313         sfx    (:,:) = 0._wp   ; 
     314         sfx_bri(:,:) = 0._wp   ;   sfx_dyn(:,:) = 0._wp  
     315         sfx_sni(:,:) = 0._wp   ;   sfx_opw(:,:) = 0._wp 
     316         sfx_bog(:,:) = 0._wp   ;   sfx_dyn(:,:) = 0._wp 
     317         sfx_bom(:,:) = 0._wp   ;   sfx_sum(:,:) = 0._wp 
     318         sfx_res(:,:) = 0._wp 
     319 
     320         wfx_snw(:,:) = 0._wp   ;   wfx_ice(:,:) = 0._wp 
     321         wfx_sni(:,:) = 0._wp   ;   wfx_opw(:,:) = 0._wp 
     322         wfx_bog(:,:) = 0._wp   ;   wfx_dyn(:,:) = 0._wp 
     323         wfx_bom(:,:) = 0._wp   ;   wfx_sum(:,:) = 0._wp 
     324         wfx_res(:,:) = 0._wp   ;   wfx_sub(:,:) = 0._wp 
     325         wfx_spr(:,:) = 0._wp   ;    
     326 
     327         hfx_in (:,:) = 0._wp   ;   hfx_out(:,:) = 0._wp 
     328         hfx_thd(:,:) = 0._wp   ;    
     329         hfx_snw(:,:) = 0._wp   ;   hfx_opw(:,:) = 0._wp 
     330         hfx_bog(:,:) = 0._wp   ;   hfx_dyn(:,:) = 0._wp 
     331         hfx_bom(:,:) = 0._wp   ;   hfx_sum(:,:) = 0._wp 
     332         hfx_res(:,:) = 0._wp   ;   hfx_sub(:,:) = 0._wp 
     333         hfx_spr(:,:) = 0._wp   ;   hfx_dif(:,:) = 0._wp  
     334         hfx_err(:,:) = 0._wp   ;   hfx_err_rem(:,:) = 0._wp 
     335 
     336         ! 
     337         fhld  (:,:)    = 0._wp  
     338         fmmflx(:,:)    = 0._wp      
     339         ! part of solar radiation transmitted through the ice 
     340         ftr_ice(:,:,:) = 0._wp 
     341 
     342         ! diags 
     343         diag_trp_vi  (:,:) = 0._wp  ; diag_trp_vs(:,:) = 0._wp  ;  diag_trp_ei(:,:) = 0._wp  ;  diag_trp_es(:,:) = 0._wp 
     344         diag_heat_dhc(:,:) = 0._wp   
     345 
    327346         ! dynamical invariants 
    328347         delta_i(:,:) = 0._wp       ;   divu_i(:,:) = 0._wp       ;   shear_i(:,:) = 0._wp 
     
    375394                          zcoef = rdt_ice /rday           !  Ice natural aging 
    376395                          oa_i(:,:,:) = oa_i(:,:,:) + a_i(:,:,:) * zcoef 
    377                           CALL lim_var_glo2eqv            ! this CALL is maybe not necessary (Martin) 
    378396         IF( ln_nicep )   CALL lim_prt_state( kt, jiindx, jjindx, 1, ' - ice thermodyn. - ' )   ! control print 
    379397                          CALL lim_itd_th( kt )           !  Remap ice categories, lateral accretion  ! 
     
    391409         !                                           ! Diagnostics and outputs  
    392410         IF (ln_limdiaout) CALL lim_diahsb 
    393 !clem # if ! defined key_iomput 
     411 
    394412                          CALL lim_wri( 1  )              ! Ice outputs  
    395 !clem # endif 
     413 
    396414         IF( kt == nit000 .AND. ln_rstart )   & 
    397415            &             CALL iom_close( numrir )        ! clem: close input ice restart file 
     
    413431       
    414432!!gm   remark, the ocean-ice stress is not saved in ice diag call above .....  find a solution!!! 
    415       ! 
    416       CALL wrk_dealloc( jpi,jpj,jpl, zalb_ice_os, zalb_ice_cs ) 
    417  
    418 #if defined key_coupled 
    419       IF ( ln_cpl .OR. ln_iceflx_ave .OR. ln_iceflx_linear ) CALL wrk_dealloc( jpi,jpj,jpl, zalb_ice) 
    420       IF ( ln_iceflx_ave .OR. ln_iceflx_linear ) & 
    421          &    CALL wrk_dealloc( jpi,jpj, ztem_ice_all, zalb_ice_all, z_qsr_ice_all, z_qns_ice_all, z_qla_ice_all, z_dqns_ice_all, z_dqla_ice_all) 
    422 #endif 
     433      CALL wrk_dealloc( jpi,jpj,jpl, zalb_ice_os, zalb_ice_cs, zalb_ice ) 
     434 
     435      IF( lk_cpl ) THEN 
     436         IF ( ln_iceflx_ave .OR. ln_iceflx_linear ) & 
     437            &    CALL wrk_dealloc( jpi,jpj, ztem_ice_all, zalb_ice_all, z_qsr_ice_all, z_qns_ice_all, z_qla_ice_all, z_dqns_ice_all, z_dqla_ice_all) 
     438      ENDIF 
    423439      ! 
    424440      IF( nn_timing == 1 )  CALL timing_stop('sbc_ice_lim') 
     
    534550!                 WRITE(numout,*) ' sst                  : ', sst_m(ji,jj) 
    535551!                 WRITE(numout,*) ' sss                  : ', sss_m(ji,jj) 
    536 !                 WRITE(numout,*) ' s_i_newice           : ', s_i_newice(ji,jj,1:jpl) 
    537552!                 WRITE(numout,*)  
    538553                  inb_alp(ialert_id) = inb_alp(ialert_id) + 1 
     
    591606               !WRITE(numout,*) ' sst       : ', sst_m(ji,jj) 
    592607               !WRITE(numout,*) ' sss       : ', sss_m(ji,jj) 
    593                !WRITE(numout,*) ' qcmif     : ', qcmif(ji,jj) 
    594                !WRITE(numout,*) ' qldif     : ', qldif(ji,jj) 
    595                !WRITE(numout,*) ' qcmif     : ', qcmif(ji,jj) / rdt_ice 
    596                !WRITE(numout,*) ' qldif     : ', qldif(ji,jj) / rdt_ice 
    597                !WRITE(numout,*) ' qfvbq     : ', qfvbq(ji,jj) 
    598                !WRITE(numout,*) ' qdtcn     : ', qdtcn(ji,jj) 
    599                !WRITE(numout,*) ' qfvbq / dt: ', qfvbq(ji,jj) / rdt_ice 
    600                !WRITE(numout,*) ' qdtcn / dt: ', qdtcn(ji,jj) / rdt_ice 
    601                !WRITE(numout,*) ' fdtcn     : ', fdtcn(ji,jj)  
    602                !WRITE(numout,*) ' fhmec     : ', fhmec(ji,jj)  
    603                !WRITE(numout,*) ' fheat_mec : ', fheat_mec(ji,jj)  
    604                !WRITE(numout,*) ' fheat_res : ', fheat_res(ji,jj)  
    605                !WRITE(numout,*) ' fhbri     : ', fhbri(ji,jj)  
    606608               ! 
    607609               !CALL lim_prt_state( kt, ji, jj, 2, '   ') 
     
    790792               WRITE(numout,*) ' - Heat / FW fluxes ' 
    791793               WRITE(numout,*) '   ~~~~~~~~~~~~~~~~ ' 
    792                WRITE(numout,*) ' emp        : ', emp      (ji,jj) 
    793                WRITE(numout,*) ' sfx        : ', sfx      (ji,jj) 
    794                WRITE(numout,*) ' sfx_thd    : ', sfx_thd(ji,jj) 
    795                WRITE(numout,*) ' sfx_bri    : ', sfx_bri  (ji,jj) 
    796                WRITE(numout,*) ' sfx_mec    : ', sfx_mec  (ji,jj) 
    797                WRITE(numout,*) ' sfx_res    : ', sfx_res(ji,jj) 
    798                WRITE(numout,*) ' fmmec      : ', fmmec    (ji,jj) 
    799                WRITE(numout,*) ' fhmec      : ', fhmec    (ji,jj) 
    800                WRITE(numout,*) ' fhbri      : ', fhbri    (ji,jj) 
    801                WRITE(numout,*) ' fheat_mec  : ', fheat_mec(ji,jj) 
     794               WRITE(numout,*) ' - Heat fluxes in and out the ice ***' 
     795               WRITE(numout,*) ' qsr_ini       : ', pfrld(ji,jj) * qsr(ji,jj) + SUM( old_a_i(ji,jj,:) * qsr_ice(ji,jj,:) ) 
     796               WRITE(numout,*) ' qns_ini       : ', pfrld(ji,jj) * qns(ji,jj) + SUM( old_a_i(ji,jj,:) * qns_ice(ji,jj,:) ) 
     797               WRITE(numout,*) 
    802798               WRITE(numout,*)  
    803799               WRITE(numout,*) ' sst        : ', sst_m(ji,jj)   
     
    829825               WRITE(numout,*) ' qsr       : ', qsr(ji,jj) 
    830826               WRITE(numout,*) ' qns       : ', qns(ji,jj) 
    831                WRITE(numout,*) ' fdtcn     : ', fdtcn(ji,jj) 
    832                WRITE(numout,*) ' qcmif     : ', qcmif(ji,jj) * r1_rdtice 
    833                WRITE(numout,*) ' qldif     : ', qldif(ji,jj) * r1_rdtice 
     827               WRITE(numout,*) 
     828               WRITE(numout,*) ' hfx_mass     : ', hfx_thd(ji,jj) + hfx_dyn(ji,jj) + hfx_snw(ji,jj) + hfx_res(ji,jj) 
     829               WRITE(numout,*) ' hfx_in       : ', hfx_in(ji,jj) 
     830               WRITE(numout,*) ' hfx_out      : ', hfx_out(ji,jj) 
     831               WRITE(numout,*) ' dhc          : ', diag_heat_dhc(ji,jj)               
     832               WRITE(numout,*) 
     833               WRITE(numout,*) ' hfx_dyn      : ', hfx_dyn(ji,jj) 
     834               WRITE(numout,*) ' hfx_thd      : ', hfx_thd(ji,jj) 
     835               WRITE(numout,*) ' hfx_res      : ', hfx_res(ji,jj) 
     836               WRITE(numout,*) ' fhtur        : ', fhtur(ji,jj)  
     837               WRITE(numout,*) ' qlead        : ', qlead(ji,jj) * r1_rdtice 
    834838               WRITE(numout,*) 
    835839               WRITE(numout,*) ' - Salt fluxes at bottom interface ***' 
    836840               WRITE(numout,*) ' emp       : ', emp    (ji,jj) 
    837                WRITE(numout,*) ' sfx_bri   : ', sfx_bri(ji,jj) 
    838841               WRITE(numout,*) ' sfx       : ', sfx    (ji,jj) 
    839842               WRITE(numout,*) ' sfx_res   : ', sfx_res(ji,jj) 
    840                WRITE(numout,*) ' sfx_mec   : ', sfx_mec(ji,jj) 
    841                WRITE(numout,*) ' - Heat fluxes at bottom interface ***' 
    842                WRITE(numout,*) ' fheat_res : ', fheat_res(ji,jj) 
     843               WRITE(numout,*) ' sfx_bri   : ', sfx_bri(ji,jj) 
     844               WRITE(numout,*) ' sfx_dyn   : ', sfx_dyn(ji,jj) 
    843845               WRITE(numout,*) 
    844846               WRITE(numout,*) ' - Momentum fluxes ' 
    845847               WRITE(numout,*) ' utau      : ', utau(ji,jj)  
    846848               WRITE(numout,*) ' vtau      : ', vtau(ji,jj) 
    847             ENDIF 
     849            ENDIF  
    848850            WRITE(numout,*) ' ' 
    849851            ! 
Note: See TracChangeset for help on using the changeset viewer.