Changeset 6507


Ignore:
Timestamp:
2016-05-03T14:28:12+02:00 (6 years ago)
Author:
timgraham
Message:

First attempt at merging in science changes from GO6 package branch at v3.6 stable (Note-namelists not yet dealt with)

Location:
branches/UKMO/dev_r6501_GO6_package_trunk/NEMOGCM/NEMO/OPA_SRC
Files:
20 edited

Legend:

Unmodified
Added
Removed
  • branches/UKMO/dev_r6501_GO6_package_trunk/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90

    r6503 r6507  
    3939   USE zdfmxl          ! mixed layer 
    4040   USE dianam          ! build name of file (routine) 
     41   USE zdftke          ! vertical physics: one-equation scheme  
    4142   USE zdfddm          ! vertical  physics: double diffusion 
    4243   USE diahth          ! thermocline diagnostics 
     
    231232      CALL iom_put( "avt" , avt                        )    ! T vert. eddy diff. coef. 
    232233      CALL iom_put( "avm" , avmu                       )    ! T vert. eddy visc. coef. 
     234      IF( lk_zdftke ) THEN    
     235         CALL iom_put( "tke"      , en                               )    ! TKE budget: Turbulent Kinetic Energy    
     236         CALL iom_put( "tke_niw"  , e_niw                            )    ! TKE budget: Near-inertial waves    
     237      ENDIF  
    233238      CALL iom_put( "avs" , fsavs(:,:,:)               )    ! S vert. eddy diff. coef. (useful only with key_zdfddm) 
     239                                                            ! Log of eddy diff coef 
     240      IF( iom_use('logavt') )   CALL iom_put( "logavt", LOG( MAX( 1.e-20_wp, avt  (:,:,:) ) ) ) 
     241      IF( iom_use('logavs') )   CALL iom_put( "logavs", LOG( MAX( 1.e-20_wp, fsavs(:,:,:) ) ) ) 
    234242 
    235243      IF( iom_use('logavt') )   CALL iom_put( "logavt", LOG( MAX( 1.e-20_wp, avt  (:,:,:) ) ) ) 
     
    974982      CALL histwrite( id_i, "sozotaux", kt, utau             , jpi*jpj    , idex )    ! i-wind stress 
    975983      CALL histwrite( id_i, "sometauy", kt, vtau             , jpi*jpj    , idex )    ! j-wind stress 
     984      IF( lk_vvl ) THEN 
     985         CALL histwrite( id_i, "vovvldep", kt, gdept_n(:,:,:), jpi*jpj*jpk, idex )!  T-cell depth        
     986         CALL histwrite( id_i, "vovvle3t", kt, e3t_n (:,:,:), jpi*jpj*jpk, idex )!  T-cell thickness   
     987      END IF 
    976988 
    977989      IF(  .NOT.ln_linssh  ) THEN              
  • branches/UKMO/dev_r6501_GO6_package_trunk/NEMOGCM/NEMO/OPA_SRC/DOM/closea.F90

    r6503 r6507  
    155155            ncsi2(4)   = 76  ;  ncsj2(4)   = 61 
    156156            ncsir(4,1) = 84  ;  ncsjr(4,1) = 59  
    157             !                                        ! ======================= 
    158          CASE ( 025 )                                ! ORCA_R025 configuration 
    159             !                                        ! ======================= 
    160             ncsnr(1)   = 1    ; ncstt(1)   = 0               ! Caspian + Aral sea 
    161             ncsi1(1)   = 1330 ; ncsj1(1)   = 645 
    162             ncsi2(1)   = 1400 ; ncsj2(1)   = 795 
     157            !                                        ! ================================ 
     158         CASE ( 025 )                                ! ORCA_R025 extended configuration 
     159            !                                        ! ================================ 
     160            ncsnr(1)   = 1    ; ncstt(1)   = 0               ! Caspian sea 
     161            ncsi1(1)   = 1330 ; ncsj1(1)   = 831 
     162            ncsi2(1)   = 1375 ; ncsj2(1)   = 981 
    163163            ncsir(1,1) = 1    ; ncsjr(1,1) = 1 
    164164            !                                         
    165             ncsnr(2)   = 1    ; ncstt(2)   = 0               ! Azov Sea  
    166             ncsi1(2)   = 1284 ; ncsj1(2)   = 722 
    167             ncsi2(2)   = 1304 ; ncsj2(2)   = 747 
     165            ncsnr(2)   = 1    ; ncstt(2)   = 0               ! Aral sea 
     166            ncsi1(2)   = 1376 ; ncsj1(2)   = 900 
     167            ncsi2(2)   = 1400 ; ncsj2(2)   = 981 
    168168            ncsir(2,1) = 1    ; ncsjr(2,1) = 1 
     169            !                                         
     170            ncsnr(3)   = 1    ; ncstt(3)   = 0               ! Azov Sea  
     171            ncsi1(3)   = 1284 ; ncsj1(3)   = 908 
     172            ncsi2(3)   = 1304 ; ncsj2(3)   = 933 
     173            ncsir(3,1) = 1    ; ncsjr(3,1) = 1 
     174            ! 
     175            ncsnr(4)   = 1    ; ncstt(4)   = 0               ! Lake Superior   
     176            ncsi1(4)   = 781  ; ncsj1(4)   = 904  
     177            ncsi2(4)   = 815  ; ncsj2(4)   = 926  
     178            ncsir(4,1) = 1    ; ncsjr(4,1) = 1  
     179            !  
     180            ncsnr(5)   = 1    ; ncstt(5)   = 0               ! Lake Michigan 
     181            ncsi1(5)   = 795  ; ncsj1(5)   = 871              
     182            ncsi2(5)   = 813  ; ncsj2(5)   = 905  
     183            ncsir(5,1) = 1    ; ncsjr(5,1) = 1  
     184            !  
     185            ncsnr(6)   = 1    ; ncstt(6)   = 0               ! Lake Huron part 1 
     186            ncsi1(6)   = 814  ; ncsj1(6)   = 882              
     187            ncsi2(6)   = 825  ; ncsj2(6)   = 905  
     188            ncsir(6,1) = 1    ; ncsjr(6,1) = 1  
     189            !  
     190            ncsnr(7)   = 1    ; ncstt(7)   = 0               ! Lake Huron part 2 
     191            ncsi1(7)   = 826  ; ncsj1(7)   = 889              
     192            ncsi2(7)   = 833  ; ncsj2(7)   = 905  
     193            ncsir(7,1) = 1    ; ncsjr(7,1) = 1  
     194            !  
     195            ncsnr(8)   = 1    ; ncstt(8)   = 0               ! Lake Erie 
     196            ncsi1(8)   = 816  ; ncsj1(8)   = 871              
     197            ncsi2(8)   = 837  ; ncsj2(8)   = 881  
     198            ncsir(8,1) = 1    ; ncsjr(8,1) = 1  
     199            !  
     200            ncsnr(9)   = 1    ; ncstt(9)   = 0               ! Lake Ontario 
     201            ncsi1(9)   = 831  ; ncsj1(9)   = 882              
     202            ncsi2(9)   = 847  ; ncsj2(9)   = 889  
     203            ncsir(9,1) = 1    ; ncsjr(9,1) = 1  
     204            !  
     205            ncsnr(10)   = 1    ; ncstt(10)   = 0               ! Lake Victoria   
     206            ncsi1(10)   = 1274 ; ncsj1(10)   = 672  
     207            ncsi2(10)   = 1289 ; ncsj2(10)   = 687  
     208            ncsir(10,1) = 1    ; ncsjr(10,1) = 1  
    169209            ! 
    170210         END SELECT 
  • branches/UKMO/dev_r6501_GO6_package_trunk/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90

    r6503 r6507  
    161161      USE ioipsl 
    162162      NAMELIST/namrun/ cn_ocerst_indir, cn_ocerst_outdir, nn_stocklist, ln_rst_list,                 & 
    163                        nn_no   , cn_exp   , cn_ocerst_in, cn_ocerst_out, ln_rstart , nn_rstctl ,     & 
     163         &             nn_no   , cn_exp    , cn_ocerst_in, cn_ocerst_out, ln_rstart , ln_rstdate, nn_rstctl,   & 
    164164         &             nn_it000, nn_itend , nn_date0    , nn_time0     , nn_leapy  , nn_istate ,     & 
    165165         &             nn_stock, nn_write , ln_mskland  , ln_clobber   , nn_chunksz, nn_euler  ,     & 
     
    196196         WRITE(numout,*) '      file prefix restart output      cn_ocerst_out= ', cn_ocerst_out 
    197197         WRITE(numout,*) '      restart output directory        cn_ocerst_outdir= ', cn_ocerst_outdir 
    198          WRITE(numout,*) '      restart logical                 ln_rstart  = ', ln_rstart 
     198         WRITE(numout,*) '      restart logical                 ln_rstart  = ' , ln_rstart 
     199         WRITE(numout,*) '      datestamping of restarts        ln_rstdate  = ', ln_rstdate 
    199200         WRITE(numout,*) '      start with forward time step    nn_euler   = ', nn_euler 
    200201         WRITE(numout,*) '      control of time step            nn_rstctl  = ', nn_rstctl 
  • branches/UKMO/dev_r6501_GO6_package_trunk/NEMOGCM/NEMO/OPA_SRC/DOM/dommsk.F90

    r6503 r6507  
    105105      INTEGER  ::   isrow                    ! index for ORCA1 starting row 
    106106      INTEGER , POINTER, DIMENSION(:,:) ::  imsk 
     107      REAL(wp) ::  zphi_drake_passage, zshlat_antarc 
    107108      REAL(wp), POINTER, DIMENSION(:,:) ::  zwf 
    108109      !! 
     
    345346      ENDIF 
    346347      ! 
     348      IF( cp_cfg == "orca" .AND. jp_cfg == 025 .AND. rn_shlat == 0.0 ) THEN    
     349         !                                              ! ORCA_R025 configuration 
     350         !                                              ! Increased lateral friction on parts of Antarctic coastline 
     351         !                                              ! for increased stability 
     352         !                                              ! NB. This only works to do this here if we have free slip  
     353         !                                              ! generally, so fmask is zero at coast points. 
     354         IF(lwp) WRITE(numout,*) 
     355         IF(lwp) WRITE(numout,*) '   orca_r025: increase friction in following regions : ' 
     356         IF(lwp) WRITE(numout,*) '      whole Antarctic coastline: partial slip shlat=1 ' 
     357 
     358         zphi_drake_passage = -58.0_wp 
     359         zshlat_antarc = 1.0_wp 
     360         zwf(:,:) = fmask(:,:,1)          
     361         DO jj = 2, jpjm1 
     362            DO ji = fs_2, fs_jpim1   ! vector opt. 
     363               IF( gphif(ji,jj) .lt. zphi_drake_passage .and. fmask(ji,jj,1) == 0._wp ) THEN 
     364                  fmask(ji,jj,:) = zshlat_antarc * MIN( 1._wp , MAX( zwf(ji+1,jj), zwf(ji,jj+1),   & 
     365                     &                                           zwf(ji-1,jj), zwf(ji,jj-1)  )  ) 
     366               ENDIF 
     367            END DO 
     368         END DO 
     369      END IF 
     370      ! 
    347371      CALL lbc_lnk( fmask, 'F', 1._wp )      ! Lateral boundary conditions on fmask 
    348372      ! 
  • branches/UKMO/dev_r6501_GO6_package_trunk/NEMOGCM/NEMO/OPA_SRC/DOM/phycst.F90

    r6503 r6507  
    6868   REAL(wp), PUBLIC ::   rcdsn    =    0.31_wp       !: thermal conductivity of snow 
    6969   REAL(wp), PUBLIC ::   cpic     = 2067.0_wp        !: specific heat for ice  
     70#if defined key_cice 
     71   REAL(wp), PUBLIC ::   lsub     =    2.835e+6_wp   !: pure ice latent heat of sublimation                   [J/kg] 
     72#else 
    7073   REAL(wp), PUBLIC ::   lsub     =    2.834e+6_wp   !: pure ice latent heat of sublimation                   [J/kg] 
     74#endif 
    7175   REAL(wp), PUBLIC ::   lfus     =    0.334e+6_wp   !: latent heat of fusion of fresh ice                    [J/kg] 
    7276   REAL(wp), PUBLIC ::   tmut     =    0.054_wp      !: decrease of seawater meltpoint with salinity 
  • branches/UKMO/dev_r6501_GO6_package_trunk/NEMOGCM/NEMO/OPA_SRC/ICB/icbclv.F90

    r6503 r6507  
    2525   USE icbutl         ! iceberg utility routines 
    2626 
     27   USE sbc_oce        ! for icesheet freshwater input variables 
     28   USE in_out_manager 
     29   USE iom 
     30 
    2731   IMPLICIT NONE 
    2832   PRIVATE 
     
    4852      ! 
    4953      REAL(wp)                        :: zcalving_used, zdist, zfact 
     54      REAL(wp)                        :: zgreenland_calving_sum, zantarctica_calving_sum 
    5055      INTEGER                         :: jn, ji, jj                    ! loop counters 
    5156      INTEGER                         :: imx                           ! temporary integer for max berg class 
     
    5964      zfact = ( (1000._wp)**3 / ( NINT(rday) * nyear_len(1) ) ) * 850._wp 
    6065      berg_grid%calving(:,:) = src_calving(:,:) * tmask_i(:,:) * zfact 
     66 
     67      IF( lk_oasis) THEN 
     68      ! ln_coupled_iceshelf_fluxes uninitialised unless lk_oasis=true 
     69      IF( ln_coupled_iceshelf_fluxes ) THEN 
     70 
     71        ! Adjust total calving rates so that sum of iceberg calving and iceshelf melting in the northern 
     72        ! and southern hemispheres equals rate of increase of mass of greenland and antarctic ice sheets 
     73        ! to preserve total freshwater conservation in coupled models without an active ice sheet model. 
     74 
     75         zgreenland_calving_sum = SUM( berg_grid%calving(:,:) * greenland_icesheet_mask(:,:) ) 
     76         IF( lk_mpp ) CALL mpp_sum( zgreenland_calving_sum ) 
     77         WHERE( greenland_icesheet_mask(:,:) == 1.0 )                                                                                 & 
     78        &    berg_grid%calving(:,:) = berg_grid%calving(:,:) * greenland_icesheet_mass_rate_of_change * rn_greenland_calving_fraction & 
     79        &                                     / ( zgreenland_calving_sum + 1.0e-10_wp ) 
     80 
     81         ! check 
     82         IF(lwp) WRITE(numout, *) 'Greenland iceberg calving climatology (kg/s) : ',zgreenland_calving_sum 
     83         zgreenland_calving_sum = SUM( berg_grid%calving(:,:) * greenland_icesheet_mask(:,:) ) 
     84         IF( lk_mpp ) CALL mpp_sum( zgreenland_calving_sum ) 
     85         IF(lwp) WRITE(numout, *) 'Greenland iceberg calving adjusted value (kg/s) : ',zgreenland_calving_sum 
     86 
     87         zantarctica_calving_sum = SUM( berg_grid%calving(:,:) * antarctica_icesheet_mask(:,:) ) 
     88         IF( lk_mpp ) CALL mpp_sum( zantarctica_calving_sum ) 
     89         WHERE( antarctica_icesheet_mask(:,:) == 1.0 )                                                                              & 
     90         berg_grid%calving(:,:) = berg_grid%calving(:,:) * antarctica_icesheet_mass_rate_of_change * rn_antarctica_calving_fraction & 
     91        &                           / ( zantarctica_calving_sum + 1.0e-10_wp ) 
     92  
     93         ! check 
     94         IF(lwp) WRITE(numout, *) 'Antarctica iceberg calving climatology (kg/s) : ',zantarctica_calving_sum 
     95         zantarctica_calving_sum = SUM( berg_grid%calving(:,:) * antarctica_icesheet_mask(:,:) ) 
     96         IF( lk_mpp ) CALL mpp_sum( zantarctica_calving_sum ) 
     97         IF(lwp) WRITE(numout, *) 'Antarctica iceberg calving adjusted value (kg/s) : ',zantarctica_calving_sum 
     98 
     99      ENDIF 
     100      ENDIF 
     101    
     102      CALL iom_put( 'berg_calve', berg_grid%calving(:,:) ) 
    61103 
    62104      ! Heat in units of W/m2, and mask (just in case) 
  • branches/UKMO/dev_r6501_GO6_package_trunk/NEMOGCM/NEMO/OPA_SRC/ICB/icbdia.F90

    r6503 r6507  
    371371      IF( .NOT. ln_bergdia )   RETURN            !!gm useless iom will control whether it is output or not 
    372372      ! 
     373      CALL iom_put( "berg_total_melt"      , berg_grid%floating_melt(:,:)   )   ! Total melt flux to ocean      [kg/m2/s] 
     374      CALL iom_put( "berg_total_heat_flux" , berg_grid%calving_hflx(:,:)    )   ! Total iceberg-ocean heat flux [W/m2] 
    373375      CALL iom_put( "berg_melt"        , berg_melt   (:,:)   )   ! Melt rate of icebergs                     [kg/m2/s] 
    374376      CALL iom_put( "berg_buoy_melt"   , buoy_melt   (:,:)   )   ! Buoyancy component of iceberg melt rate   [kg/m2/s] 
  • branches/UKMO/dev_r6501_GO6_package_trunk/NEMOGCM/NEMO/OPA_SRC/ICB/icbrst.F90

    r6503 r6507  
    158158      INTEGER ::   jn   ! dummy loop index 
    159159      INTEGER ::   ix_dim, iy_dim, ik_dim, in_dim 
    160       CHARACTER(len=256)     :: cl_path 
    161       CHARACTER(len=256)     :: cl_filename 
     160      INTEGER             ::   iyear, imonth, iday 
     161      REAL (wp)           ::   zsec 
     162      CHARACTER(len=256)  :: cl_path 
     163      CHARACTER(len=256)  :: cl_filename 
     164      CHARACTER(LEN=20)   ::   clkt     ! ocean time-step deine as a character 
    162165      TYPE(iceberg), POINTER :: this 
    163166      TYPE(point)  , POINTER :: pt 
     
    168171      IF( cl_path(LEN_TRIM(cl_path):) /= '/' ) cl_path = TRIM(cl_path) // '/' 
    169172      IF( lk_mpp ) THEN 
    170          WRITE(cl_filename,'(A,"_icebergs_",I8.8,"_restart_",I4.4,".nc")') TRIM(cexper), kt, narea-1 
     173         WRITE(cl_filename,'(A,"_icebergs_",A,"_restart_",I4.4,".nc")') TRIM(cexper), TRIM(ADJUSTL(clkt)), narea-1 
    171174      ELSE 
    172          WRITE(cl_filename,'(A,"_icebergs_",I8.8,"_restart.nc")') TRIM(cexper), kt 
     175         WRITE(cl_filename,'(A,"_icebergs_",A,"_restart.nc")') TRIM(cexper), TRIM(ADJUSTL(clkt)) 
    173176      ENDIF 
    174177      IF (nn_verbose_level >= 0) WRITE(numout,'(2a)') 'icebergs, write_restart: creating ',TRIM(cl_path)//TRIM(cl_filename) 
  • branches/UKMO/dev_r6501_GO6_package_trunk/NEMOGCM/NEMO/OPA_SRC/ICB/icbthm.F90

    r6503 r6507  
    1818   USE dom_oce        ! NEMO domain 
    1919   USE in_out_manager ! NEMO IO routines, numout in particular 
     20   USE iom 
    2021   USE lib_mpp        ! NEMO MPI routines, ctl_stop in particular 
    2122   USE phycst         ! NEMO physical constants 
     
    160161            zmelt    = ( zdM - ( zdMbitsE - zdMbitsM ) ) * z1_dt   ! kg/s 
    161162            berg_grid%floating_melt(ii,ij) = berg_grid%floating_melt(ii,ij) + zmelt    * z1_e1e2    ! kg/m2/s 
    162             zheat = zmelt * pt%heat_density              ! kg/s x J/kg = J/s 
     163!            zheat = zmelt * pt%heat_density              ! kg/s x J/kg = J/s 
     164            zheat = zmelt * lfus                           !rma kg/s x J/kg (latent heat of fusion) = J/s 
    163165            berg_grid%calving_hflx (ii,ij) = berg_grid%calving_hflx (ii,ij) + zheat    * z1_e1e2    ! W/m2 
    164166            CALL icb_dia_melt( ii, ij, zMnew, zheat, this%mass_scaling,       & 
     
    208210      IF(.NOT. ln_passive_mode ) THEN 
    209211         emp (:,:) = emp (:,:) - berg_grid%floating_melt(:,:) 
    210 !!       qns (:,:) = qns (:,:) + berg_grid%calving_hflx (:,:)  !!gm heat flux not yet properly coded ==>> need it, SOLVE that! 
     212         qns (:,:) = qns (:,:) - berg_grid%calving_hflx (:,:)   
    211213      ENDIF 
    212214      ! 
  • branches/UKMO/dev_r6501_GO6_package_trunk/NEMOGCM/NEMO/OPA_SRC/IOM/in_out_manager.F90

    r6503 r6507  
    2929   CHARACTER(lc) ::   cn_ocerst_outdir !: restart output directory 
    3030   LOGICAL       ::   ln_rstart        !: start from (F) rest or (T) a restart file 
     31   LOGICAL       ::   ln_rstdate       !: datestamping of restarts 
    3132   LOGICAL       ::   ln_rst_list      !: output restarts at list of times (T) or by frequency (F) 
    3233   INTEGER       ::   nn_no            !: job number 
  • branches/UKMO/dev_r6501_GO6_package_trunk/NEMOGCM/NEMO/OPA_SRC/IOM/restart.F90

    r6503 r6507  
    5757      !!---------------------------------------------------------------------- 
    5858      INTEGER, INTENT(in) ::   kt     ! ocean time-step 
     59      INTEGER             ::   iyear, imonth, iday 
     60      REAL (wp)           ::   zsec 
    5961      !! 
    6062      CHARACTER(LEN=20)   ::   clkt     ! ocean time-step deine as a character 
     
    8486      IF( kt == nitrst - 1 .OR. nstock == 1 .OR. ( kt == nitend .AND. .NOT. lrst_oce ) ) THEN 
    8587         IF( nitrst <= nitend .AND. nitrst > 0 ) THEN  
    86             ! beware of the format used to write kt (default is i8.8, that should be large enough...) 
    87             IF( nitrst > 999999999 ) THEN   ;   WRITE(clkt, *       ) nitrst 
    88             ELSE                            ;   WRITE(clkt, '(i8.8)') nitrst 
     88            IF ( ln_rstdate ) THEN 
     89               CALL ju2ymds( fjulday + rdttra(1) / rday, iyear, imonth, iday, zsec )            
     90               WRITE(clkt, '(i4.4,2i2.2)') iyear, imonth, iday 
     91            ELSE 
     92               ! beware of the format used to write kt (default is i8.8, that should be large enough...) 
     93               IF( nitrst > 999999999 ) THEN   ;   WRITE(clkt, *       ) nitrst 
     94               ELSE                            ;   WRITE(clkt, '(i8.8)') nitrst 
     95               ENDIF 
    8996            ENDIF 
    9097            ! create the file 
     
    155162      IF (ln_diurnal) CALL iom_rstput( kt, nitrst, numrow, 'Dsst', x_dsst    )   
    156163 
     164                     IF( lk_oasis) THEN 
     165                     ! ln_coupled_iceshelf_fluxes uninitialised unless lk_oasis=true 
     166                     IF( ln_coupled_iceshelf_fluxes ) THEN 
     167                        CALL iom_rstput( kt, nitrst, numrow, 'greenland_icesheet_mass', greenland_icesheet_mass ) 
     168                        CALL iom_rstput( kt, nitrst, numrow, 'greenland_icesheet_timelapsed', greenland_icesheet_timelapsed ) 
     169                        CALL iom_rstput( kt, nitrst, numrow, 'greenland_icesheet_mass_roc', greenland_icesheet_mass_rate_of_change ) 
     170                        CALL iom_rstput( kt, nitrst, numrow, 'antarctica_icesheet_mass', antarctica_icesheet_mass ) 
     171                        CALL iom_rstput( kt, nitrst, numrow, 'antarctica_icesheet_timelapsed', antarctica_icesheet_timelapsed ) 
     172                        CALL iom_rstput( kt, nitrst, numrow, 'antarctica_icesheet_mass_roc', antarctica_icesheet_mass_rate_of_change ) 
     173                     ENDIF 
     174                     ENDIF 
     175 
    157176      IF( kt == nitrst ) THEN 
    158177         CALL iom_close( numrow )     ! close the restart file (only at last time step) 
     
    252271      ENDIF 
    253272      ! 
     273      IF( iom_varid( numror, 'greenland_icesheet_mass', ldstop = .FALSE. ) > 0 )   THEN 
     274         CALL iom_get( numror, 'greenland_icesheet_mass', greenland_icesheet_mass ) 
     275         CALL iom_get( numror, 'greenland_icesheet_timelapsed', greenland_icesheet_timelapsed ) 
     276         CALL iom_get( numror, 'greenland_icesheet_mass_roc', greenland_icesheet_mass_rate_of_change ) 
     277      ELSE 
     278         greenland_icesheet_mass = 0.0  
     279         greenland_icesheet_mass_rate_of_change = 0.0  
     280         greenland_icesheet_timelapsed = 0.0 
     281      ENDIF 
     282      IF( iom_varid( numror, 'antarctica_icesheet_mass', ldstop = .FALSE. ) > 0 )   THEN 
     283         CALL iom_get( numror, 'antarctica_icesheet_mass', antarctica_icesheet_mass ) 
     284         CALL iom_get( numror, 'antarctica_icesheet_timelapsed', antarctica_icesheet_timelapsed ) 
     285         CALL iom_get( numror, 'antarctica_icesheet_mass_roc', antarctica_icesheet_mass_rate_of_change ) 
     286      ELSE 
     287         antarctica_icesheet_mass = 0.0  
     288         antarctica_icesheet_mass_rate_of_change = 0.0  
     289         antarctica_icesheet_timelapsed = 0.0 
     290      ENDIF 
     291 
    254292      IF( neuler == 0 ) THEN                                  ! Euler restart (neuler=0) 
    255293         tsb  (:,:,:,:) = tsn  (:,:,:,:)                             ! all before fields set to now values 
  • branches/UKMO/dev_r6501_GO6_package_trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_ice.F90

    r6503 r6507  
    102102   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   fr_iu              !: ice fraction at NEMO U point 
    103103   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   fr_iv              !: ice fraction at NEMO V point 
    104     
     104   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   sstfrz             !: sea surface freezing temperature 
     105   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   tsfc_ice           !: sea-ice surface skin temperature (on categories) 
     106   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   kn_ice             !: sea-ice surface layer thermal conductivity (on cats) 
     107 
    105108   ! variables used in the coupled interface 
    106109   INTEGER , PUBLIC, PARAMETER ::   jpl = ncat 
     
    153156 
    154157#if defined key_cice 
    155       ALLOCATE( qla_ice(jpi,jpj,1)    , qlw_ice(jpi,jpj,1)    , qsr_ice(jpi,jpj,1)    , & 
     158      ALLOCATE( qla_ice(jpi,jpj,ncat) , qlw_ice(jpi,jpj,1)    , qsr_ice(jpi,jpj,1)    , & 
    156159                wndi_ice(jpi,jpj)     , tatm_ice(jpi,jpj)     , qatm_ice(jpi,jpj)     , & 
    157160                wndj_ice(jpi,jpj)     , nfrzmlt(jpi,jpj)      , ss_iou(jpi,jpj)       , & 
    158161                ss_iov(jpi,jpj)       , fr_iu(jpi,jpj)        , fr_iv(jpi,jpj)        , & 
    159162                a_i(jpi,jpj,ncat)     , topmelt(jpi,jpj,ncat) , botmelt(jpi,jpj,ncat) , & 
    160                 STAT= ierr(1) ) 
    161       IF( ln_cpl )   ALLOCATE( u_ice(jpi,jpj)        , fr1_i0(jpi,jpj)       , tn_ice (jpi,jpj,1)    , & 
     163                sstfrz(jpi,jpj)       , STAT= ierr(1) ) 
     164   ! Alex West: Allocating tn_ice with 5 categories.  When NEMO is used with CICE, this variable 
     165   ! represents top layer ice temperature, which is multi-category. 
     166      IF( ln_cpl )   ALLOCATE( u_ice(jpi,jpj)        , fr1_i0(jpi,jpj)       , tn_ice (jpi,jpj,jpl)  , & 
    162167         &                     v_ice(jpi,jpj)        , fr2_i0(jpi,jpj)       , alb_ice(jpi,jpj,1)    , & 
    163168         &                     emp_ice(jpi,jpj)      , qns_ice(jpi,jpj,1)    , dqns_ice(jpi,jpj,1)   , & 
    164          &                     STAT= ierr(2) ) 
     169         &                     a_p(jpi,jpj,jpl)      , ht_p(jpi,jpj,jpl)     , tsfc_ice(jpi,jpj,jpl) , & 
     170         &                     kn_ice(jpi,jpj,jpl) ,    STAT=ierr(2) ) 
    165171       
    166172#endif 
  • branches/UKMO/dev_r6501_GO6_package_trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_oce.F90

    r6503 r6507  
    124124#endif 
    125125   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: xcplmask          !: coupling mask for ln_mixcpl (warning: allocated in sbccpl) 
     126   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   greenland_icesheet_mass_array, greenland_icesheet_mask 
     127   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   antarctica_icesheet_mass_array, antarctica_icesheet_mask 
    126128 
    127129   !!---------------------------------------------------------------------- 
     
    136138   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   e3t_m     !: mean (nn_fsbc time-step) sea surface layer thickness       [m] 
    137139   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   frq_m     !: mean (nn_fsbc time-step) fraction of solar net radiation absorbed in the 1st T level [-] 
     140    
     141   !!---------------------------------------------------------------------- 
     142   !!  Surface scalars of total ice sheet mass for Greenland and Antarctica,  
     143   !! passed from atmosphere to be converted to dvol and hence a freshwater  
     144   !! flux  by using old values. New values are saved in the dump, to become 
     145   !! old values next coupling timestep. Freshwater fluxes split between  
     146   !! sub iceshelf melting and iceberg calving, scalled to flux per second 
     147   !!---------------------------------------------------------------------- 
     148    
     149   REAL(wp), PUBLIC  :: greenland_icesheet_mass, greenland_icesheet_mass_rate_of_change, greenland_icesheet_timelapsed  
     150   REAL(wp), PUBLIC  :: antarctica_icesheet_mass, antarctica_icesheet_mass_rate_of_change, antarctica_icesheet_timelapsed 
     151 
     152   ! sbccpl namelist parameters associated with icesheet freshwater input code. Included here rather than in sbccpl.F90 to  
     153   ! avoid circular dependencies. 
     154   LOGICAL, PUBLIC     ::   ln_coupled_iceshelf_fluxes     ! If true use rate of change of mass of Greenland and Antarctic icesheets to set the  
     155                                                           ! combined magnitude of the iceberg calving and iceshelf melting freshwater fluxes. 
     156   REAL(wp), PUBLIC    ::   rn_greenland_calving_fraction  ! Set fraction of total freshwater flux for iceberg calving - remainder goes to iceshelf melting. 
     157   REAL(wp), PUBLIC    ::   rn_antarctica_calving_fraction ! Set fraction of total freshwater flux for iceberg calving - remainder goes to iceshelf melting. 
     158   REAL(wp), PUBLIC    ::   rn_iceshelf_fluxes_tolerance   ! Absolute tolerance for detecting differences in icesheet masses.  
    138159 
    139160   !! * Substitutions 
     
    171192         &      ssu_m  (jpi,jpj) , sst_m(jpi,jpj) , frq_m(jpi,jpj) ,      & 
    172193         &      ssv_m  (jpi,jpj) , sss_m(jpi,jpj) , ssh_m(jpi,jpj) , STAT=ierr(4) ) 
     194      ALLOCATE( greenland_icesheet_mass_array(jpi,jpj) , antarctica_icesheet_mass_array(jpi,jpj) ) 
     195      ALLOCATE( greenland_icesheet_mask(jpi,jpj) , antarctica_icesheet_mask(jpi,jpj) ) 
    173196         ! 
    174197      ALLOCATE( e3t_m(jpi,jpj) , STAT=ierr(5) ) 
  • branches/UKMO/dev_r6501_GO6_package_trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90

    r6503 r6507  
    106106   INTEGER, PARAMETER ::   jpr_e3t1st = 41   ! first T level thickness  
    107107   INTEGER, PARAMETER ::   jpr_fraqsr = 42   ! fraction of solar net radiation absorbed in the first ocean level 
    108    INTEGER, PARAMETER ::   jprcv      = 42   ! total number of fields received 
     108   INTEGER, PARAMETER ::   jpr_ts_ice = 43   ! skin temperature of sea-ice (used for melt-ponds) 
     109   INTEGER, PARAMETER ::   jpr_grnm   = 44   ! Greenland ice mass 
     110   INTEGER, PARAMETER ::   jpr_antm   = 45   ! Antarctic ice mass 
     111   INTEGER, PARAMETER ::   jprcv      = 45   ! total number of fields received 
    109112 
    110113   INTEGER, PARAMETER ::   jps_fice   =  1   ! ice fraction sent to the atmosphere 
     
    136139   INTEGER, PARAMETER ::   jps_e3t1st = 27   ! first level depth (vvl) 
    137140   INTEGER, PARAMETER ::   jps_fraqsr = 28   ! fraction of solar net radiation absorbed in the first ocean level 
    138    INTEGER, PARAMETER ::   jpsnd      = 28   ! total number of fields sended 
     141   INTEGER, PARAMETER ::   jps_a_p    = 29   ! meltpond fraction   
     142   INTEGER, PARAMETER ::   jps_ht_p   = 30   ! meltpond depth (m)  
     143   INTEGER, PARAMETER ::   jps_kice   = 31   ! ice surface layer thermal conductivity 
     144   INTEGER, PARAMETER ::   jps_sstfrz = 32   ! sea-surface freezing temperature 
     145   INTEGER, PARAMETER ::   jps_fice1  = 33   ! first-order ice concentration (for time-travelling ice coupling) 
     146   INTEGER, PARAMETER ::   jpsnd      = 33   ! total number of fields sended 
    139147 
    140148   !                                  !!** namelist namsbc_cpl ** 
     
    147155   END TYPE FLD_C 
    148156   !                                   ! Send to the atmosphere   
    149    TYPE(FLD_C) ::   sn_snd_temp, sn_snd_alb, sn_snd_thick, sn_snd_crt, sn_snd_co2                         
     157   TYPE(FLD_C) ::   sn_snd_temp, sn_snd_alb, sn_snd_thick, sn_snd_crt, sn_snd_co2, sn_snd_cond, sn_snd_mpnd, sn_snd_sstfrz, sn_snd_thick1 
     158 
    150159   !                                   ! Received from the atmosphere 
    151160   TYPE(FLD_C) ::   sn_rcv_w10m, sn_rcv_taumod, sn_rcv_tau, sn_rcv_dqnsdt, sn_rcv_qsr, sn_rcv_qns, sn_rcv_emp, sn_rcv_rnf 
    152    TYPE(FLD_C) ::   sn_rcv_cal, sn_rcv_iceflx, sn_rcv_co2                         
     161   TYPE(FLD_C) ::   sn_rcv_cal, sn_rcv_iceflx, sn_rcv_co2, sn_rcv_ts_ice, sn_rcv_grnm, sn_rcv_antm 
    153162   !                                   ! Other namelist parameters 
    154163   INTEGER     ::   nn_cplmodel           ! Maximum number of models to/from which NEMO is potentialy sending/receiving data 
     
    214223      REAL(wp), POINTER, DIMENSION(:,:) ::   zacs, zaos 
    215224      !! 
    216       NAMELIST/namsbc_cpl/  sn_snd_temp, sn_snd_alb   , sn_snd_thick, sn_snd_crt   , sn_snd_co2,      & 
    217          &                  sn_rcv_w10m, sn_rcv_taumod, sn_rcv_tau  , sn_rcv_dqnsdt, sn_rcv_qsr,      & 
    218          &                  sn_rcv_qns , sn_rcv_emp   , sn_rcv_rnf  , sn_rcv_cal   , sn_rcv_iceflx,   & 
    219          &                  sn_rcv_co2 , nn_cplmodel  , ln_usecplmask 
     225      NAMELIST/namsbc_cpl/  sn_snd_temp, sn_snd_alb   , sn_snd_thick , sn_snd_crt   , sn_snd_co2,     & 
     226         &                  sn_snd_cond, sn_snd_mpnd  , sn_snd_sstfrz, sn_snd_thick1,                 & 
     227         &                  sn_rcv_w10m, sn_rcv_taumod, sn_rcv_tau   , sn_rcv_dqnsdt, sn_rcv_qsr,     & 
     228         &                  sn_rcv_qns , sn_rcv_emp   , sn_rcv_rnf   , sn_rcv_cal   , sn_rcv_iceflx,  & 
     229         &                  sn_rcv_co2 , sn_rcv_grnm  , sn_rcv_antm  , sn_rcv_ts_ice, nn_cplmodel  ,  & 
     230         &                  ln_usecplmask, ln_coupled_iceshelf_fluxes, rn_greenland_calving_fraction, & 
     231         &                  rn_antarctica_calving_fraction, rn_iceshelf_fluxes_tolerance 
    220232      !!--------------------------------------------------------------------- 
    221233      ! 
     
    256268         WRITE(numout,*)'      runoffs                         = ', TRIM(sn_rcv_rnf%cldes   ), ' (', TRIM(sn_rcv_rnf%clcat   ), ')' 
    257269         WRITE(numout,*)'      calving                         = ', TRIM(sn_rcv_cal%cldes   ), ' (', TRIM(sn_rcv_cal%clcat   ), ')' 
     270         WRITE(numout,*)'      Greenland ice mass              = ', TRIM(sn_rcv_grnm%cldes  ), ' (', TRIM(sn_rcv_grnm%clcat  ), ')' 
     271         WRITE(numout,*)'      Antarctica ice mass             = ', TRIM(sn_rcv_antm%cldes  ), ' (', TRIM(sn_rcv_antm%clcat  ), ')' 
    258272         WRITE(numout,*)'      sea ice heat fluxes             = ', TRIM(sn_rcv_iceflx%cldes), ' (', TRIM(sn_rcv_iceflx%clcat), ')' 
    259273         WRITE(numout,*)'      atm co2                         = ', TRIM(sn_rcv_co2%cldes   ), ' (', TRIM(sn_rcv_co2%clcat   ), ')' 
     
    267281         WRITE(numout,*)'                      - mesh          = ', sn_snd_crt%clvgrd 
    268282         WRITE(numout,*)'      oce co2 flux                    = ', TRIM(sn_snd_co2%cldes   ), ' (', TRIM(sn_snd_co2%clcat   ), ')' 
     283         WRITE(numout,*)'      ice effective conductivity      = ', TRIM(sn_snd_cond%cldes   ), ' (', TRIM(sn_snd_cond%clcat   ), ')' 
     284         WRITE(numout,*)'      meltponds fraction & depth      = ', TRIM(sn_snd_mpnd%cldes  ), ' (', TRIM(sn_snd_mpnd%clcat   ), ')' 
     285         WRITE(numout,*)'      sea surface freezing temp       = ', TRIM(sn_snd_sstfrz%cldes   ), ' (', TRIM(sn_snd_sstfrz%clcat   ), ')' 
     286 
    269287         WRITE(numout,*)'  nn_cplmodel                         = ', nn_cplmodel 
    270288         WRITE(numout,*)'  ln_usecplmask                       = ', ln_usecplmask 
     289         WRITE(numout,*)'  ln_coupled_iceshelf_fluxes          = ', ln_coupled_iceshelf_fluxes 
     290         WRITE(numout,*)'  rn_greenland_calving_fraction       = ', rn_greenland_calving_fraction 
     291         WRITE(numout,*)'  rn_antarctica_calving_fraction      = ', rn_antarctica_calving_fraction 
     292         WRITE(numout,*)'  rn_iceshelf_fluxes_tolerance        = ', rn_iceshelf_fluxes_tolerance 
    271293      ENDIF 
    272294 
     
    381403      srcv(jpr_snow)%clname = 'OTotSnow'      ! Snow = solid precipitation 
    382404      srcv(jpr_tevp)%clname = 'OTotEvap'      ! total evaporation (over oce + ice sublimation) 
    383       srcv(jpr_ievp)%clname = 'OIceEvap'      ! evaporation over ice = sublimation 
     405      srcv(jpr_ievp)%clname = 'OIceEvp'      ! evaporation over ice = sublimation 
    384406      srcv(jpr_sbpr)%clname = 'OSubMPre'      ! sublimation - liquid precipitation - solid precipitation  
    385407      srcv(jpr_semp)%clname = 'OISubMSn'      ! ice solid water budget = sublimation - solid precipitation 
     
    395417      END SELECT 
    396418      ! 
     419      !Set the number of categories for coupling of sublimation 
     420      IF ( TRIM( sn_rcv_emp%clcat ) == 'yes' ) srcv(jpr_ievp)%nct = jpl 
     421      ! 
    397422      !                                                      ! ------------------------- ! 
    398423      !                                                      !     Runoffs & Calving     !    
     
    408433      ! 
    409434      srcv(jpr_cal   )%clname = 'OCalving'   ;   IF( TRIM( sn_rcv_cal%cldes ) == 'coupled' )   srcv(jpr_cal)%laction = .TRUE. 
     435      srcv(jpr_grnm  )%clname = 'OGrnmass'   ;   IF( TRIM( sn_rcv_grnm%cldes ) == 'coupled' )   srcv(jpr_grnm)%laction = .TRUE. 
     436      srcv(jpr_antm  )%clname = 'OAntmass'   ;   IF( TRIM( sn_rcv_antm%cldes ) == 'coupled' )   srcv(jpr_antm)%laction = .TRUE. 
    410437      ! 
    411438      !                                                      ! ------------------------- ! 
     
    481508         srcv(jpr_topm:jpr_botm)%laction = .TRUE. 
    482509      ENDIF 
     510       
     511#if defined key_cice && ! defined key_cice4 
     512      !                                                      ! ----------------------------- ! 
     513      !                                                      !  sea-ice skin temperature     !    
     514      !                                                      !  used in meltpond scheme      ! 
     515      !                                                      !  May be calculated in Atm     ! 
     516      !                                                      ! ----------------------------- ! 
     517      srcv(jpr_ts_ice)%clname = 'OTsfIce' 
     518      IF ( TRIM( sn_rcv_ts_ice%cldes ) == 'ice' ) srcv(jpr_ts_ice)%laction = .TRUE. 
     519      IF ( TRIM( sn_rcv_ts_ice%clcat ) == 'yes' ) srcv(jpr_ts_ice)%nct = jpl 
     520      !TODO: Should there be a consistency check here? 
     521#endif 
     522 
    483523      !                                                      ! ------------------------------- ! 
    484524      !                                                      !   OPA-SAS coupling - rcv by opa !    
     
    598638      !                                                      ! ------------------------- ! 
    599639      ssnd(jps_toce)%clname = 'O_SSTSST' 
    600       ssnd(jps_tice)%clname = 'O_TepIce' 
     640      ssnd(jps_tice)%clname = 'OTepIce' 
    601641      ssnd(jps_tmix)%clname = 'O_TepMix' 
    602642      SELECT CASE( TRIM( sn_snd_temp%cldes ) ) 
    603643      CASE( 'none'                                 )       ! nothing to do 
    604644      CASE( 'oce only'                             )   ;   ssnd( jps_toce )%laction = .TRUE. 
    605       CASE( 'oce and ice' , 'weighted oce and ice' ) 
     645      CASE( 'oce and ice' , 'weighted oce and ice' , 'oce and weighted ice') 
    606646         ssnd( (/jps_toce, jps_tice/) )%laction = .TRUE. 
    607647         IF ( TRIM( sn_snd_temp%clcat ) == 'yes' )  ssnd(jps_tice)%nct = jpl 
     
    637677      ssnd(jps_hice)%clname = 'OIceTck' 
    638678      ssnd(jps_hsnw)%clname = 'OSnwTck' 
     679      ssnd(jps_a_p)%clname  = 'OPndFrc' 
     680      ssnd(jps_ht_p)%clname = 'OPndTck' 
     681      ssnd(jps_fice1)%clname = 'OIceFrd' 
    639682      IF( k_ice /= 0 ) THEN 
    640683         ssnd(jps_fice)%laction = .TRUE.                  ! if ice treated in the ocean (even in climato case) 
     684         ssnd(jps_fice1)%laction = .TRUE.                 ! First-order regridded ice concentration, to be used 
     685                                                     ! in producing atmos-to-ice fluxes 
    641686! Currently no namelist entry to determine sending of multi-category ice fraction so use the thickness entry for now 
    642687         IF ( TRIM( sn_snd_thick%clcat ) == 'yes' ) ssnd(jps_fice)%nct = jpl 
     
    655700      CASE default   ;   CALL ctl_stop( 'sbc_cpl_init: wrong definition of sn_snd_thick%cldes' ) 
    656701      END SELECT 
     702 
     703      !                                                      ! ------------------------- ! 
     704      !                                                      ! Ice Meltponds             ! 
     705      !                                                      ! ------------------------- ! 
     706#if defined key_cice && ! defined key_cice4 
     707      ! Meltponds only CICE5  
     708      ssnd(jps_a_p)%clname = 'OPndFrc'    
     709      ssnd(jps_ht_p)%clname = 'OPndTck'    
     710      SELECT CASE ( TRIM( sn_snd_mpnd%cldes ) ) 
     711      CASE ( 'none' ) 
     712         ssnd(jps_a_p)%laction = .FALSE. 
     713         ssnd(jps_ht_p)%laction = .FALSE. 
     714      CASE ( 'ice only' )  
     715         ssnd(jps_a_p)%laction = .TRUE. 
     716         ssnd(jps_ht_p)%laction = .TRUE. 
     717         IF ( TRIM( sn_snd_mpnd%clcat ) == 'yes' ) THEN 
     718            ssnd(jps_a_p)%nct = jpl 
     719            ssnd(jps_ht_p)%nct = jpl 
     720         ELSE 
     721            IF ( jpl > 1 ) THEN 
     722               CALL ctl_stop( 'sbc_cpl_init: use weighted ice option for sn_snd_mpnd%cldes if not exchanging category fields' ) 
     723            ENDIF 
     724         ENDIF 
     725      CASE ( 'weighted ice' )  
     726         ssnd(jps_a_p)%laction = .TRUE. 
     727         ssnd(jps_ht_p)%laction = .TRUE. 
     728         IF ( TRIM( sn_snd_mpnd%clcat ) == 'yes' ) THEN 
     729            ssnd(jps_a_p)%nct = jpl  
     730            ssnd(jps_ht_p)%nct = jpl  
     731         ENDIF 
     732      CASE default   ;   CALL ctl_stop( 'sbc_cpl_init: wrong definition of sn_snd_mpnd%cldes' ) 
     733      END SELECT 
     734#else 
     735      IF( TRIM( sn_snd_mpnd%cldes /= 'none' ) THEN 
     736         CALL ctl_stop('Meltponds can only be used with CICEv5') 
     737      ENDIF 
     738#endif 
    657739 
    658740      !                                                      ! ------------------------- ! 
     
    687769      !                                                      ! ------------------------- ! 
    688770      ssnd(jps_co2)%clname = 'O_CO2FLX' ;  IF( TRIM(sn_snd_co2%cldes) == 'coupled' )    ssnd(jps_co2 )%laction = .TRUE. 
     771      ! 
     772       
     773      !                                                      ! ------------------------- ! 
     774      !                                                      ! Sea surface freezing temp ! 
     775      !                                                      ! ------------------------- ! 
     776      ssnd(jps_sstfrz)%clname = 'O_SSTFrz' ; IF( TRIM(sn_snd_sstfrz%cldes) == 'coupled' )  ssnd(jps_sstfrz)%laction = .TRUE. 
     777      ! 
     778      !                                                      ! ------------------------- ! 
     779      !                                                      !    Ice conductivity       ! 
     780      !                                                      ! ------------------------- ! 
     781      ! Note that ultimately we will move to passing an ocean effective conductivity as well so there 
     782      ! will be some changes to the parts of the code which currently relate only to ice conductivity 
     783      ssnd(jps_kice )%clname = 'OIceKn' 
     784      SELECT CASE ( TRIM( sn_snd_cond%cldes ) ) 
     785      CASE ( 'none' ) 
     786         ssnd(jps_kice)%laction = .FALSE. 
     787      CASE ( 'ice only' ) 
     788         ssnd(jps_kice)%laction = .TRUE. 
     789         IF ( TRIM( sn_snd_cond%clcat ) == 'yes' ) THEN 
     790            ssnd(jps_kice)%nct = jpl 
     791         ELSE 
     792            IF ( jpl > 1 ) THEN 
     793               CALL ctl_stop( 'sbc_cpl_init: use weighted ice option for sn_snd_cond%cldes if not exchanging category fields' ) 
     794            ENDIF 
     795         ENDIF 
     796      CASE ( 'weighted ice' ) 
     797         ssnd(jps_kice)%laction = .TRUE. 
     798         IF ( TRIM( sn_snd_cond%clcat ) == 'yes' ) ssnd(jps_kice)%nct = jpl 
     799      CASE default   ;   CALL ctl_stop( 'sbc_cpl_init: wrong definition of sn_snd_cond%cldes' ) 
     800      END SELECT 
     801      ! 
     802       
    689803 
    690804      !                                                      ! ------------------------------- ! 
     
    783897      ncpl_qsr_freq = 86400 / ncpl_qsr_freq 
    784898 
    785       CALL wrk_dealloc( jpi,jpj,   zacs, zaos ) 
     899      IF( ln_coupled_iceshelf_fluxes ) THEN 
     900          ! Crude masks to separate the Antarctic and Greenland icesheets. Obviously something 
     901          ! more complicated could be done if required. 
     902          greenland_icesheet_mask = 0.0 
     903          WHERE( gphit >= 0.0 ) greenland_icesheet_mask = 1.0 
     904          antarctica_icesheet_mask = 0.0 
     905          WHERE( gphit < 0.0 ) antarctica_icesheet_mask = 1.0 
     906 
     907          ! initialise other variables 
     908          greenland_icesheet_mass_array(:,:) = 0.0 
     909          antarctica_icesheet_mass_array(:,:) = 0.0 
     910 
     911          IF( .not. ln_rstart ) THEN 
     912             greenland_icesheet_mass = 0.0  
     913             greenland_icesheet_mass_rate_of_change = 0.0  
     914             greenland_icesheet_timelapsed = 0.0 
     915             antarctica_icesheet_mass = 0.0  
     916             antarctica_icesheet_mass_rate_of_change = 0.0  
     917             antarctica_icesheet_timelapsed = 0.0 
     918          ENDIF 
     919 
     920      ENDIF 
     921 
     922      CALL wrk_dealloc( jpi,jpj, zacs, zaos ) 
    786923      ! 
    787924      IF( nn_timing == 1 )   CALL timing_stop('sbc_cpl_init') 
     
    841978      !! 
    842979      LOGICAL  ::   llnewtx, llnewtau      ! update wind stress components and module?? 
    843       INTEGER  ::   ji, jj, jn             ! dummy loop indices 
     980      INTEGER  ::   ji, jj, jl, jn         ! dummy loop indices 
    844981      INTEGER  ::   isec                   ! number of seconds since nit000 (assuming rdt did not change since nit000) 
    845982      REAL(wp) ::   zcumulneg, zcumulpos   ! temporary scalars      
     983      REAL(wp) ::   zgreenland_icesheet_mass_in, zantarctica_icesheet_mass_in 
     984      REAL(wp) ::   zgreenland_icesheet_mass_b, zantarctica_icesheet_mass_b 
     985      REAL(wp) ::   zmask_sum, zepsilon       
    846986      REAL(wp) ::   zcoef                  ! temporary scalar 
    847987      REAL(wp) ::   zrhoa  = 1.22          ! Air density kg/m3 
     
    8941034               IF( srcv(jpr_otx2)%laction ) THEN 
    8951035                  CALL rot_rep( frcv(jpr_otx2)%z3(:,:,1), frcv(jpr_oty2)%z3(:,:,1), srcv(jpr_otx2)%clgrid, 'en->j', zty )    
    896                ELSE   
     1036               ELSE 
    8971037                  CALL rot_rep( frcv(jpr_otx1)%z3(:,:,1), frcv(jpr_oty1)%z3(:,:,1), srcv(jpr_otx1)%clgrid, 'en->j', zty )   
    8981038               ENDIF 
     
    9911131#endif 
    9921132 
     1133#if defined key_cice && ! defined key_cice4 
     1134      !  ! Sea ice surface skin temp: 
     1135      IF( srcv(jpr_ts_ice)%laction ) THEN 
     1136        DO jl = 1, jpl 
     1137          DO jj = 1, jpj 
     1138            DO ji = 1, jpi 
     1139              IF (frcv(jpr_ts_ice)%z3(ji,jj,jl) > 0.0) THEN 
     1140                tsfc_ice(ji,jj,jl) = 0.0 
     1141              ELSE IF (frcv(jpr_ts_ice)%z3(ji,jj,jl) < -60.0) THEN 
     1142                tsfc_ice(ji,jj,jl) = -60.0 
     1143              ELSE 
     1144                tsfc_ice(ji,jj,jl) = frcv(jpr_ts_ice)%z3(ji,jj,jl) 
     1145              ENDIF 
     1146            END DO 
     1147          END DO 
     1148        END DO 
     1149      ENDIF 
     1150#endif 
     1151 
    9931152      !  Fields received by SAS when OASIS coupling 
    9941153      !  (arrays no more filled at sbcssm stage) 
     
    11051264         ! 
    11061265      ENDIF 
     1266       
     1267      !                                                        ! land ice masses : Greenland 
     1268      zepsilon = rn_iceshelf_fluxes_tolerance 
     1269 
     1270      IF( srcv(jpr_grnm)%laction ) THEN 
     1271         greenland_icesheet_mass_array(:,:) = frcv(jpr_grnm)%z3(:,:,1) 
     1272         ! take average over ocean points of input array to avoid cumulative error over time 
     1273         zgreenland_icesheet_mass_in = SUM( greenland_icesheet_mass_array(:,:) * tmask(:,:,1) ) 
     1274         IF(lk_mpp) CALL mpp_sum( zgreenland_icesheet_mass_in ) 
     1275         zmask_sum = SUM( tmask(:,:,1) ) 
     1276         IF(lk_mpp) CALL mpp_sum( zmask_sum )  
     1277         zgreenland_icesheet_mass_in = zgreenland_icesheet_mass_in / zmask_sum 
     1278         greenland_icesheet_timelapsed = greenland_icesheet_timelapsed + rdt          
     1279         IF( ABS( zgreenland_icesheet_mass_in - greenland_icesheet_mass ) > zepsilon ) THEN 
     1280            zgreenland_icesheet_mass_b = greenland_icesheet_mass 
     1281             
     1282            ! Only update the mass if it has increased 
     1283            IF ( (zgreenland_icesheet_mass_in - greenland_icesheet_mass) > 0.0 ) THEN 
     1284               greenland_icesheet_mass = zgreenland_icesheet_mass_in 
     1285            ENDIF 
     1286             
     1287            IF( zgreenland_icesheet_mass_b /= 0.0 ) & 
     1288           &     greenland_icesheet_mass_rate_of_change = ( greenland_icesheet_mass - zgreenland_icesheet_mass_b ) / greenland_icesheet_timelapsed  
     1289            greenland_icesheet_timelapsed = 0.0_wp        
     1290         ENDIF 
     1291         IF(lwp) WRITE(numout,*) 'Greenland icesheet mass (kg) read in is ', zgreenland_icesheet_mass_in 
     1292         IF(lwp) WRITE(numout,*) 'Greenland icesheet mass (kg) used is    ', greenland_icesheet_mass 
     1293         IF(lwp) WRITE(numout,*) 'Greenland icesheet mass rate of change (kg/s) is ', greenland_icesheet_mass_rate_of_change 
     1294         IF(lwp) WRITE(numout,*) 'Greenland icesheet seconds lapsed since last change is ', greenland_icesheet_timelapsed 
     1295      ENDIF 
     1296 
     1297      !                                                        ! land ice masses : Antarctica 
     1298      IF( srcv(jpr_antm)%laction ) THEN 
     1299         antarctica_icesheet_mass_array(:,:) = frcv(jpr_antm)%z3(:,:,1) 
     1300         ! take average over ocean points of input array to avoid cumulative error from rounding errors over time 
     1301         zantarctica_icesheet_mass_in = SUM( antarctica_icesheet_mass_array(:,:) * tmask(:,:,1) ) 
     1302         IF(lk_mpp) CALL mpp_sum( zantarctica_icesheet_mass_in ) 
     1303         zmask_sum = SUM( tmask(:,:,1) ) 
     1304         IF(lk_mpp) CALL mpp_sum( zmask_sum )  
     1305         zantarctica_icesheet_mass_in = zantarctica_icesheet_mass_in / zmask_sum 
     1306         antarctica_icesheet_timelapsed = antarctica_icesheet_timelapsed + rdt          
     1307         IF( ABS( zantarctica_icesheet_mass_in - antarctica_icesheet_mass ) > zepsilon ) THEN 
     1308            zantarctica_icesheet_mass_b = antarctica_icesheet_mass 
     1309             
     1310            ! Only update the mass if it has increased 
     1311            IF ( (zantarctica_icesheet_mass_in - antarctica_icesheet_mass) > 0.0 ) THEN 
     1312               antarctica_icesheet_mass = zantarctica_icesheet_mass_in 
     1313            END IF 
     1314             
     1315            IF( zantarctica_icesheet_mass_b /= 0.0 ) & 
     1316          &      antarctica_icesheet_mass_rate_of_change = ( antarctica_icesheet_mass - zantarctica_icesheet_mass_b ) / antarctica_icesheet_timelapsed  
     1317            antarctica_icesheet_timelapsed = 0.0_wp        
     1318         ENDIF 
     1319         IF(lwp) WRITE(numout,*) 'Antarctica icesheet mass (kg) read in is ', zantarctica_icesheet_mass_in 
     1320         IF(lwp) WRITE(numout,*) 'Antarctica icesheet mass (kg) used is    ', antarctica_icesheet_mass 
     1321         IF(lwp) WRITE(numout,*) 'Antarctica icesheet mass rate of change (kg/s) is ', antarctica_icesheet_mass_rate_of_change 
     1322         IF(lwp) WRITE(numout,*) 'Antarctica icesheet seconds lapsed since last change is ', antarctica_icesheet_timelapsed 
     1323      ENDIF 
     1324 
    11071325      ! 
    11081326      CALL wrk_dealloc( jpi,jpj,   ztx, zty, zmsk, zemp, zqns, zqsr ) 
     
    14001618         ztprecip(:,:) = frcv(jpr_rain)%z3(:,:,1) + zsprecip(:,:)  ! May need to ensure positive here 
    14011619         zemp_tot(:,:) = frcv(jpr_tevp)%z3(:,:,1) - ztprecip(:,:) 
     1620#if defined key_cice 
     1621         IF ( TRIM(sn_rcv_emp%clcat) == 'yes' ) THEN 
     1622            ! zemp_ice is the sum of frcv(jpr_ievp)%z3(:,:,1) over all layers - snow 
     1623            zemp_ice(:,:) = - frcv(jpr_snow)%z3(:,:,1) 
     1624            DO jl=1,jpl 
     1625               zemp_ice(:,:   ) = zemp_ice(:,:) + frcv(jpr_ievp)%z3(:,:,jl) 
     1626            ENDDO 
     1627            ! latent heat coupled for each category in CICE 
     1628            qla_ice(:,:,1:jpl) = - frcv(jpr_ievp)%z3(:,:,1:jpl) * lsub 
     1629         ELSE 
     1630            ! If CICE has multicategories it still expects coupling fields for 
     1631            ! each even if we treat as a single field 
     1632            ! The latent heat flux is split between the ice categories according 
     1633            ! to the fraction of the ice in each category 
     1634            zemp_ice(:,:) = frcv(jpr_ievp)%z3(:,:,1) - frcv(jpr_snow)%z3(:,:,1) 
     1635            WHERE ( zicefr(:,:) /= 0._wp )  
     1636               ztmp(:,:) = 1./zicefr(:,:) 
     1637            ELSEWHERE  
     1638               ztmp(:,:) = 0.e0 
     1639            END WHERE   
     1640            DO jl=1,jpl 
     1641               qla_ice(:,:,jl) = - a_i(:,:,jl) * ztmp(:,:) * frcv(jpr_ievp)%z3(:,:,1) * lsub  
     1642            END DO 
     1643            WHERE ( zicefr(:,:) == 0._wp )  qla_ice(:,:,1) = -frcv(jpr_ievp)%z3(:,:,1) * lsub  
     1644         ENDIF 
     1645#else          
    14021646         zemp_ice(:,:) = frcv(jpr_ievp)%z3(:,:,1) - frcv(jpr_snow)%z3(:,:,1) 
     1647#endif                   
    14031648            CALL iom_put( 'rain'         , frcv(jpr_rain)%z3(:,:,1)              )   ! liquid precipitation  
    14041649         IF( iom_use('hflx_rain_cea') )   & 
     
    17952040               CASE default                  ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_temp%clcat' ) 
    17962041               END SELECT 
     2042            CASE( 'oce and weighted ice' )   ;   ztmp1(:,:) =   tsn(:,:,1,jp_tem) + rt0  
     2043               SELECT CASE( sn_snd_temp%clcat ) 
     2044               CASE( 'yes' )    
     2045                  ztmp3(:,:,1:jpl) = tn_ice(:,:,1:jpl) * a_i(:,:,1:jpl) 
     2046               CASE( 'no' ) 
     2047                  ztmp3(:,:,:) = 0.0 
     2048                  DO jl=1,jpl 
     2049                     ztmp3(:,:,1) = ztmp3(:,:,1) + tn_ice(:,:,jl) * a_i(:,:,jl) 
     2050                  ENDDO 
     2051               CASE default                  ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_temp%clcat' ) 
     2052               END SELECT 
    17972053            CASE( 'mixed oce-ice'        )    
    17982054               ztmp1(:,:) = ( ztmp1(:,:) + rt0 ) * zfr_l(:,:)  
     
    18652121         END SELECT 
    18662122         IF( ssnd(jps_fice)%laction )   CALL cpl_snd( jps_fice, isec, ztmp3, info ) 
     2123      ENDIF 
     2124       
     2125      ! Send ice fraction field (first order interpolation), for weighting UM fluxes to be passed to NEMO 
     2126      IF (ssnd(jps_fice1)%laction) THEN 
     2127         SELECT CASE (sn_snd_thick1%clcat) 
     2128         CASE( 'yes' )   ;   ztmp3(:,:,1:jpl) =  a_i(:,:,1:jpl) 
     2129         CASE( 'no' )    ;   ztmp3(:,:,1) = fr_i(:,:) 
     2130         CASE default    ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_thick1%clcat' ) 
     2131      END SELECT 
     2132         CALL cpl_snd (jps_fice1, isec, ztmp3, info) 
    18672133      ENDIF 
    18682134       
     
    19102176         IF( ssnd(jps_hsnw)%laction )   CALL cpl_snd( jps_hsnw, isec, ztmp4, info ) 
    19112177      ENDIF 
     2178      ! 
     2179      ! Send meltpond fields  
     2180      IF( ssnd(jps_a_p)%laction .OR. ssnd(jps_ht_p)%laction ) THEN 
     2181         SELECT CASE( sn_snd_mpnd%cldes)  
     2182         CASE( 'weighted ice' )  
     2183            SELECT CASE( sn_snd_mpnd%clcat )  
     2184            CASE( 'yes' )  
     2185               ztmp3(:,:,1:jpl) =  a_p(:,:,1:jpl) * a_i(:,:,1:jpl)  
     2186               ztmp4(:,:,1:jpl) =  ht_p(:,:,1:jpl) * a_i(:,:,1:jpl)  
     2187            CASE( 'no' )  
     2188               ztmp3(:,:,:) = 0.0  
     2189               ztmp4(:,:,:) = 0.0  
     2190               DO jl=1,jpl  
     2191                 ztmp3(:,:,1) = ztmp3(:,:,1) + a_p(:,:,jpl) * a_i(:,:,jpl)  
     2192                 ztmp4(:,:,1) = ztmp4(:,:,1) + ht_p(:,:,jpl) * a_i(:,:,jpl)  
     2193               ENDDO  
     2194            CASE default    ;   CALL ctl_stop( 'sbc_cpl_mpd: wrong definition of sn_snd_mpnd%clcat' )  
     2195            END SELECT  
     2196         CASE( 'ice only' )     
     2197            ztmp3(:,:,1:jpl) = a_p(:,:,1:jpl)  
     2198            ztmp4(:,:,1:jpl) = ht_p(:,:,1:jpl)  
     2199         END SELECT  
     2200         IF( ssnd(jps_a_p)%laction )   CALL cpl_snd( jps_a_p, isec, ztmp3, info )     
     2201         IF( ssnd(jps_ht_p)%laction )   CALL cpl_snd( jps_ht_p, isec, ztmp4, info )     
     2202         ! 
     2203         ! Send ice effective conductivity 
     2204         SELECT CASE( sn_snd_cond%cldes) 
     2205         CASE( 'weighted ice' )    
     2206            SELECT CASE( sn_snd_cond%clcat ) 
     2207            CASE( 'yes' )    
     2208               ztmp3(:,:,1:jpl) =  kn_ice(:,:,1:jpl) * a_i(:,:,1:jpl) 
     2209            CASE( 'no' ) 
     2210               ztmp3(:,:,:) = 0.0 
     2211               DO jl=1,jpl 
     2212                 ztmp3(:,:,1) = ztmp3(:,:,1) + kn_ice(:,:,jl) * a_i(:,:,jl) 
     2213               ENDDO 
     2214            CASE default                  ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_cond%clcat' ) 
     2215            END SELECT 
     2216         CASE( 'ice only' )    
     2217           ztmp3(:,:,1:jpl) = kn_ice(:,:,1:jpl) 
     2218         END SELECT 
     2219         IF( ssnd(jps_kice)%laction )   CALL cpl_snd( jps_kice, isec, ztmp3, info ) 
     2220      ENDIF 
     2221      ! 
    19122222      ! 
    19132223#if defined key_cpl_carbon_cycle 
     
    20892399      IF( ssnd(jps_rnf   )%laction )  CALL cpl_snd( jps_rnf   , isec, RESHAPE ( rnf , (/jpi,jpj,1/) ), info ) 
    20902400      IF( ssnd(jps_taum  )%laction )  CALL cpl_snd( jps_taum  , isec, RESHAPE ( taum, (/jpi,jpj,1/) ), info ) 
    2091  
     2401       
     2402      ztmp1(:,:) = sstfrz(:,:) + rt0 
     2403      IF( ssnd(jps_sstfrz)%laction )  CALL cpl_snd( jps_sstfrz, isec, RESHAPE ( ztmp1, (/jpi,jpj,1/) ), info ) 
     2404      ! 
    20922405      CALL wrk_dealloc( jpi,jpj,       zfr_l, ztmp1, ztmp2, zotx1, zoty1, zotz1, zitx1, zity1, zitz1 ) 
    20932406      CALL wrk_dealloc( jpi,jpj,jpl,   ztmp3, ztmp4 ) 
  • branches/UKMO/dev_r6501_GO6_package_trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90

    r6503 r6507  
    1313   USE dom_oce         ! ocean space and time domain 
    1414   USE domvvl 
    15    USE phycst, only : rcp, rau0, r1_rau0, rhosn, rhoic 
     15   USE eosbn2, only : eos_fzp ! Function to calculate freezing point of seawater 
     16   USE phycst, only : rcp, rau0, r1_rau0, rhosn, rhoic, rt0 
    1617   USE in_out_manager  ! I/O manager 
    1718   USE iom, ONLY : iom_put,iom_use              ! I/O manager library !!Joakim edit 
     
    3536   USE ice_gather_scatter 
    3637   USE ice_calendar, only: dt 
     38# if defined key_cice4 
    3739   USE ice_state, only: aice,aicen,uvel,vvel,vsno,vsnon,vice,vicen 
    38 # if defined key_cice4 
    3940   USE ice_flux, only: strax,stray,strocnx,strocny,frain,fsnow,  & 
    4041                strocnxT,strocnyT,                               &  
     
    4344                flatn_f,fsurfn_f,fcondtopn_f,                    & 
    4445                uatm,vatm,wind,fsw,flw,Tair,potT,Qa,rhoa,zlvl,   & 
    45                 swvdr,swvdf,swidr,swidf 
     46                swvdr,swvdf,swidr,swidf,Tf 
    4647   USE ice_therm_vertical, only: calc_Tsfc 
    4748#else 
     49   USE ice_state, only: aice,aicen,uvel,nt_hpnd,trcrn,vvel,vsno,& 
     50                vsnon,vice,vicen,nt_Tsfc 
    4851   USE ice_flux, only: strax,stray,strocnx,strocny,frain,fsnow,  & 
    4952                strocnxT,strocnyT,                               &  
    50                 sst,sss,uocn,vocn,ss_tltx,ss_tlty,fsalt_ai,     & 
    51                 fresh_ai,fhocn_ai,fswthru_ai,frzmlt,          & 
     53                sst,sss,uocn,vocn,ss_tltx,ss_tlty,fsalt_ai,      & 
     54                fresh_ai,fhocn_ai,fswthru_ai,frzmlt,             & 
    5255                flatn_f,fsurfn_f,fcondtopn_f,                    & 
    5356                uatm,vatm,wind,fsw,flw,Tair,potT,Qa,rhoa,zlvl,   & 
    54                 swvdr,swvdf,swidr,swidf 
    55    USE ice_therm_shared, only: calc_Tsfc 
     57                swvdr,swvdf,swidr,swidf,Tf,                      & 
     58                !! When using NEMO with CICE, this change requires use of  
     59                !! one of the following two CICE branches: 
     60                !! - at CICE5.0,   hadax/r1015_GSI8_with_GSI7 
     61                !! - at CICE5.1.2, hadax/vn5.1.2_GSI8 
     62                keffn_top,Tn_top 
     63 
     64   USE ice_therm_shared, only: calc_Tsfc, heat_capacity 
     65   USE ice_shortwave, only: apeffn 
    5666#endif 
    5767   USE ice_forcing, only: frcvdr,frcvdf,frcidr,frcidf 
     
    162172      REAL(wp), DIMENSION(:,:), POINTER :: ztmp1, ztmp2 
    163173      REAL(wp) ::   zcoefu, zcoefv, zcoeff            ! local scalar 
     174      REAL(wp), DIMENSION(:,:,:), POINTER :: ztfrz3d 
    164175      INTEGER  ::   ji, jj, jl, jk                    ! dummy loop indices 
    165176      !!--------------------------------------------------------------------- 
     
    173184      ji_off = INT ( (jpiglo - nx_global) / 2 ) 
    174185      jj_off = INT ( (jpjglo - ny_global) / 2 ) 
    175  
    176 #if defined key_nemocice_decomp 
    177       ! Pass initial SST from NEMO to CICE so ice is initialised correctly if 
    178       ! there is no restart file. 
    179       ! Values from a CICE restart file would overwrite this 
    180       IF ( .NOT. ln_rstart ) THEN     
    181          CALL nemo2cice( tsn(:,:,1,jp_tem) , sst , 'T' , 1.)  
    182       ENDIF   
    183 #endif 
    184186 
    185187! Initialize CICE 
     
    199201 
    200202! allocate sbc_ice and sbc_cice arrays 
    201       IF( sbc_ice_alloc()      /= 0 )   CALL ctl_stop( 'STOP', 'sbc_ice_cice_alloc : unable to allocate arrays' ) 
     203      IF( sbc_ice_alloc()      /= 0 )   CALL ctl_stop( 'STOP', 'sbc_ice_alloc : unable to allocate arrays' ) 
    202204      IF( sbc_ice_cice_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'sbc_ice_cice_alloc : unable to allocate cice arrays' ) 
    203205 
    204 ! Ensure ocean temperatures are nowhere below freezing if not a NEMO restart 
     206      ! Ensure that no temperature points are below freezing if not a NEMO restart 
    205207      IF( .NOT. ln_rstart ) THEN 
    206          tsn(:,:,:,jp_tem) = MAX (tsn(:,:,:,jp_tem),Tocnfrz) 
     208         CALL wrk_alloc( jpi,jpj,jpk, ztfrz3d )  
     209         DO jk=1,jpk 
     210             CALL eos_fzp( tsn(:,:,jk,jp_sal), ztfrz3d(:,:,jk), fsdept_n(:,:,jk) ) 
     211         ENDDO 
     212         tsn(:,:,:,jp_tem) = MAX( tsn(:,:,:,jp_tem), ztfrz3d ) 
    207213         tsb(:,:,:,jp_tem) = tsn(:,:,:,jp_tem) 
    208       ENDIF 
    209  
    210       fr_iu(:,:)=0.0 
    211       fr_iv(:,:)=0.0 
     214         CALL wrk_dealloc( jpi,jpj,jpk, ztfrz3d )  
     215 
     216#if defined key_nemocice_decomp 
     217         ! Pass initial SST from NEMO to CICE so ice is initialised correctly if 
     218         ! there is no restart file. 
     219         ! Values from a CICE restart file would overwrite this 
     220         CALL nemo2cice( tsn(:,:,1,jp_tem) , sst , 'T' , 1.)  
     221#endif 
     222 
     223      ENDIF 
     224 
     225      ! calculate surface freezing temperature and send to CICE 
     226      CALL  eos_fzp(sss_m(:,:), sstfrz(:,:), fsdept_n(:,:,1))  
     227      CALL nemo2cice(sstfrz,Tf, 'T', 1. ) 
    212228 
    213229      CALL cice2nemo(aice,fr_i, 'T', 1. ) 
     
    220236! T point to U point 
    221237! T point to V point 
     238      fr_iu(:,:)=0.0 
     239      fr_iv(:,:)=0.0 
    222240      DO jj=1,jpjm1 
    223241         DO ji=1,jpim1 
     
    345363         CALL nemo2cice(ztmp,stray,'F', -1. ) 
    346364 
     365 
     366! Alex West: From configuration GSI8 onwards, when NEMO is used with CICE in 
     367! HadGEM3 the 'time-travelling ice' coupling approach is used, whereby  
     368! atmosphere-ice fluxes are passed as pseudo-local values, formed by dividing 
     369! gridbox mean fluxes in the UM by future ice concentration obtained through   
     370! OASIS.  This allows for a much more realistic apportionment of energy through 
     371! the ice - and conserves energy. 
     372! Therefore the fluxes are now divided by ice concentration in the coupled 
     373! formulation (jp_purecpl) as well as for jp_flx.  This NEMO branch should only 
     374! be used at UM10.2 onwards (unless an explicit GSI8 UM branch is included), at 
     375! which point the GSI8 UM changes were committed. 
     376 
    347377! Surface downward latent heat flux (CI_5) 
    348378         IF (ksbc == jp_flx) THEN 
     
    350380               ztmpn(:,:,jl)=qla_ice(:,:,1)*a_i(:,:,jl) 
    351381            ENDDO 
     382         ELSE IF (ksbc == jp_purecpl) THEN 
     383            DO jl=1,ncat 
     384               ztmpn(:,:,jl)=qla_ice(:,:,jl)*a_i(:,:,jl) 
     385            ENDDO 
    352386         ELSE 
    353 ! emp_ice is set in sbc_cpl_ice_flx as sublimation-snow 
    354             qla_ice(:,:,1)= - ( emp_ice(:,:)+sprecip(:,:) ) * Lsub 
    355 ! End of temporary code 
    356             DO jj=1,jpj 
    357                DO ji=1,jpi 
    358                   IF (fr_i(ji,jj).eq.0.0) THEN 
    359                      DO jl=1,ncat 
    360                         ztmpn(ji,jj,jl)=0.0 
    361                      ENDDO 
    362                      ! This will then be conserved in CICE 
    363                      ztmpn(ji,jj,1)=qla_ice(ji,jj,1) 
    364                   ELSE 
    365                      DO jl=1,ncat 
    366                         ztmpn(ji,jj,jl)=qla_ice(ji,jj,1)*a_i(ji,jj,jl)/fr_i(ji,jj) 
    367                      ENDDO 
    368                   ENDIF 
    369                ENDDO 
    370             ENDDO 
     387           !In coupled mode - qla_ice calculated in sbc_cpl for each category 
     388           ztmpn(:,:,1:ncat)=qla_ice(:,:,1:ncat) 
    371389         ENDIF 
     390 
    372391         DO jl=1,ncat 
    373392            CALL nemo2cice(ztmpn(:,:,jl),flatn_f(:,:,jl,:),'T', 1. ) 
     
    375394! GBM conductive flux through ice (CI_6) 
    376395!  Convert to GBM 
    377             IF (ksbc == jp_flx) THEN 
     396            IF (ksbc == jp_flx .OR. ksbc == jp_purecpl) THEN 
    378397               ztmp(:,:) = botmelt(:,:,jl)*a_i(:,:,jl) 
    379398            ELSE 
     
    384403! GBM surface heat flux (CI_7) 
    385404!  Convert to GBM 
    386             IF (ksbc == jp_flx) THEN 
     405            IF (ksbc == jp_flx .OR. ksbc == jp_purecpl) THEN 
    387406               ztmp(:,:) = (topmelt(:,:,jl)+botmelt(:,:,jl))*a_i(:,:,jl)  
    388407            ELSE 
     
    456475      CALL nemo2cice(sss_m,sss,'T', 1. ) 
    457476 
     477      IF( ksbc == jp_purecpl ) THEN 
     478! Sea ice surface skin temperature 
     479         DO jl=1,ncat 
     480           CALL nemo2cice(tsfc_ice(:,:,jl), trcrn(:,:,nt_tsfc,jl,:),'T',1.) 
     481         ENDDO  
     482      ENDIF 
     483 
    458484! x comp and y comp of surface ocean current 
    459485! U point to F point 
     
    730756         CALL cice2nemo(vicen(:,:,jl,:),ht_i(:,:,jl),'T', 1. ) 
    731757      ENDDO 
     758 
     759#if ! defined key_cice4 
     760! Meltpond fraction and depth 
     761      DO jl = 1,ncat 
     762         CALL cice2nemo(apeffn(:,:,jl,:),a_p(:,:,jl),'T', 1. ) 
     763         CALL cice2nemo(trcrn(:,:,nt_hpnd,jl,:),ht_p(:,:,jl),'T', 1. ) 
     764      ENDDO 
     765#endif 
     766 
     767 
     768! If using multilayers thermodynamics in CICE then get top layer temperature 
     769! and effective conductivity        
     770!! When using NEMO with CICE, this change requires use of  
     771!! one of the following two CICE branches: 
     772!! - at CICE5.0,   hadax/r1015_GSI8_with_GSI7 
     773!! - at CICE5.1.2, hadax/vn5.1.2_GSI8 
     774      IF (heat_capacity) THEN 
     775         DO jl = 1,ncat 
     776            CALL cice2nemo(Tn_top(:,:,jl,:),tn_ice(:,:,jl),'T', 1. ) 
     777            CALL cice2nemo(keffn_top(:,:,jl,:),kn_ice(:,:,jl),'T', 1. ) 
     778         ENDDO 
     779! Convert surface temperature to Kelvin 
     780         tn_ice(:,:,:)=tn_ice(:,:,:)+rt0 
     781      ELSE 
     782         tn_ice(:,:,:) = 0.0 
     783         kn_ice(:,:,:) = 0.0 
     784      ENDIF        
     785 
    732786      ! 
    733787      IF( nn_timing == 1 )  CALL timing_stop('cice_sbc_hadgam') 
  • branches/UKMO/dev_r6501_GO6_package_trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcisf.F90

    r6503 r6507  
    133133            CALL fld_read ( kt, nn_fsbc, sf_rnfisf   ) 
    134134            fwfisf(:,:) = - sf_rnfisf(1)%fnow(:,:,1)         ! fwf  flux from the isf (fwfisf <0 mean melting)  
     135 
     136            IF( lk_oasis) THEN 
     137            ! ln_coupled_iceshelf_fluxes uninitialised unless lk_oasis=true 
     138            IF( ln_coupled_iceshelf_fluxes ) THEN 
     139 
     140              ! Adjust total iceshelf melt rates so that sum of iceberg calving and iceshelf melting in the northern 
     141              ! and southern hemispheres equals rate of increase of mass of greenland and antarctic ice sheets 
     142              ! to preserve total freshwater conservation in coupled models without an active ice sheet model. 
     143 
     144               zgreenland_fwfisf_sum = SUM( fwfisf(:,:) * e1t(:,:) * e2t(:,:) * greenland_icesheet_mask(:,:) ) 
     145               IF( lk_mpp ) CALL mpp_sum( zgreenland_fwfisf_sum ) 
     146               ! use ABS function because we need to preserve the sign of fwfisf 
     147               WHERE( greenland_icesheet_mask(:,:) == 1.0 )                                                                  & 
     148              &    fwfisf(:,:) = fwfisf(:,:)  * ABS( greenland_icesheet_mass_rate_of_change * (1.0-rn_greenland_calving_fraction) & 
     149              &                           / ( zgreenland_fwfisf_sum + 1.0e-10_wp ) ) 
     150 
     151               ! check 
     152               IF(lwp) WRITE(numout, *) 'Greenland iceshelf melting climatology (kg/s) : ',zgreenland_fwfisf_sum 
     153               zgreenland_fwfisf_sum = SUM( fwfisf(:,:) * e1t(:,:) * e2t(:,:) * greenland_icesheet_mask(:,:) ) 
     154               IF( lk_mpp ) CALL mpp_sum( zgreenland_fwfisf_sum ) 
     155               IF(lwp) WRITE(numout, *) 'Greenland iceshelf melting adjusted value (kg/s) : ',zgreenland_fwfisf_sum 
     156 
     157               zantarctica_fwfisf_sum = SUM( fwfisf(:,:) * e1t(:,:) * e2t(:,:) * antarctica_icesheet_mask(:,:) ) 
     158               IF( lk_mpp ) CALL mpp_sum( zantarctica_fwfisf_sum ) 
     159               ! use ABS function because we need to preserve the sign of fwfisf 
     160               WHERE( antarctica_icesheet_mask(:,:) == 1.0 ) & 
     161              &    fwfisf(:,:) = fwfisf(:,:)  * ABS( antarctica_icesheet_mass_rate_of_change * (1.0-rn_antarctica_calving_fraction) & 
     162              &                           / ( zantarctica_fwfisf_sum + 1.0e-10_wp ) ) 
     163       
     164               ! check 
     165               IF(lwp) WRITE(numout, *) 'Antarctica iceshelf melting climatology (kg/s) : ',zantarctica_fwfisf_sum 
     166               zantarctica_fwfisf_sum = SUM( fwfisf(:,:) * e1t(:,:) * e2t(:,:) * antarctica_icesheet_mask(:,:) ) 
     167               IF( lk_mpp ) CALL mpp_sum( zantarctica_fwfisf_sum ) 
     168               IF(lwp) WRITE(numout, *) 'Antarctica iceshelf melting adjusted value (kg/s) : ',zantarctica_fwfisf_sum 
     169 
     170            ENDIF 
     171            ENDIF 
     172 
    135173            qisf(:,:)   = fwfisf(:,:) * rlfusisf             ! heat flux 
    136174            stbl(:,:)   = soce 
     
    139177            CALL fld_read ( kt, nn_fsbc, sf_fwfisf   ) 
    140178            fwfisf(:,:) = - sf_fwfisf(1)%fnow(:,:,1)           ! fwf  flux from the isf (fwfisf <0 mean melting) 
     179 
     180            IF( lk_oasis) THEN 
     181            ! ln_coupled_iceshelf_fluxes uninitialised unless lk_oasis=true 
     182            IF( ln_coupled_iceshelf_fluxes ) THEN 
     183 
     184              ! Adjust total iceshelf melt rates so that sum of iceberg calving and iceshelf melting in the northern 
     185              ! and southern hemispheres equals rate of increase of mass of greenland and antarctic ice sheets 
     186              ! to preserve total freshwater conservation in coupled models without an active ice sheet model. 
     187 
     188               zgreenland_fwfisf_sum = SUM( fwfisf(:,:) * e1t(:,:) * e2t(:,:) * greenland_icesheet_mask(:,:) ) 
     189               IF( lk_mpp ) CALL mpp_sum( zgreenland_fwfisf_sum ) 
     190               ! use ABS function because we need to preserve the sign of fwfisf 
     191               WHERE( greenland_icesheet_mask(:,:) == 1.0 )                                                                  & 
     192              &    fwfisf(:,:) = fwfisf(:,:)  * ABS( greenland_icesheet_mass_rate_of_change * (1.0-rn_greenland_calving_fraction) & 
     193              &                           / ( zgreenland_fwfisf_sum + 1.0e-10_wp ) ) 
     194 
     195               ! check 
     196               IF(lwp) WRITE(numout, *) 'Greenland iceshelf melting climatology (kg/s) : ',zgreenland_fwfisf_sum 
     197               zgreenland_fwfisf_sum = SUM( fwfisf(:,:) * e1t(:,:) * e2t(:,:) * greenland_icesheet_mask(:,:) ) 
     198               IF( lk_mpp ) CALL mpp_sum( zgreenland_fwfisf_sum ) 
     199               IF(lwp) WRITE(numout, *) 'Greenland iceshelf melting adjusted value (kg/s) : ',zgreenland_fwfisf_sum 
     200 
     201               zantarctica_fwfisf_sum = SUM( fwfisf(:,:) * e1t(:,:) * e2t(:,:) * antarctica_icesheet_mask(:,:) ) 
     202               IF( lk_mpp ) CALL mpp_sum( zantarctica_fwfisf_sum ) 
     203               ! use ABS function because we need to preserve the sign of fwfisf 
     204               WHERE( antarctica_icesheet_mask(:,:) == 1.0 ) & 
     205              &    fwfisf(:,:) = fwfisf(:,:)  * ABS( antarctica_icesheet_mass_rate_of_change * (1.0-rn_antarctica_calving_fraction) & 
     206              &                           / ( zantarctica_fwfisf_sum + 1.0e-10_wp ) ) 
     207       
     208               ! check 
     209               IF(lwp) WRITE(numout, *) 'Antarctica iceshelf melting climatology (kg/s) : ',zantarctica_fwfisf_sum 
     210               zantarctica_fwfisf_sum = SUM( fwfisf(:,:) * e1t(:,:) * e2t(:,:) * antarctica_icesheet_mask(:,:) ) 
     211               IF( lk_mpp ) CALL mpp_sum( zantarctica_fwfisf_sum ) 
     212               IF(lwp) WRITE(numout, *) 'Antarctica iceshelf melting adjusted value (kg/s) : ',zantarctica_fwfisf_sum 
     213 
     214            ENDIF 
     215            ENDIF 
     216 
    141217            qisf(:,:)   = fwfisf(:,:) * rlfusisf               ! heat flux 
    142218            stbl(:,:)   = soce 
     
    155231         CALL eos_fzp( stbl(:,:), zt_frz(:,:), zdep(:,:) ) 
    156232          
     233         ! if apply only on the trend and not as a volume flux (rdivisf = 0), fwfisf have to be set to 0 now 
     234         fwfisf(:,:) = rdivisf * fwfisf(:,:)          
     235  
    157236         risf_tsc(:,:,jp_tem) = qisf(:,:) * r1_rau0_rcp - fwfisf(:,:) * zt_frz(:,:) * r1_rau0 ! 
    158237         risf_tsc(:,:,jp_sal) = 0.0_wp 
  • branches/UKMO/dev_r6501_GO6_package_trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90

    r6503 r6507  
    189189         fwfisf  (:,:)   = 0.0_wp ; fwfisf_b  (:,:)   = 0.0_wp 
    190190         risf_tsc(:,:,:) = 0.0_wp ; risf_tsc_b(:,:,:) = 0.0_wp 
     191         rdivisf       = 0.0_wp 
    191192      END IF 
    192193      IF( nn_ice == 0 .AND. nn_components /= jp_iam_opa )   fr_i(:,:) = 0._wp    ! no ice in the domain, ice fraction is always zero 
  • branches/UKMO/dev_r6501_GO6_package_trunk/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfmxl.F90

    r6503 r6507  
    1818   USE phycst          ! physical constants 
    1919   USE iom             ! I/O library 
     20   USE eosbn2          ! for zdf_mxl_zint 
    2021   USE lib_mpp         ! MPP library 
    2122   USE wrk_nemo        ! work arrays 
     
    3233   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hmlp    !: mixed layer depth  (rho=rho0+zdcrit) [m] 
    3334   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hmlpt   !: depth of the last T-point inside the mixed layer [m] 
     35   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hmld_zint  !: vertically-interpolated mixed layer depth   [m]  
     36   LOGICAL, PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)    :: ll_found   ! Is T_b to be found by interpolation ?  
     37   LOGICAL, PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  :: ll_belowml ! Flag points below mixed layer when ll_found=F 
    3438 
    3539   REAL(wp), PUBLIC ::   rho_c = 0.01_wp    !: density criterion for mixed layer depth 
     
    4953      zdf_mxl_alloc = 0      ! set to zero if no array to be allocated 
    5054      IF( .NOT. ALLOCATED( nmln ) ) THEN 
    51          ALLOCATE( nmln(jpi,jpj), hmld(jpi,jpj), hmlp(jpi,jpj), hmlpt(jpi,jpj), STAT= zdf_mxl_alloc ) 
     55         ALLOCATE( nmln(jpi,jpj), hmld(jpi,jpj), hmlp(jpi,jpj), hmlpt(jpi,jpj), hmld_zint(jpi,jpj),       & 
     56        &          ll_found(jpi,jpj), ll_belowml(jpi,jpj,jpk), STAT= zdf_mxl_alloc ) 
    5257         ! 
    5358         IF( lk_mpp             )   CALL mpp_sum ( zdf_mxl_alloc ) 
     
    145150         END IF 
    146151      ENDIF 
    147        
     152      ! 
     153      ! Vertically-interpolated mixed-layer depth diagnostic 
     154      IF( iom_use( "mldzint" ) ) THEN 
     155         CALL zdf_mxl_zint( kt ) 
     156         CALL iom_put( "mldzint" , hmld_zint ) 
     157      ENDIF 
     158      ! 
    148159      IF(ln_ctl)   CALL prt_ctl( tab2d_1=REAL(nmln,wp), clinfo1=' nmln : ', tab2d_2=hmlp, clinfo2=' hmlp : ', ovlap=1 ) 
    149160      ! 
     
    153164      ! 
    154165   END SUBROUTINE zdf_mxl 
     166 
     167   SUBROUTINE zdf_mxl_zint( kt )  
     168      !!----------------------------------------------------------------------------------  
     169      !!                    ***  ROUTINE zdf_mxl_zint  ***  
     170      !                                                                         
     171      !   Calculate vertically-interpolated mixed layer depth diagnostic.  
     172      !             
     173      !   This routine can calculate the mixed layer depth diagnostic suggested by 
     174      !   Kara et al, 2000, JGR, 105, 16803, but is more general and can calculate 
     175      !   vertically-interpolated mixed-layer depth diagnostics with other parameter 
     176      !   settings set in the namzdf_mldzint namelist.   
     177      !  
     178      !   If mld_type=1 the mixed layer depth is calculated as the depth at which the   
     179      !   density has increased by an amount equivalent to a temperature difference of   
     180      !   0.8C at the surface.  
     181      !  
     182      !   For other values of mld_type the mixed layer is calculated as the depth at   
     183      !   which the temperature differs by 0.8C from the surface temperature.   
     184      !                                                                         
     185      !   David Acreman, Daley Calvert                                       
     186      !  
     187      !!-----------------------------------------------------------------------------------  
     188 
     189      INTEGER, INTENT(in) ::   kt   ! ocean time-step index 
     190      ! 
     191      ! Local variables 
     192      INTEGER, POINTER, DIMENSION(:,:) :: ikmt          ! number of active tracer levels  
     193      INTEGER, POINTER, DIMENSION(:,:) :: ik_ref        ! index of reference level  
     194      INTEGER, POINTER, DIMENSION(:,:) :: ik_iso        ! index of last uniform temp level  
     195      REAL, POINTER, DIMENSION(:,:,:)  :: zT            ! Temperature or density  
     196      REAL, POINTER, DIMENSION(:,:)    :: ppzdep        ! depth for use in calculating d(rho)  
     197      REAL, POINTER, DIMENSION(:,:)    :: zT_ref        ! reference temperature  
     198      REAL    :: zT_b                                   ! base temperature  
     199      REAL, POINTER, DIMENSION(:,:,:)  :: zdTdz         ! gradient of zT  
     200      REAL, POINTER, DIMENSION(:,:,:)  :: zmoddT        ! Absolute temperature difference  
     201      REAL    :: zdz                                    ! depth difference  
     202      REAL    :: zdT                                    ! temperature difference  
     203      REAL, POINTER, DIMENSION(:,:)    :: zdelta_T      ! difference critereon  
     204      REAL, POINTER, DIMENSION(:,:)    :: zRHO1, zRHO2  ! Densities  
     205      INTEGER :: ji, jj, jk                             ! loop counter  
     206      INTEGER :: ios 
     207 
     208      NAMELIST/namzdf_mldzint/ nn_mld_type, rn_zref, rn_dT_crit, rn_iso_frac 
     209 
     210      !!-------------------------------------------------------------------------------------  
     211      !   
     212      CALL wrk_alloc( jpi, jpj, ikmt, ik_ref, ik_iso)  
     213      CALL wrk_alloc( jpi, jpj, ppzdep, zT_ref, zdelta_T, zRHO1, zRHO2 )  
     214      CALL wrk_alloc( jpi, jpj, jpk, zT, zdTdz, zmoddT )  
     215  
     216      IF( kt == nit000 ) THEN 
     217         REWIND( numnam_ref )              ! Namelist namzdf_mldzint in reference namelist  
     218         READ  ( numnam_ref, namzdf_mldzint, IOSTAT = ios, ERR = 901) 
     219901      IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_mldzint in reference namelist', lwp ) 
     220 
     221         REWIND( numnam_cfg )              ! Namelist namzdf_mldzint in configuration namelist  
     222         READ  ( numnam_cfg, namzdf_mldzint, IOSTAT = ios, ERR = 902 ) 
     223902      IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_mldzint in configuration namelist', lwp ) 
     224         IF(lwm) WRITE ( numond, namzdf_mldzint ) 
     225 
     226         WRITE(numout,*) '===== Vertically-interpolated mixed layer =====' 
     227         WRITE(numout,*) 'nn_mld_type = ',nn_mld_type 
     228         WRITE(numout,*) 'rn_zref = ',rn_zref 
     229         WRITE(numout,*) 'rn_dT_crit = ',rn_dT_crit 
     230         WRITE(numout,*) 'rn_iso_frac = ',rn_iso_frac 
     231         WRITE(numout,*) '===============================================' 
     232      ENDIF 
     233  
     234      ! Set the mixed layer depth criterion at each grid point  
     235      IF (nn_mld_type == 1) THEN                                          
     236         ppzdep(:,:)=0.0  
     237         call eos ( tsn(:,:,1,:), ppzdep(:,:), zRHO1(:,:) )  
     238! Use zT temporarily as a copy of tsn with rn_dT_crit added to SST  
     239! [assumes number of tracers less than number of vertical levels]  
     240         zT(:,:,1:jpts)=tsn(:,:,1,1:jpts)  
     241         zT(:,:,jp_tem)=zT(:,:,1)+rn_dT_crit  
     242         CALL eos( zT(:,:,1:jpts), ppzdep(:,:), zRHO2(:,:) )  
     243         zdelta_T(:,:) = abs( zRHO1(:,:) - zRHO2(:,:) ) * rau0  
     244         ! RHO from eos (2d version) doesn't calculate north or east halo:  
     245         CALL lbc_lnk( zdelta_T, 'T', 1. )  
     246         zT(:,:,:) = rhop(:,:,:)  
     247      ELSE  
     248         zdelta_T(:,:) = rn_dT_crit                       
     249         zT(:,:,:) = tsn(:,:,:,jp_tem)                            
     250      END IF  
     251 
     252      ! Calculate the gradient of zT and absolute difference for use later  
     253      DO jk = 1 ,jpk-2  
     254         zdTdz(:,:,jk)  =    ( zT(:,:,jk+1) - zT(:,:,jk) ) / e3w_n(:,:,jk+1)  
     255         zmoddT(:,:,jk) = abs( zT(:,:,jk+1) - zT(:,:,jk) )  
     256      END DO  
     257 
     258      ! Find density/temperature at the reference level (Kara et al use 10m).           
     259      ! ik_ref is the index of the box centre immediately above or at the reference level  
     260      ! Find rn_zref in the array of model level depths and find the ref     
     261      ! density/temperature by linear interpolation.                                    
     262      DO jk = jpkm1, 2, -1  
     263         WHERE ( gdept_n(:,:,jk) > rn_zref )  
     264           ik_ref(:,:) = jk - 1  
     265           zT_ref(:,:) = zT(:,:,jk-1) + zdTdz(:,:,jk-1) * ( rn_zref - gdept_n(:,:,jk-1) )  
     266         END WHERE  
     267      END DO  
     268 
     269      ! If the first grid box centre is below the reference level then use the  
     270      ! top model level to get zT_ref  
     271      WHERE ( gdept_n(:,:,1) > rn_zref )   
     272         zT_ref = zT(:,:,1)  
     273         ik_ref = 1  
     274      END WHERE  
     275 
     276      ! The number of active tracer levels is 1 less than the number of active w levels  
     277      ikmt(:,:) = mbathy(:,:) - 1  
     278 
     279      ! Search for a uniform density/temperature region where adjacent levels           
     280      ! differ by less than rn_iso_frac * deltaT.                                       
     281      ! ik_iso is the index of the last level in the uniform layer   
     282      ! ll_found indicates whether the mixed layer depth can be found by interpolation  
     283      ik_iso(:,:)   = ik_ref(:,:)  
     284      ll_found(:,:) = .false.  
     285      DO jj = 1, nlcj  
     286         DO ji = 1, nlci  
     287!CDIR NOVECTOR  
     288            DO jk = ik_ref(ji,jj), ikmt(ji,jj)-1  
     289               IF ( zmoddT(ji,jj,jk) > ( rn_iso_frac * zdelta_T(ji,jj) ) ) THEN  
     290                  ik_iso(ji,jj)   = jk  
     291                  ll_found(ji,jj) = ( zmoddT(ji,jj,jk) > zdelta_T(ji,jj) )  
     292                  EXIT  
     293               END IF  
     294            END DO  
     295         END DO  
     296      END DO  
     297 
     298      ! Use linear interpolation to find depth of mixed layer base where possible  
     299      hmld_zint(:,:) = rn_zref  
     300      DO jj = 1, jpj  
     301         DO ji = 1, jpi  
     302            IF (ll_found(ji,jj) .and. tmask(ji,jj,1) == 1.0) THEN  
     303               zdz =  abs( zdelta_T(ji,jj) / zdTdz(ji,jj,ik_iso(ji,jj)) )  
     304               hmld_zint(ji,jj) = gdept_n(ji,jj,ik_iso(ji,jj)) + zdz  
     305            END IF  
     306         END DO  
     307      END DO  
     308 
     309      ! If ll_found = .false. then calculate MLD using difference of zdelta_T     
     310      ! from the reference density/temperature  
     311  
     312! Prevent this section from working on land points  
     313      WHERE ( tmask(:,:,1) /= 1.0 )  
     314         ll_found = .true.  
     315      END WHERE  
     316  
     317      DO jk=1, jpk  
     318         ll_belowml(:,:,jk) = abs( zT(:,:,jk) - zT_ref(:,:) ) >= zdelta_T(:,:)   
     319      END DO  
     320  
     321! Set default value where interpolation cannot be used (ll_found=false)   
     322      DO jj = 1, jpj  
     323         DO ji = 1, jpi  
     324            IF ( .not. ll_found(ji,jj) )  hmld_zint(ji,jj) = gdept_n(ji,jj,ikmt(ji,jj))  
     325         END DO  
     326      END DO  
     327 
     328      DO jj = 1, jpj  
     329         DO ji = 1, jpi  
     330!CDIR NOVECTOR  
     331            DO jk = ik_ref(ji,jj)+1, ikmt(ji,jj)  
     332               IF ( ll_found(ji,jj) ) EXIT  
     333               IF ( ll_belowml(ji,jj,jk) ) THEN                 
     334                  zT_b = zT_ref(ji,jj) + zdelta_T(ji,jj) * SIGN(1.0, zdTdz(ji,jj,jk-1) )  
     335                  zdT  = zT_b - zT(ji,jj,jk-1)                                       
     336                  zdz  = zdT / zdTdz(ji,jj,jk-1)                                        
     337                  hmld_zint(ji,jj) = gdept_n(ji,jj,jk-1) + zdz  
     338                  EXIT                                                    
     339               END IF  
     340            END DO  
     341         END DO  
     342      END DO  
     343 
     344      hmld_zint(:,:) = hmld_zint(:,:)*tmask(:,:,1)  
     345      !   
     346      CALL wrk_dealloc( jpi, jpj, ikmt, ik_ref, ik_iso)  
     347      CALL wrk_dealloc( jpi, jpj, ppzdep, zT_ref, zdelta_T, zRHO1, zRHO2 )  
     348      CALL wrk_dealloc( jpi,jpj, jpk, zT, zdTdz, zmoddT )  
     349      !  
     350   END SUBROUTINE zdf_mxl_zint 
     351 
    155352 
    156353   !!====================================================================== 
  • branches/UKMO/dev_r6501_GO6_package_trunk/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90

    r6503 r6507  
    8181   INTEGER  ::   nn_htau   ! type of tke profile of penetration (=0/1) 
    8282   REAL(wp) ::   rn_efr    ! fraction of TKE surface value which penetrates in the ocean 
     83   REAL(wp) ::   rn_c      ! fraction of TKE added within the mixed layer by nn_etau 
    8384   LOGICAL  ::   ln_lc     ! Langmuir cells (LC) as a source term of TKE or not 
    8485   REAL(wp) ::   rn_lc     ! coef to compute vertical velocity of Langmuir cells 
     
    8687   REAL(wp) ::   ri_cri    ! critic Richardson number (deduced from rn_ediff and rn_ediss values) 
    8788   REAL(wp) ::   rmxl_min  ! minimum mixing length value (deduced from rn_ediff and rn_emin values)  [m] 
     89   REAL(wp) ::   rhtau                     ! coefficient to relate MLD to htau when nn_htau == 2 
    8890   REAL(wp) ::   rhftau_add = 1.e-3_wp     ! add offset   applied to HF part of taum  (nn_etau=3) 
    8991   REAL(wp) ::   rhftau_scl = 1.0_wp       ! scale factor applied to HF part of taum  (nn_etau=3) 
    9092 
    9193   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   htau           ! depth of tke penetration (nn_htau) 
     94   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   efr            ! surface boundary condition for nn_etau = 4 
    9295   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   dissl          ! now mixing lenght of dissipation 
    9396   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   apdlr          ! now mixing lenght of dissipation 
     
    112115      !!---------------------------------------------------------------------- 
    113116      ALLOCATE(                                                                    & 
     117         &      efr  (jpi,jpj)     ,                                               &       
    114118#if defined key_c1d 
    115119         &      e_dis(jpi,jpj,jpk) , e_mix(jpi,jpj,jpk) ,                          & 
     
    446450      !                            !  TKE due to surface and internal wave breaking 
    447451      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
     452      IF( nn_htau == 2 ) THEN           !* mixed-layer depth dependant length scale 
     453         DO jj = 2, jpjm1 
     454            DO ji = fs_2, fs_jpim1   ! vector opt. 
     455               htau(ji,jj) = rhtau * hmlp(ji,jj) 
     456            END DO 
     457         END DO 
     458      ENDIF 
     459#if defined key_iomput 
     460      ! 
     461      CALL iom_put( "htau", htau(:,:) )  ! Check htau (even if constant in time) 
     462#endif 
     463      ! 
    448464!!gm BUG : in the exp  remove the depth of ssh !!! 
    449        
    450465       
    451466      IF( nn_etau == 1 ) THEN           !* penetration below the mixed layer (rn_efr fraction) 
     
    477492                  en(ji,jj,jk) = en(ji,jj,jk) + zbbrau * zdif * EXP( -gdepw_n(ji,jj,jk) / htau(ji,jj) )   & 
    478493                     &                        * ( 1._wp - fr_i(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1) 
     494               END DO 
     495            END DO 
     496         END DO 
     497      ELSEIF( nn_etau == 4 ) THEN       !* column integral independant of htau (rn_efr must be scaled up) 
     498         IF( nn_htau == 2 ) THEN        ! efr dependant on time-varying htau  
     499            DO jj = 2, jpjm1 
     500               DO ji = fs_2, fs_jpim1   ! vector opt. 
     501                  efr(ji,jj) = rn_efr / ( htau(ji,jj) * ( 1._wp - EXP( -bathy(ji,jj) / htau(ji,jj) ) ) ) 
     502               END DO 
     503            END DO 
     504         ENDIF 
     505         DO jk = 2, jpkm1 
     506            DO jj = 2, jpjm1 
     507               DO ji = fs_2, fs_jpim1   ! vector opt. 
     508                  en(ji,jj,jk) = en(ji,jj,jk) + efr(ji,jj) * en(ji,jj,1) * EXP( -gdepw_n(ji,jj,jk) / htau(ji,jj) )   & 
     509                     &                                                   * ( 1._wp - fr_i(ji,jj) )  * tmask(ji,jj,jk) 
    479510               END DO 
    480511            END DO 
     
    723754      !!---------------------------------------------------------------------- 
    724755      INTEGER ::   ji, jj, jk   ! dummy loop indices 
    725       INTEGER ::   ios 
     756      INTEGER ::   ios, ierr 
    726757      !! 
    727758      NAMELIST/namzdf_tke/ rn_ediff, rn_ediss , rn_ebb , rn_emin  ,   & 
    728759         &                 rn_emin0, rn_bshear, nn_mxl , ln_mxl0  ,   & 
    729760         &                 rn_mxl0 , nn_pdl   , ln_lc  , rn_lc    ,   & 
    730          &                 nn_etau , nn_htau  , rn_efr    
     761         &                 nn_etau , nn_htau  , rn_efr , rn_c    
    731762      !!---------------------------------------------------------------------- 
    732763      ! 
     
    774805         WRITE(numout,*) '      flag for computation of exp. tke profile    nn_htau   = ', nn_htau 
    775806         WRITE(numout,*) '      fraction of en which pene. the thermocline  rn_efr    = ', rn_efr 
     807         WRITE(numout,*) '      fraction of TKE added within the mixed layer by nn_etau rn_c    = ', rn_c 
    776808         WRITE(numout,*) 
    777809         WRITE(numout,*) '      critical Richardson nb with your parameters  ri_cri = ', ri_cri 
     
    784816      IF( nn_mxl  < 0   .OR.  nn_mxl  > 3 )   CALL ctl_stop( 'bad flag: nn_mxl is  0, 1 or 2 ' ) 
    785817      IF( nn_pdl  < 0   .OR.  nn_pdl  > 1 )   CALL ctl_stop( 'bad flag: nn_pdl is  0 or 1    ' ) 
    786       IF( nn_htau < 0   .OR.  nn_htau > 1 )   CALL ctl_stop( 'bad flag: nn_htau is 0, 1 or 2 ' ) 
     818      IF( nn_htau < 0  .OR.  nn_htau > 5 )   CALL ctl_stop( 'bad flag: nn_htau is 0 to 5    ' ) 
    787819      IF( nn_etau == 3 .AND. .NOT. ln_cpl )   CALL ctl_stop( 'nn_etau == 3 : HF taum only known in coupled mode' ) 
    788820 
     
    792824      ENDIF 
    793825       
    794       IF( nn_etau == 2  )   CALL zdf_mxl( nit000 )      ! Initialization of nmln  
     826      IF( nn_etau == 2  ) THEN 
     827          ierr = zdf_mxl_alloc() 
     828          nmln(:,:) = nlb10           ! Initialization of nmln 
     829      ENDIF 
     830 
     831      IF( nn_etau /= 0 .and. nn_htau == 2 ) THEN 
     832          ierr = zdf_mxl_alloc() 
     833          nmln(:,:) = nlb10           ! Initialization of nmln 
     834      ENDIF 
    795835 
    796836      !                               !* depth of penetration of surface tke 
    797837      IF( nn_etau /= 0 ) THEN       
     838         htau(:,:) = 0._wp 
    798839         SELECT CASE( nn_htau )             ! Choice of the depth of penetration 
    799840         CASE( 0 )                                 ! constant depth penetration (here 10 meters) 
     
    801842         CASE( 1 )                                 ! F(latitude) : 0.5m to 30m poleward of 40 degrees 
    802843            htau(:,:) = MAX(  0.5_wp, MIN( 30._wp, 45._wp* ABS( SIN( rpi/180._wp * gphit(:,:) ) ) )   )             
     844         CASE( 2 )                                 ! fraction of depth-integrated TKE within mixed-layer 
     845            rhtau = -1._wp / LOG( 1._wp - rn_c ) 
     846         CASE( 3 )                                 ! F(latitude) : 0.5m to 15m poleward of 20 degrees 
     847            htau(:,:) = MAX(  0.5_wp, MIN( 15._wp, 45._wp* ABS( SIN( rpi/180._wp * gphit(:,:) ) ) )   ) 
     848         CASE( 4 )                                 ! F(latitude) : 0.5m to 10m/30m poleward of 13/40 degrees north/south 
     849            DO jj = 2, jpjm1 
     850               DO ji = fs_2, fs_jpim1   ! vector opt. 
     851                  IF( gphit(ji,jj) <= 0._wp ) THEN 
     852                     htau(ji,jj) = MAX(  0.5_wp, MIN( 30._wp, 45._wp* ABS( SIN( rpi/180._wp * gphit(ji,jj) ) ) )   ) 
     853                  ELSE 
     854                     htau(ji,jj) = MAX(  0.5_wp, MIN( 10._wp, 45._wp* ABS( SIN( rpi/180._wp * gphit(ji,jj) ) ) )   ) 
     855                  ENDIF 
     856               END DO 
     857            END DO 
     858         CASE ( 5 )                                ! F(latitude) : 0.5m to 10m poleward of 13 degrees north/south, 
     859            DO jj = 2, jpjm1                       !               10m to 30m between 30/45 degrees south 
     860               DO ji = fs_2, fs_jpim1   ! vector opt. 
     861                  IF( gphit(ji,jj) <= -30._wp ) THEN 
     862                     htau(ji,jj) = MAX(  10._wp, MIN( 30._wp, 55._wp* ABS( SIN( rpi/120._wp * ( gphit(ji,jj) + 23._wp ) ) ) )   ) 
     863                  ELSE 
     864                     htau(ji,jj) = MAX(  0.5_wp, MIN( 10._wp, 45._wp* ABS( SIN( rpi/180._wp * gphit(ji,jj) ) ) )   ) 
     865                  ENDIF 
     866               END DO 
     867            END DO 
    803868         END SELECT 
     869         ! 
     870         IF( nn_etau == 4 .AND. nn_htau /= 2 ) THEN            ! efr dependant on constant htau 
     871            DO jj = 2, jpjm1 
     872               DO ji = fs_2, fs_jpim1   ! vector opt. 
     873                  efr(ji,jj) = rn_efr / ( htau(ji,jj) * ( 1._wp - EXP( -bathy(ji,jj) / htau(ji,jj) ) ) ) 
     874               END DO 
     875            END DO 
     876         ENDIF 
    804877      ENDIF 
    805878      !                               !* set vertical eddy coef. to the background value 
  • branches/UKMO/dev_r6501_GO6_package_trunk/NEMOGCM/NEMO/OPA_SRC/step_oce.F90

    r6503 r6507  
    8080 
    8181   USE diawri           ! Standard run outputs             (dia_wri routine) 
     82   USE diaprod          ! Product diagnostics              (dia_prod routine) 
    8283   USE diaptr           ! poleward transports              (dia_ptr routine) 
    8384   USE diadct           ! sections transports              (dia_dct routine) 
Note: See TracChangeset for help on using the changeset viewer.