!!---------------------------------------------------------------------- !! *** flx_core.h90 *** !!---------------------------------------------------------------------- !! flx : update surface thermohaline fluxes from the NCAR data set !! read in a NetCDF file !!---------------------------------------------------------------------- !! * Modules used C A U T I O N already defined in flxmod.F90 !! * Shared module variables REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: & gsss, & !: SSS mean on nfbulk ocean time step gu , gv , & !: (un,vn)(jk=1) mean over nfbulk ocean time step ! ! these variables are used for output in diawri qlw_oce , & !: long wave flx over ocean qla_oce , & !: latent heat flx over ocean qsb_oce , & !: sensible heat flx over ocean qlw_ice , & !: long wave flx over ice qsb_ice !: sensible heat flx over ice !! * Module variables INTEGER, PARAMETER :: jpfile = 8 ! maximum number of files to read CHARACTER (LEN=32), DIMENSION (jpfile) :: clvarname CHARACTER (LEN=50), DIMENSION (jpfile) :: clname INTEGER :: isnap INTEGER, DIMENSION(jpfile) :: & numflxall, & ! logical units for surface fluxes data nrecflx , nrecflx2 ! first and second record to be read in flux file REAL(wp), DIMENSION(jpfile) :: & freqh ! Frequency for each forcing file (hours) ! ! a negative value of -12 corresponds to monthly REAL(wp), DIMENSION(jpi,jpj,jpfile,2) :: & flxdta !: set of NCAR 6hourly/daily/monthly fluxes !!---------------------------------------------------------------------- !! History : !! 9.0 ! 04-08 (U. Schweckendiek) Original code !! ! 05-04 (L. Brodeau, A.M. Treguier) severals modifications: !! ! - new bulk routine for efficiency !! ! - WINDS ARE NOW ASSUMED TO BE AT T POINTS in input files !! ! - file names and file characteristics in namelist !! ! - Implement reading of 6-hourly fields !! ! 06-02 (S. Masson, G. Madec) IOM read + NEMO "style" !!---------------------------------------------------------------------- !! OPA 9.0 , LOCEAN-IPSL (2006) !! $Header$ !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) !!---------------------------------------------------------------------- CONTAINS SUBROUTINE flx( kt ) !!--------------------------------------------------------------------- !! *** ROUTINE flx *** !! !! ** Purpose : provide the thermohaline fluxes (heat and freshwater) !! and momentum fluxes (tau) !! to the ocean at each time step. !! !! ** Method : Read NCAR data in a NetCDF file !! (file names and frequency of inputs specified in namelist) !! precipitation: total (rain+snow) !! precipitation: snow only !! u10,v10 -> scalar wind at 10m in m/s - ON 'T' GRID POINTS!!! !! solar radiation (short wave) in W/m2 !! thermal radiation (long wave) in W/m2 !! specific humidity in % !! temperature at 10m in degrees K !! !! ** Action : !! call flx_blk_albedo to compute ocean and ice albedo !! Calculates forcing fluxes to input into ice model ! or to be used directly for the ocean in ocesbc. !! !! ** Outputs !! COMMENTS TO BE ADDED -units to be verified! !! taux : zonal wind stress on "u" points (N/m2) !! tauy : meridional wind stress on "v" points (N/m2) !! qsr_oce : Solar flux over the ocean (W/m2) !! qnsr_oce: longwave flux over the ocean (W/m2) !! qsr_ice : solar flux over the ice (W/m2) !! qnsr_ice: longwave flux over the ice (W/m2) !! qla_ice : latent flux over sea-ice (W/m2) !! dqns_ice: total heat fluxes sensitivity over ice (dQ/dT) !! dqla_ice: latent heat flux sensitivity over ice (dQla/dT) !! fr1_i0 !! fr2_i0 !! tprecip !! sprecip !! evap !! !! caution : now, in the opa global model, the net upward water flux is !! ------- with mm/day unit,i.e. Kg/m2/s. !!--------------------------------------------------------------------- !! * modules used USE iom USE par_oce USE flx_oce USE blk_oce ! bulk variable USE taumod USE bulk USE phycst USE lbclnk USE dtatem ! FOR Ce = F(SST(levitus)): !! * arguments INTEGER, INTENT( in ) :: kt ! ocean time step !! * Local declarations INTEGER , PARAMETER :: jpmaxtime = 366*4 ! maximum time steps in file will be for a 6 hourly field ! and a leap year (necessary variable for flinopen??) INTEGER :: irectot, irecflx INTEGER :: ihour ! current hour of the day at which we read the forcings INTEGER :: & imois, imois2, itime, & ! temporary integers i15 , iman , & ! " " ifpr , jfpr , & ! frequency of index for print in prihre jj , ji ! Loop indices REAL(wp) :: & zadatrjb, & ! date (noninteger day) at which we read the forcings zsecond, & ! temporary scalars zxy , & ! scalar for temporal interpolation zcof ! scalar REAL(wp) :: em , aa REAL(wp) :: zind1, zind2, zind3 REAL(wp), DIMENSION(jpi,jpj, jpfile) :: & flxnow ! input flux values at current time step REAL(wp), DIMENSION(jpi,jpj) :: & sstk , t04 REAL(wp), DIMENSION(jpi,jpj) :: & qtune , & ! artifical field for tuning q dqlw_ice , & ! long wave heat flx sensitivity over ice dqsb_ice , & ! sensible heat flx sensitivity over ice alb_ice , & ! albedo of ice alb_oce_os, & ! albedo of the ocean under overcast sky alb_ice_os, & ! albedo of the ice under overcast sky alb_ice_cs, & ! albedo of ice under clear sky alb_oce_cs, & ! albedo of the ocean under clear sky zsstlev, & ! SST from Levitus in K zsst, & ! nfbulk : mean SST zsss, & ! nfbulk : mean tn_ice(:,:,1) zut, & ! nfbulk : mean U at T-point zvt, & ! nfbulk : mean V at T-point dUnormt, & ! scalar wind (norm) on T points tauxt, tauyt, & ! wind stress computed at T-point qsatw, & ! specific humidity at zSST qsat, & ! specific humidity at zsss Ch, & ! coefficient for sensible heat flux Ce, & ! coefficient for latent heat flux Cd ! coefficient for wind stress REAL(wp), PARAMETER :: & rhoa = 1.22, & ! air density cp = 1000.5, & ! specific heat of air Lv = 2.5e6, & ! latent heat of vaporization Ls = 2.839e6, & ! latent heat of sublimation Stef = 5.67e-8, & ! Stefan Boltzmann constant Cice = 1.63e-3 ! transfert coefficient over ice REAL(wp), DIMENSION(jpi,jpj) :: & catm1 ! NAMELIST /namcore/ clname, clvarname, freqh !!--------------------------------------------------------------------- !! =============================================== !! !! PART A: READING FLUX FILES WHEN NECESSARY !! !! =============================================== !! !--- calculation for monthly data i15 = INT( 2 * FLOAT( nday ) / ( FLOAT( nobis(nmonth) ) + 0.5 ) ) iman = 12 imois = nmonth + i15 - 1 IF( imois == 0 ) imois = iman imois2 = nmonth !--- calculation for 6-hourly data ! we need the fraction of day, measured in hours at the date of forcings. ! This is the time step before the date we are calculating, zadatrj ! "adatrj" (real time in days (noninteger)) zadatrjb = adatrj - rdt / rday ihour = INT( ( zadatrjb - INT(zadatrjb) ) * 24 ) ! ! ------------------------ ! IF( kt == nit000 ) THEN ! first call kt=nit000 ! ! ! ------------------------ ! !--- Initializes default values of file names, frequency of forcings ! and variable to be read in the files, before reading namelist clname(1) = 'precip_core.nc' ; freqh(1) = -12 ; clvarname(1) = 'precip' ! monthly clname(2) = 'u10_core.nc' ; freqh(2) = 24 ; clvarname(2) = 'u_10' ! 6 hourly clname(3) = 'v10_core.nc' ; freqh(3) = 24 ; clvarname(3) = 'v_10' ! 6 hourly clname(4) = 'q10_core.nc' ; freqh(4) = 24 ; clvarname(4) = 'q_10' ! 6 hourly clname(5) = 'tot_solar_core.nc'; freqh(5) = 24 ; clvarname(5) = 'tot_solar' ! daily clname(6) = 'therm_rad_core.nc'; freqh(6) = 24 ; clvarname(6) = 'therm_rad' ! daily clname(7) = 'temp_10m_core.nc' ; freqh(7) = 24 ; clvarname(7) = 't_10' ! 6 hourly clname(8) = 'snow_core.nc' ; freqh(8) = -12 ; clvarname(8) = 'snow' ! monthly !--- Read Namelist namcore : OMIP/CORE REWIND ( numnam ) READ ( numnam, namcore ) IF(lwp) THEN WRITE(numout,*)' ' WRITE(numout,*)' flxmod/flx_core : global CORE fields in NetCDF format' WRITE(numout,*)' ~~~~~~~~~~~~~~~ ' WRITE(numout,*) ' list of files and frequency (hour), or monthly (-12) ' DO ji = 1, jpfile WRITE(numout,*) trim(clname(ji)), ' frequency:', freqh(ji) END DO ENDIF !--- Initialization to zero qsr_oce (:,:) = 0.e0 qsb_oce (:,:) = 0.e0 qla_oce (:,:) = 0.e0 qlw_oce (:,:) = 0.e0 qnsr_oce(:,:) = 0.e0 qsr_ice (:,:) = 0.e0 qnsr_ice(:,:) = 0.e0 qla_ice (:,:) = 0.e0 qlw_ice (:,:) = 0.e0 qsb_ice (:,:) = 0.e0 dqns_ice(:,:) = 0.e0 dqla_ice(:,:) = 0.e0 tprecip (:,:) = 0.e0 sprecip (:,:) = 0.e0 evap (:,:) = 0.e0 flxnow (:,:,:) = 1.e-15 !RB flxdta (:,:,:,:) = 1.e-15 !RB nrecflx (:) = 0 ! switch for reading flux data for each file nrecflx2(:) = 0 ! switch for reading flux data for each file !--- Open the files of the list DO ji = 1, jpfile CALL iom_open( clname(ji), numflxall(ji) ) END DO ENDIF ! ! ------------------------ ! ! ! Read files if required ! ! ! ------------------------ ! !--- read data if necessary checks each file in turn DO ji = 1, jpfile IF ( freqh(ji) .GT.0 .AND. freqh(ji) .LT. 24 ) THEN !--- Case of 6-hourly flux data !--- calculate current snapshot from hour of day isnap = ihour / INT( freqh(ji) ) + 1 !--- reading is necessary when nrecflx(ji) differs from isnap IF( nrecflx(ji) /= isnap ) THEN nrecflx(ji) = isnap irecflx = (nday_year-1)*24 / freqh(ji) + isnap irectot = 365*24 / freqh(ji) ! this variable is not used by iom_get IF(lwp) WRITE (numout,*)' CORE flux 6-h record :',irecflx, ' var:',trim(clvarname(ji)) CALL iom_get( numflxall(ji), jpdom_data, clvarname(ji), flxdta(:,:,ji,1), irecflx ) ENDIF ELSE IF ( freqh(ji) .EQ. 24 ) THEN !--- Case of daily flux data !--- reading is necessary when nrecflx(ji) differs from nday IF( nrecflx(ji) /= nday ) THEN nrecflx(ji) = nday !! remember present read day irecflx = nday_year !! irectot = 365 ! this variable is not used by iom_get IF(lwp) WRITE (numout,*)'DAILY CORE flux record :',irecflx, ' var:',trim(clvarname(ji)) CALL iom_get( numflxall(ji), jpdom_data, clvarname(ji), flxdta(:,:,ji,1), irecflx ) ENDIF ELSE IF ( freqh(ji) .EQ. -12 ) THEN !--- Case monthly data from CORE ! we read two months all the time although we ! could only read one and swap the arrays IF( kt == nit000 .OR. imois /= nrecflx(ji) ) THEN ! calendar computation ! nrecflx number of the first file record used in the simulation ! nrecflx2 number of the last file record nrecflx(ji) = imois nrecflx2(ji) = nrecflx(ji)+1 nrecflx(ji) = MOD( nrecflx(ji), iman ) IF( nrecflx(ji) == 0 ) nrecflx(ji) = iman nrecflx2(ji) = MOD( nrecflx2(ji), iman ) IF ( nrecflx2(ji) == 0 ) nrecflx2(ji) = iman irectot = 12 ! this variable is not used by iom_get IF(lwp) WRITE(numout,*) 'MONTHLY CORE flux record 1:', nrecflx(ji), & & ' rec 2:', nrecflx2(ji), ' var:', trim(clvarname(ji)) CALL iom_get( numflxall(ji), jpdom_data, clvarname(ji), flxdta(:,:,ji,1), nrecflx (ji) ) CALL iom_get( numflxall(ji), jpdom_data, clvarname(ji), flxdta(:,:,ji,2), nrecflx2(ji) ) ENDIF ENDIF END DO ! end of loop over forcing files to verify if reading is necessary !--- Printout if required ! IF( MOD( kt - 1 , nfbulk ) == 0 ) THEN ! printout at the initial time step only (unless you want to debug!!) IF( kt == nit000 ) THEN IF(lwp) THEN WRITE(numout,*) WRITE(numout,*) ' read daily and monthly CORE fluxes: ok' WRITE(numout,*) ifpr = INT(jpi/8) ; jfpr = INT(jpj/10) DO ji = 1, jpfile WRITE(numout,*) trim(clvarname(ji)),' day: ',ndastp CALL prihre(flxdta(1,1,ji,1),jpi,jpj,1,jpi,ifpr,1,jpj,jfpr,0.,numout) WRITE(numout,*) END DO ENDIF ENDIF ! ! ------------------------ ! ! ! Interpolates in time ! ! ! ------------------------ ! ! N.B. For now, only monthly values are interpolated, and they ! are interpolated to the current day, not to the time step. ! DO ji = 1, jpfile IF( freqh(ji) .EQ. -12 ) THEN zxy = REAL(nday) / REAL( nobis(imois) ) + 0.5 - i15 !!! <== Caution nobis hard coded !!! flxnow(:,:,ji) = ( ( 1. - zxy) * flxdta(:,:,ji,1) + zxy * flxdta(:,:,ji,2) ) ELSE flxnow(:,:,ji) = flxdta(:,:,ji,1) ENDIF END DO ! JMM : add vatm needed in tracer routines ==> wind module ??? vatm(:,:) = SQRT( flxnow(:,:,2) * flxnow(:,:,2) + flxnow(:,:,3) * flxnow(:,:,3) ) !! =============================================== !! !! PART B: CORE BULK CALCULATION !! !! =============================================== !! ! for other forcing cases this is done in modules bulk.F90 and flxblk ! ! ------------------------ ! ! ! Bulk initialisation ! ! ! ------------------------ ! ! code part from bulk.F90 : IF( kt == nit000) THEN !--- computation of rdtbs2 ===> !gm is it really usefull ???? IF( nacc == 1 ) THEN rdtbs2 = nfbulk * rdtmin * 0.5 ELSE rdtbs2 = nfbulk * rdt * 0.5 ENDIF IF( .NOT.ln_rstart ) THEN zcof = REAL( nfbulk - 1, wp ) gsst(:,:) = zcof * ( tn(:,:,1) + rt0 ) gsss(:,:) = zcof * tn_ice(:,:) gu (:,:) = zcof * un(:,:,1) gv (:,:) = zcof * vn(:,:,1) ENDIF ENDIF ! ----------------------------------------------------------------------------- ! ! 0 Mean SST, ice temperature and ocean currents over nfbulk pdt ! ! ----------------------------------------------------------------------------- ! ! gsst(:,:) = gsst(:,:) + (tn(:,:,1) + rt0 ) gsss(:,:) = gsss(:,:) + tn_ice(:,:) gu (:,:) = gu (:,:) + un (:,:,1) gv (:,:) = gv (:,:) + vn (:,:,1) IF( MOD( kt - 1 , nfbulk ) == 0 ) THEN ! zsst(:,:) = gsst(:,:) / REAL( nfbulk, wp ) * tmask(:,:,1) ! mean sst in K zsss(:,:) = gsss(:,:) / REAL( nfbulk ) * tmask(:,:,1) ! mean tn_ice in K where ( zsst(:,:) .eq. 0.e0 ) zsst(:,:) = rt0 !lb vilain !!??? where ( zsss(:,:) .eq. 0.e0 ) zsss(:,:) = rt0 !lb // !!gm better coded: ! Caution set to rt0 over land, not 0 Kelvin (otherwise floating point exception in bulk computation) ! zcof = 1. / REAL( nfbulk, wp ) ! zsst(:,:) = gsst(:,:) * zcof * tmask(:,:,1) + rt0 * (1.- tmask(:,:,1) ! mean sst in K ! zsss(:,:) = gsss(:,:) * zcof * tmask(:,:,1) + rt0 * (1.- tmask(:,:,1) ! mean tn_ice in K !!gm zut(:,:) = 0.e0 !!gm not necessary but at least first and last column zvt(:,:) = 0.e0 ! lb ! Interpolation of surface current at T-point, zut and zvt : DO ji=2, jpi zut(ji,:) = 0.5*(gu(ji-1,:) + gu(ji,:)) / REAL( nfbulk ) END DO ! DO jj=2, jpj zvt(:,jj) = 0.5*(gv(:,jj-1) + gv(:,jj)) / REAL( nfbulk ) END DO !!gm better coded ! zcof = 0.5 / REAL( nfbulk, wp ) ! DO jj = 2, jpjm1 ! DO ji = fs_2, fs_jpim1 ! vector opt. ! zut(ji,jj) = ( gu(ji-1,jj ) + gu(ji,jj) ) * zcof ! zvt(ji,jj) = ( gv(ji ,jj-1) + gv(ji,jj) ) * zcof ! END DO ! END DO !!gm CALL lbc_lnk( zut(:,:), 'T', -1. ) CALL lbc_lnk( zvt(:,:), 'T', -1. ) ! ----------------------------------------------------------------------------- ! ! I Radiative FLUXES ! ! ----------------------------------------------------------------------------- ! !--- Ocean & Ice Albedo alb_oce_os(:,:) = 0. ; alb_oce_cs(:,:) = 0. !gm is it necessary ??? alb_ice_os(:,:) = 0. ; alb_ice_cs(:,:) = 0. CALL flx_blk_albedo( alb_ice_os, alb_oce_os, alb_ice_cs, alb_oce_cs ) !--- Radiative fluxes over ocean qsr_oce(:,:) = ( 1. - 0.066 ) * flxnow(:,:,5) ! Solar Radiation qlw_oce(:,:) = flxnow(:,:,6) - Stef*zsst(:,:)*zsst(:,:)*zsst(:,:)*zsst(:,:) ! Long Waves (Infra-red) !--- Radiative fluxes over ice qsr_ice(:,:) = ( 1. - alb_ice_cs(:,:) ) * flxnow(:,:,5) ! Solar Radiation qlw_ice(:,:) = 0.95 * ( flxnow(:,:,6) - Stef*zsss(:,:)*zsss(:,:)*zsss(:,:)*zsss(:,:) ) ! Long Waves (Infra-red) !-------------------------------------------------------------------------------! ! ----------------------------------------------------------------------------- ! ! II Turbulent FLUXES ! ! ----------------------------------------------------------------------------- ! ! ! scalar wind ( = | U10m - SSU | ) ! It is important to take into account the sea surface courant ! lb ! Now, wind components are provided on T-points within the netcdf input file. ! Thus, the wind module is computded at T-points taking into account the sea !--- Module of relative wind dUnormt(:,:) = sqrt( (flxnow(:,:,2) - zut(:,:))*(flxnow(:,:,2) - zut(:,:)) & & + (flxnow(:,:,3) - zvt(:,:))*(flxnow(:,:,3) - zvt(:,:)) ) !!gm more efficient coding: ! DO jj = 1, jpi ! DO ji = 1, jpi ! zzu = flxnow(ji,jj,2) - zut(ji,jj) ! zzv = flxnow(ji,jj,3) - zvt(ji,jj) ! dUnormt(ji,jj) = SQRT( zzu*zzu + zzv*zzv ) ! END DO ! END DO !!gm ! lb ! !--- specific humidity at temp SST qsatw(:,:) = 0.98*640380*exp(-5107.4/zsst(:,:))/rhoa ! !--- specific humidity at temp tn_ice qsat(:,:) = 11637800*exp(-5897.8/zsss(:,:))/rhoa ! !--- NCAR Bulk formulae, computation of Cd, Ch, Ce at T-point (NCAR Bulk formulae) CALL TURB_CORE( 10., zsst(:,:), flxnow(:,:,7), qsatw(:,:), flxnow(:,:,4), & & dUnormt(:,:), Cd(:,:), Ch(:,:), Ce(:,:) ) ! II.1) Momentum over ocean and ice ! --------------------------------- ! Tau_x at T-point tauxt(:,:) = rhoa*dUnormt(:,:)*( (1. - freeze(:,:))*Cd(:,:)*(flxnow(:,:,2) - zut(:,:)) & & + freeze(:,:)*Cice*flxnow(:,:,2) ) !lb correct pour glace ! Tau_y at T-point tauyt(:,:) = rhoa*dUnormt(:,:)*( (1. - freeze(:,:))*Cd(:,:)*(flxnow(:,:,3) - zvt(:,:)) & & + freeze(:,:)*Cice*flxnow(:,:,3) ) ! CALL lbc_lnk( tauxt(:,:), 'T', -1. ) !!gm seems not decessary at this point.... CALL lbc_lnk( tauyt(:,:), 'T', -1. ) ! !Tau_x at U-point DO ji = 1, jpim1 taux(ji,:) = 0.5*(tauxt(ji,:) + tauxt(ji+1,:)) END DO ! ! Tau_y at V-point DO jj = 1, jpjm1 tauy(:,jj) = 0.5*(tauyt(:,jj) + tauyt(:,jj+1)) END DO !!gm ==> at this point, a lbc is required, no? ! ! lb : should we do this here? !!gm yes should we remove that??? tauxg(:,:) = taux(:,:) ! Save components in tauyg(:,:) = tauy(:,:) ! geographical ref on U grid ! ! ! II.2) Turbulent fluxes over ocean ! --------------------------------- ! ! Sensible Heat : qsb_oce(:,:) = rhoa * cp * Ch(:,:) * ( zsst(:,:) - flxnow(:,:,7) ) * dUnormt(:,:) ! ! Latent Heat : evap(:,:) = rhoa * Ce(:,:) * ( qsatw(:,:) - flxnow(:,:,4) ) * dUnormt(:,:) qla_oce(:,:) = Lv * evap(:,:) ! ! ! II.3) Turbulent fluxes over ice ! ------------------------------- ! ! Sensible Heat : qsb_ice(:,:) = rhoa*cp*Cice*(zsss(:,:) - flxnow(:,:,7))*dUnormt(:,:) !lb use ! !lb dUnormt??? ! Latent Heat : !lb or rather Unormt? qla_ice(:,:) = Ls*rhoa*Cice*(qsat(:,:) - flxnow(:,:,4))*dUnormt(:,:) ! !------------------------------------------------------------------------------------- ! ----------------------------------------------------------------------------- ! ! III Total FLUXES ! ! ----------------------------------------------------------------------------- ! ! ! III.1) Downward Non Solar flux over ocean ! ----------------------------------------- qnsr_oce(:,:) = qlw_oce(:,:) - qsb_oce(:,:) - qla_oce(:,:) ! ! III.1) Downward Non Solar flux over ice ! --------------------------------------- qnsr_ice(:,:) = qlw_ice(:,:) - qsb_ice(:,:) - qla_ice(:,:) ! !-------------------------------------------------------------------------------! ! 6. TOTAL NON SOLAR SENSITIVITY dqlw_ice(:,:)= 4.0*0.95*Stef*zsss(:,:)*zsss(:,:)*zsss(:,:) ! d qla_ice/ d zsss dqla_ice(:,:) = -Ls*Cice*0.98*11637800/(rhoa*rhoa) & * (-5897.8)/(zsss(:,:)*zsss(:,:)) & * exp(-5897.8/zsss(:,:)) * dUnormt(:,:) ! d qsb_ice/ d zsss dqsb_ice(:,:) = rhoa * cp * Cice * dUnormt(:,:) dqns_ice(:,:) = - ( dqlw_ice(:,:) + dqsb_ice(:,:) + dqla_ice(:,:) ) !-------------------------------------------------------------------- ! FRACTION of net shortwave radiation which is not absorbed in the ! thin surface layer and penetrates inside the ice cover ! ( Maykut and Untersteiner, 1971 ; Elbert and Curry, 1993 ) catm1(:,:) = 1.0 - 0.3 ! flxnow(:,:,8) fr1_i0(:,:) = & (0.18 * catm1(:,:) + 0.35 * flxnow(:,:,8)) fr2_i0(:,:) = & (0.82 * catm1(:,:) + 0.65 * flxnow(:,:,8)) ! Total PRECIPITATION (kg/m2/s) tprecip(:,:) = flxnow(:,:,1) ! rename precipitation for freshwater budget calculations: watm(:,:) = flxnow(:,:,1)*1000 ! ! SNOW PRECIPITATION (kg/m2/s) sprecip(:,:) = flxnow(:,:,8) !--------------------------------------------------------------------- CALL lbc_lnk( taux (:,:) , 'U', -1. ) CALL lbc_lnk( tauy (:,:) , 'V', -1. ) CALL lbc_lnk( qsr_oce (:,:) , 'T', 1. ) CALL lbc_lnk( qnsr_oce(:,:) , 'T', 1. ) CALL lbc_lnk( qsr_ice (:,:) , 'T', 1. ) CALL lbc_lnk( qnsr_ice(:,:) , 'T', 1. ) CALL lbc_lnk( qla_ice (:,:) , 'T', 1. ) CALL lbc_lnk( dqns_ice(:,:) , 'T', 1. ) CALL lbc_lnk( dqla_ice(:,:) , 'T', 1. ) CALL lbc_lnk( fr1_i0 (:,:) , 'T', 1. ) CALL lbc_lnk( fr2_i0 (:,:) , 'T', 1. ) CALL lbc_lnk( tprecip (:,:) , 'T', 1. ) CALL lbc_lnk( sprecip (:,:) , 'T', 1. ) CALL lbc_lnk( evap (:,:) , 'T', 1. ) !! !! !! NEVER mask the windstress!! qsr_oce (:,:) = qsr_oce (:,:)*tmask(:,:,1) qnsr_oce(:,:) = qnsr_oce(:,:)*tmask(:,:,1) qsr_ice (:,:) = qsr_ice (:,:)*tmask(:,:,1) qnsr_ice(:,:) = qnsr_ice(:,:)*tmask(:,:,1) qla_ice (:,:) = qla_ice (:,:)*tmask(:,:,1) dqns_ice(:,:) = dqns_ice(:,:)*tmask(:,:,1) dqla_ice(:,:) = dqla_ice(:,:)*tmask(:,:,1) fr1_i0 (:,:) = fr1_i0 (:,:)*tmask(:,:,1) fr2_i0 (:,:) = fr2_i0 (:,:)*tmask(:,:,1) tprecip (:,:) = tprecip (:,:)*tmask(:,:,1) sprecip (:,:) = sprecip (:,:)*tmask(:,:,1) evap (:,:) = evap (:,:)*tmask(:,:,1) gsst(:,:) = 0.e0 gsss(:,:) = 0.e0 gu (:,:) = 0.e0 gv (:,:) = 0.e0 ! Degugging print (A.M. Treguier) ! write (55) taux, tauy, qnsr_oce, qsb_ice, qnsr_ice, qla_ice, dqns_ice & ! & , dqla_ice, fr1_i0, fr2_i0, qlw_oce, qla_oce, qsb_oce, evap ! write(numout,*) 'written 14 arrays on unit fort.55' ENDIF ! ------------------- ! ! Last call kt=nitend ! ! ------------------- ! ! Closing of the numflx file (required in mpp) IF( kt == nitend ) THEN DO ji=1, jpfile CALL iom_close(numflxall(ji)) END DO ENDIF END SUBROUTINE flx SUBROUTINE TURB_CORE( zzu, T_0, T_a, q_sat, q_a, & & dU , C_d, C_h , C_e ) !!---------------------------------------------------------------------- !! *** ROUTINE turb_core *** !! !! ** Purpose : Computes turbulent transfert coefficients of surface !! fluxes according to Large & Yeager (2004) !! !! ** Method : I N E R T I A L D I S S I P A T I O N M E T H O D !! Momentum, Latent and sensible heat exchange coefficients !! Caution: this procedure should only be used in cases when air !! temperature (T_air), air specific humidity (q_air) and wind (dU) !! are provided at the same height 'zzu'! !! !! References : !! Large & Yeager, 2004 : ??? !! History : !! ! XX-XX (??? ) Original code !! 9.0 ! 05-08 (L. Brodeau) Optimisation !!---------------------------------------------------------------------- !! * Arguments REAL(wp), INTENT( in ) :: & zzu ! height for all given atmospheric variables [m] REAL(wp), INTENT( in ), DIMENSION(jpi,jpj) :: & T_0, & ! sea surface temperature [Kelvin] T_a, & ! potential air temperature at zzu [Kelvin] q_sat, & ! sea surface specific humidity at zzu [kg/kg] q_a, & ! specific air humidity at zzu [kg/kg] dU ! relative wind module |U(zzu)-U(0)| at zzu [m/s] REAL(wp), INTENT(inout), DIMENSION(jpi,jpj) :: & C_d, & ! transfer coefficient for momentum (tau) C_h, & ! transfer coefficient for sensible heat (Q_sens) C_e ! tansfert coefficient for evaporation (Q_lat) !! * Local declarations REAL(wp), DIMENSION(jpi,jpj) :: & dU10, & ! dU [m/s] dT, & ! air/sea temperature differeence [K] dq, & ! air/sea humidity difference [K] Cd_n10, & ! 10m neutral drag coefficient Ce_n10, & ! 10m neutral latent coefficient Ch_n10, & ! 10m neutral sensible coefficient Cd, & ! drag coefficient Ce, & ! latent coefficient Ch, & ! sensible coefficient sqrt_Cd_n10, & ! root square of Cd_n10 sqrt_Cd, & ! root square of Cd T_vpot, & ! virtual potential temperature [K] T_star, & ! turbulent scale of tem. fluct. q_star, & ! turbulent humidity of temp. fluct. U_star, & ! turb. scale of velocity fluct. L, & ! Monin-Obukov length [m] zeta, & ! stability parameter at height zzu X2, X, & psi_m, & psi_h, & U_n10, & ! neutral wind velocity at 10m [m] xlogt INTEGER :: jkk,ji,jj,ii,ij,ilocu(2), jk ! doomy loop for iterations INTEGER, PARAMETER :: nit = 3 ! number of iterations REAL(wp), DIMENSION(jpi,jpj) :: dbg1,dbg2,dbg3,dbg4 REAL :: zXmax INTEGER, DIMENSION(jpi,jpj) :: & & stab, & ! stability test integer & stabit ! stability within iterative loop REAL(wp), PARAMETER :: & & pi = 3.14159, & ! Pi & grav = 9.8, & ! gravity & kappa = 0.4 ! von Karman's constant !!---------------------------------------------------------------------- !! !! ------------- !! S T A R T !! ------------- !! !! I.1 ) Preliminary stuffs !! ------------------------ !! !! Air/sea differences !! ------------------- dU10 = max(0.5, dU) ! we don't want to fall under 0.5 m/s dT = T_a - T_0 ! assuming that T_a is allready the potential temp. at zzu dq = q_a - q_sat !! !! Virtual potential temperature !! ----------------------------- T_vpot = T_a*(1. + 0.608*q_a) !! !! !! I.2 ) Computing Neutral Drag Coefficient !! ---------------------------------------- Cd_n10 = 1E-3 * ( 2.7/dU10 + 0.142 + dU10/13.09 ) ! \\ L & Y eq. (6a) sqrt_Cd_n10 = sqrt(Cd_n10) !! Ce_n10 = 1E-3 * ( 34.6 * sqrt_Cd_n10 ) ! \\ L & Y eq. (6b) !! !! First guess of stabilitty : stab = 0.5 + sign(0.5,dT) ! stable : stab = 1 ; unstable : stab = 0 Ch_n10 = 1E-3*sqrt_Cd_n10*(18*stab + 32.7*(1-stab)) ! \\ L & Y eq. (6c), (6d) !! !! Initializing transfert coefficients with their first guess neutral equivalents : Cd = Cd_n10 ; Ce = Ce_n10 ; Ch = Ch_n10 !! !! !! !! II ) Now starting iteration loop (IDM) !! ---------------------------------------- !! DO jk=1, nit !! sqrt_Cd = sqrt(Cd) !! !! Turbulent scales : !! ------------------ U_star = sqrt_Cd*dU10 ! \\ L & Y eq. (7a) T_star = Ch/sqrt_Cd*dT ! \\ L & Y eq. (7b) q_star = Ce/sqrt_Cd*dq ! \\ L & Y eq. (7c) !! !! Estimate the Monin-Obukov length : !! ---------------------------------- !dbg1 = T_star/T_vpot !dbg2 = q_star/(q_a + 1./0.608) L = (U_star**2)/( kappa*grav*(T_star/T_vpot + q_star/(q_a + 1./0.608)) ) !! !! Stability parameters : !! ---------------------- zeta = zzu/L zeta = sign( min(abs(zeta),10.0), zeta ) !! !! Psis, L & Y eq. (8c), (8d), (8e) : !! ---------------------------------- X2 = sqrt(abs(1. - 16.*zeta)) ; X2 = max(X2 , 1.0) ; X = sqrt(X2) !! stabit = 0.5 + sign(0.5,zeta) !! ! dbg1 = -5*zeta*stabit ! dbg2 = 2*log((1. + X)/2) ! dbg3 = log((1. + X2)/2) ! dbg4 = 2*atan(X) psi_m = -5*zeta*stabit & ! Stable & + (1 - stabit)*(2*log((1. + X)/2) + log((1. + X2)/2) - 2*atan(X) + pi/2) ! Unstable !! psi_h = -5*zeta*stabit & ! Stable & + (1 - stabit)*(2*log( (1. + X2)/2 )) ! Unstable !! !! !! Shifting the wind speed to 10m and neutral stability : !! ------------------------------------------------------ U_n10 = dU10*1./(1. + sqrt(Cd_n10)/kappa*(log(zzu/10.) - psi_m)) ! \\ L & Y eq. (9a) !! !! !! Updating the neutral 10m transfer coefficients : !! ------------------------------------------------ Cd_n10 = 1E-3 * (2.7/U_n10 + 0.142 + U_n10/13.09) ! \\ L & Y eq. (6a) sqrt_Cd_n10 = sqrt(Cd_n10) !! Ce_n10 = 1E-3 * (34.6 * sqrt_Cd_n10) ! \\ L & Y eq. (6b) !! stab = 0.5 + sign(0.5,zeta) Ch_n10 = 1E-3*sqrt_Cd_n10*(18.*stab + 32.7*(1-stab)) ! \\ L & Y eq. (6c), (6d) !! !! !! Shifting the neutral 10m transfer coefficients to ( zzu , zeta ) : !! -------------------------------------------------------------------- !! Problem here, formulation used within L & Y differs from the one provided !! in their fortran code (only for Ce and Ch) !! Cd = Cd_n10/(1. + sqrt_Cd_n10/kappa*(log(zzu/10) - psi_m))**2 ! \\ L & Y eq. (10a) !! xlogt = log(zzu/10) - psi_h !? Ch = Ch_n10*sqrt(Cd/Cd_n10)/(1. + Ch_n10/(kappa*sqrt_Cd_n10)*xlogt) Ch = Ch_n10/( 1. + Ch_n10*xlogt/kappa/sqrt_Cd_n10 )**2 ! \\ L & Y eq. (10b) !! !? Ce = Ce_n10*sqrt(Cd/Cd_n10)/(1. + Ce_n10/(kappa*sqrt_Cd_n10)*xlogt) Ce = Ce_n10/( 1. + Ce_n10*xlogt/kappa/sqrt_Cd_n10 )**2 ! \\ L & Y eq. (10c) !! !! END DO !! C_d = Cd C_h = Ch C_e = Ce !! END SUBROUTINE TURB_CORE SUBROUTINE flx_blk_albedo( palb , palcn , palbp , palcnp ) !!---------------------------------------------------------------------- !! *** ROUTINE flx_blk_albedo *** !! !! ** Purpose : Computation of the albedo of the snow/ice system !! as well as the ocean one !! !! ** Method : - Computation of the albedo of snow or ice (choose the !! right one by a large number of tests !! - Computation of the albedo of the ocean !! !! References : !! Shine and Hendersson-Sellers 1985, JGR, 90(D1), 2243-2250. !! !! History : !! 8.0 ! 01-04 (LIM 1.0) !! 8.5 ! 03-07 (C. Ethe, G. Madec) Optimization (old name:shine) !!---------------------------------------------------------------------- !! * Modules used USE ice ! ??? !! * Arguments REAL(wp), DIMENSION(jpi,jpj), INTENT(out) :: & palb , & ! albedo of ice under overcast sky palcn , & ! albedo of ocean under overcast sky palbp , & ! albedo of ice under clear sky palcnp ! albedo of ocean under clear sky !! * Local variables INTEGER :: & ji, jj ! dummy loop indices REAL(wp) :: & c1 = 0.05 , & ! constants values c2 = 0.1 , & albice = 0.50 , & ! albedo of melting ice in the arctic and antarctic (Shine & Hendersson-Sellers) cgren = 0.06 , & ! correction of the snow or ice albedo to take into account ! effects of cloudiness (Grenfell & Perovich, 1984) alphd = 0.80 , & ! coefficients for linear interpolation used to compute alphdi = 0.72 , & ! albedo between two extremes values (Pyane, 1972) alphc = 0.65 , & zmue = 0.4 , & ! cosine of local solar altitude zzero = 0.0 , & zone = 1.0 REAL(wp) :: & zmue14 , & ! zmue**1.4 zalbpsnm , & ! albedo of ice under clear sky when snow is melting zalbpsnf , & ! albedo of ice under clear sky when snow is freezing zalbpsn , & ! albedo of snow/ice system when ice is coverd by snow zalbpic , & ! albedo of snow/ice system when ice is free of snow zithsn , & ! = 1 for hsn >= 0 ( ice is cov. by snow ) ; = 0 otherwise (ice is free of snow) zitmlsn , & ! = 1 freezinz snow (sist >=rt0_snow) ; = 0 melting snow (sist c1 zihsc2 ! = 1 hsn >= c2 ; = 0 hsn < c2 REAL(wp), DIMENSION(jpi,jpj) :: & zalbfz , & ! ( = alphdi for freezing ice ; = albice for melting ice ) zficeth ! function of ice thickness LOGICAL , DIMENSION(jpi,jpj) :: & llmask !! to be included for without seaice !! REAL(wp), DIMENSION(jpi,jpj) :: & !: !! sist , & !: Sea-Ice Surface Temperature (Kelvin ) !! hsnif , & !: Snow thickness !! hicif !: Ice thickness !!--------------------------------------------------------------------- !------------------------- ! Computation of zficeth !-------------------------- llmask = (hsnif == 0.0) .AND. ( sist >= rt0_ice ) WHERE ( llmask ) ! ice free of snow and melts zalbfz = albice ELSEWHERE zalbfz = alphdi END WHERE DO jj = 1, jpj DO ji = 1, jpi IF( hicif(ji,jj) > 1.5 ) THEN zficeth(ji,jj) = zalbfz(ji,jj) ELSEIF( hicif(ji,jj) > 1.0 .AND. hicif(ji,jj) <= 1.5 ) THEN zficeth(ji,jj) = 0.472 + 2.0 * ( zalbfz(ji,jj) - 0.472 ) * ( hicif(ji,jj) - 1.0 ) ELSEIF( hicif(ji,jj) > 0.05 .AND. hicif(ji,jj) <= 1.0 ) THEN zficeth(ji,jj) = 0.2467 + 0.7049 * hicif(ji,jj) & & - 0.8608 * hicif(ji,jj) * hicif(ji,jj) & & + 0.3812 * hicif(ji,jj) * hicif(ji,jj) * hicif (ji,jj) ELSE zficeth(ji,jj) = 0.1 + 3.6 * hicif(ji,jj) ENDIF END DO END DO !----------------------------------------------- ! Computation of the snow/ice albedo system !-------------------------- --------------------- ! Albedo of snow-ice for clear sky. !----------------------------------------------- DO jj = 1, jpj DO ji = 1, jpi ! Case of ice covered by snow. ! melting snow zihsc1 = 1.0 - MAX ( zzero , SIGN ( zone , - ( hsnif(ji,jj) - c1 ) ) ) zalbpsnm = ( 1.0 - zihsc1 ) & * ( zficeth(ji,jj) + hsnif(ji,jj) * ( alphd - zficeth(ji,jj) ) / c1 ) & & + zihsc1 * alphd ! freezing snow zihsc2 = MAX ( zzero , SIGN ( zone , hsnif(ji,jj) - c2 ) ) zalbpsnf = ( 1.0 - zihsc2 ) * & ( albice + hsnif(ji,jj) * ( alphc - albice ) / c2 ) & & + zihsc2 * alphc zitmlsn = MAX ( zzero , SIGN ( zone , sist(ji,jj) - rt0_snow ) ) zalbpsn = zitmlsn * zalbpsnf + ( 1.0 - zitmlsn ) * zalbpsnm ! Case of ice free of snow. zalbpic = zficeth(ji,jj) ! albedo of the system zithsn = 1.0 - MAX ( zzero , SIGN ( zone , - hsnif(ji,jj) ) ) palbp(ji,jj) = zithsn * zalbpsn + ( 1.0 - zithsn ) * zalbpic END DO END DO ! Albedo of snow-ice for overcast sky. !---------------------------------------------- palb(:,:) = palbp(:,:) + cgren !-------------------------------------------- ! Computation of the albedo of the ocean !-------------------------- ----------------- ! Parameterization of Briegled and Ramanathan, 1982 zmue14 = zmue**1.4 palcnp(:,:) = 0.05 / ( 1.1 * zmue14 + 0.15 ) ! Parameterization of Kondratyev, 1969 and Payne, 1972 palcn(:,:) = 0.06 END SUBROUTINE flx_blk_albedo