New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Changeset 6491 for branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO – NEMO

Ignore:
Timestamp:
2016-04-21T18:15:17+02:00 (8 years ago)
Author:
davestorkey
Message:

Commit remaining changes to UKMO/r5518_GO6_package branch from following branches:
UKMO/dev_r5021_nn_etau_revision@6238
UKMO/dev_r5107_mld_zint@5534
UKMO/dev_r5107_eorca025_closea@6390
UKMO/restart_datestamp@5539
UKMO/icebergs_latent_heat@5821
UKMO/icebergs_restart_single_file_corrected@6480
UKMO/product_diagnostics@5971
UKMO/antarctic_partial_slip@5961
Custom merge into /branches/UKMO/dev_r5518_GO6_package/NEMOGCM: r5961 cf. r5958 of /branches/UKMO/antarctic_partial_slip/NEMOGCM@6490

Custom merge into /branches/UKMO/dev_r5518_GO6_package/NEMOGCM: r6349 cf. r5962 of /branches/UKMO/product_diagnostics/NEMOGCM@6490

Custom merge into /branches/UKMO/dev_r5518_GO6_package/NEMOGCM: r6480 cf. r6479 of /branches/UKMO/icebergs_restart_single_file_corrected/NEMOGCM@6490

Custom merge into /branches/UKMO/dev_r5518_GO6_package/NEMOGCM: r5986 cf. r5852 of /branches/UKMO/icebergs_restart_single_file/NEMOGCM@6490

Custom merge into /branches/UKMO/dev_r5518_GO6_package/NEMOGCM: r5821 cf. r5808 of /branches/UKMO/icebergs_latent_heat/NEMOGCM@6490

Location:
branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC
Files:
17 edited
1 copied

Legend:

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

    r6487 r6491  
    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 
     
    242243      CALL iom_put( "avt" , avt                        )    ! T vert. eddy diff. coef. 
    243244      CALL iom_put( "avm" , avmu                       )    ! T vert. eddy visc. coef. 
     245      IF( lk_zdftke ) THEN    
     246         CALL iom_put( "tke"      , en                               )    ! TKE budget: Turbulent Kinetic Energy    
     247         CALL iom_put( "tke_niw"  , e_niw                            )    ! TKE budget: Near-inertial waves    
     248      ENDIF  
    244249      CALL iom_put( "avs" , fsavs(:,:,:)               )    ! S vert. eddy diff. coef. (useful only with key_zdfddm) 
    245250 
  • branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/DOM/closea.F90

    r6486 r6491  
    3636   PUBLIC clo_bat      ! routine called in domzgr module 
    3737 
    38    INTEGER, PUBLIC, PARAMETER          ::   jpncs   = 4      !: number of closed sea 
     38   INTEGER, PUBLIC, PARAMETER          ::   jpncs   = 10      !: number of closed sea 
    3939   INTEGER, PUBLIC, DIMENSION(jpncs)   ::   ncstt            !: Type of closed sea 
    4040   INTEGER, PUBLIC, DIMENSION(jpncs)   ::   ncsi1, ncsj1     !: south-west closed sea limits (i,j) 
     
    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_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90

    r6486 r6491  
    136136      USE ioipsl 
    137137      NAMELIST/namrun/ cn_ocerst_indir, cn_ocerst_outdir, nn_stocklist, ln_rst_list,               & 
    138          &             nn_no   , cn_exp    , cn_ocerst_in, cn_ocerst_out, ln_rstart , nn_rstctl,   & 
     138         &             nn_no   , cn_exp    , cn_ocerst_in, cn_ocerst_out, ln_rstart , ln_rstdate, nn_rstctl,   & 
    139139         &             nn_it000, nn_itend  , nn_date0    , nn_leapy     , nn_istate , nn_stock ,   & 
    140140         &             nn_write, ln_dimgnnn, ln_mskland  , ln_cfmeta    , ln_clobber, nn_chunksz, nn_euler 
     
    173173         WRITE(numout,*) '      file prefix restart output      cn_ocerst_out= ', cn_ocerst_out 
    174174         WRITE(numout,*) '      restart output directory        cn_ocerst_outdir= ', cn_ocerst_outdir 
    175          WRITE(numout,*) '      restart logical                 ln_rstart  = ', ln_rstart 
     175         WRITE(numout,*) '      restart logical                 ln_rstart  = ' , ln_rstart 
     176         WRITE(numout,*) '      datestamping of restarts        ln_rstdate  = ', ln_rstdate 
    176177         WRITE(numout,*) '      start with forward time step    nn_euler   = ', nn_euler 
    177178         WRITE(numout,*) '      control of time step            nn_rstctl  = ', nn_rstctl 
  • branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/DOM/dommsk.F90

    r6487 r6491  
    136136      INTEGER  ::   isrow                    ! index for ORCA1 starting row 
    137137      INTEGER , POINTER, DIMENSION(:,:) ::  imsk 
     138      REAL(wp) ::  zphi_drake_passage, zshlat_antarc 
    138139      REAL(wp), POINTER, DIMENSION(:,:) ::  zwf 
    139140      !! 
     
    445446      ENDIF 
    446447      ! 
     448      IF( cp_cfg == "orca" .AND. jp_cfg == 025 .AND. rn_shlat == 0.0 ) THEN    
     449         !                                              ! ORCA_R025 configuration 
     450         !                                              ! Increased lateral friction on parts of Antarctic coastline 
     451         !                                              ! for increased stability 
     452         !                                              ! NB. This only works to do this here if we have free slip  
     453         !                                              ! generally, so fmask is zero at coast points. 
     454         IF(lwp) WRITE(numout,*) 
     455         IF(lwp) WRITE(numout,*) '   orca_r025: increase friction in following regions : ' 
     456         IF(lwp) WRITE(numout,*) '      whole Antarctic coastline: partial slip shlat=1 ' 
     457 
     458         zphi_drake_passage = -58.0_wp 
     459         zshlat_antarc = 1.0_wp 
     460         zwf(:,:) = fmask(:,:,1)          
     461         DO jj = 2, jpjm1 
     462            DO ji = fs_2, fs_jpim1   ! vector opt. 
     463               IF( gphif(ji,jj) .lt. zphi_drake_passage .and. fmask(ji,jj,1) == 0._wp ) THEN 
     464                  fmask(ji,jj,:) = zshlat_antarc * MIN( 1._wp , MAX( zwf(ji+1,jj), zwf(ji,jj+1),   & 
     465                     &                                           zwf(ji-1,jj), zwf(ji,jj-1)  )  ) 
     466               ENDIF 
     467            END DO 
     468         END DO 
     469      END IF 
     470      ! 
    447471      CALL lbc_lnk( fmask, 'F', 1._wp )      ! Lateral boundary conditions on fmask 
    448472 
  • branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90

    r6486 r6491  
    594594         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~   - interpolate scale factors and compute depths for next time step' 
    595595      ENDIF 
     596 
     597      ! Write outputs 
     598      ! ============= 
     599      CALL iom_put(     "e3t" , fse3t_n  (:,:,:) ) 
     600      CALL iom_put(     "e3u" , fse3u_n  (:,:,:) ) 
     601      CALL iom_put(     "e3v" , fse3v_n  (:,:,:) ) 
     602      CALL iom_put(     "e3w" , fse3w_n  (:,:,:) ) 
     603      CALL iom_put( "tpt_dep" , fsde3w_n (:,:,:) ) 
     604      IF( iom_use("e3tdef") )   & 
     605         CALL iom_put( "e3tdef"  , ( ( fse3t_n(:,:,:) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100 * tmask(:,:,:) ) ** 2 ) 
     606 
    596607      ! 
    597608      ! Time filter and swap of scale factors 
     
    666677      END DO 
    667678 
    668       ! Write outputs 
    669       ! ============= 
    670       CALL iom_put(     "e3t" , fse3t_n  (:,:,:) ) 
    671       CALL iom_put(     "e3u" , fse3u_n  (:,:,:) ) 
    672       CALL iom_put(     "e3v" , fse3v_n  (:,:,:) ) 
    673       CALL iom_put(     "e3w" , fse3w_n  (:,:,:) ) 
    674       CALL iom_put( "tpt_dep" , fsde3w_n (:,:,:) ) 
    675       IF( iom_use("e3tdef") )   & 
    676          CALL iom_put( "e3tdef"  , ( ( fse3t_n(:,:,:) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100 * tmask(:,:,:) ) ** 2 ) 
    677679 
    678680      ! write restart file 
  • branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/ICB/icbdia.F90

    r6486 r6491  
    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_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/ICB/icbrst.F90

    r6486 r6491  
    1212   !!            -    !                            Currently needs a fixed processor 
    1313   !!            -    !                            layout between restarts 
     14   !!            -    !  2015-11  Dave Storkey     Convert icb_rst_read to use IOM so can 
     15   !!                                              read single restart files 
    1416   !!---------------------------------------------------------------------- 
    1517   !!---------------------------------------------------------------------- 
     
    1820   !!---------------------------------------------------------------------- 
    1921   USE par_oce        ! NEMO parameters 
     22   USE phycst         ! for rday 
    2023   USE dom_oce        ! NEMO domain 
    2124   USE in_out_manager ! NEMO IO routines 
     25   USE ioipsl, ONLY : ju2ymds    ! for calendar 
    2226   USE lib_mpp        ! NEMO MPI library, lk_mpp in particular 
    2327   USE netcdf         ! netcdf routines for IO 
     28   USE iom 
    2429   USE icb_oce        ! define iceberg arrays 
    2530   USE icbutl         ! iceberg utility routines 
     
    5762      INTEGER                      ::   idim, ivar, iatt 
    5863      INTEGER                      ::   jn, iunlim_dim, ibergs_in_file 
    59       INTEGER                      ::   iclass 
    60       INTEGER, DIMENSION(1)        ::   istrt, ilngth, idata 
    61       INTEGER, DIMENSION(2)        ::   istrt2, ilngth2 
    62       INTEGER, DIMENSION(nkounts)  ::   idata2 
    63       REAL(wp), DIMENSION(1)       ::   zdata                                         ! need 1d array to read in with 
    64                                                                                             ! start and count arrays 
     64      INTEGER                      ::   ii,ij,iclass 
     65      REAL(wp), DIMENSION(nkounts) ::   zdata       
    6566      LOGICAL                      ::   ll_found_restart 
    6667      CHARACTER(len=256)           ::   cl_path 
     
    7172      !!---------------------------------------------------------------------- 
    7273 
    73       ! Find a restart file. Assume iceberg restarts in same directory as ocean restarts.  
     74      ! Find a restart file. Assume iceberg restarts in same directory as ocean restarts 
     75      ! and are called TRIM(cn_ocerst)//'_icebergs' 
    7476      cl_path = TRIM(cn_ocerst_indir) 
    7577      IF( cl_path(LEN_TRIM(cl_path):) /= '/' ) cl_path = TRIM(cl_path) // '/' 
    76       cl_filename = ' ' 
    77       IF ( lk_mpp ) THEN 
    78          cl_filename = ' ' 
    79          WRITE( cl_filename, '("restart_icebergs_",I4.4,".nc")' ) narea-1 
    80          INQUIRE( file=TRIM(cl_path)//TRIM(cl_filename), exist=ll_found_restart ) 
    81       ELSE 
    82          cl_filename = 'restart_icebergs.nc' 
    83          INQUIRE( file=TRIM(cl_path)//TRIM(cl_filename), exist=ll_found_restart ) 
    84       ENDIF 
    85  
    86       IF ( .NOT. ll_found_restart) THEN                     ! only do the following if a file was found 
    87          CALL ctl_stop('icebergs: no restart file found') 
    88       ENDIF 
    89  
    90       IF (nn_verbose_level >= 0 .AND. lwp)  & 
    91          WRITE(numout,'(2a)') 'icebergs, read_restart_bergs: found restart file = ',TRIM(cl_path)//TRIM(cl_filename) 
    92  
    93       nret = NF90_OPEN(TRIM(cl_path)//TRIM(cl_filename), NF90_NOWRITE, ncid) 
    94       IF (nret .ne. NF90_NOERR) CALL ctl_stop('icebergs, read_restart_bergs: nf_open failed') 
    95  
    96       nret = nf90_inquire(ncid, idim, ivar, iatt, iunlim_dim) 
    97       IF (nret .ne. NF90_NOERR) CALL ctl_stop('icebergs, read_restart_bergs: nf_inquire failed') 
    98  
    99       IF( iunlim_dim .NE. -1) THEN 
    100  
    101          nret = nf90_inquire_dimension(ncid, iunlim_dim, cl_dname, ibergs_in_file) 
    102          IF (nret .ne. NF90_NOERR) CALL ctl_stop('icebergs, read_restart_bergs: nf_inq_dimlen failed') 
    103  
    104          nret = NF90_INQ_VARID(ncid, 'number', numberid) 
    105          nret = NF90_INQ_VARID(ncid, 'mass_scaling', nscaling_id) 
    106          nret = NF90_INQ_VARID(ncid, 'xi', nxid) 
    107          nret = NF90_INQ_VARID(ncid, 'yj', nyid) 
    108          nret = NF90_INQ_VARID(ncid, 'lon', nlonid) 
    109          nret = NF90_INQ_VARID(ncid, 'lat', nlatid) 
    110          nret = NF90_INQ_VARID(ncid, 'uvel', nuvelid) 
    111          nret = NF90_INQ_VARID(ncid, 'vvel', nvvelid) 
    112          nret = NF90_INQ_VARID(ncid, 'mass', nmassid) 
    113          nret = NF90_INQ_VARID(ncid, 'thickness', nthicknessid) 
    114          nret = NF90_INQ_VARID(ncid, 'width', nwidthid) 
    115          nret = NF90_INQ_VARID(ncid, 'length', nlengthid) 
    116          nret = NF90_INQ_VARID(ncid, 'year', nyearid) 
    117          nret = NF90_INQ_VARID(ncid, 'day', ndayid) 
    118          nret = NF90_INQ_VARID(ncid, 'mass_of_bits', nmass_of_bits_id) 
    119          nret = NF90_INQ_VARID(ncid, 'heat_density', nheat_density_id) 
    120  
    121          ilngth(1) = 1 
    122          istrt2(1) = 1 
    123          ilngth2(1) = nkounts 
    124          ilngth2(2) = 1 
    125          DO jn=1, ibergs_in_file 
    126  
    127             istrt(1) = jn 
    128             istrt2(2) = jn 
    129  
    130             nret = NF90_GET_VAR(ncid, numberid, idata2, istrt2, ilngth2 ) 
    131             localberg%number(:) = idata2(:) 
    132  
    133             nret = NF90_GET_VAR(ncid, nscaling_id, zdata, istrt, ilngth ) 
    134             localberg%mass_scaling = zdata(1) 
    135  
    136             nret = NF90_GET_VAR(ncid, nlonid, zdata, istrt, ilngth) 
    137             localpt%lon = zdata(1) 
    138             nret = NF90_GET_VAR(ncid, nlatid, zdata, istrt, ilngth) 
    139             localpt%lat = zdata(1) 
    140             IF (nn_verbose_level >= 2 .AND. lwp) THEN 
    141                WRITE(numout,'(a,i5,a,2f10.4,a,i5)') 'icebergs, read_restart_bergs: berg ',jn,' is at ', & 
    142                                               localpt%lon,localpt%lat,' on PE ',narea-1 
     78      cl_filename = TRIM(cn_ocerst_in)//'_icebergs' 
     79      CALL iom_open( TRIM(cl_path)//cl_filename, ncid ) 
     80 
     81      IF( iom_file(ncid)%iduld .GE. 0) THEN 
     82 
     83         ibergs_in_file = iom_file(ncid)%lenuld 
     84         DO jn = 1,ibergs_in_file 
     85 
     86            ! iom_get treats the unlimited dimension as time. Here the unlimited dimension  
     87            ! is the iceberg index, but we can still use the ktime keyword to get the iceberg we want.  
     88 
     89            CALL iom_get( ncid, 'xi'     ,localpt%xi  , ktime=jn ) 
     90            CALL iom_get( ncid, 'yj'     ,localpt%yj  , ktime=jn ) 
     91 
     92            ii = INT( localpt%xi + 0.5 ) 
     93            ij = INT( localpt%yj + 0.5 ) 
     94            ! Only proceed if this iceberg is on the local processor (excluding halos). 
     95            IF ( ii .GE. nldi+nimpp-1 .AND. ii .LE. nlei+nimpp-1 .AND. & 
     96           &     ij .GE. nldj+njmpp-1 .AND. ij .LE. nlej+njmpp-1 ) THEN            
     97 
     98               CALL iom_get( ncid, jpdom_unknown, 'number'       , zdata(:) , ktime=jn, kstart=(/1/), kcount=(/nkounts/) ) 
     99               localberg%number(:) = INT(zdata(:)) 
     100               CALL iom_get( ncid, 'mass_scaling' , localberg%mass_scaling, ktime=jn ) 
     101               CALL iom_get( ncid, 'lon'          , localpt%lon           , ktime=jn ) 
     102               CALL iom_get( ncid, 'lat'          , localpt%lat           , ktime=jn ) 
     103               CALL iom_get( ncid, 'uvel'         , localpt%uvel          , ktime=jn ) 
     104               CALL iom_get( ncid, 'vvel'         , localpt%vvel          , ktime=jn ) 
     105               CALL iom_get( ncid, 'mass'         , localpt%mass          , ktime=jn ) 
     106               CALL iom_get( ncid, 'thickness'    , localpt%thickness     , ktime=jn ) 
     107               CALL iom_get( ncid, 'width'        , localpt%width         , ktime=jn ) 
     108               CALL iom_get( ncid, 'length'       , localpt%length        , ktime=jn ) 
     109               CALL iom_get( ncid, 'year'         , zdata(1)              , ktime=jn ) 
     110               localpt%year = INT(zdata(1)) 
     111               CALL iom_get( ncid, 'day'          , localpt%day           , ktime=jn ) 
     112               CALL iom_get( ncid, 'mass_of_bits' , localpt%mass_of_bits  , ktime=jn ) 
     113               CALL iom_get( ncid, 'heat_density' , localpt%heat_density  , ktime=jn ) 
     114 
     115               ! 
     116               CALL icb_utl_add( localberg, localpt ) 
     117 
    143118            ENDIF 
    144             nret = NF90_GET_VAR(ncid, nxid, zdata, istrt, ilngth) 
    145             localpt%xi = zdata(1) 
    146             nret = NF90_GET_VAR(ncid, nyid, zdata, istrt, ilngth) 
    147             localpt%yj = zdata(1) 
    148             nret = NF90_GET_VAR(ncid, nuvelid, zdata, istrt, ilngth ) 
    149             localpt%uvel = zdata(1) 
    150             nret = NF90_GET_VAR(ncid, nvvelid, zdata, istrt, ilngth ) 
    151             localpt%vvel = zdata(1) 
    152             nret = NF90_GET_VAR(ncid, nmassid, zdata, istrt, ilngth ) 
    153             localpt%mass = zdata(1) 
    154             nret = NF90_GET_VAR(ncid, nthicknessid, zdata, istrt, ilngth ) 
    155             localpt%thickness = zdata(1) 
    156             nret = NF90_GET_VAR(ncid, nwidthid, zdata, istrt, ilngth ) 
    157             localpt%width = zdata(1) 
    158             nret = NF90_GET_VAR(ncid, nlengthid, zdata, istrt, ilngth ) 
    159             localpt%length = zdata(1) 
    160             nret = NF90_GET_VAR(ncid, nyearid, idata, istrt, ilngth ) 
    161             localpt%year = idata(1) 
    162             nret = NF90_GET_VAR(ncid, ndayid, zdata, istrt, ilngth ) 
    163             localpt%day = zdata(1) 
    164             nret = NF90_GET_VAR(ncid, nmass_of_bits_id, zdata, istrt, ilngth ) 
    165             localpt%mass_of_bits = zdata(1) 
    166             nret = NF90_GET_VAR(ncid, nheat_density_id, zdata, istrt, ilngth ) 
    167             localpt%heat_density = zdata(1) 
    168             ! 
    169             CALL icb_utl_add( localberg, localpt ) 
     119 
    170120         END DO 
    171          ! 
    172       ENDIF 
    173  
    174       nret = NF90_INQ_DIMID( ncid, 'c', nc_dim ) 
    175       IF (nret .ne. NF90_NOERR) CALL ctl_stop('icebergs, read_restart: nf_inq_dimid c failed') 
    176  
    177       nret = NF90_INQUIRE_DIMENSION( ncid, nc_dim, cl_dname, iclass ) 
    178       IF (nret .ne. NF90_NOERR) CALL ctl_stop('icebergs, read_restart: nf_inquire_dimension failed') 
    179  
    180       nret = NF90_INQ_VARID(ncid, 'kount'       , nkountid) 
    181       nret = NF90_INQ_VARID(ncid, 'calving'     , ncalvid) 
    182       nret = NF90_INQ_VARID(ncid, 'calving_hflx', ncalvhid) 
    183       nret = NF90_INQ_VARID(ncid, 'stored_ice'  , nsiceid) 
    184       nret = NF90_INQ_VARID(ncid, 'stored_heat' , nsheatid) 
    185  
    186       nstrt3(1) = 1 
    187       nstrt3(2) = 1 
    188       nlngth3(1) = jpi 
    189       nlngth3(2) = jpj 
    190       nlngth3(3) = 1 
    191  
    192       DO jn = 1, iclass 
    193          nstrt3(3) = jn 
    194          nret      = NF90_GET_VAR( ncid, nsiceid , griddata, nstrt3, nlngth3 ) 
    195          berg_grid%stored_ice(:,:,jn) = griddata(:,:,1) 
    196       END DO 
    197  
    198       nret = NF90_GET_VAR( ncid, ncalvid , src_calving          (:,:) ) 
    199       nret = NF90_GET_VAR( ncid, ncalvhid, src_calving_hflx     (:,:) ) 
    200       nret = NF90_GET_VAR( ncid, nsheatid, berg_grid%stored_heat(:,:) ) 
    201       nret = NF90_GET_VAR( ncid, nkountid, idata2(:) ) 
    202       num_bergs(:) = idata2(:) 
    203  
    204       ! Finish up 
    205       nret = NF90_CLOSE(ncid) 
    206       IF (nret .ne. NF90_NOERR) CALL ctl_stop('icebergs, read_restart: nf_close failed') 
     121 
     122      ENDIF  
     123 
     124      ! Gridded variables 
     125      CALL iom_get( ncid, jpdom_autoglo,    'calving'     , src_calving  ) 
     126      CALL iom_get( ncid, jpdom_autoglo,    'calving_hflx', src_calving_hflx  ) 
     127      CALL iom_get( ncid, jpdom_autoglo,    'stored_heat' , berg_grid%stored_heat  ) 
     128      CALL iom_get( ncid, jpdom_autoglo_xy, 'stored_ice'  , berg_grid%stored_ice, kstart=(/1,1,1/), kcount=(/1,1,nclasses/) ) 
     129       
     130      CALL iom_get( ncid, jpdom_unknown, 'kount' , zdata(:) ) 
     131      num_bergs(:) = INT(zdata(:)) 
    207132 
    208133      ! Sanity check 
     
    211136         WRITE(numout,'(2(a,i5))') 'icebergs, read_restart_bergs: # bergs =',jn,' on PE',narea-1 
    212137      IF( lk_mpp ) THEN 
    213          CALL mpp_sum(ibergs_in_file) 
     138         ! Only mpp_sum ibergs_in_file if we are reading from multiple restart files.  
     139         IF( INDEX(iom_file(ncid)%name,'icebergs.nc' ) .EQ. 0 ) CALL mpp_sum(ibergs_in_file) 
    214140         CALL mpp_sum(jn) 
    215141      ENDIF 
     
    217143         &                                    ' bergs in the restart file and', jn,' bergs have been read' 
    218144      ! 
     145      ! Finish up 
     146      CALL iom_close( ncid ) 
     147      ! 
    219148      IF( lwp .and. nn_verbose_level >= 0)  WRITE(numout,'(a)') 'icebergs, read_restart_bergs: completed' 
    220149      ! 
     
    231160      INTEGER ::   jn   ! dummy loop index 
    232161      INTEGER ::   ix_dim, iy_dim, ik_dim, in_dim 
    233       CHARACTER(len=256)     :: cl_path 
    234       CHARACTER(len=256)     :: cl_filename 
     162      INTEGER             ::   iyear, imonth, iday 
     163      REAL (wp)           ::   zsec 
     164      CHARACTER(len=256)  :: cl_path 
     165      CHARACTER(len=256)  :: cl_filename 
     166      CHARACTER(LEN=20)   ::   clkt     ! ocean time-step deine as a character 
    235167      TYPE(iceberg), POINTER :: this 
    236168      TYPE(point)  , POINTER :: pt 
     
    240172      cl_path = TRIM(cn_ocerst_outdir) 
    241173      IF( cl_path(LEN_TRIM(cl_path):) /= '/' ) cl_path = TRIM(cl_path) // '/' 
     174      IF ( ln_rstdate ) THEN 
     175         CALL ju2ymds( fjulday + rdttra(1) / rday, iyear, imonth, iday, zsec )            
     176         WRITE(clkt, '(i4.4,2i2.2)') iyear, imonth, iday 
     177      ELSE 
     178         IF( kt > 999999999 ) THEN   ;   WRITE(clkt, *       ) kt 
     179         ELSE                        ;   WRITE(clkt, '(i8.8)') kt 
     180         ENDIF 
     181      ENDIF 
    242182      IF( lk_mpp ) THEN 
    243          WRITE(cl_filename,'(A,"_icebergs_",I8.8,"_restart_",I4.4,".nc")') TRIM(cexper), kt, narea-1 
     183         WRITE(cl_filename,'(A,"_icebergs_",A,"_restart_",I4.4,".nc")') TRIM(cexper), TRIM(ADJUSTL(clkt)), narea-1 
    244184      ELSE 
    245          WRITE(cl_filename,'(A,"_icebergs_",I8.8,"_restart.nc")') TRIM(cexper), kt 
     185         WRITE(cl_filename,'(A,"_icebergs_",A,"_restart.nc")') TRIM(cexper), TRIM(ADJUSTL(clkt)) 
    246186      ENDIF 
    247187      IF (nn_verbose_level >= 0) WRITE(numout,'(2a)') 'icebergs, write_restart: creating ',TRIM(cl_path)//TRIM(cl_filename) 
  • branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/ICB/icbthm.F90

    r6486 r6491  
    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_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/IOM/in_out_manager.F90

    r6486 r6491  
    3030   CHARACTER(lc) ::   cn_ocerst_outdir !: restart output directory 
    3131   LOGICAL       ::   ln_rstart        !: start from (F) rest or (T) a restart file 
     32   LOGICAL       ::   ln_rstdate       !: datestamping of restarts 
    3233   LOGICAL       ::   ln_rst_list      !: output restarts at list of times (T) or by frequency (F) 
    3334   INTEGER       ::   nn_no            !: job number 
  • branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/IOM/iom.F90

    r6487 r6491  
    710710      CHARACTER(LEN=256)             ::   clname      ! file name 
    711711      CHARACTER(LEN=1)               ::   clrankpv, cldmspc      !  
     712      LOGICAL                        ::   ll_depth_spec ! T => if kstart, kcount present then *only* use values for 3rd spatial dimension. 
    712713      !--------------------------------------------------------------------- 
    713714      ! 
     
    722723      IF( PRESENT(kcount) .AND. (.NOT. PRESENT(kstart)) ) CALL ctl_stop(trim(clinfo), 'kcount present needs kstart present') 
    723724      IF( PRESENT(kstart) .AND. (.NOT. PRESENT(kcount)) ) CALL ctl_stop(trim(clinfo), 'kstart present needs kcount present') 
    724       IF( PRESENT(kstart) .AND. idom /= jpdom_unknown   ) CALL ctl_stop(trim(clinfo), 'kstart present needs kdom = jpdom_unknown') 
     725      IF( PRESENT(kstart) .AND. idom /= jpdom_unknown .AND.  idom /= jpdom_autoglo_xy  ) & 
     726     &           CALL ctl_stop(trim(clinfo), 'kstart present needs kdom = jpdom_unknown or kdom = jpdom_autoglo_xy') 
    725727 
    726728      luse_jattr = .false. 
     
    755757         ! update idom definition... 
    756758         ! Identify the domain in case of jpdom_auto(glo/dta) definition 
     759         IF( idom == jpdom_autoglo_xy ) THEN 
     760            ll_depth_spec = .TRUE. 
     761            idom = jpdom_autoglo 
     762         ELSE 
     763            ll_depth_spec = .FALSE. 
     764         ENDIF 
    757765         IF( idom == jpdom_autoglo .OR. idom == jpdom_autodta ) THEN             
    758766            IF( idom == jpdom_autoglo ) THEN   ;   idom = jpdom_global  
     
    808816         istart(idmspc+1) = itime 
    809817 
    810          IF(              PRESENT(kstart)      ) THEN ; istart(1:idmspc) = kstart(1:idmspc) ; icnt(1:idmspc) = kcount(1:idmspc) 
     818         IF( PRESENT(kstart) .AND. .NOT. ll_depth_spec ) THEN ; istart(1:idmspc) = kstart(1:idmspc) ; icnt(1:idmspc) = kcount(1:idmspc) 
    811819         ELSE 
    812             IF(           idom == jpdom_unknown ) THEN                                       ; icnt(1:idmspc) = idimsz(1:idmspc) 
     820            IF(           idom == jpdom_unknown ) THEN                                                ; icnt(1:idmspc) = idimsz(1:idmspc) 
    813821            ELSE  
    814822               IF( .NOT. PRESENT(pv_r1d) ) THEN   !   not a 1D array 
     
    833841                  ENDIF 
    834842                  IF( PRESENT(pv_r3d) ) THEN 
    835                      IF( idom == jpdom_data ) THEN   ; icnt(3) = jpkdta 
    836                      ELSE                            ; icnt(3) = jpk 
     843                     IF( idom == jpdom_data ) THEN                                  ; icnt(3) = jpkdta 
     844                     ELSE IF( ll_depth_spec .AND. PRESENT(kstart) ) THEN            ; istart(3) = kstart(3); icnt(3) = kcount(3) 
     845                     ELSE                                                           ; icnt(3) = jpk 
    837846                     ENDIF 
    838847                  ENDIF 
  • branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/IOM/iom_def.F90

    r6486 r6491  
    2626   INTEGER, PARAMETER, PUBLIC ::   jpdom_unknown       = 7   !: No dimension checking 
    2727   INTEGER, PARAMETER, PUBLIC ::   jpdom_autoglo       = 8   !:  
    28    INTEGER, PARAMETER, PUBLIC ::   jpdom_autodta       = 9   !:  
     28   INTEGER, PARAMETER, PUBLIC ::   jpdom_autoglo_xy    = 9   !: Automatically set horizontal dimensions only 
     29   INTEGER, PARAMETER, PUBLIC ::   jpdom_autodta       = 10  !:  
    2930 
    3031   INTEGER, PARAMETER, PUBLIC ::   jpioipsl    = 100      !: Use ioipsl (fliocom only) library 
     
    5758      INTEGER                                   ::   nvars    !: number of identified varibles in the file 
    5859      INTEGER                                   ::   iduld    !: id of the unlimited dimension 
     60      INTEGER                                   ::   lenuld   !: length of the unlimited dimension (number of records in file) 
    5961      INTEGER                                   ::   irec     !: writing record position   
    6062      CHARACTER(LEN=32)                         ::   uldname  !: name of the unlimited dimension 
  • branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/IOM/iom_nf90.F90

    r6486 r6491  
    154154         CALL iom_nf90_check(NF90_Inquire(if90id, unlimitedDimId = iom_file(kiomid)%iduld), clinfo) 
    155155         IF ( iom_file(kiomid)%iduld .GE. 0 ) THEN 
    156            CALL iom_nf90_check(NF90_Inquire_Dimension(if90id, iom_file(kiomid)%iduld,   & 
    157         &                                               name = iom_file(kiomid)%uldname), clinfo) 
     156           CALL iom_nf90_check(NF90_Inquire_Dimension(if90id, iom_file(kiomid)%iduld,     &  
     157        &                                               name = iom_file(kiomid)%uldname,  & 
     158        &                                               len  = iom_file(kiomid)%lenuld ), clinfo ) 
    158159         ENDIF 
    159160         IF(lwp) WRITE(numout,*) '                   ---> '//TRIM(cdname)//' OK' 
  • branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/IOM/restart.F90

    r6488 r6491  
    2121   USE in_out_manager  ! I/O manager 
    2222   USE iom             ! I/O module 
     23   USE ioipsl, ONLY : ju2ymds    ! for calendar 
    2324   USE eosbn2          ! equation of state            (eos bn2 routine) 
    2425   USE trdmxl_oce      ! ocean active mixed layer tracers trends variables 
     
    5556      !!---------------------------------------------------------------------- 
    5657      INTEGER, INTENT(in) ::   kt     ! ocean time-step 
     58      INTEGER             ::   iyear, imonth, iday 
     59      REAL (wp)           ::   zsec 
    5760      !! 
    5861      CHARACTER(LEN=20)   ::   clkt     ! ocean time-step deine as a character 
    5962      CHARACTER(LEN=50)   ::   clname   ! ocean output restart file name 
    60       CHARACTER(lc)       ::   clpath   ! full path to ocean output restart file 
     63      CHARACTER(LEN=150)  ::   clpath   ! full path to ocean output restart file 
    6164      !!---------------------------------------------------------------------- 
    6265      ! 
     
    8285      IF( kt == nitrst - 1 .OR. nstock == 1 .OR. ( kt == nitend .AND. .NOT. lrst_oce ) ) THEN 
    8386         IF( nitrst <= nitend .AND. nitrst > 0 ) THEN  
    84             ! beware of the format used to write kt (default is i8.8, that should be large enough...) 
    85             IF( nitrst > 999999999 ) THEN   ;   WRITE(clkt, *       ) nitrst 
    86             ELSE                            ;   WRITE(clkt, '(i8.8)') nitrst 
     87            IF ( ln_rstdate ) THEN 
     88               CALL ju2ymds( fjulday + rdttra(1) / rday, iyear, imonth, iday, zsec )            
     89               WRITE(clkt, '(i4.4,2i2.2)') iyear, imonth, iday 
     90            ELSE 
     91               ! beware of the format used to write kt (default is i8.8, that should be large enough...) 
     92               IF( nitrst > 999999999 ) THEN   ;   WRITE(clkt, *       ) nitrst 
     93               ELSE                            ;   WRITE(clkt, '(i8.8)') nitrst 
     94               ENDIF 
    8795            ENDIF 
    8896            ! create the file 
  • branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfmxl.F90

    r6487 r6491  
    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 
     
    3334   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hmlp    !: mixed layer depth  (rho=rho0+zdcrit) [m] 
    3435   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hmlpt   !: mixed layer depth at t-points        [m] 
     36   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hmld_zint  !: vertically-interpolated mixed layer depth   [m]  
     37   LOGICAL, PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)    :: ll_found   ! Is T_b to be found by interpolation ?  
     38   LOGICAL, PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  :: ll_belowml ! Flag points below mixed layer when ll_found=F 
    3539 
    3640   REAL(wp), PUBLIC ::   rho_c = 0.01_wp    !: density criterion for mixed layer depth 
    3741   REAL(wp)         ::   avt_c = 5.e-4_wp   ! Kz criterion for the turbocline depth 
     42 
     43   ! Namelist variables for  namzdf_mldzint 
     44   INTEGER          :: nn_mld_type         ! mixed layer type             
     45   REAL(wp)         :: rn_zref            ! depth of initial T_ref 
     46   REAL(wp)         :: rn_dT_crit          ! Critical temp diff  
     47   REAL(wp)         :: rn_iso_frac         ! Fraction of rn_dT_crit used  
    3848 
    3949   !! * Substitutions 
     
    5262      zdf_mxl_alloc = 0      ! set to zero if no array to be allocated 
    5363      IF( .NOT. ALLOCATED( nmln ) ) THEN 
    54          ALLOCATE( nmln(jpi,jpj), hmld(jpi,jpj), hmlp(jpi,jpj), hmlpt(jpi,jpj), STAT= zdf_mxl_alloc ) 
     64         ALLOCATE( nmln(jpi,jpj), hmld(jpi,jpj), hmlp(jpi,jpj), hmlpt(jpi,jpj), hmld_zint(jpi,jpj),       & 
     65        &          ll_found(jpi,jpj), ll_belowml(jpi,jpj,jpk), STAT= zdf_mxl_alloc ) 
    5566         ! 
    5667         IF( lk_mpp             )   CALL mpp_sum ( zdf_mxl_alloc ) 
     
    90101      CALL wrk_alloc( jpi,jpj, imld ) 
    91102 
    92       IF( kt == nit000 ) THEN 
     103      IF( kt <= nit000 ) THEN 
    93104         IF(lwp) WRITE(numout,*) 
    94105         IF(lwp) WRITE(numout,*) 'zdf_mxl : mixed layer depth' 
     
    134145      END DO 
    135146      IF( .NOT.lk_offline ) THEN            ! no need to output in offline mode 
     147      IF( kt >= nit000 ) THEN               ! workaround for calls before SOMETHING reads the XIOS namelist 
    136148         CALL iom_put( "mldr10_1", hmlp )   ! mixed layer depth 
    137149         CALL iom_put( "mldkz5"  , hmld )   ! turbocline depth 
    138150      ENDIF 
     151      ENDIF 
    139152       
     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 
    140159      IF(ln_ctl)   CALL prt_ctl( tab2d_1=REAL(nmln,wp), clinfo1=' nmln : ', tab2d_2=hmlp, clinfo2=' hmlp : ', ovlap=1 ) 
    141160      ! 
     
    145164      ! 
    146165   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) ) / fse3w(:,:,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 ( fsdept(:,:,jk) > rn_zref )  
     264           ik_ref(:,:) = jk - 1  
     265           zT_ref(:,:) = zT(:,:,jk-1) + zdTdz(:,:,jk-1) * ( rn_zref - fsdept(:,:,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 ( fsdept(:,:,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) = fsdept(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) = fsdept(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) = fsdept(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 
    147351 
    148352   !!====================================================================== 
  • branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90

    r6487 r6491  
    8383   INTEGER  ::   nn_htau   ! type of tke profile of penetration (=0/1) 
    8484   REAL(wp) ::   rn_efr    ! fraction of TKE surface value which penetrates in the ocean 
     85   REAL(wp) ::   rn_c      ! fraction of TKE added within the mixed layer by nn_etau 
    8586   LOGICAL  ::   ln_lc     ! Langmuir cells (LC) as a source term of TKE or not 
    8687   REAL(wp) ::   rn_lc     ! coef to compute vertical velocity of Langmuir cells 
     
    8889   REAL(wp) ::   ri_cri    ! critic Richardson number (deduced from rn_ediff and rn_ediss values) 
    8990   REAL(wp) ::   rmxl_min  ! minimum mixing length value (deduced from rn_ediff and rn_emin values)  [m] 
     91   REAL(wp) ::   rhtau                     ! coefficient to relate MLD to htau when nn_htau == 2 
    9092   REAL(wp) ::   rhftau_add = 1.e-3_wp     ! add offset   applied to HF part of taum  (nn_etau=3) 
    9193   REAL(wp) ::   rhftau_scl = 1.0_wp       ! scale factor applied to HF part of taum  (nn_etau=3) 
    9294 
    9395   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   htau           ! depth of tke penetration (nn_htau) 
     96   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   e_niw          !: TKE budget- near-inertial waves term 
     97   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   efr            ! surface boundary condition for nn_etau = 4 
    9498   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   dissl          ! now mixing lenght of dissipation 
    9599#if defined key_c1d 
     
    114118      !!---------------------------------------------------------------------- 
    115119      ALLOCATE(                                                                    & 
     120         &      efr  (jpi,jpj)     , e_niw(jpi,jpj,jpk) ,                         &       
    116121#if defined key_c1d 
    117122         &      e_dis(jpi,jpj,jpk) , e_mix(jpi,jpj,jpk) ,                          & 
     
    420425      END DO 
    421426 
     427      !                                 ! Save TKE prior to nn_etau addition   
     428      e_niw(:,:,:) = en(:,:,:)   
     429      !   
    422430      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    423431      !                            !  TKE due to surface and internal wave breaking 
    424432      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
     433      IF( nn_htau == 2 ) THEN           !* mixed-layer depth dependant length scale 
     434         DO jj = 2, jpjm1 
     435            DO ji = fs_2, fs_jpim1   ! vector opt. 
     436               htau(ji,jj) = rhtau * hmlp(ji,jj) 
     437            END DO 
     438         END DO 
     439      ENDIF 
     440#if defined key_iomput 
     441      ! 
     442      CALL iom_put( "htau", htau(:,:) )  ! Check htau (even if constant in time) 
     443#endif 
     444      ! 
    425445      IF( nn_etau == 1 ) THEN           !* penetration below the mixed layer (rn_efr fraction) 
    426446         DO jk = 2, jpkm1 
     
    457477            END DO 
    458478         END DO 
     479      ELSEIF( nn_etau == 4 ) THEN       !* column integral independant of htau (rn_efr must be scaled up) 
     480         IF( nn_htau == 2 ) THEN        ! efr dependant on time-varying htau  
     481            DO jj = 2, jpjm1 
     482               DO ji = fs_2, fs_jpim1   ! vector opt. 
     483                  efr(ji,jj) = rn_efr / ( htau(ji,jj) * ( 1._wp - EXP( -bathy(ji,jj) / htau(ji,jj) ) ) ) 
     484               END DO 
     485            END DO 
     486         ENDIF 
     487         DO jk = 2, jpkm1 
     488            DO jj = 2, jpjm1 
     489               DO ji = fs_2, fs_jpim1   ! vector opt. 
     490                  en(ji,jj,jk) = en(ji,jj,jk) + efr(ji,jj) * en(ji,jj,1) * EXP( -fsdepw(ji,jj,jk) / htau(ji,jj) )   & 
     491                     &                                                   * ( 1._wp - fr_i(ji,jj) )  * tmask(ji,jj,jk) 
     492               END DO 
     493            END DO 
     494         END DO 
    459495      ENDIF 
    460496      CALL lbc_lnk( en, 'W', 1. )      ! Lateral boundary conditions (sign unchanged) 
     497      ! 
     498      DO jk = 2, jpkm1                             ! TKE budget: near-inertial waves term   
     499         DO jj = 2, jpjm1   
     500            DO ji = fs_2, fs_jpim1   ! vector opt.   
     501               e_niw(ji,jj,jk) = en(ji,jj,jk) - e_niw(ji,jj,jk)   
     502            END DO   
     503         END DO   
     504      END DO   
     505      !   
     506      CALL lbc_lnk( e_niw, 'W', 1. )   
    461507      ! 
    462508      CALL wrk_dealloc( jpi,jpj, imlc )    ! integer 
     
    722768         &                 rn_emin0, rn_bshear, nn_mxl , ln_mxl0  ,   & 
    723769         &                 rn_mxl0 , nn_pdl   , ln_lc  , rn_lc    ,   & 
    724          &                 nn_etau , nn_htau  , rn_efr    
    725       !!---------------------------------------------------------------------- 
    726       ! 
     770         &                 nn_etau , nn_htau  , rn_efr , rn_c    
     771      !!---------------------------------------------------------------------- 
     772 
    727773      REWIND( numnam_ref )              ! Namelist namzdf_tke in reference namelist : Turbulent Kinetic Energy 
    728774      READ  ( numnam_ref, namzdf_tke, IOSTAT = ios, ERR = 901) 
     
    757803         WRITE(numout,*) '      flag for computation of exp. tke profile    nn_htau   = ', nn_htau 
    758804         WRITE(numout,*) '      fraction of en which pene. the thermocline  rn_efr    = ', rn_efr 
     805         WRITE(numout,*) '      fraction of TKE added within the mixed layer by nn_etau rn_c    = ', rn_c 
    759806         WRITE(numout,*) 
    760807         WRITE(numout,*) '      critical Richardson nb with your parameters  ri_cri = ', ri_cri 
     
    767814      IF( nn_mxl  < 0   .OR.  nn_mxl  > 3 )   CALL ctl_stop( 'bad flag: nn_mxl is  0, 1 or 2 ' ) 
    768815      IF( nn_pdl  < 0   .OR.  nn_pdl  > 1 )   CALL ctl_stop( 'bad flag: nn_pdl is  0 or 1    ' ) 
    769       IF( nn_htau < 0   .OR.  nn_htau > 1 )   CALL ctl_stop( 'bad flag: nn_htau is 0, 1 or 2 ' ) 
     816      IF( nn_htau < 0  .OR.  nn_htau > 5 )   CALL ctl_stop( 'bad flag: nn_htau is 0 to 5    ' ) 
    770817      IF( nn_etau == 3 .AND. .NOT. ln_cpl )   CALL ctl_stop( 'nn_etau == 3 : HF taum only known in coupled mode' ) 
    771818 
     
    780827      ENDIF 
    781828 
     829      IF( nn_etau /= 0 .and. nn_htau == 2 ) THEN 
     830          ierr = zdf_mxl_alloc() 
     831          nmln(:,:) = nlb10           ! Initialization of nmln 
     832      ENDIF 
     833 
    782834      !                               !* depth of penetration of surface tke 
    783835      IF( nn_etau /= 0 ) THEN       
     836         htau(:,:) = 0._wp 
    784837         SELECT CASE( nn_htau )             ! Choice of the depth of penetration 
    785838         CASE( 0 )                                 ! constant depth penetration (here 10 meters) 
     
    787840         CASE( 1 )                                 ! F(latitude) : 0.5m to 30m poleward of 40 degrees 
    788841            htau(:,:) = MAX(  0.5_wp, MIN( 30._wp, 45._wp* ABS( SIN( rpi/180._wp * gphit(:,:) ) ) )   )             
     842         CASE( 2 )                                 ! fraction of depth-integrated TKE within mixed-layer 
     843            rhtau = -1._wp / LOG( 1._wp - rn_c ) 
     844         CASE( 3 )                                 ! F(latitude) : 0.5m to 15m poleward of 20 degrees 
     845            htau(:,:) = MAX(  0.5_wp, MIN( 15._wp, 45._wp* ABS( SIN( rpi/180._wp * gphit(:,:) ) ) )   ) 
     846         CASE( 4 )                                 ! F(latitude) : 0.5m to 10m/30m poleward of 13/40 degrees north/south 
     847            DO jj = 2, jpjm1 
     848               DO ji = fs_2, fs_jpim1   ! vector opt. 
     849                  IF( gphit(ji,jj) <= 0._wp ) THEN 
     850                     htau(ji,jj) = MAX(  0.5_wp, MIN( 30._wp, 45._wp* ABS( SIN( rpi/180._wp * gphit(ji,jj) ) ) )   ) 
     851                  ELSE 
     852                     htau(ji,jj) = MAX(  0.5_wp, MIN( 10._wp, 45._wp* ABS( SIN( rpi/180._wp * gphit(ji,jj) ) ) )   ) 
     853                  ENDIF 
     854               END DO 
     855            END DO 
     856         CASE ( 5 )                                ! F(latitude) : 0.5m to 10m poleward of 13 degrees north/south, 
     857            DO jj = 2, jpjm1                       !               10m to 30m between 30/45 degrees south 
     858               DO ji = fs_2, fs_jpim1   ! vector opt. 
     859                  IF( gphit(ji,jj) <= -30._wp ) THEN 
     860                     htau(ji,jj) = MAX(  10._wp, MIN( 30._wp, 55._wp* ABS( SIN( rpi/120._wp * ( gphit(ji,jj) + 23._wp ) ) ) )   ) 
     861                  ELSE 
     862                     htau(ji,jj) = MAX(  0.5_wp, MIN( 10._wp, 45._wp* ABS( SIN( rpi/180._wp * gphit(ji,jj) ) ) )   ) 
     863                  ENDIF 
     864               END DO 
     865            END DO 
    789866         END SELECT 
     867         ! 
     868         IF( nn_etau == 4 .AND. nn_htau /= 2 ) THEN            ! efr dependant on constant htau 
     869            DO jj = 2, jpjm1 
     870               DO ji = fs_2, fs_jpim1   ! vector opt. 
     871                  efr(ji,jj) = rn_efr / ( htau(ji,jj) * ( 1._wp - EXP( -bathy(ji,jj) / htau(ji,jj) ) ) ) 
     872               END DO 
     873            END DO 
     874         ENDIF 
    790875      ENDIF 
    791876      !                               !* set vertical eddy coef. to the background value 
  • branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/step.F90

    r6487 r6491  
    234234      IF( lk_diaar5  )      CALL dia_ar5( kstp )         ! ar5 diag 
    235235      IF( lk_diaharm )      CALL dia_harm( kstp )        ! Tidal harmonic analysis 
     236                            CALL dia_prod( kstp )        ! ocean model: product diagnostics 
    236237                            CALL dia_wri( kstp )         ! ocean model: outputs 
    237238      ! 
  • branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/step_oce.F90

    r6487 r6491  
    8888 
    8989   USE diawri           ! Standard run outputs             (dia_wri routine) 
     90   USE diaprod          ! Product diagnostics              (dia_prod routine) 
    9091   USE diaptr           ! poleward transports              (dia_ptr routine) 
    9192   USE diadct           ! sections transports              (dia_dct routine) 
Note: See TracChangeset for help on using the changeset viewer.