Changeset 6507
- Timestamp:
- 2016-05-03T14:28:12+02:00 (8 years ago)
- Location:
- branches/UKMO/dev_r6501_GO6_package_trunk/NEMOGCM/NEMO/OPA_SRC
- Files:
-
- 20 edited
-
DIA/diawri.F90 (modified) (3 diffs)
-
DOM/closea.F90 (modified) (1 diff)
-
DOM/domain.F90 (modified) (2 diffs)
-
DOM/dommsk.F90 (modified) (2 diffs)
-
DOM/phycst.F90 (modified) (1 diff)
-
ICB/icbclv.F90 (modified) (3 diffs)
-
ICB/icbdia.F90 (modified) (1 diff)
-
ICB/icbrst.F90 (modified) (2 diffs)
-
ICB/icbthm.F90 (modified) (3 diffs)
-
IOM/in_out_manager.F90 (modified) (1 diff)
-
IOM/restart.F90 (modified) (4 diffs)
-
SBC/sbc_ice.F90 (modified) (2 diffs)
-
SBC/sbc_oce.F90 (modified) (3 diffs)
-
SBC/sbccpl.F90 (modified) (24 diffs)
-
SBC/sbcice_cice.F90 (modified) (13 diffs)
-
SBC/sbcisf.F90 (modified) (3 diffs)
-
SBC/sbcmod.F90 (modified) (1 diff)
-
ZDF/zdfmxl.F90 (modified) (5 diffs)
-
ZDF/zdftke.F90 (modified) (10 diffs)
-
step_oce.F90 (modified) (1 diff)
Legend:
- Unmodified
- Added
- Removed
-
branches/UKMO/dev_r6501_GO6_package_trunk/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90
r6503 r6507 39 39 USE zdfmxl ! mixed layer 40 40 USE dianam ! build name of file (routine) 41 USE zdftke ! vertical physics: one-equation scheme 41 42 USE zdfddm ! vertical physics: double diffusion 42 43 USE diahth ! thermocline diagnostics … … 231 232 CALL iom_put( "avt" , avt ) ! T vert. eddy diff. coef. 232 233 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 233 238 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(:,:,:) ) ) ) 234 242 235 243 IF( iom_use('logavt') ) CALL iom_put( "logavt", LOG( MAX( 1.e-20_wp, avt (:,:,:) ) ) ) … … 974 982 CALL histwrite( id_i, "sozotaux", kt, utau , jpi*jpj , idex ) ! i-wind stress 975 983 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 976 988 977 989 IF( .NOT.ln_linssh ) THEN -
branches/UKMO/dev_r6501_GO6_package_trunk/NEMOGCM/NEMO/OPA_SRC/DOM/closea.F90
r6503 r6507 155 155 ncsi2(4) = 76 ; ncsj2(4) = 61 156 156 ncsir(4,1) = 84 ; ncsjr(4,1) = 59 157 ! ! ======================= 158 CASE ( 025 ) ! ORCA_R025 configuration159 ! ! ======================= 160 ncsnr(1) = 1 ; ncstt(1) = 0 ! Caspian + Aralsea161 ncsi1(1) = 1330 ; ncsj1(1) = 645162 ncsi2(1) = 1 400 ; ncsj2(1) = 795157 ! ! ================================ 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 163 163 ncsir(1,1) = 1 ; ncsjr(1,1) = 1 164 164 ! 165 ncsnr(2) = 1 ; ncstt(2) = 0 ! A zov Sea166 ncsi1(2) = 1 284 ; ncsj1(2) = 722167 ncsi2(2) = 1 304 ; ncsj2(2) = 747165 ncsnr(2) = 1 ; ncstt(2) = 0 ! Aral sea 166 ncsi1(2) = 1376 ; ncsj1(2) = 900 167 ncsi2(2) = 1400 ; ncsj2(2) = 981 168 168 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 169 209 ! 170 210 END SELECT -
branches/UKMO/dev_r6501_GO6_package_trunk/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90
r6503 r6507 161 161 USE ioipsl 162 162 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, & 164 164 & nn_it000, nn_itend , nn_date0 , nn_time0 , nn_leapy , nn_istate , & 165 165 & nn_stock, nn_write , ln_mskland , ln_clobber , nn_chunksz, nn_euler , & … … 196 196 WRITE(numout,*) ' file prefix restart output cn_ocerst_out= ', cn_ocerst_out 197 197 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 199 200 WRITE(numout,*) ' start with forward time step nn_euler = ', nn_euler 200 201 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 105 105 INTEGER :: isrow ! index for ORCA1 starting row 106 106 INTEGER , POINTER, DIMENSION(:,:) :: imsk 107 REAL(wp) :: zphi_drake_passage, zshlat_antarc 107 108 REAL(wp), POINTER, DIMENSION(:,:) :: zwf 108 109 !! … … 345 346 ENDIF 346 347 ! 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 ! 347 371 CALL lbc_lnk( fmask, 'F', 1._wp ) ! Lateral boundary conditions on fmask 348 372 ! -
branches/UKMO/dev_r6501_GO6_package_trunk/NEMOGCM/NEMO/OPA_SRC/DOM/phycst.F90
r6503 r6507 68 68 REAL(wp), PUBLIC :: rcdsn = 0.31_wp !: thermal conductivity of snow 69 69 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 70 73 REAL(wp), PUBLIC :: lsub = 2.834e+6_wp !: pure ice latent heat of sublimation [J/kg] 74 #endif 71 75 REAL(wp), PUBLIC :: lfus = 0.334e+6_wp !: latent heat of fusion of fresh ice [J/kg] 72 76 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 25 25 USE icbutl ! iceberg utility routines 26 26 27 USE sbc_oce ! for icesheet freshwater input variables 28 USE in_out_manager 29 USE iom 30 27 31 IMPLICIT NONE 28 32 PRIVATE … … 48 52 ! 49 53 REAL(wp) :: zcalving_used, zdist, zfact 54 REAL(wp) :: zgreenland_calving_sum, zantarctica_calving_sum 50 55 INTEGER :: jn, ji, jj ! loop counters 51 56 INTEGER :: imx ! temporary integer for max berg class … … 59 64 zfact = ( (1000._wp)**3 / ( NINT(rday) * nyear_len(1) ) ) * 850._wp 60 65 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(:,:) ) 61 103 62 104 ! 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 371 371 IF( .NOT. ln_bergdia ) RETURN !!gm useless iom will control whether it is output or not 372 372 ! 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] 373 375 CALL iom_put( "berg_melt" , berg_melt (:,:) ) ! Melt rate of icebergs [kg/m2/s] 374 376 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 158 158 INTEGER :: jn ! dummy loop index 159 159 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 162 165 TYPE(iceberg), POINTER :: this 163 166 TYPE(point) , POINTER :: pt … … 168 171 IF( cl_path(LEN_TRIM(cl_path):) /= '/' ) cl_path = TRIM(cl_path) // '/' 169 172 IF( lk_mpp ) THEN 170 WRITE(cl_filename,'(A,"_icebergs_", I8.8,"_restart_",I4.4,".nc")') TRIM(cexper), kt, narea-1173 WRITE(cl_filename,'(A,"_icebergs_",A,"_restart_",I4.4,".nc")') TRIM(cexper), TRIM(ADJUSTL(clkt)), narea-1 171 174 ELSE 172 WRITE(cl_filename,'(A,"_icebergs_", I8.8,"_restart.nc")') TRIM(cexper), kt175 WRITE(cl_filename,'(A,"_icebergs_",A,"_restart.nc")') TRIM(cexper), TRIM(ADJUSTL(clkt)) 173 176 ENDIF 174 177 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 18 18 USE dom_oce ! NEMO domain 19 19 USE in_out_manager ! NEMO IO routines, numout in particular 20 USE iom 20 21 USE lib_mpp ! NEMO MPI routines, ctl_stop in particular 21 22 USE phycst ! NEMO physical constants … … 160 161 zmelt = ( zdM - ( zdMbitsE - zdMbitsM ) ) * z1_dt ! kg/s 161 162 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 163 165 berg_grid%calving_hflx (ii,ij) = berg_grid%calving_hflx (ii,ij) + zheat * z1_e1e2 ! W/m2 164 166 CALL icb_dia_melt( ii, ij, zMnew, zheat, this%mass_scaling, & … … 208 210 IF(.NOT. ln_passive_mode ) THEN 209 211 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 (:,:) 211 213 ENDIF 212 214 ! -
branches/UKMO/dev_r6501_GO6_package_trunk/NEMOGCM/NEMO/OPA_SRC/IOM/in_out_manager.F90
r6503 r6507 29 29 CHARACTER(lc) :: cn_ocerst_outdir !: restart output directory 30 30 LOGICAL :: ln_rstart !: start from (F) rest or (T) a restart file 31 LOGICAL :: ln_rstdate !: datestamping of restarts 31 32 LOGICAL :: ln_rst_list !: output restarts at list of times (T) or by frequency (F) 32 33 INTEGER :: nn_no !: job number -
branches/UKMO/dev_r6501_GO6_package_trunk/NEMOGCM/NEMO/OPA_SRC/IOM/restart.F90
r6503 r6507 57 57 !!---------------------------------------------------------------------- 58 58 INTEGER, INTENT(in) :: kt ! ocean time-step 59 INTEGER :: iyear, imonth, iday 60 REAL (wp) :: zsec 59 61 !! 60 62 CHARACTER(LEN=20) :: clkt ! ocean time-step deine as a character … … 84 86 IF( kt == nitrst - 1 .OR. nstock == 1 .OR. ( kt == nitend .AND. .NOT. lrst_oce ) ) THEN 85 87 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 89 96 ENDIF 90 97 ! create the file … … 155 162 IF (ln_diurnal) CALL iom_rstput( kt, nitrst, numrow, 'Dsst', x_dsst ) 156 163 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 157 176 IF( kt == nitrst ) THEN 158 177 CALL iom_close( numrow ) ! close the restart file (only at last time step) … … 252 271 ENDIF 253 272 ! 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 254 292 IF( neuler == 0 ) THEN ! Euler restart (neuler=0) 255 293 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 102 102 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: fr_iu !: ice fraction at NEMO U point 103 103 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 105 108 ! variables used in the coupled interface 106 109 INTEGER , PUBLIC, PARAMETER :: jpl = ncat … … 153 156 154 157 #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) , & 156 159 wndi_ice(jpi,jpj) , tatm_ice(jpi,jpj) , qatm_ice(jpi,jpj) , & 157 160 wndj_ice(jpi,jpj) , nfrzmlt(jpi,jpj) , ss_iou(jpi,jpj) , & 158 161 ss_iov(jpi,jpj) , fr_iu(jpi,jpj) , fr_iv(jpi,jpj) , & 159 162 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) , & 162 167 & v_ice(jpi,jpj) , fr2_i0(jpi,jpj) , alb_ice(jpi,jpj,1) , & 163 168 & 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) ) 165 171 166 172 #endif -
branches/UKMO/dev_r6501_GO6_package_trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_oce.F90
r6503 r6507 124 124 #endif 125 125 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 126 128 127 129 !!---------------------------------------------------------------------- … … 136 138 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: e3t_m !: mean (nn_fsbc time-step) sea surface layer thickness [m] 137 139 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. 138 159 139 160 !! * Substitutions … … 171 192 & ssu_m (jpi,jpj) , sst_m(jpi,jpj) , frq_m(jpi,jpj) , & 172 193 & 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) ) 173 196 ! 174 197 ALLOCATE( e3t_m(jpi,jpj) , STAT=ierr(5) ) -
branches/UKMO/dev_r6501_GO6_package_trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90
r6503 r6507 106 106 INTEGER, PARAMETER :: jpr_e3t1st = 41 ! first T level thickness 107 107 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 109 112 110 113 INTEGER, PARAMETER :: jps_fice = 1 ! ice fraction sent to the atmosphere … … 136 139 INTEGER, PARAMETER :: jps_e3t1st = 27 ! first level depth (vvl) 137 140 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 139 147 140 148 ! !!** namelist namsbc_cpl ** … … 147 155 END TYPE FLD_C 148 156 ! ! 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 150 159 ! ! Received from the atmosphere 151 160 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 153 162 ! ! Other namelist parameters 154 163 INTEGER :: nn_cplmodel ! Maximum number of models to/from which NEMO is potentialy sending/receiving data … … 214 223 REAL(wp), POINTER, DIMENSION(:,:) :: zacs, zaos 215 224 !! 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 220 232 !!--------------------------------------------------------------------- 221 233 ! … … 256 268 WRITE(numout,*)' runoffs = ', TRIM(sn_rcv_rnf%cldes ), ' (', TRIM(sn_rcv_rnf%clcat ), ')' 257 269 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 ), ')' 258 272 WRITE(numout,*)' sea ice heat fluxes = ', TRIM(sn_rcv_iceflx%cldes), ' (', TRIM(sn_rcv_iceflx%clcat), ')' 259 273 WRITE(numout,*)' atm co2 = ', TRIM(sn_rcv_co2%cldes ), ' (', TRIM(sn_rcv_co2%clcat ), ')' … … 267 281 WRITE(numout,*)' - mesh = ', sn_snd_crt%clvgrd 268 282 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 269 287 WRITE(numout,*)' nn_cplmodel = ', nn_cplmodel 270 288 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 271 293 ENDIF 272 294 … … 381 403 srcv(jpr_snow)%clname = 'OTotSnow' ! Snow = solid precipitation 382 404 srcv(jpr_tevp)%clname = 'OTotEvap' ! total evaporation (over oce + ice sublimation) 383 srcv(jpr_ievp)%clname = 'OIceEv ap' ! evaporation over ice = sublimation405 srcv(jpr_ievp)%clname = 'OIceEvp' ! evaporation over ice = sublimation 384 406 srcv(jpr_sbpr)%clname = 'OSubMPre' ! sublimation - liquid precipitation - solid precipitation 385 407 srcv(jpr_semp)%clname = 'OISubMSn' ! ice solid water budget = sublimation - solid precipitation … … 395 417 END SELECT 396 418 ! 419 !Set the number of categories for coupling of sublimation 420 IF ( TRIM( sn_rcv_emp%clcat ) == 'yes' ) srcv(jpr_ievp)%nct = jpl 421 ! 397 422 ! ! ------------------------- ! 398 423 ! ! Runoffs & Calving ! … … 408 433 ! 409 434 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. 410 437 ! 411 438 ! ! ------------------------- ! … … 481 508 srcv(jpr_topm:jpr_botm)%laction = .TRUE. 482 509 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 483 523 ! ! ------------------------------- ! 484 524 ! ! OPA-SAS coupling - rcv by opa ! … … 598 638 ! ! ------------------------- ! 599 639 ssnd(jps_toce)%clname = 'O_SSTSST' 600 ssnd(jps_tice)%clname = 'O _TepIce'640 ssnd(jps_tice)%clname = 'OTepIce' 601 641 ssnd(jps_tmix)%clname = 'O_TepMix' 602 642 SELECT CASE( TRIM( sn_snd_temp%cldes ) ) 603 643 CASE( 'none' ) ! nothing to do 604 644 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') 606 646 ssnd( (/jps_toce, jps_tice/) )%laction = .TRUE. 607 647 IF ( TRIM( sn_snd_temp%clcat ) == 'yes' ) ssnd(jps_tice)%nct = jpl … … 637 677 ssnd(jps_hice)%clname = 'OIceTck' 638 678 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' 639 682 IF( k_ice /= 0 ) THEN 640 683 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 641 686 ! Currently no namelist entry to determine sending of multi-category ice fraction so use the thickness entry for now 642 687 IF ( TRIM( sn_snd_thick%clcat ) == 'yes' ) ssnd(jps_fice)%nct = jpl … … 655 700 CASE default ; CALL ctl_stop( 'sbc_cpl_init: wrong definition of sn_snd_thick%cldes' ) 656 701 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 657 739 658 740 ! ! ------------------------- ! … … 687 769 ! ! ------------------------- ! 688 770 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 689 803 690 804 ! ! ------------------------------- ! … … 783 897 ncpl_qsr_freq = 86400 / ncpl_qsr_freq 784 898 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 ) 786 923 ! 787 924 IF( nn_timing == 1 ) CALL timing_stop('sbc_cpl_init') … … 841 978 !! 842 979 LOGICAL :: llnewtx, llnewtau ! update wind stress components and module?? 843 INTEGER :: ji, jj, j n! dummy loop indices980 INTEGER :: ji, jj, jl, jn ! dummy loop indices 844 981 INTEGER :: isec ! number of seconds since nit000 (assuming rdt did not change since nit000) 845 982 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 846 986 REAL(wp) :: zcoef ! temporary scalar 847 987 REAL(wp) :: zrhoa = 1.22 ! Air density kg/m3 … … 894 1034 IF( srcv(jpr_otx2)%laction ) THEN 895 1035 CALL rot_rep( frcv(jpr_otx2)%z3(:,:,1), frcv(jpr_oty2)%z3(:,:,1), srcv(jpr_otx2)%clgrid, 'en->j', zty ) 896 ELSE 1036 ELSE 897 1037 CALL rot_rep( frcv(jpr_otx1)%z3(:,:,1), frcv(jpr_oty1)%z3(:,:,1), srcv(jpr_otx1)%clgrid, 'en->j', zty ) 898 1038 ENDIF … … 991 1131 #endif 992 1132 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 993 1152 ! Fields received by SAS when OASIS coupling 994 1153 ! (arrays no more filled at sbcssm stage) … … 1105 1264 ! 1106 1265 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 1107 1325 ! 1108 1326 CALL wrk_dealloc( jpi,jpj, ztx, zty, zmsk, zemp, zqns, zqsr ) … … 1400 1618 ztprecip(:,:) = frcv(jpr_rain)%z3(:,:,1) + zsprecip(:,:) ! May need to ensure positive here 1401 1619 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 1402 1646 zemp_ice(:,:) = frcv(jpr_ievp)%z3(:,:,1) - frcv(jpr_snow)%z3(:,:,1) 1647 #endif 1403 1648 CALL iom_put( 'rain' , frcv(jpr_rain)%z3(:,:,1) ) ! liquid precipitation 1404 1649 IF( iom_use('hflx_rain_cea') ) & … … 1795 2040 CASE default ; CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_temp%clcat' ) 1796 2041 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 1797 2053 CASE( 'mixed oce-ice' ) 1798 2054 ztmp1(:,:) = ( ztmp1(:,:) + rt0 ) * zfr_l(:,:) … … 1865 2121 END SELECT 1866 2122 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) 1867 2133 ENDIF 1868 2134 … … 1910 2176 IF( ssnd(jps_hsnw)%laction ) CALL cpl_snd( jps_hsnw, isec, ztmp4, info ) 1911 2177 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 ! 1912 2222 ! 1913 2223 #if defined key_cpl_carbon_cycle … … 2089 2399 IF( ssnd(jps_rnf )%laction ) CALL cpl_snd( jps_rnf , isec, RESHAPE ( rnf , (/jpi,jpj,1/) ), info ) 2090 2400 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 ! 2092 2405 CALL wrk_dealloc( jpi,jpj, zfr_l, ztmp1, ztmp2, zotx1, zoty1, zotz1, zitx1, zity1, zitz1 ) 2093 2406 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 13 13 USE dom_oce ! ocean space and time domain 14 14 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 16 17 USE in_out_manager ! I/O manager 17 18 USE iom, ONLY : iom_put,iom_use ! I/O manager library !!Joakim edit … … 35 36 USE ice_gather_scatter 36 37 USE ice_calendar, only: dt 38 # if defined key_cice4 37 39 USE ice_state, only: aice,aicen,uvel,vvel,vsno,vsnon,vice,vicen 38 # if defined key_cice439 40 USE ice_flux, only: strax,stray,strocnx,strocny,frain,fsnow, & 40 41 strocnxT,strocnyT, & … … 43 44 flatn_f,fsurfn_f,fcondtopn_f, & 44 45 uatm,vatm,wind,fsw,flw,Tair,potT,Qa,rhoa,zlvl, & 45 swvdr,swvdf,swidr,swidf 46 swvdr,swvdf,swidr,swidf,Tf 46 47 USE ice_therm_vertical, only: calc_Tsfc 47 48 #else 49 USE ice_state, only: aice,aicen,uvel,nt_hpnd,trcrn,vvel,vsno,& 50 vsnon,vice,vicen,nt_Tsfc 48 51 USE ice_flux, only: strax,stray,strocnx,strocny,frain,fsnow, & 49 52 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, & 52 55 flatn_f,fsurfn_f,fcondtopn_f, & 53 56 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 56 66 #endif 57 67 USE ice_forcing, only: frcvdr,frcvdf,frcidr,frcidf … … 162 172 REAL(wp), DIMENSION(:,:), POINTER :: ztmp1, ztmp2 163 173 REAL(wp) :: zcoefu, zcoefv, zcoeff ! local scalar 174 REAL(wp), DIMENSION(:,:,:), POINTER :: ztfrz3d 164 175 INTEGER :: ji, jj, jl, jk ! dummy loop indices 165 176 !!--------------------------------------------------------------------- … … 173 184 ji_off = INT ( (jpiglo - nx_global) / 2 ) 174 185 jj_off = INT ( (jpjglo - ny_global) / 2 ) 175 176 #if defined key_nemocice_decomp177 ! Pass initial SST from NEMO to CICE so ice is initialised correctly if178 ! there is no restart file.179 ! Values from a CICE restart file would overwrite this180 IF ( .NOT. ln_rstart ) THEN181 CALL nemo2cice( tsn(:,:,1,jp_tem) , sst , 'T' , 1.)182 ENDIF183 #endif184 186 185 187 ! Initialize CICE … … 199 201 200 202 ! 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' ) 202 204 IF( sbc_ice_cice_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'sbc_ice_cice_alloc : unable to allocate cice arrays' ) 203 205 204 ! Ensure ocean temperatures are nowhere below freezing if not a NEMO restart206 ! Ensure that no temperature points are below freezing if not a NEMO restart 205 207 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 ) 207 213 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. ) 212 228 213 229 CALL cice2nemo(aice,fr_i, 'T', 1. ) … … 220 236 ! T point to U point 221 237 ! T point to V point 238 fr_iu(:,:)=0.0 239 fr_iv(:,:)=0.0 222 240 DO jj=1,jpjm1 223 241 DO ji=1,jpim1 … … 345 363 CALL nemo2cice(ztmp,stray,'F', -1. ) 346 364 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 347 377 ! Surface downward latent heat flux (CI_5) 348 378 IF (ksbc == jp_flx) THEN … … 350 380 ztmpn(:,:,jl)=qla_ice(:,:,1)*a_i(:,:,jl) 351 381 ENDDO 382 ELSE IF (ksbc == jp_purecpl) THEN 383 DO jl=1,ncat 384 ztmpn(:,:,jl)=qla_ice(:,:,jl)*a_i(:,:,jl) 385 ENDDO 352 386 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) 371 389 ENDIF 390 372 391 DO jl=1,ncat 373 392 CALL nemo2cice(ztmpn(:,:,jl),flatn_f(:,:,jl,:),'T', 1. ) … … 375 394 ! GBM conductive flux through ice (CI_6) 376 395 ! Convert to GBM 377 IF (ksbc == jp_flx ) THEN396 IF (ksbc == jp_flx .OR. ksbc == jp_purecpl) THEN 378 397 ztmp(:,:) = botmelt(:,:,jl)*a_i(:,:,jl) 379 398 ELSE … … 384 403 ! GBM surface heat flux (CI_7) 385 404 ! Convert to GBM 386 IF (ksbc == jp_flx ) THEN405 IF (ksbc == jp_flx .OR. ksbc == jp_purecpl) THEN 387 406 ztmp(:,:) = (topmelt(:,:,jl)+botmelt(:,:,jl))*a_i(:,:,jl) 388 407 ELSE … … 456 475 CALL nemo2cice(sss_m,sss,'T', 1. ) 457 476 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 458 484 ! x comp and y comp of surface ocean current 459 485 ! U point to F point … … 730 756 CALL cice2nemo(vicen(:,:,jl,:),ht_i(:,:,jl),'T', 1. ) 731 757 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 732 786 ! 733 787 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 133 133 CALL fld_read ( kt, nn_fsbc, sf_rnfisf ) 134 134 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 135 173 qisf(:,:) = fwfisf(:,:) * rlfusisf ! heat flux 136 174 stbl(:,:) = soce … … 139 177 CALL fld_read ( kt, nn_fsbc, sf_fwfisf ) 140 178 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 141 217 qisf(:,:) = fwfisf(:,:) * rlfusisf ! heat flux 142 218 stbl(:,:) = soce … … 155 231 CALL eos_fzp( stbl(:,:), zt_frz(:,:), zdep(:,:) ) 156 232 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 157 236 risf_tsc(:,:,jp_tem) = qisf(:,:) * r1_rau0_rcp - fwfisf(:,:) * zt_frz(:,:) * r1_rau0 ! 158 237 risf_tsc(:,:,jp_sal) = 0.0_wp -
branches/UKMO/dev_r6501_GO6_package_trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90
r6503 r6507 189 189 fwfisf (:,:) = 0.0_wp ; fwfisf_b (:,:) = 0.0_wp 190 190 risf_tsc(:,:,:) = 0.0_wp ; risf_tsc_b(:,:,:) = 0.0_wp 191 rdivisf = 0.0_wp 191 192 END IF 192 193 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 18 18 USE phycst ! physical constants 19 19 USE iom ! I/O library 20 USE eosbn2 ! for zdf_mxl_zint 20 21 USE lib_mpp ! MPP library 21 22 USE wrk_nemo ! work arrays … … 32 33 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hmlp !: mixed layer depth (rho=rho0+zdcrit) [m] 33 34 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 34 38 35 39 REAL(wp), PUBLIC :: rho_c = 0.01_wp !: density criterion for mixed layer depth … … 49 53 zdf_mxl_alloc = 0 ! set to zero if no array to be allocated 50 54 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 ) 52 57 ! 53 58 IF( lk_mpp ) CALL mpp_sum ( zdf_mxl_alloc ) … … 145 150 END IF 146 151 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 ! 148 159 IF(ln_ctl) CALL prt_ctl( tab2d_1=REAL(nmln,wp), clinfo1=' nmln : ', tab2d_2=hmlp, clinfo2=' hmlp : ', ovlap=1 ) 149 160 ! … … 153 164 ! 154 165 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) 219 901 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 ) 223 902 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 155 352 156 353 !!====================================================================== -
branches/UKMO/dev_r6501_GO6_package_trunk/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90
r6503 r6507 81 81 INTEGER :: nn_htau ! type of tke profile of penetration (=0/1) 82 82 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 83 84 LOGICAL :: ln_lc ! Langmuir cells (LC) as a source term of TKE or not 84 85 REAL(wp) :: rn_lc ! coef to compute vertical velocity of Langmuir cells … … 86 87 REAL(wp) :: ri_cri ! critic Richardson number (deduced from rn_ediff and rn_ediss values) 87 88 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 88 90 REAL(wp) :: rhftau_add = 1.e-3_wp ! add offset applied to HF part of taum (nn_etau=3) 89 91 REAL(wp) :: rhftau_scl = 1.0_wp ! scale factor applied to HF part of taum (nn_etau=3) 90 92 91 93 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 92 95 REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: dissl ! now mixing lenght of dissipation 93 96 REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: apdlr ! now mixing lenght of dissipation … … 112 115 !!---------------------------------------------------------------------- 113 116 ALLOCATE( & 117 & efr (jpi,jpj) , & 114 118 #if defined key_c1d 115 119 & e_dis(jpi,jpj,jpk) , e_mix(jpi,jpj,jpk) , & … … 446 450 ! ! TKE due to surface and internal wave breaking 447 451 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 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 ! 448 464 !!gm BUG : in the exp remove the depth of ssh !!! 449 450 465 451 466 IF( nn_etau == 1 ) THEN !* penetration below the mixed layer (rn_efr fraction) … … 477 492 en(ji,jj,jk) = en(ji,jj,jk) + zbbrau * zdif * EXP( -gdepw_n(ji,jj,jk) / htau(ji,jj) ) & 478 493 & * ( 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) 479 510 END DO 480 511 END DO … … 723 754 !!---------------------------------------------------------------------- 724 755 INTEGER :: ji, jj, jk ! dummy loop indices 725 INTEGER :: ios 756 INTEGER :: ios, ierr 726 757 !! 727 758 NAMELIST/namzdf_tke/ rn_ediff, rn_ediss , rn_ebb , rn_emin , & 728 759 & rn_emin0, rn_bshear, nn_mxl , ln_mxl0 , & 729 760 & rn_mxl0 , nn_pdl , ln_lc , rn_lc , & 730 & nn_etau , nn_htau , rn_efr 761 & nn_etau , nn_htau , rn_efr , rn_c 731 762 !!---------------------------------------------------------------------- 732 763 ! … … 774 805 WRITE(numout,*) ' flag for computation of exp. tke profile nn_htau = ', nn_htau 775 806 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 776 808 WRITE(numout,*) 777 809 WRITE(numout,*) ' critical Richardson nb with your parameters ri_cri = ', ri_cri … … 784 816 IF( nn_mxl < 0 .OR. nn_mxl > 3 ) CALL ctl_stop( 'bad flag: nn_mxl is 0, 1 or 2 ' ) 785 817 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 ' ) 787 819 IF( nn_etau == 3 .AND. .NOT. ln_cpl ) CALL ctl_stop( 'nn_etau == 3 : HF taum only known in coupled mode' ) 788 820 … … 792 824 ENDIF 793 825 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 795 835 796 836 ! !* depth of penetration of surface tke 797 837 IF( nn_etau /= 0 ) THEN 838 htau(:,:) = 0._wp 798 839 SELECT CASE( nn_htau ) ! Choice of the depth of penetration 799 840 CASE( 0 ) ! constant depth penetration (here 10 meters) … … 801 842 CASE( 1 ) ! F(latitude) : 0.5m to 30m poleward of 40 degrees 802 843 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 803 868 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 804 877 ENDIF 805 878 ! !* set vertical eddy coef. to the background value -
branches/UKMO/dev_r6501_GO6_package_trunk/NEMOGCM/NEMO/OPA_SRC/step_oce.F90
r6503 r6507 80 80 81 81 USE diawri ! Standard run outputs (dia_wri routine) 82 USE diaprod ! Product diagnostics (dia_prod routine) 82 83 USE diaptr ! poleward transports (dia_ptr routine) 83 84 USE diadct ! sections transports (dia_dct routine)
Note: See TracChangeset
for help on using the changeset viewer.
