Changeset 3894
- Timestamp:
- 2013-04-27T17:46:22+02:00 (11 years ago)
- Location:
- branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM
- Files:
-
- 24 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/CONFIG/ORCA2_LIM/EXP00/namelist
r3893 r3894 877 877 ln_glo_trd = .FALSE. ! (T) global domain averaged diag for T, T^2, KE, and PE 878 878 ln_dyn_trd = .FALSE. ! (T) 3D momentum trend output 879 ln_dyn_m ld= .FALSE. ! (T) 2D momentum trends averaged over the mixed layer879 ln_dyn_mxl = .FALSE. ! (T) 2D momentum trends averaged over the mixed layer 880 880 ln_vor_trd = .FALSE. ! (T) 2D barotropic vorticity trends 881 881 ln_KE_trd = .FALSE. ! (T) 3D Kinetic Energy trends 882 882 ln_PE_trd = .FALSE. ! (T) 3D Potential Energy trends 883 883 ln_tra_trd = .FALSE. ! (T) 3D tracer trend output 884 ln_tra_m ld= .FALSE. ! (T) 2D tracer trends averaged over the mixed layer884 ln_tra_mxl = .FALSE. ! (T) 2D tracer trends averaged over the mixed layer 885 885 nn_trd = 365 ! print frequency (ln_glo_trd=T) (unit=time step) 886 886 / … … 894 894 cn_trdrst_in = "restart_mld" ! suffix of ocean restart name (input) 895 895 cn_trdrst_out = "restart_mld" ! suffix of ocean restart name (output) 896 ln_trdm ld_restart = .false. ! restart for ML diagnostics897 ln_trdm ld_instant = .false. ! flag to diagnose trends of instantantaneous or mean ML T/S896 ln_trdmxl_restart = .false. ! restart for ML diagnostics 897 ln_trdmxl_instant = .false. ! flag to diagnose trends of instantantaneous or mean ML T/S 898 898 / 899 899 !----------------------------------------------------------------------- -
branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/LIM_SRC_2/limistate_2.F90
r3625 r3894 71 71 IF( .NOT. ln_limini ) THEN 72 72 73 tfu(:,:) = tfreez( tsn(:,:,1,jp_sal) ) * tmask(:,:,1) ! freezing/melting point of sea water [Celcius]73 tfu(:,:) = eos_fzp( tsn(:,:,1,jp_sal) ) * tmask(:,:,1) ! freezing/melting point of sea water [Celcius] 74 74 75 75 DO jj = 1, jpj -
branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/LIM_SRC_3/limistate.F90
r3625 r3894 94 94 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ' 95 95 96 t_bo(:,:) = tfreez( tsn(:,:,1,jp_sal) ) * tmask(:,:,1) ! freezing/melting point of sea water [Celcius]96 t_bo(:,:) = eos_fzp( tsn(:,:,1,jp_sal) ) * tmask(:,:,1) ! freezing/melting point of sea water [Celcius] 97 97 98 98 DO jj = 1, jpj ! ice if sst <= t-freez + ttest -
branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OFF_SRC/dtadyn.F90
r3827 r3894 243 243 tsn(:,:,:,jp_sal) = sf_dyn(jf_sal)%fnow(:,:,:) * tmask(:,:,:) ! salinity 244 244 ! 245 CALL eos ( tsn, rhd, rhop )! In any case, we need rhop246 CALL zdf_mxl( kt ) 245 CALL eos( tsn, fsdept(:,:,:), rhd, rhop ) ! In any case, we need rhop 246 CALL zdf_mxl( kt ) ! In any case, we need mxl 247 247 ! 248 248 avt(:,:,:) = sf_dyn(jf_avt)%fnow(:,:,:) * tmask(:,:,:) ! vertical diffusive coefficient … … 546 546 !!--------------------------------------------------------------------- 547 547 #if defined key_ldfslp && ! defined key_c1d 548 CALL eos( pts, rhd, rhop ) ! Time-filtered in situ density549 CALL bn2( pts, rn2 )! before Brunt-Vaisala frequency548 CALL eos( pts, fsdept(:,:,:), rhd, rhop ) ! Time-filtered in situ density 549 CALL eos( tsn, rab_n, rn2 ) ! before Brunt-Vaisala frequency 550 550 IF( ln_zps ) & 551 551 & CALL zps_hde( kt, jpts, pts, gtsu, gtsv, rhd, gru, grv ) ! Partial steps: before Horizontal DErivative -
branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90
r3785 r3894 679 679 !!---------------------------------------------------------------------- 680 680 681 ! freezing point calculation taken from oc_fz_pt (but calculated for all depths)681 ! freezing point 682 682 ! used to prevent the applied increments taking the temperature below the local freezing point 683 684 683 #if defined key_cice 685 684 fzptnz(:,:,:) = -1.8_wp 686 685 #else 687 DO jk = 1, jpk 688 DO jj = 1, jpj 689 DO ji = 1, jpk 690 fzptnz (ji,jj,jk) = ( -0.0575_wp + 1.710523e-3_wp * SQRT( tsn(ji,jj,jk,jp_sal) ) & 691 - 2.154996e-4_wp * tsn(ji,jj,jk,jp_sal) ) * tsn(ji,jj,jk,jp_sal) & 692 - 7.53e-4_wp * fsdepw(ji,jj,jk) ! (pressure in dbar) 693 END DO 694 END DO 695 END DO 686 DO jk = 1, jpk 687 fzptnz(:,:,jk) =eos_fzp( tsn(:,:,jk,jp_sal), fsdept(:,:,jk) ) 688 END DO 696 689 #endif 697 690 … … 709 702 IF(lwp) THEN 710 703 WRITE(numout,*) 711 WRITE(numout,*) 'tra_asm_inc : Tracer IAU at time step = ', & 712 & kt,' with IAU weight = ', wgtiau(it) 704 WRITE(numout,*) 'tra_asm_inc : Tracer IAU at time step = ', kt,' with IAU weight = ', wgtiau(it) 713 705 WRITE(numout,*) '~~~~~~~~~~~~' 714 706 ENDIF … … 758 750 IF (ln_temnofreeze) THEN 759 751 ! Do not apply negative increments if the temperature will fall below freezing 760 WHERE(t_bkginc(:,:,:) > 0.0_wp .OR. & 761 & tsn(:,:,:,jp_tem) + t_bkginc(:,:,:) > fzptnz(:,:,:) ) 752 WHERE( t_bkginc(:,:,:) > 0.0_wp .OR. tsn(:,:,:,jp_tem) + t_bkginc(:,:,:) > fzptnz(:,:,:) ) 762 753 tsn(:,:,:,jp_tem) = t_bkg(:,:,:) + t_bkginc(:,:,:) 763 754 END WHERE … … 768 759 ! Do not apply negative increments if the salinity will fall below a specified 769 760 ! minimum value salfixmin 770 WHERE(s_bkginc(:,:,:) > 0.0_wp .OR. & 771 & tsn(:,:,:,jp_sal) + s_bkginc(:,:,:) > salfixmin ) 761 WHERE( s_bkginc(:,:,:) > 0.0_wp .OR. tsn(:,:,:,jp_sal) + s_bkginc(:,:,:) > salfixmin ) 772 762 tsn(:,:,:,jp_sal) = s_bkg(:,:,:) + s_bkginc(:,:,:) 773 763 END WHERE … … 778 768 tsb(:,:,:,:) = tsn(:,:,:,:) ! Update before fields 779 769 780 CALL eos( tsb, rhd, rhop ) ! Before potential and in situ densities770 CALL eos( tsb, fsdept(:,:,:), rhd, rhop ) ! Before potential and in situ densities 781 771 782 772 IF( ln_zps .AND. .NOT. lk_c1d ) & … … 786 776 787 777 #if defined key_zdfkpp 788 CALL eos( tsn, rhd ) ! Compute rhd778 CALL eos( tsn, fsdept(:,:,:), rhd ) ! Compute rhd 789 779 #endif 790 780 -
branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/C1D/step_c1d.F90
r3680 r3894 70 70 ! Ocean physics update (ua, va, ta, sa used as workspace) 71 71 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 72 CALL bn2( tsb, rn2b ) ! before Brunt-Vaisala frequency 73 CALL bn2( tsn, rn2 ) ! now Brunt-Vaisala frequency 72 ! ! thermal/haline expansion coef. & Brunt-Vaisala frequency 73 CALL eos( tsb, rab_b, rn2b ) ! before 74 CALL eos( tsn, rab_n, rn2 ) ! now 75 ! 74 76 ! VERTICAL PHYSICS 75 77 CALL zdf_bfr( kstp ) ! bottom friction … … 125 127 CALL tra_nxt ( kstp ) ! tracer fields at next time step 126 128 IF( ln_zdfnpc ) CALL tra_npc ( kstp ) ! applied non penetrative convective adjustment on (t,s) 127 CALL eos( tsb, rhd, rhop ) ! now (swap=before) in situ density for dynhpg module 129 ! now (swap=before) in situ density for dynhpg module 130 CALL eos( tsb, fsdept(:,:,:), rhd, rhop ) 128 131 129 132 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 130 133 ! Dynamics (ta, sa used as workspace) 131 134 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 132 133 135 ua(:,:,:) = 0.e0 ! set dynamics trends to zero 136 va(:,:,:) = 0.e0 134 137 135 136 137 138 CALL dyn_cor_c1d( kstp ) ! vorticity term including Coriolis 139 CALL dyn_zdf ( kstp ) ! vertical diffusion 140 CALL dyn_nxt_c1d( kstp ) ! lateral velocity at next time step 138 141 139 142 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> … … 141 144 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 142 145 CALL stp_ctl( kstp, indic ) 143 IF( kstp == nit000 144 IF( lrst_oce 146 IF( kstp == nit000 ) CALL iom_close( numror ) ! close input ocean restart file 147 IF( lrst_oce ) CALL rst_write ( kstp ) ! write output ocean restart file 145 148 ! 146 149 END SUBROUTINE stp_c1d -
branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/DIA/diaar5.F90
r3294 r3894 96 96 97 97 ! 98 ztsn(:,:,:,jp_tem) = tsn(:,:,:,jp_tem) 98 ztsn(:,:,:,jp_tem) = tsn(:,:,:,jp_tem) ! thermosteric ssh 99 99 ztsn(:,:,:,jp_sal) = sn0(:,:,:) 100 CALL eos( ztsn, zrhd )! now in situ density using initial salinity101 ! 102 zbotpres(:,:) = 0._wp ! no atmospheric surface pressure, levitating sea-ice100 CALL eos( ztsn, fsdept(:,:,:), zrhd ) ! now in situ density using initial salinity 101 ! 102 zbotpres(:,:) = 0._wp ! no atmospheric surface pressure, levitating sea-ice 103 103 DO jk = 1, jpkm1 104 104 zbotpres(:,:) = zbotpres(:,:) + fse3t(:,:,jk) * zrhd(:,:,jk) … … 112 112 113 113 ! ! steric sea surface height 114 CALL eos( tsn, zrhd, zrhop )! now in situ and potential density114 CALL eos( tsn, fsdept(:,:,:), zrhd, zrhop ) ! now in situ and potential density 115 115 zrhop(:,:,jpk) = 0._wp 116 116 CALL iom_put( 'rhop', zrhop ) 117 117 ! 118 zbotpres(:,:) = 0._wp ! no atmospheric surface pressure, levitating sea-ice118 zbotpres(:,:) = 0._wp ! no atmospheric surface pressure, levitating sea-ice 119 119 DO jk = 1, jpkm1 120 120 zbotpres(:,:) = zbotpres(:,:) + fse3t(:,:,jk) * zrhd(:,:,jk) -
branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/DOM/istate.F90
r3893 r3894 122 122 ENDIF 123 123 ! 124 CALL eos( tsb, rhd, rhop ) ! before potential and in situ densities124 CALL eos( tsb, fsdept(:,:,:), rhd, rhop ) ! before potential and in situ densities 125 125 #if ! defined key_c1d 126 126 IF( ln_zps ) CALL zps_hde( nit000, jpts, tsb, gtsu, gtsv, & ! zps: before hor. gradient -
branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/IOM/restart.F90
r3878 r3894 209 209 CALL iom_get( numror, jpdom_autoglo, 'rhop' , rhop ) ! now potential density 210 210 ELSE 211 CALL eos ( tsn, rhd, rhop )211 CALL eos( tsn, fsdept(:,:,:), rhd, rhop ) 212 212 ENDIF 213 213 #if defined key_zdfkpp … … 215 215 CALL iom_get( numror, jpdom_autoglo, 'rhd' , rhd ) ! now in situ density anomaly 216 216 ELSE 217 CALL eos( tsn, rhd ) ! compute rhd217 CALL eos( tsn, fsdept(:,:,:), rhd ) ! compute rhd 218 218 ENDIF 219 219 #endif -
branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_if.F90
r3625 r3894 98 98 ! ( d rho / dt ) / ( d rho / ds ) ( s = 34, t = -1.8 ) 99 99 100 fr_i(:,:) = tfreez( sss_m ) * tmask(:,:,1) ! sea surface freezing temperature [Celcius]100 fr_i(:,:) = eos_fzp( sss_m ) * tmask(:,:,1) ! sea surface freezing temperature [Celcius] 101 101 #if defined key_coupled 102 102 a_i(:,:,1) = fr_i(:,:) -
branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90
r3625 r3894 122 122 v_oce(:,:) = ssv_m(:,:) ! (C-grid dynamics : U- & V-points as the ocean) 123 123 ! 124 t_bo(:,:) = tfreez( sss_m ) + rt0 ! masked sea surface freezing temperature [Kelvin]124 t_bo(:,:) = eos_fzp( sss_m ) + rt0 ! masked sea surface freezing temperature [Kelvin] 125 125 ! ! (set to rt0 over land) 126 126 CALL albedo_ice( t_su, ht_i, ht_s, zalb_ice_cs, zalb_ice_os ) ! ... ice albedo -
branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim_2.F90
r3680 r3894 140 140 141 141 ! ... masked sea surface freezing temperature [Kelvin] (set to rt0 over land) 142 tfu(:,:) = tfreez( sss_m ) + rt0142 tfu(:,:) = eos_fzp( sss_m ) + rt0 143 143 144 144 zsist (:,:,1) = sist (:,:) + rt0 * ( 1. - tmask(:,:,1) ) -
branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/TRA/eosbn2.F90
r3893 r3894 15 15 !! - ! 2002-11 (G. Madec, A. Bozec) partial step, eos_insitu_2d 16 16 !! - ! 2003-08 (G. Madec) F90, free form 17 !! 3.0 ! 2006-08 (G. Madec) add tfreez function 17 !! 3.0 ! 2006-08 (G. Madec) add tfreez function (now eos_fzp function) 18 18 !! 3.3 ! 2010-05 (C. Ethe, G. Madec) merge TRC-TRA 19 !! - ! 2010-10 (G. Nurser, G. Madec) add eos_abused in ldfslp20 !! 3.5 ! 2012-03 (F. Roquet, G. Madec) add p en_ddt_dds used in trdpen19 !! - ! 2010-10 (G. Nurser, G. Madec) add alpha/beta used in ldfslp 20 !! 3.5 ! 2012-03 (F. Roquet, G. Madec) add primitive of alpha and beta used in PE computation 21 21 !! - ! 2012-05 (F. Roquet) add Vallis and original JM95 equation of state 22 !! - ! 201 2-05 (F. Roquet) add eos_alpha_beta to be used in dynhpg22 !! - ! 2013-04 (F. Roquet, G. Madec) add eos_rab, change bn2 computation and reorganize the module 23 23 !!---------------------------------------------------------------------- 24 24 25 25 !!---------------------------------------------------------------------- 26 26 !! eos : generic interface of the equation of state 27 !! eos_insitu : density anomaly (rho/rho0 - 1)28 !! eos_insitu_pot : density and surface referenced potential density anomalies29 !! eos_insitu_2d : density anomaly for 2d fields30 !! eos_bn2 : Brunt-Vaisala frequency 31 !! eos_ rab : partial derivative of density anomaly with respect to T and S32 !! tfreez : surface freezing temperature27 !! eos_rhd_pot : density anomaly (rho/rho0 - 1) and surface referenced potential density (T-point) 28 !! eos_rab_bn2 : partial derivative of density anomaly with respect to T and S (T-point) 29 !! and Brunt-Vaisala frequency (w-point) 30 !! eos_bn2 : Brunt-Vaisala frequency (w-point) computed from alpha & beta 31 !! eos_fzp : freezing temperature 32 !! eos_pen : primitive of partial derivative of density anomaly with respect to T and S (T-point) 33 33 !! eos_init : set eos parameters (namelist) 34 34 !!---------------------------------------------------------------------- 35 USE oce ! ocean variables 35 36 USE dom_oce ! ocean space and time domain 36 37 USE phycst ! physical constants 38 USE ldfslp ! lateral physics: iso-neutral slopes 37 39 USE zdfddm ! vertical physics: double diffusion 40 ! 38 41 USE in_out_manager ! I/O manager 39 42 USE lib_mpp ! MPP library … … 41 44 USE wrk_nemo ! Memory Allocation 42 45 USE timing ! Timing 46 USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 43 47 44 48 IMPLICIT NONE 45 49 PRIVATE 46 50 47 ! !! * Interface48 51 INTERFACE eos 49 MODULE PROCEDURE eos_insitu, eos_insitu_pot, eos_insitu_2d, & 50 & eos_ab, eos_insitu_ab 52 MODULE PROCEDURE eos_rhd_pot, eos_rab_bn2 51 53 END INTERFACE 52 INTERFACE bn2 53 MODULE PROCEDURE eos_bn2, eos_bn2_ab 54 END INTERFACE 55 56 PUBLIC eos ! called by step, istate, tranpc and zpsgrd modules 57 PUBLIC eos_init ! called by istate module 58 PUBLIC bn2 ! called by step module, update alpha/beta and compute bn2 59 PUBLIC eos_pen ! used for pe diagnostics in trdpen 60 PUBLIC tfreez ! called by sbcice_... modules 61 62 ! !!* Namelist (nameos) * 63 INTEGER , PUBLIC :: nn_eos = 0 !: = -1/0/1 type of eq. of state and Brunt-Vaisala frequ. 64 65 REAL(wp) :: rn_alph0 = 1.67e-4 !: thermal expansion coeff. (simplified equation of state) 66 REAL(wp) :: rn_beta0 = 7.8e-4 !: saline expansion coeff. (simplified equation of state) 67 REAL(wp) :: rn_gamm0 = 4.4e-6 !: depth expansion coeff. (simplified equation of state) 68 REAL(wp) :: rn_cabbel = 3.0e-2 !: cabbeling coeff. (simplified equation of state) 69 REAL(wp) :: rn_thermo = 1.1e-4 !: thermobaric coeff. (simplified equation of state) 70 REAL(wp) :: offset = 0 !: offset applied on Vallis density anomaly, so that rho(10,35,0)=1027 54 55 PUBLIC eos ! called by step, istate, tranpc and zpsgrd modules 56 PUBLIC eos_bn2 ! not use outside eosbn2 module for the moment 57 PUBLIC eos_fzp ! called by sbcice_... modules 58 PUBLIC eos_pen ! used for pe diagnostics in trdpen module 59 PUBLIC eos_init ! called by istate module 60 61 ! !!* Namelist (nameos) * 62 INTEGER , PUBLIC :: nn_eos = 0 !: = -1/0/1 type of equation of state 63 64 ! !!! simplified eos coefficients 65 REAL(wp) :: rn_alph0 = 1.67e-4 ! thermal expansion coeff. 66 REAL(wp) :: rn_beta0 = 7.8e-4 ! saline expansion coeff. 67 REAL(wp) :: rn_gamm0 = 4.4e-6 ! depth expansion coeff. 68 REAL(wp) :: rn_cabbel = 3.0e-2 ! cabbeling coeff. 69 REAL(wp) :: rn_thermo = 1.1e-4 ! thermobaric coeff. 70 REAL(wp) :: offset = 0 ! offset so that rho(10,35,0)=1027kg/m3 71 71 72 72 !! * Substitutions … … 74 74 # include "vectopt_loop_substitute.h90" 75 75 !!---------------------------------------------------------------------- 76 !! NEMO/OPA 3. 3 , NEMO Consortium (2010)76 !! NEMO/OPA 3.5 , NEMO Consortium (2013) 77 77 !! $Id$ 78 78 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) … … 80 80 CONTAINS 81 81 82 SUBROUTINE eos_insitu( pts, prd ) 83 !!---------------------------------------------------------------------- 84 !! *** ROUTINE eos_insitu *** 85 !! 86 !! ** Purpose : Compute the in situ density (ratio rho/rau0) from 87 !! potential temperature and salinity using an equation of state 88 !! defined through the namelist parameter nn_eos. 89 !! 90 !! ** Method : 3 cases: 91 !! nn_eos = -1 : Jackett and McDougall (1995) equation of state. 92 !! Check value: rho = 1041.83267 kg/m**3 for p=3000 dbar, 93 !! t = 3 deg celcius, s=35.5 psu 94 !! nn_eos = 0 : standard NEMO equation of state, based on 95 !! Jackett and McDougall (1995) equation of state. 96 !! the in situ density is computed directly as a function of 97 !! potential temperature relative to the surface (the opa t 98 !! variable), salt and pressure (assuming no pressure variation 99 !! along geopotential surfaces, i.e. the pressure p in decibars 100 !! is approximated by the depth in meters. 101 !! prd(t,s,z) = ( rho(t,s,z) - rau0 ) / rau0 102 !! with depth z meters 103 !! potential temperature t deg celsius 104 !! salinity s psu 105 !! reference volumic mass rau0 kg/m**3 106 !! in situ volumic mass rho kg/m**3 107 !! in situ density anomalie prd no units 108 !! Check value: rho = 1060.93298 kg/m**3 for z=10000 dbar, 109 !! t = 40 deg celcius, s=40 psu 110 !! nn_eos = 1 : Vallis simplified equation of state (Vallis book, 2006) 111 !! prd(t,s,z) = ( rho(t,s,p) - rau0 ) / rau0 112 !! = -alpha*(1+cabbel*(T-T0)+thermo*z)*(T-T0) + beta*(S-S0) + gamma*z 82 SUBROUTINE eos_rhd_pot( pts, pdep, prd, prhop ) 83 !!---------------------------------------------------------------------- 84 !! *** ROUTINE eos_rhd_pot *** 85 !! 86 !! ** Purpose : Compute the density anomaly (rho/rho0-1) from potential 87 !! temperature and salinity and optionally the density referenced 88 !! to the surface using an equation of state defined through 89 !! nn_eos namelist parameter. 90 !! 91 !! ** Method : prd(t,s,z) = ( rho(t,s,z) - rau0 ) / rau0 92 !! with prd in situ density anomalie no units 93 !! t potential temperature Celsius 94 !! s salinity psu 95 !! z depth meters 96 !! rho in situ density kg/m^3 97 !! rau0 reference density kg/m^3 98 !! 99 !! nn_eos = -1 : Jackett and McDougall (1995) equation of state is used for rho(t,s,z). 100 !! Check value: rho = 1041.83267 kg/m^3 for p=3000 dbar, t=3 Celcius, s=35.5 psu 101 !! 102 !! nn_eos = 0 : standard NEMO equation of state, based also on Jackett and McDougall (1995) 103 !! but with slightly different coefficients. 104 !! Check value: rho = 1060.93298 kg/m^3 for z=10000 dbar, t=40 Celcius, s=40 psu 105 !! 106 !! nn_eos = 1 : simplified equation of state (inspired Vallis 2006) 107 !! prd(t,s,z) = -alpha*(1+cabbel*(T-T0)+thermo*z)*(T-T0) + beta*(S-S0) + gamma*z 113 108 !! linear case function of T only: rn_alpha<>0, other coefficients = 0 114 109 !! linear eos function of T and S: rn_alpha and rn_beta<>0, other coefficients=0 115 !! Vallis (2006) equation: use default values of coefficients 116 !! Note that no boundary condition problem occurs in this routine 117 !! as pts are defined over the whole domain. 118 !! 119 !! ** Action : compute prd , the in situ density (no units) 110 !! Vallis like equation: use default values of coefficients 111 !! 112 !! ** Action : prd density anomaly (no units) 120 113 !! 121 114 !! References : Jackett and McDougall, J. Atmos. Ocean. Tech., 1995 122 115 !! Vallis, Atmospheric and Oceanic Fluid Dynamics, 2006 123 116 !!---------------------------------------------------------------------- 124 !!125 REAL(wp), DIMENSION(:,:,: ,:), INTENT(in ) :: pts ! 1 : potential temperature [Celcius]126 ! ! 2 : salinity [psu]127 REAL(wp), DIMENSION(:,:,:) , INTENT( out) :: prd ! in situ density [-]128 ! !117 REAL(wp), DIMENSION(:,:,:,:), INTENT(in ) :: pts ! pot. temperature and salinity [Celcius,psu] 118 REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: pdep ! depth of pts data [m] 119 REAL(wp), DIMENSION(:,:,:) , INTENT( out) :: prd ! density anomaly [-] 120 REAL(wp), DIMENSION(:,:,:) , INTENT( out), OPTIONAL :: prhop ! potential density (ref. to z=0) [kg/m3] 121 ! 129 122 INTEGER :: ji, jj, jk ! dummy loop indices 123 INTEGER :: ikm1 ! local integer 130 124 REAL(wp) :: zt , zs , zh , zsr ! local scalars 131 125 REAL(wp) :: zn1, zn2, zn3, zn4 ! - - … … 133 127 REAL(wp) :: zb1, zb2, zb3, zb ! - - 134 128 REAL(wp) :: zc1, zc2, zc3, zc ! - - 135 REAL(wp), POINTER, DIMENSION(:,:,:) :: zws136 ! !----------------------------------------------------------------------137 138 ! 139 IF( nn_timing == 1 ) CALL timing_start('eos')129 !!---------------------------------------------------------------------- 130 ! 131 IF( nn_timing == 1 ) CALL timing_start('eos_rhd_pot') 132 ! 133 ikm1 = MAX( SIZE( pts, 3 ) - 1 , 1 ) ! =1 for 2D calculation, jpkm1 otherwise 140 134 ! 141 135 SELECT CASE( nn_eos ) … … 143 137 CASE( -1 ) !== Jackett and McDougall (1995) formulation ==! 144 138 ! 145 CALL wrk_alloc( jpi, jpj, jpk, zws ) 146 ! 147 !CDIR NOVERRCHK 148 zws(:,:,:) = SQRT( ABS( pts(:,:,:,jp_sal) ) ) 149 ! 150 DO jk = 1, jpkm1 139 DO jk = 1, ikm1 151 140 DO jj = 1, jpj 152 141 DO ji = 1, jpi 153 zt = pts(ji,jj,jk,jp_tem)154 zs = pts(ji,jj,jk,jp_sal)155 zh = fsdept(ji,jj,jk)! depth156 zsr = zws (ji,jj,jk)! square root salinity142 zt = pts (ji,jj,jk,jp_tem) 143 zs = pts (ji,jj,jk,jp_sal) 144 zh = pdep(ji,jj,jk) ! depth 145 zsr = SQRT( ABS( pts(ji,jj,jk,jp_sal) ) ) ! square root salinity 157 146 ! 158 147 ! compute volumic mass pure water at atm pressure … … 182 171 ! 183 172 ! masked in situ density anomaly. Important: rau0=1035, should be 1027! 184 prd(ji,jj,jk) = ( zn * ( 1.0_wp + zh / ( ( za * zh + zb ) * zh + zc ) ) / rau0 & 185 & - 1._wp ) * tmask(ji,jj,jk) 173 prd(ji,jj,jk) = ( zn * ( 1.0_wp + zh / ( ( za * zh + zb ) * zh + zc ) ) * r1_rau0 - 1._wp ) * tmask(ji,jj,jk) 174 ! 175 IF( PRESENT(prhop) ) THEN ! potential density referenced at the surface 176 prhop(ji,jj,jk) = zn * tmask(ji,jj,jk) 177 ENDIF 186 178 END DO 187 179 END DO 188 180 END DO 189 181 ! 190 CALL wrk_dealloc( jpi, jpj, jpk, zws )191 !192 182 CASE( 0 ) !== modified Jackett and McDougall (1995) formulation ==! 193 183 ! 194 CALL wrk_alloc( jpi, jpj, jpk, zws ) 195 ! 196 !CDIR NOVERRCHK 197 zws(:,:,:) = SQRT( ABS( pts(:,:,:,jp_sal) ) ) 198 ! 199 DO jk = 1, jpkm1 184 DO jk = 1, ikm1 200 185 DO jj = 1, jpj 201 186 DO ji = 1, jpi 202 187 zt = pts (ji,jj,jk,jp_tem) 203 188 zs = pts (ji,jj,jk,jp_sal) 204 zh = fsdept(ji,jj,jk) ! depth205 zsr= zws (ji,jj,jk)! square root salinity189 zh = fsdept(ji,jj,jk) ! depth 190 zsr= SQRT( ABS( pts(ji,jj,jk,jp_sal) ) ) ! square root salinity 206 191 ! 207 192 ! compute volumic mass pure water at atm pressure … … 230 215 zc = ( zc3*zsr + zc2 )*zs + zc1 231 216 ! 232 ! masked in situ density anomaly 233 prd(ji,jj,jk) = ( zn * ( 1.0_wp + zh / ( ( za * zh + zb ) * zh + zc ) ) / rau0 & 234 & - 1._wp ) * tmask(ji,jj,jk) 217 ! ! density anomaly (masked) 218 prd(ji,jj,jk) = ( zn * ( 1.0_wp + zh / ( ( za * zh + zb ) * zh + zc ) ) * r1_rau0 - 1._wp ) * tmask(ji,jj,jk) 219 ! 220 IF( PRESENT(prhop) ) THEN ! potential density referenced at the surface 221 prhop(ji,jj,jk) = zn * tmask(ji,jj,jk) 222 ENDIF 235 223 END DO 236 224 END DO 237 225 END DO 238 226 ! 239 CALL wrk_dealloc( jpi, jpj, jpk, zws ) 240 ! 241 CASE( 1 ) !== Vallis (2006) simplified EOS ==! 242 DO jk = 1, jpkm1 227 CASE( 1 ) !== simplified EOS ==! 228 DO jk = 1, ikm1 243 229 DO jj = 1, jpj 244 230 DO ji = 1, jpi 245 zt = pts (ji,jj,jk,jp_tem) - 10._wp 246 zs = pts (ji,jj,jk,jp_sal) - 35._wp 247 zh = fsdept(ji,jj,jk) ! depth 248 ! masked in situ density anomaly 249 prd(ji,jj,jk) = ( offset - rn_alph0 * ( 1._wp + rn_cabbel*zt + rn_thermo*zh ) *zt & & 250 & + rn_beta0 * zs + rn_gamm0 * zh ) * tmask(ji,jj,jk) 231 zt = pts (ji,jj,jk,jp_tem) - 10._wp 232 zs = pts (ji,jj,jk,jp_sal) - 35._wp 233 zh = pdep(ji,jj,jk) 234 ! ! density anomaly (masked) 235 prd(ji,jj,jk) = ( offset - rn_alph0 * ( 1._wp + rn_cabbel*zt + rn_thermo*zh ) * zt & 236 & + rn_beta0 * zs + rn_gamm0 * zh ) * tmask(ji,jj,jk) 237 ! 238 IF( PRESENT(prhop) ) THEN ! potential density referenced at z=0 239 prhop(ji,jj,jk) = rau0 * ( 1._wp & 240 & + offset - rn_alph0 * ( 1._wp + rn_cabbel*zt ) * zt & 241 & + rn_beta0 * zs ) * tmask(ji,jj,jk) 242 ENDIF 251 243 END DO 252 244 END DO … … 255 247 END SELECT 256 248 ! 257 IF(ln_ctl) CALL prt_ctl( tab3d_1=prd, clinfo1=' eos : ', ovlap=1, kdim=jpk ) 258 ! 259 IF( nn_timing == 1 ) CALL timing_stop('eos') 260 ! 261 END SUBROUTINE eos_insitu 262 263 264 SUBROUTINE eos_insitu_pot( pts, prd, prhop ) 265 !!---------------------------------------------------------------------- 266 !! *** ROUTINE eos_insitu_pot *** 267 !! 268 !! ** Purpose : Compute the in situ density (ratio (rho-rau0)/rau0) and the 269 !! potential volumic mass (Kg/m3) from potential temperature and 270 !! salinity fields using an equation of state defined through the 271 !! namelist parameter nn_eos. 272 !! 273 !! ** Method : Same as for eos_insitu 274 !! 275 !! ** Action : - prd , the in situ density (no units) 276 !! - prhop, the potential volumic mass (Kg/m3) 277 !! 278 !! References : Vallis, Atmospheric and Oceanic Fluid Dynamics, 2006 279 !! Jackett and McDougall, J. Atmos. Ocean. Tech., 1995 280 !! Brown and Campana, Mon. Weather Rev., 1978 281 !!---------------------------------------------------------------------- 282 !! 283 REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in ) :: pts ! 1 : potential temperature [Celcius] 284 ! ! 2 : salinity [psu] 285 REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT( out) :: prd ! in situ density [-] 286 REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT( out) :: prhop ! potential density (surface referenced) 249 IF(ln_ctl) CALL prt_ctl( tab3d_1=prd, clinfo1=' eos_rhd_pot : ', ovlap=1, kdim=jpk ) 250 ! 251 IF( nn_timing == 1 ) CALL timing_stop('eos_rhd_pot') 252 ! 253 END SUBROUTINE eos_rhd_pot 254 255 256 SUBROUTINE eos_rab_bn2( pts, pab, pn2 ) 257 !!---------------------------------------------------------------------- 258 !! *** ROUTINE eos_ab *** 259 !! 260 !! ** Purpose : Calculates the thermal/haline expansion terms at T-points 261 !! and optionally the Brunt-Vaisala frequency square 262 !! 263 !! ** Method : - calculates alpha and beta at T-points 264 !! * nn_eos = -1 : Jackett and McDougall (1995) equation of state 265 !! * nn_eos = 0 : modified Jackett and McDougall (1995) equation of state 266 !! * nn_eos = 1 : Vallis equation of state (Vallis 2006) 267 !! alpha and beta ratios are returned as 3-D arrays. 268 !! - calculates N^2 through a call to eos_bn2 (if pn2 argument present) 269 !! 270 !! ** Action : pab partial derivative of density anomaly with respect to T & S at T-points 271 !! pn2 square of the Brunt-Vaisala frequency at w-point [1/s2] 272 !!---------------------------------------------------------------------- 273 REAL(wp), DIMENSION(:,:,:,:), INTENT(in ) :: pts ! pot. temperature & salinity [Celcius,psu] 274 REAL(wp), DIMENSION(:,:,:,:), INTENT( out) :: pab ! alpha and beta [1/Celcius, 1/psu] 275 REAL(wp), DIMENSION(:,:,:) , INTENT( out), OPTIONAL :: pn2 ! Brunt-Vaisala frequency**2 [1/s2] 287 276 ! 288 277 INTEGER :: ji, jj, jk ! dummy loop indices … … 292 281 REAL(wp) :: zb1, zb2, zb3, zb ! - - 293 282 REAL(wp) :: zc1, zc2, zc3, zc ! - - 294 REAL(wp), POINTER, DIMENSION(:,:,:) :: zws 295 !!---------------------------------------------------------------------- 296 ! 297 IF( nn_timing == 1 ) CALL timing_start('eos-p') 298 ! 299 CALL wrk_alloc( jpi, jpj, jpk, zws ) 300 ! 301 SELECT CASE ( nn_eos ) 302 ! 303 CASE( -1 ) !== Jackett and McDougall (1995) formulation ==! 304 !CDIR NOVERRCHK 305 zws(:,:,:) = SQRT( ABS( pts(:,:,:,jp_sal) ) ) 283 REAL(wp) :: zm_inv, zdm, zdn ! - - 284 REAL(wp) :: zdza, zdzb, zdzc ! - - 285 !!---------------------------------------------------------------------- 286 ! 287 IF ( nn_timing == 1 ) CALL timing_start('eos_rab_bn2') 288 ! 289 SELECT CASE ( nn_eos ) !++ partial derivative of density anomaly with respect to T & S ++! 290 ! 291 CASE ( -1 ) !== Jackett and McDougall (1995) formulation ==! 306 292 ! 307 293 DO jk = 1, jpkm1 … … 310 296 zt = pts (ji,jj,jk,jp_tem) 311 297 zs = pts (ji,jj,jk,jp_sal) 312 zh = fsdept(ji,jj,jk) ! depth313 zsr= zws (ji,jj,jk)! square root salinity298 zh = fsdept(ji,jj,jk) ! depth 299 zsr= SQRT( ABS( pts(ji,jj,jk,jp_sal) ) ) ! square root salinity 314 300 ! 315 301 ! compute volumic mass pure water at atm pressure … … 333 319 zb = ( zb3*zsr + zb2 ) *zs + zb1 334 320 ! 335 zc3 = (-4.619924e-3_wp*zt+9.085835e-2_wp ) *zt+3.886640_wp336 zc2 = ( ( -5.084188e-4_wp*zt+6.283263e-2_wp) *zt-3.101089_wp ) *zt+528.4855_wp337 321 zc1 = ( ( (-4.190253e-4_wp*zt+9.648704e-2_wp)*zt-17.06103_wp )*zt+1444.304_wp)*zt+196593.3_wp 322 zc2 = ( ( -5.084188e-4_wp*zt+6.283263e-2_wp) *zt-3.101089_wp ) *zt+528.4855_wp 323 zc3 = (-4.619924e-3_wp*zt+9.085835e-2_wp ) *zt+3.886640_wp 338 324 zc = ( zc3*zsr + zc2 )*zs + zc1 339 325 ! 340 ! masked in situ density anomaly. Important: rau0=1035, should be 1027! 341 prd(ji,jj,jk) = ( zn * ( 1.0_wp + zh / ( ( za * zh + zb ) * zh + zc ) ) / rau0 & 342 & - 1._wp ) * tmask(ji,jj,jk) 343 prhop(ji,jj,jk) = zn * tmask(ji,jj,jk) 326 zm_inv = 1._wp / ( ( za*zh + zb )*zh + zc ) 327 ! ! alpha 328 zdza = (2.78936e-8_wp*zt-1.202016e-6_wp) + (12.414646e-11_wp*zt+6.128773e-9_wp ) * zs 329 zdzb = (4.118662e-7_wp*zt-1.847318e-4_wp)*zs + & 330 & ( 5.869245e-6_wp*zt-5.969284e-4_wp)*zt+2.212276e-2_wp 331 zdzc = ( (-9.239848e-3_wp*zt+9.085835e-2_wp)*zsr + & 332 & (-15.252564e-4_wp*zt+12.566526e-2_wp)*zt-3.101089_wp ) * zs + & 333 & ((-16.761012e-4_wp*zt+28.946112e-2_wp)*zt-34.12206_wp )*zt+1444.304_wp 334 zdn = ((-3.3092e-6_wp*zt+1.0227e-4_wp)*zsr & 335 & +(((21.55e-9_wp*zt-24.7401e-7_wp)*zt+15.2876e-5_wp)*zt-4.0899e-3_wp))*zs & 336 & +((((32.68166e-9_wp*zt-4.480332e-6_wp)*zt+3.005055e-4_wp)*zt & 337 & -18.19058e-3_wp)*zt+6.793952e-2_wp) 338 zdm = (zdza*zh+zdzb)*zh+zdzc 339 ! 340 pab(ji,jj,jk,jp_tem) = - ( zdn * ( 1._wp + zh * zm_inv ) - zn * zh * zdm * zm_inv * zm_inv ) & 341 & * r1_rau0 * tmask(ji,jj,jk) 342 ! 343 zdza = za2 ! beta 344 zdzb = 2.220399e-4_wp * zsr + zb2 345 zdzc = 1.5_wp * zc3 * zsr + zc2 346 zdn = 9.6628e-4_wp * zs + 1.5_wp * zn3 * zsr + zn2 347 zdm = ( zdza * zh + zdzb ) * zh + zdzc 348 pab(ji,jj,jk,jp_sal) = ( zdn * ( 1._wp + zh * zm_inv ) - zn * zh * zdm * zm_inv * zm_inv ) & 349 & * r1_rau0 * tmask(ji,jj,jk) 344 350 END DO 345 351 END DO 346 352 END DO 347 353 ! 348 CASE( 0 ) !== modified Jackett and McDougall (1995) formulation ==! 349 !CDIR NOVERRCHK 350 zws(:,:,:) = SQRT( ABS( pts(:,:,:,jp_sal) ) ) 354 CASE ( 0 ) !== modified Jackett and McDougall (1995) formulation ==! 351 355 ! 352 356 DO jk = 1, jpkm1 … … 355 359 zt = pts (ji,jj,jk,jp_tem) 356 360 zs = pts (ji,jj,jk,jp_sal) 357 zh = fsdept(ji,jj,jk) ! depth358 zsr= zws (ji,jj,jk)! square root salinity361 zh = fsdept(ji,jj,jk) ! depth 362 zsr= SQRT( ABS( pts(ji,jj,jk,jp_sal) ) ) ! square root salinity 359 363 ! 360 364 ! compute volumic mass pure water at atm pressure … … 383 387 zc = ( zc3*zsr + zc2 )*zs + zc1 384 388 ! 385 ! masked in situ density anomaly386 prd(ji,jj,jk) = ( zn * ( 1.0_wp + zh / ( ( za * zh + zb ) * zh + zc ) ) / rau0 &387 & - 1._wp ) * tmask(ji,jj,jk)388 prhop(ji,jj,jk) = zn * tmask(ji,jj,jk)389 END DO390 END DO391 END DO392 !393 CASE( 1 ) !== Vallis (2006) simplified EOS ==!394 DO jk = 1, jpkm1395 DO jj = 1, jpj396 DO ji = 1, jpi397 zt = pts (ji,jj,jk,jp_tem) - 10._wp398 zs = pts (ji,jj,jk,jp_sal) - 35._wp399 zh = fsdept(ji,jj,jk) ! depth400 ! masked in situ density anomaly401 prd(ji,jj,jk) = ( offset - rn_alph0 * ( 1._wp + rn_cabbel*zt + rn_thermo*zh ) *zt & &402 & + rn_beta0 * zs + rn_gamm0 * zh ) * tmask(ji,jj,jk)403 ! masked in situ density anomaly404 prhop(ji,jj,jk) = ( 1.0_wp + offset - rn_alph0 * ( 1._wp + rn_cabbel*zt ) *zt + rn_beta0 *zs ) &405 & * rau0 * tmask(ji,jj,jk)406 !407 END DO408 END DO409 END DO410 !411 END SELECT412 !413 IF(ln_ctl) CALL prt_ctl( tab3d_1=prd, clinfo1=' eos-p: ', tab3d_2=prhop, clinfo2=' pot : ', ovlap=1, kdim=jpk )414 !415 CALL wrk_dealloc( jpi, jpj, jpk, zws )416 !417 IF( nn_timing == 1 ) CALL timing_stop('eos-p')418 !419 END SUBROUTINE eos_insitu_pot420 421 422 SUBROUTINE eos_insitu_2d( pts, pdep, prd )423 !!----------------------------------------------------------------------424 !! *** ROUTINE eos_insitu_2d ***425 !!426 !! ** Purpose : Compute the in situ density (ratio rho/rau0) from427 !! potential temperature and salinity using an equation of state428 !! defined through the namelist parameter nn_eos. * 2D field case429 !!430 !! ** Method : Same as for eos_insitu431 !!432 !! ** Action : - prd , the in situ density (no units)433 !!434 !! References : Vallis, Atmospheric and Oceanic Fluid Dynamics, 2006435 !! Jackett and McDougall, J. Atmos. Ocean. Tech., 1995436 !!----------------------------------------------------------------------437 !!438 REAL(wp), DIMENSION(jpi,jpj,jpts), INTENT(in ) :: pts ! 1 : potential temperature [Celcius]439 ! ! 2 : salinity [psu]440 REAL(wp), DIMENSION(jpi,jpj) , INTENT(in ) :: pdep ! depth [m]441 REAL(wp), DIMENSION(jpi,jpj) , INTENT( out) :: prd ! in situ density442 !!443 INTEGER :: ji, jj ! dummy loop indices444 REAL(wp) :: zt , zs , zh , zsr ! local scalars445 REAL(wp) :: zn1, zn2, zn3, zn4 ! - -446 REAL(wp) :: zn, za1, za2, za ! - -447 REAL(wp) :: zb1, zb2, zb3, zb ! - -448 REAL(wp) :: zc1, zc2, zc3, zc ! - -449 REAL(wp), POINTER, DIMENSION(:,:) :: zws450 !!----------------------------------------------------------------------451 !452 IF( nn_timing == 1 ) CALL timing_start('eos2d')453 !454 CALL wrk_alloc( jpi, jpj, zws )455 !456 457 prd(:,:) = 0._wp458 459 SELECT CASE( nn_eos )460 !461 CASE( -1 ) !== Jackett and McDougall (1995) formulation ==!462 !463 !CDIR NOVERRCHK464 DO jj = 1, jpjm1465 !CDIR NOVERRCHK466 DO ji = 1, fs_jpim1 ! vector opt.467 zws(ji,jj) = SQRT( ABS( pts(ji,jj,jp_sal) ) )468 END DO469 END DO470 DO jj = 1, jpjm1471 DO ji = 1, fs_jpim1 ! vector opt.472 zt = pts (ji,jj,jp_tem) ! interpolated T473 zs = pts (ji,jj,jp_sal) ! interpolated S474 zsr = zws (ji,jj) ! square root of interpolated S475 zh = pdep (ji,jj) ! depth at the partial step level476 !477 ! compute volumic mass pure water at atm pressure478 zn1= ( ( ( ( 6.536332e-9_wp *zt - 1.120083e-6_wp )*zt + 1.001685e-4_wp )*zt &479 & -9.095290e-3_wp )*zt + 6.793952e-2_wp )*zt + 999.842594_wp480 ! seawater volumic mass atm pressure481 zn2= ( ( ( 5.3875e-9_wp*zt-8.2467e-7_wp ) *zt+7.6438e-5_wp ) *zt &482 & -4.0899e-3_wp ) *zt+0.824493_wp483 zn3= ( -1.6546e-6_wp*zt+1.0227e-4_wp ) *zt-5.72466e-3_wp484 zn4= 4.8314e-4_wp485 zn = ( zn4*zs + zn3*zsr + zn2 ) *zs + zn1486 !487 ! add the compression terms488 za1= ( 1.394680e-8_wp*zt-1.202016e-6_wp ) *zt+2.102898e-5_wp489 za2 = ( 6.207323e-11_wp*zt+6.128773e-9_wp ) *zt-2.040237e-7_wp490 za = za1 + za2 * zs491 !492 zb1 = ( ( 1.956415e-6_wp*zt-2.984642e-4_wp ) *zt+2.212276e-2_wp ) *zt +2.186519_wp493 zb2 = ( 2.059331e-7_wp*zt-1.847318e-4_wp ) *zt+6.704388e-3494 zb3 = 1.480266e-4_wp495 zb = ( zb3*zsr + zb2 ) *zs + zb1496 !497 zc3 = (-4.619924e-3_wp*zt+9.085835e-2_wp ) *zt+3.886640_wp498 zc2 = ( ( -5.084188e-4_wp*zt+6.283263e-2_wp) *zt-3.101089_wp ) *zt+528.4855_wp499 zc1 = ( ( (-4.190253e-4_wp*zt+9.648704e-2_wp)*zt-17.06103_wp )*zt+1444.304_wp)*zt+196593.3_wp500 zc = ( zc3*zsr + zc2 )*zs + zc1501 !502 ! masked in situ density anomaly503 prd(ji,jj) = ( zn * ( 1.0_wp + zh / ( ( za * zh + zb ) * zh + zc ) ) / rau0 &504 & - 1._wp ) * tmask(ji,jj,1)505 END DO506 END DO507 !508 CASE( 0 ) !== modified Jackett and McDougall (1995) formulation ==!509 !510 !CDIR NOVERRCHK511 DO jj = 1, jpjm1512 !CDIR NOVERRCHK513 DO ji = 1, fs_jpim1 ! vector opt.514 zws(ji,jj) = SQRT( ABS( pts(ji,jj,jp_sal) ) )515 END DO516 END DO517 DO jj = 1, jpjm1518 DO ji = 1, fs_jpim1 ! vector opt.519 zt = pts (ji,jj,jp_tem) ! interpolated T520 zs = pts (ji,jj,jp_sal) ! interpolated S521 zsr = zws (ji,jj) ! square root of interpolated S522 zh = pdep (ji,jj) ! depth at the partial step level523 !524 ! compute volumic mass pure water at atm pressure525 zn1= ( ( ( ( 6.536332e-9_wp *zt - 1.120083e-6_wp )*zt + 1.001685e-4_wp )*zt &526 & -9.095290e-3_wp )*zt + 6.793952e-2_wp )*zt + 999.842594_wp527 ! seawater volumic mass atm pressure528 zn2= ( ( ( 5.3875e-9_wp*zt-8.2467e-7_wp ) *zt+7.6438e-5_wp ) *zt &529 & -4.0899e-3_wp ) *zt+0.824493_wp530 zn3= ( -1.6546e-6_wp*zt+1.0227e-4_wp ) *zt-5.72466e-3_wp531 zn4= 4.8314e-4_wp532 zn= ( zn4*zs + zn3*zsr + zn2 ) *zs + zn1533 !534 ! add the compression terms535 za2= ( -3.508914e-8_wp*zt-1.248266e-8_wp ) *zt-2.595994e-6_wp536 za1= ( 1.296821e-6_wp*zt-5.782165e-9_wp ) *zt+1.045941e-4_wp537 za = za1 + za2 * zs538 !539 zb1= ( (-5.939910e-6_wp*zt-2.512549e-3_wp ) *zt+0.1028859_wp ) *zt + 3.721788_wp540 zb2= (7.267926e-5_wp*zt-2.598241e-3_wp ) *zt-0.1571896_wp541 zb3= 2.042967e-2_wp542 zb = ( zb3*zsr + zb2 ) *zs + zb1543 !544 zc1= ( ( (-1.361629e-4_wp*zt-1.852732e-2_wp ) *zt-30.41638_wp ) *zt + 2098.925_wp ) *zt+190925.6_wp545 zc2= ( ( 2.326469e-3_wp*zt+1.553190_wp) *zt-65.00517_wp ) *zt+1044.077_wp546 zc3= (-0.1909078_wp*zt+7.390729_wp ) *zt-55.87545_wp547 zc = ( zc3*zsr + zc2 )*zs + zc1548 !549 ! masked in situ density anomaly550 prd(ji,jj) = ( zn * ( 1.0_wp + zh / ( ( za * zh + zb ) * zh + zc ) ) / rau0 &551 & - 1._wp ) * tmask(ji,jj,1)552 END DO553 END DO554 !555 CASE( 1 ) !== Vallis (2006) simplified EOS ==!556 DO jj = 1, jpjm1557 DO ji = 1, fs_jpim1 ! vector opt.558 zt = pts (ji,jj,jp_tem) - 10._wp ! interpolated T559 zs = pts (ji,jj,jp_sal) - 35._wp ! interpolated S560 zh = pdep (ji,jj) ! depth at the partial step level561 ! masked in situ density anomaly562 prd(ji,jj) = ( offset - rn_alph0 * ( 1._wp + rn_cabbel*zt + rn_thermo*zh ) *zt & &563 & + rn_beta0 * zs + rn_gamm0 * zh ) * tmask(ji,jj,1)564 !565 END DO566 END DO567 !568 END SELECT569 570 IF(ln_ctl) CALL prt_ctl( tab2d_1=prd, clinfo1=' eos2d: ' )571 !572 CALL wrk_dealloc( jpi, jpj, zws )573 !574 IF( nn_timing == 1 ) CALL timing_stop('eos2d')575 !576 END SUBROUTINE eos_insitu_2d577 578 579 580 SUBROUTINE eos_ab( pts, pab )581 !!----------------------------------------------------------------------582 !! *** ROUTINE eos_ab ***583 !!584 !! ** Purpose : Calculates the in situ thermal/haline expansion terms at T-points585 !!586 !! ** Method : calculates alpha and beta at T-points587 !! * nn_eos = -1 : Jackett and McDougall (1995) equation of state588 !! * nn_eos = 0 : modified Jackett and McDougall (1995) equation of state589 !! * nn_eos = 1 : Vallis equation of state (Vallis 2006)590 !! alpha and beta ratios are returned as 3-D arrays.591 !!592 !! ** Action : - pab : partial derivative of in situ density with respect to T & S at T-points593 !!----------------------------------------------------------------------594 REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in ) :: pts ! pot. temperature & salinity595 REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT( out) :: pab ! alpha and beta596 !597 INTEGER :: ji, jj, jk ! dummy loop indices598 REAL(wp) :: zt , zs , zh , zsr ! local scalars599 REAL(wp) :: zn1, zn2, zn3, zn4 ! - -600 REAL(wp) :: zn, za1, za2, za ! - -601 REAL(wp) :: zb1, zb2, zb3, zb ! - -602 REAL(wp) :: zc1, zc2, zc3, zc ! - -603 REAL(wp) :: zm_inv, zdm, zdn ! - -604 REAL(wp) :: zdza, zdzb, zdzc ! - -605 REAL(wp), POINTER, DIMENSION(:,:,:) :: zws606 !!----------------------------------------------------------------------607 !608 IF ( nn_timing == 1 ) CALL timing_start('eos_ab')609 !610 CALL wrk_alloc( jpi, jpj, jpk, zws )611 !612 SELECT CASE ( nn_eos )613 !614 CASE ( -1 ) ! Jackett and McDougall (1995) formulation615 !616 zws(:,:,:) = SQRT( ABS( pts(:,:,:,jp_sal) ) )617 DO jk = 1, jpkm1618 DO jj = 1, jpj619 DO ji = 1, jpi620 zt = pts (ji,jj,jk,jp_tem)621 zs = pts (ji,jj,jk,jp_sal)622 zh = fsdept(ji,jj,jk) ! depth623 zsr= zws (ji,jj,jk) ! square root salinity624 !625 ! compute volumic mass pure water at atm pressure626 zn1= ( ( ( ( 6.536332e-9_wp *zt - 1.120083e-6_wp )*zt + 1.001685e-4_wp )*zt &627 & -9.095290e-3_wp )*zt + 6.793952e-2_wp )*zt + 999.842594_wp628 ! seawater volumic mass atm pressure629 zn2= ( ( ( 5.3875e-9_wp*zt-8.2467e-7_wp ) *zt+7.6438e-5_wp ) *zt &630 & -4.0899e-3_wp ) *zt+0.824493_wp631 zn3= ( -1.6546e-6_wp*zt+1.0227e-4_wp ) *zt-5.72466e-3_wp632 zn4= 4.8314e-4_wp633 zn = ( zn4*zs + zn3*zsr + zn2 ) *zs + zn1634 !635 ! add the compression terms636 za1= ( 1.394680e-8_wp*zt-1.202016e-6_wp ) *zt+2.102898e-5_wp637 za2 = ( 6.207323e-11_wp*zt+6.128773e-9_wp ) *zt-2.040237e-7_wp638 za = za1 + za2 * zs639 !640 zb1 = ( ( 1.956415e-6_wp*zt-2.984642e-4_wp ) *zt+2.212276e-2_wp ) *zt +2.186519_wp641 zb2 = ( 2.059331e-7_wp*zt-1.847318e-4_wp ) *zt+6.704388e-3642 zb3 = 1.480266e-4_wp643 zb = ( zb3*zsr + zb2 ) *zs + zb1644 !645 zc1 = ( ( (-4.190253e-4_wp*zt+9.648704e-2_wp)*zt-17.06103_wp )*zt+1444.304_wp)*zt+196593.3_wp646 zc2 = ( ( -5.084188e-4_wp*zt+6.283263e-2_wp) *zt-3.101089_wp ) *zt+528.4855_wp647 zc3 = (-4.619924e-3_wp*zt+9.085835e-2_wp ) *zt+3.886640_wp648 zc = ( zc3*zsr + zc2 )*zs + zc1649 !650 389 zm_inv = 1._wp / ( ( za*zh + zb )*zh + zc ) 651 ! alpha 652 zdza = (2.78936e-8_wp*zt-1.202016e-6_wp) + (12.414646e-11_wp*zt+6.128773e-9_wp ) * zs 653 zdzb = (4.118662e-7_wp*zt-1.847318e-4_wp)*zs + & 654 & ( 5.869245e-6_wp*zt-5.969284e-4_wp)*zt+2.212276e-2_wp 655 zdzc = ( (-9.239848e-3_wp*zt+9.085835e-2_wp)*zsr + & 656 & (-15.252564e-4_wp*zt+12.566526e-2_wp)*zt-3.101089_wp ) * zs + & 657 & ((-16.761012e-4_wp*zt+28.946112e-2_wp)*zt-34.12206_wp )*zt+1444.304_wp 658 659 zdn = ((-3.3092e-6_wp*zt+1.0227e-4_wp)*zsr & 660 & +(((21.55e-9_wp*zt-24.7401e-7_wp)*zt+15.2876e-5_wp)*zt-4.0899e-3_wp))*zs & 661 & +((((32.68166e-9_wp*zt-4.480332e-6_wp)*zt+3.005055e-4_wp)*zt & 662 & -18.19058e-3_wp)*zt+6.793952e-2_wp) 663 zdm = (zdza*zh+zdzb)*zh+zdzc 664 ! 665 pab(ji,jj,jk,jp_tem) = - ( & 666 & zdn * ( 1._wp + zh * zm_inv ) - zn * zh * zdm * zm_inv * zm_inv ) & 667 & / rau0 * tmask(ji,jj,jk) 668 ! 669 ! beta 670 zdza = za2 671 zdzb = 2.220399e-4_wp*zsr+zb2 672 zdzc = 1.5_wp*zc3*zsr+zc2 673 zdn = 9.6628e-4_wp*zs+1.5_wp*zn3*zsr+zn2 674 zdm = (zdza*zh+zdzb)*zh+zdzc 675 pab(ji,jj,jk,jp_sal) = ( & 676 & zdn * ( 1._wp + zh * zm_inv ) - zn * zh * zdm * zm_inv * zm_inv ) & 677 & / rau0 * tmask(ji,jj,jk) 678 ! 679 END DO 680 END DO 681 END DO 682 ! 683 CASE ( 0 ) ! modified Jackett and McDougall (1995) formulation 684 ! 685 zws(:,:,:) = SQRT( ABS( pts(:,:,:,jp_sal) ) ) 686 DO jk = 1, jpkm1 687 DO jj = 1, jpj 688 DO ji = 1, jpi 689 zt = pts (ji,jj,jk,jp_tem) 690 zs = pts (ji,jj,jk,jp_sal) 691 zh = fsdept(ji,jj,jk) ! depth 692 zsr= zws (ji,jj,jk) ! square root salinity 693 ! 694 ! compute volumic mass pure water at atm pressure 695 zn1= ( ( ( ( 6.536332e-9_wp *zt - 1.120083e-6_wp )*zt + 1.001685e-4_wp )*zt & 696 & -9.095290e-3_wp )*zt + 6.793952e-2_wp )*zt + 999.842594_wp 697 ! seawater volumic mass atm pressure 698 zn2= ( ( ( 5.3875e-9_wp*zt-8.2467e-7_wp ) *zt+7.6438e-5_wp ) *zt & 699 & -4.0899e-3_wp ) *zt+0.824493_wp 700 zn3= ( -1.6546e-6_wp*zt+1.0227e-4_wp ) *zt-5.72466e-3_wp 701 zn4= 4.8314e-4_wp 702 zn= ( zn4*zs + zn3*zsr + zn2 ) *zs + zn1 703 ! 704 ! add the compression terms 705 za2= ( -3.508914e-8_wp*zt-1.248266e-8_wp ) *zt-2.595994e-6_wp 706 za1= ( 1.296821e-6_wp*zt-5.782165e-9_wp ) *zt+1.045941e-4_wp 707 za = za1 + za2 * zs 708 ! 709 zb1= ( (-5.939910e-6_wp*zt-2.512549e-3_wp ) *zt+0.1028859_wp ) *zt + 3.721788_wp 710 zb2= (7.267926e-5_wp*zt-2.598241e-3_wp ) *zt-0.1571896_wp 711 zb3= 2.042967e-2_wp 712 zb = ( zb3*zsr + zb2 ) *zs + zb1 713 ! 714 zc1= ( ( (-1.361629e-4_wp*zt-1.852732e-2_wp ) *zt-30.41638_wp ) *zt + 2098.925_wp ) *zt+190925.6_wp 715 zc2= ( ( 2.326469e-3_wp*zt+1.553190_wp) *zt-65.00517_wp ) *zt+1044.077_wp 716 zc3= (-0.1909078_wp*zt+7.390729_wp ) *zt-55.87545_wp 717 zc = ( zc3*zsr + zc2 )*zs + zc1 718 ! 719 zm_inv = 1._wp / ( ( za*zh + zb )*zh + zc ) 720 ! alpha 390 ! ! alpha 721 391 zdza = (2.593642e-6_wp*zt-5.782165e-9_wp) + (-7.017828e-8_wp*zt-1.248266e-8_wp) * zs 722 392 zdzb = (14.535852e-5_wp*zt-2.598241e-3)*zs+(-17.81973e-6_wp*zt-5.025098e-3_wp)*zt+0.1028859_wp … … 727 397 & +((((32.68166e-9_wp*zt-4.480332e-6_wp)*zt+3.005055e-4_wp)*zt-18.19058e-3_wp)*zt+6.793952e-2_wp) 728 398 zdm = (zdza*zh+zdzb)*zh+zdzc 729 ! 730 pab(ji,jj,jk,jp_tem) = - ( & 731 & zdn * ( 1._wp + zh * zm_inv ) - zn * zh * zdm * zm_inv * zm_inv ) & 732 & / rau0 * tmask(ji,jj,jk) 733 ! beta 734 zdza = za2 399 pab(ji,jj,jk,jp_tem) = - ( zdn * ( 1._wp + zh * zm_inv ) - zn * zh * zdm * zm_inv * zm_inv ) & 400 & * r1_rau0 * tmask(ji,jj,jk) 401 ! 402 zdza = za2 ! beta 735 403 zdzb = 3.0644505e-2_wp*zsr+zb2 736 404 zdzc = 1.5_wp*zc3*zsr+zc2 737 405 zdn = 9.6628e-4_wp*zs+1.5_wp*zn3*zsr+zn2 738 406 zdm = (zdza*zh+zdzb)*zh+zdzc 739 pab(ji,jj,jk,jp_sal) = ( & 740 & zdn * ( 1._wp + zh * zm_inv ) - zn * zh * zdm * zm_inv * zm_inv ) & 741 & / rau0 * tmask(ji,jj,jk) 742 ! 407 pab(ji,jj,jk,jp_sal) = ( zdn * ( 1._wp + zh * zm_inv ) - zn * zh * zdm * zm_inv * zm_inv ) & 408 & * r1_rau0 * tmask(ji,jj,jk) 743 409 END DO 744 410 END DO 745 411 END DO 746 412 ! 747 CASE( 1 ) !== Vallis (2006)simplified EOS ==!413 CASE( 1 ) !== simplified EOS ==! 748 414 DO jk = 1, jpkm1 749 415 DO jj = 1, jpj 750 416 DO ji = 1, jpi 751 zt = pts (ji,jj,jk,jp_tem) - 10._wp ! potentialtemperature anomaly (t-T0)752 zh = fsdept(ji,jj,jk) ! depth in meters at t-point753 ! 417 zt = pts (ji,jj,jk,jp_tem) - 10._wp ! pot. temperature anomaly (t-T0) 418 zh = fsdept(ji,jj,jk) ! depth in meters at t-point 419 ! ! alpha 754 420 pab(ji,jj,jk,jp_tem) = rn_alph0 * ( 1._wp + 0.5_wp*rn_cabbel*zt + rn_thermo*zh ) * tmask(ji,jj,jk) ! alpha 421 ! ! beta = constant: set one for all during the initialisation 755 422 END DO 756 423 END DO 757 424 END DO 758 pab(:,:,:,jp_sal) = rn_beta0 * tmask(:,:,:) ! beta759 !760 CASE DEFAULT761 IF(lwp) WRITE(numout,cform_err)762 IF(lwp) WRITE(numout,*) ' bad flag value for nn_eos = ', nn_eos763 nstop = nstop + 1764 425 ! 765 426 END SELECT 766 427 ! 767 CALL wrk_dealloc( jpi, jpj, jpk, zws ) 768 ! 769 IF( nn_timing == 1 ) CALL timing_stop('eos_ab') 770 ! 771 END SUBROUTINE eos_ab 772 773 774 775 SUBROUTINE eos_insitu_ab( pts, prd, pab ) 776 !!---------------------------------------------------------------------- 777 !! *** ROUTINE eos_insitu_ab *** 778 !! 779 !! ** Purpose : Calculates the in situ thermal/haline expansion terms 780 !! and the insitu density anomaly at T-points 781 !! 782 !! ** Method : calculates rhd, alpha, beta at T-points 783 !! * nn_eos = -1 : Jackett and McDougall (1995) equation of state 784 !! * nn_eos = 0 : modified Jackett and McDougall (1995) equation of state 785 !! * nn_eos = 1 : Vallis equation of state (Vallis 2006) 786 !! alpha and beta ratios are returned as 3-D arrays. 787 !! 788 !! ** Action : - pab : thermal and haline expansion ratios at T-points 789 !! - prd , the density anomaly (no units) 428 ! !++ Brunt-Vaisala frequency (only for 3D computation) ++! 429 IF( PRESENT( pn2 ) ) CALL eos_bn2( pts, pab, pn2 ) ! N^2 430 ! 431 IF( nn_timing == 1 ) CALL timing_stop('eos_rab_bn2') 432 ! 433 END SUBROUTINE eos_rab_bn2 434 435 436 SUBROUTINE eos_bn2( pts, pab, pn2 ) 437 !!---------------------------------------------------------------------- 438 !! *** ROUTINE eos_bn2 *** 439 !! 440 !! ** Purpose : Compute the local Brunt-Vaisala frequency at the 441 !! time-step of the input arguments 442 !! 443 !! ** Method : pn2 = grav * (alpha dk[T] + beta dk[S] ) / e3w 444 !! where alpha and beta have been vertically interpolated from T- to w-point. 445 !! N.B. N^2 is set one for all to zero at jk=1 in istate module. 446 !! 447 !! ** Action : pn2 : square of the brunt-vaisala frequency at w-point 790 448 !! 791 449 !! References : Vallis, Atmospheric and Oceanic Fluid Dynamics, 2006 792 450 !! Jackett and McDougall, J. Atmos. Ocean. Tech., 1995 793 !!---------------------------------------------------------------------- 794 REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in ) :: pts ! pot. temperature & salinity [Celcius,psu] 795 REAL(wp), DIMENSION(jpi,jpj,jpk) , INTENT( out) :: prd ! density anomaly [-] 796 REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT( out) :: pab ! partial derivative of density anomaly 797 ! ! with respect to T and S [1/Celcius,1/psu] 798 ! 799 INTEGER :: ji, jj, jk ! dummy loop indices 800 REAL(wp) :: zt , zs , zh , zsr ! local scalars 801 REAL(wp) :: zn1, zn2, zn3, zn4 ! - - 802 REAL(wp) :: zn, za1, za2, za ! - - 803 REAL(wp) :: zb1, zb2, zb3, zb ! - - 804 REAL(wp) :: zc1, zc2, zc3, zc ! - - 805 REAL(wp) :: zm_inv, zdm, zdn ! - - 806 REAL(wp) :: zdza, zdzb, zdzc ! - - 807 REAL(wp), POINTER, DIMENSION(:,:,:) :: zws 808 !!---------------------------------------------------------------------- 809 ! 810 IF ( nn_timing == 1 ) CALL timing_start('eos_insitu_ab') 811 ! 812 CALL wrk_alloc( jpi, jpj, jpk, zws ) 813 ! 814 SELECT CASE ( nn_eos ) 815 ! 816 CASE ( -1 ) ! Jackett and McDougall (1995) formulation 817 ! 818 zws(:,:,:) = SQRT( ABS( pts(:,:,:,jp_sal) ) ) 819 DO jk = 1, jpkm1 820 DO jj = 1, jpj 821 DO ji = 1, jpi 822 zt = pts (ji,jj,jk,jp_tem) 823 zs = pts (ji,jj,jk,jp_sal) 824 zh = fsdept(ji,jj,jk) ! depth 825 zsr= zws (ji,jj,jk) ! square root salinity 826 ! 827 ! compute volumic mass pure water at atm pressure 828 zn1= ( ( ( ( 6.536332e-9_wp *zt - 1.120083e-6_wp )*zt + 1.001685e-4_wp )*zt & 829 & -9.095290e-3_wp )*zt + 6.793952e-2_wp )*zt + 999.842594_wp 830 ! seawater volumic mass atm pressure 831 zn2= ( ( ( 5.3875e-9_wp*zt-8.2467e-7_wp ) *zt+7.6438e-5_wp ) *zt & 832 & -4.0899e-3_wp ) *zt+0.824493_wp 833 zn3= ( -1.6546e-6_wp*zt+1.0227e-4_wp ) *zt-5.72466e-3_wp 834 zn4= 4.8314e-4_wp 835 zn = ( zn4*zs + zn3*zsr + zn2 ) *zs + zn1 836 ! 837 ! add the compression terms 838 za1= ( 1.394680e-8_wp*zt-1.202016e-6_wp ) *zt+2.102898e-5_wp 839 za2 = ( 6.207323e-11_wp*zt+6.128773e-9_wp ) *zt-2.040237e-7_wp 840 za = za1 + za2 * zs 841 ! 842 zb1 = ( ( 1.956415e-6_wp*zt-2.984642e-4_wp ) *zt+2.212276e-2_wp ) *zt +2.186519_wp 843 zb2 = ( 2.059331e-7_wp*zt-1.847318e-4_wp ) *zt+6.704388e-3 844 zb3 = 1.480266e-4_wp 845 zb = ( zb3*zsr + zb2 ) *zs + zb1 846 ! 847 zc1 = ( ( (-4.190253e-4_wp*zt+9.648704e-2_wp)*zt-17.06103_wp )*zt+1444.304_wp)*zt+196593.3_wp 848 zc2 = ( ( -5.084188e-4_wp*zt+6.283263e-2_wp) *zt-3.101089_wp ) *zt+528.4855_wp 849 zc3 = (-4.619924e-3_wp*zt+9.085835e-2_wp ) *zt+3.886640_wp 850 zc = ( zc3*zsr + zc2 )*zs + zc1 851 ! 852 zm_inv = 1._wp / ( ( za*zh + zb )*zh + zc ) 853 ! 854 ! masked in situ density anomaly. Important: rau0=1035, should be 1027! 855 prd(ji,jj,jk) = ( zn * ( 1.0_wp + zh * zm_inv ) / rau0 - 1._wp ) * tmask(ji,jj,jk) 856 ! 857 ! alpha 858 zdza = (2.78936e-8_wp*zt-1.202016e-6_wp) + (12.414646e-11_wp*zt+6.128773e-9_wp ) * zs 859 zdzb = (4.118662e-7_wp*zt-1.847318e-4_wp)*zs + & 860 & ( 5.869245e-6_wp*zt-5.969284e-4_wp)*zt+2.212276e-2_wp 861 zdzc = ( (-9.239848e-3_wp*zt+9.085835e-2_wp)*zsr + & 862 & (-15.252564e-4_wp*zt+12.566526e-2_wp)*zt-3.101089_wp ) * zs + & 863 & ((-16.761012e-4_wp*zt+28.946112e-2_wp)*zt-34.12206_wp )*zt+1444.304_wp 864 865 zdn = ((-3.3092e-6_wp*zt+1.0227e-4_wp)*zsr & 866 & +(((21.55e-9_wp*zt-24.7401e-7_wp)*zt+15.2876e-5_wp)*zt-4.0899e-3_wp))*zs & 867 & +((((32.68166e-9_wp*zt-4.480332e-6_wp)*zt+3.005055e-4_wp)*zt & 868 & -18.19058e-3_wp)*zt+6.793952e-2_wp) 869 zdm = (zdza*zh+zdzb)*zh+zdzc 870 ! 871 pab(ji,jj,jk,jp_tem) = - ( & 872 & zdn * ( 1._wp + zh * zm_inv ) - zn * zh * zdm * zm_inv * zm_inv ) & 873 & / rau0 * tmask(ji,jj,jk) 874 ! 875 ! beta 876 zdza = za2 877 zdzb = 2.220399e-4_wp*zsr+zb2 878 zdzc = 1.5_wp*zc3*zsr+zc2 879 zdn = 9.6628e-4_wp*zs+1.5_wp*zn3*zsr+zn2 880 zdm = (zdza*zh+zdzb)*zh+zdzc 881 pab(ji,jj,jk,jp_sal) = ( & 882 & zdn * ( 1._wp + zh * zm_inv ) - zn * zh * zdm * zm_inv * zm_inv ) & 883 & / rau0 * tmask(ji,jj,jk) 884 ! 885 END DO 886 END DO 887 END DO 888 ! 889 CASE ( 0 ) ! modified Jackett and McDougall (1995) formulation 890 ! 891 zws(:,:,:) = SQRT( ABS( pts(:,:,:,jp_sal) ) ) 892 DO jk = 1, jpkm1 893 DO jj = 1, jpj 894 DO ji = 1, jpi 895 zt = pts (ji,jj,jk,jp_tem) 896 zs = pts (ji,jj,jk,jp_sal) 897 zh = fsdept(ji,jj,jk) ! depth 898 zsr= zws (ji,jj,jk) ! square root salinity 899 ! 900 ! compute volumic mass pure water at atm pressure 901 zn1= ( ( ( ( 6.536332e-9_wp *zt - 1.120083e-6_wp )*zt + 1.001685e-4_wp )*zt & 902 & -9.095290e-3_wp )*zt + 6.793952e-2_wp )*zt + 999.842594_wp 903 ! seawater volumic mass atm pressure 904 zn2= ( ( ( 5.3875e-9_wp*zt-8.2467e-7_wp ) *zt+7.6438e-5_wp ) *zt & 905 & -4.0899e-3_wp ) *zt+0.824493_wp 906 zn3= ( -1.6546e-6_wp*zt+1.0227e-4_wp ) *zt-5.72466e-3_wp 907 zn4= 4.8314e-4_wp 908 zn= ( zn4*zs + zn3*zsr + zn2 ) *zs + zn1 909 ! 910 ! add the compression terms 911 za2= ( -3.508914e-8_wp*zt-1.248266e-8_wp ) *zt-2.595994e-6_wp 912 za1= ( 1.296821e-6_wp*zt-5.782165e-9_wp ) *zt+1.045941e-4_wp 913 za = za1 + za2 * zs 914 ! 915 zb1= ( (-5.939910e-6_wp*zt-2.512549e-3_wp ) *zt+0.1028859_wp ) *zt + 3.721788_wp 916 zb2= (7.267926e-5_wp*zt-2.598241e-3_wp ) *zt-0.1571896_wp 917 zb3= 2.042967e-2_wp 918 zb = ( zb3*zsr + zb2 ) *zs + zb1 919 ! 920 zc1= ( ( (-1.361629e-4_wp*zt-1.852732e-2_wp ) *zt-30.41638_wp ) *zt + 2098.925_wp ) *zt+190925.6_wp 921 zc2= ( ( 2.326469e-3_wp*zt+1.553190_wp) *zt-65.00517_wp ) *zt+1044.077_wp 922 zc3= (-0.1909078_wp*zt+7.390729_wp ) *zt-55.87545_wp 923 zc = ( zc3*zsr + zc2 )*zs + zc1 924 ! 925 zm_inv = 1._wp / ( ( za*zh + zb )*zh + zc ) 926 ! 927 ! masked in situ density anomaly (rho/rho0 - 1). Important: rau0=1035, should be 1027! 928 prd(ji,jj,jk) = ( zn * ( 1.0_wp + zh * zm_inv ) * r1_rau0 - 1._wp ) * tmask(ji,jj,jk) 929 ! 930 ! ! alpha 931 zdza = (2.593642e-6_wp*zt-5.782165e-9_wp) + (-7.017828e-8_wp*zt-1.248266e-8_wp) * zs 932 zdzb = (14.535852e-5_wp*zt-2.598241e-3)*zs+(-17.81973e-6_wp*zt-5.025098e-3_wp)*zt+0.1028859_wp 933 zdzc = ((-0.3818156_wp*zt+7.390729_wp)*zsr+((6.979407e-3_wp*zt+3.10638_wp)*zt-65.00517_wp))*zs & 934 & +(((-5.446516e-4_wp*zt-5.558196e-2_wp)*zt-60.83276_wp)*zt+2098.925_wp) 935 zdn = ((-3.3092e-6_wp*zt+1.0227e-4_wp)*zsr & 936 & +(((21.55e-9_wp*zt-24.7401e-7_wp)*zt+15.2876e-5_wp)*zt-4.0899e-3_wp))*zs & 937 & +((((32.68166e-9_wp*zt-4.480332e-6_wp)*zt+3.005055e-4_wp)*zt-18.19058e-3_wp)*zt+6.793952e-2_wp) 938 zdm = (zdza*zh+zdzb)*zh+zdzc 939 ! 940 pab(ji,jj,jk,jp_tem) = - ( zdn * ( 1._wp + zh * zm_inv ) - zn * zh * zdm * zm_inv * zm_inv ) & 941 & * r1_rau0 * tmask(ji,jj,jk) 942 ! 943 ! ! beta 944 zdza = za2 945 zdzb = 3.0644505e-2_wp*zsr+zb2 946 zdzc = 1.5_wp*zc3*zsr+zc2 947 zdn = 9.6628e-4_wp*zs+1.5_wp*zn3*zsr+zn2 948 zdm = (zdza*zh+zdzb)*zh+zdzc 949 pab(ji,jj,jk,jp_sal) = ( zdn * ( 1._wp + zh * zm_inv ) - zn * zh * zdm * zm_inv * zm_inv ) & 950 & * r1_rau0 * tmask(ji,jj,jk) 951 END DO 952 END DO 953 END DO 954 ! 955 CASE( 1 ) !== Vallis (2006) simplified EOS ==! 956 DO jk = 1, jpkm1 957 DO jj = 1, jpj 958 DO ji = 1, jpi 959 zt = pts(ji,jj,jk,jp_tem) - 10._wp ! potential temperature (t-T0) 960 zh = fsdept(ji,jj,jk) ! depth in meters at t-point 961 ! masked in situ density anomaly 962 prd(ji,jj,jk) = ( offset - rn_alph0 * ( 1._wp + rn_cabbel*zt + rn_thermo*zh ) *zt & & 963 & + rn_beta0 * zs + rn_gamm0 * zh ) * tmask(ji,jj,jk) 964 ! alpha 965 pab(ji,jj,jk,jp_tem) = rn_alph0 * ( 1._wp + 0.5_wp*rn_cabbel*zt + rn_thermo*zh ) & 966 & * tmask(ji,jj,jk) ! alpha 967 ! 968 END DO 969 END DO 970 END DO 971 ! beta 972 pab (:,:,:,jp_sal) = rn_beta0 * tmask(:,:,:) ! beta 973 ! 974 CASE DEFAULT 975 IF(lwp) WRITE(numout,cform_err) 976 IF(lwp) WRITE(numout,*) ' bad flag value for nn_eos = ', nn_eos 977 nstop = nstop + 1 978 ! 979 END SELECT 980 ! 981 CALL wrk_dealloc( jpi, jpj, jpk, zws ) 982 ! 983 IF( nn_timing == 1 ) CALL timing_stop('eos_insitu_ab') 984 ! 985 END SUBROUTINE eos_insitu_ab 986 451 !! McDougall, J. Phys. Oceano., 1987 452 !!---------------------------------------------------------------------- 453 REAL(wp), DIMENSION(:,:,:,:), INTENT(in ) :: pts ! pot. temperature and salinity [Celcius,psu] 454 REAL(wp), DIMENSION(:,:,:,:), INTENT(in ) :: pab ! partial derivative of in situ 455 ! ! density with respect to T and S [1/Celcius,1/psu] 456 REAL(wp), DIMENSION(:,:,:) , INTENT( out) :: pn2 ! Brunt-Vaisala frequency [1/s^2] 457 !! 458 INTEGER :: ji, jj, jk ! dummy loop indices 459 REAL(wp) :: zaw, zbw, zrw ! local scalars 460 !!---------------------------------------------------------------------- 461 ! 462 IF( nn_timing == 1 ) CALL timing_start('eos_bn2') 463 ! 464 ! pn2 : interior points only (2=< jk =< jpkm1 ) 465 ! -------------------------- 466 ! 467 DO jk = 2, jpkm1 468 DO jj = 1, jpj 469 DO ji = 1, jpi 470 zrw = ( fsdepw(ji,jj,jk ) - fsdept(ji,jj,jk) ) & 471 & / ( fsdept(ji,jj,jk-1) - fsdept(ji,jj,jk) ) 472 ! 473 zaw = ( pab(ji,jj,jk,jp_tem) * (1. - zrw) + pab(ji,jj,jk-1,jp_tem) * zrw ) * tmask(ji,jj,jk) 474 zbw = ( pab(ji,jj,jk,jp_sal) * (1. - zrw) + pab(ji,jj,jk-1,jp_sal) * zrw ) * tmask(ji,jj,jk) 475 ! 476 pn2(ji,jj,jk) = grav * ( zaw * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) ) & 477 & - zbw * ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) ) ) / fse3w(ji,jj,jk) 478 END DO 479 END DO 480 END DO 481 ! 482 IF(ln_ctl) CALL prt_ctl( tab3d_1=pn2, clinfo1=' bn2 : ', ovlap=1, kdim=jpk ) 483 ! 484 IF( nn_timing == 1 ) CALL timing_stop('eos_bn2') 485 ! 486 END SUBROUTINE eos_bn2 487 488 489 FUNCTION eos_fzp( psal, pdep ) RESULT( ptf ) 490 !!---------------------------------------------------------------------- 491 !! *** ROUTINE eos_fzp *** 492 !! 493 !! ** Purpose : Compute the freezing point temperature [Celcius] 494 !! 495 !! ** Method : UNESCO freezing point (ptf) in Celcius is given by 496 !! ptf(t,z) = (-.0575+1.710523e-3*sqrt(abs(s))-2.154996e-4*s)*s - 7.53e-4*z 497 !! checkvalue: tf=-2.588567 Celsius for s=40psu, z=500m 498 !! 499 !! Reference : UNESCO tech. papers in the marine science no. 28. 1978 500 !!---------------------------------------------------------------------- 501 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: psal ! salinity [psu] 502 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ), OPTIONAL :: pdep ! depth [m] 503 ! Leave result array automatic rather than making explicitly allocated 504 REAL(wp), DIMENSION(jpi,jpj) :: ptf ! freezing temperature [Celcius] 505 !!---------------------------------------------------------------------- 506 ! 507 ptf(:,:) = ( - 0.0575_wp + 1.710523e-3_wp * SQRT( psal(:,:) ) & 508 & - 2.154996e-4_wp * psal(:,:) ) * psal(:,:) 509 ! 510 IF( PRESENT( pdep ) ) ptf(:,:) = ptf(:,:) - 7.53e-4 * pdep(:,:) 511 ! 512 END FUNCTION eos_fzp 987 513 988 514 … … 1027 553 REAL(wp) :: zddelta, zdAp ! - - 1028 554 REAL(wp) :: zdlogBp, zdCp, zdP ! - - 1029 REAL(wp), POINTER, DIMENSION(:,:,:) :: zws1030 555 !!---------------------------------------------------------------------- 1031 556 ! 1032 557 IF ( nn_timing == 1 ) CALL timing_start('eos_pen') 1033 558 ! 1034 CALL wrk_alloc( jpi, jpj, jpk, zws )1035 !1036 559 SELECT CASE ( nn_eos ) 1037 560 ! 1038 561 CASE ( -1 ) ! Jackett and McDougall (1995) formulation 1039 562 ! 1040 zws(:,:,:) = SQRT( ABS( pts(:,:,:,jp_sal) ) )1041 563 DO jk = 1, jpkm1 1042 564 DO jj = 1, jpj … … 1044 566 zt = pts (ji,jj,jk,jp_tem) 1045 567 zs = pts (ji,jj,jk,jp_sal) 1046 zh = fsdept(ji,jj,jk) ! depth (Important: zh must be different from 0)1047 zsr= zws (ji,jj,jk)! square root salinity568 zh = fsdept(ji,jj,jk) ! depth (Important: zh must be different from 0) 569 zsr= SQRT( ABS( pts(ji,jj,jk,jp_sal) ) ) ! square root salinity 1048 570 ! 1049 571 ! compute volumic mass pure water at atm pressure … … 1100 622 ! 1101 623 zm_inv = 1._wp / ( ( za*zh + zb )*zh + zc ) 1102 zsgn = SIGN( 1. , zdelta ) 1103 zsdelta = SQRT( zsgn * zdelta ) 624 zsdelta = SQRT( ABS( zdelta ) ) 1104 625 zAp = - zb / za / zsdelta 1105 626 zlogBp = - LOG( zc * zm_inv ) 1106 627 zCp = zh * zsdelta / (2._wp*zc+zh*zb) 1107 IF( zsgn > 0. ) THEN ; zatCp = ATAN(zCp); 1108 ELSE ; zatCp = ATANH(zCp) 1109 ENDIF 628 zsgn = SIGN( 1._wp , zdelta ) 629 zatCp = 0.5_wp * ( ATAN( zCp ) * (1._wp + zsgn) + ATANH( zCp ) * (1._wp - zsgn) ) 1110 630 zP = ( 0.5_wp*zlogBp/za + zAp * zatCp ) / zh 1111 631 ! 1112 ! computepotential energy1113 pen(ji,jj,jk) = ( zn - rau0 + zn * zP ) /rau0 * tmask(ji,jj,jk)1114 ! 1115 ! alphaPE632 ! ! potential energy 633 pen(ji,jj,jk) = ( zn - rau0 + zn * zP ) * r1_rau0 * tmask(ji,jj,jk) 634 ! 635 ! ! alphaPE 1116 636 zddelta = ( 2._wp * (za*zdzc+zc*zdza) - zb*zdzb ) / zdelta 1117 637 zdAp = zAp * (zdzb/zb - zdza/za - zddelta) … … 1121 641 & + zAp*zdCp / ( 1._wp + zsgn*zCp*zCp ) ) / zh 1122 642 ! 1123 pab_pe(ji,jj,jk,jp_tem) = - ( zdn * (1._wp + zP ) + zn * zdP ) /rau0 * tmask(ji,jj,jk)1124 ! dds1125 zdza = za2 643 pab_pe(ji,jj,jk,jp_tem) = - ( zdn * (1._wp + zP ) + zn * zdP ) * r1_rau0 * tmask(ji,jj,jk) 644 ! 645 zdza = za2 ! dds 1126 646 zdzb = 2.220399e-4_wp*zsr+zb2 1127 647 zdzc = 1.5_wp*zc3*zsr+zc2 1128 648 zdn = 9.6628e-4_wp*zs+1.5_wp*zn3*zsr+zn2 1129 649 zdm = (zdza*zh+zdzb)*zh+zdzc 1130 ! betaPE650 ! ! betaPE 1131 651 zddelta = ( 2._wp * (za*zdzc+zc*zdza) - zb*zdzb ) / zdelta 1132 652 zdAp = zAp * (zdzb/zb - zdza/za - zddelta) … … 1136 656 & zAp*zdCp/(1._wp+zsgn*zCp*zCp) ) / zh 1137 657 ! 1138 pab_pe(ji,jj,jk,jp_sal) = ( zdn * (1._wp + zP ) + zn * zdP ) & 1139 & / rau0 * tmask(ji,jj,jk) 1140 ! 658 pab_pe(ji,jj,jk,jp_sal) = ( zdn * (1._wp + zP ) + zn * zdP ) * r1_rau0 * tmask(ji,jj,jk) 1141 659 END DO 1142 660 END DO … … 1145 663 CASE ( 0 ) ! modified Jackett and McDougall (1995) formulation 1146 664 ! 1147 zws(:,:,:) = SQRT( ABS( pts(:,:,:,jp_sal) ) )1148 665 DO jk = 1, jpkm1 1149 666 DO jj = 1, jpj 1150 667 DO ji = 1, jpi 1151 !1152 668 zt = pts (ji,jj,jk,jp_tem) 1153 669 zs = pts (ji,jj,jk,jp_sal) 1154 zh = fsdept(ji,jj,jk) ! depth (Important: zh must be different from 0)1155 zsr= zws (ji,jj,jk)! square root salinity670 zh = fsdept(ji,jj,jk) ! depth (Important: zh must be different from 0) 671 zsr= SQRT( ABS( pts(ji,jj,jk,jp_sal) ) ) ! square root salinity 1156 672 ! 1157 673 ! compute volumic mass pure water at atm pressure … … 1182 698 zdelta = 4._wp * za * zc - zb * zb 1183 699 ! 1184 ! dds 1185 zdza = za2 700 zdza = za2 ! dds 1186 701 zdzb = 3.0644505e-2_wp*zsr+zb2 1187 702 zdzc = 1.5_wp*zc3*zsr+zc2 … … 1189 704 zdm = (zdza*zh+zdzb)*zh+zdzc 1190 705 ! stability condition for the calculation (add a small perturbation on S if needed) 1191 IF( ABS(zdelta) < 1.e-5 .OR. ABS(za) <1.e-9 ) THEN ;!stability condition for the calculation706 IF( ABS(zdelta) < 1.e-5 .OR. ABS(za) < 1.e-9 ) THEN !stability condition for the calculation 1192 707 zs = zs + 1.e-3_wp 1193 708 zsr= zsr + 0.5e-3_wp/zsr … … 1199 714 ! 1200 715 zdelta = 4._wp * za * zc - zb * zb 1201 !1202 716 ENDIF 1203 717 ! 1204 718 zm_inv = 1._wp / ( ( za*zh + zb )*zh + zc ) 1205 zsgn = SIGN( 1. , zdelta ) 1206 zsdelta = SQRT( zsgn * zdelta ) 719 zsdelta = SQRT( ABS( zdelta ) ) 1207 720 zAp = - zb / za / zsdelta 1208 721 zlogBp = - LOG( zc * zm_inv ) 1209 722 zCp = zh * zsdelta / (2._wp*zc+zh*zb) 1210 IF( zsgn > 0. ) THEN ; zatCp = ATAN(zCp); 1211 ELSE ; zatCp = ATANH(zCp) 1212 ENDIF 723 zsgn = SIGN( 1._wp , zdelta ) 724 zatCp = 0.5_wp * ( ATAN( zCp ) * (1._wp + zsgn) + ATANH( zCp ) * (1._wp - zsgn) ) 1213 725 zP = ( 0.5_wp*zlogBp/za + zAp * zatCp ) / zh 1214 726 ! 1215 ! computepotential energy1216 pen(ji,jj,jk) = ( zn - rau0 + zn * zP ) /rau0 * tmask(ji,jj,jk)1217 ! 1218 ! betaPE727 ! ! potential energy 728 pen(ji,jj,jk) = ( zn - rau0 + zn * zP ) * r1_rau0 * tmask(ji,jj,jk) 729 ! 730 ! ! betaPE 1219 731 zddelta = ( 2._wp * (za*zdzc+zc*zdza) - zb*zdzb ) / zdelta 1220 732 zdAp = zAp * (zdzb/zb - zdza/za - zddelta) … … 1224 736 & zAp*zdCp/(1._wp+zsgn*zCp*zCp) ) / zh 1225 737 ! 1226 pab_pe(ji,jj,jk,jp_sal) = ( zdn * (1._wp + zP ) + zn * zdP ) & 1227 & / rau0 * tmask(ji,jj,jk) 1228 ! 1229 ! ddt 738 pab_pe(ji,jj,jk,jp_sal) = ( zdn * (1._wp + zP ) + zn * zdP ) * r1_rau0 * tmask(ji,jj,jk) 739 ! 740 ! ! ddt 1230 741 zdza = (2.593642e-6_wp*zt-5.782165e-9_wp) + (-7.017828e-8_wp*zt-1.248266e-8_wp) * zs 1231 742 zdzb = (14.535852e-5_wp*zt-2.598241e-3)*zs+(-17.81973e-6_wp*zt-5.025098e-3_wp)*zt+0.1028859_wp … … 1236 747 & +((((32.68166e-9_wp*zt-4.480332e-6_wp)*zt+3.005055e-4_wp)*zt-18.19058e-3_wp)*zt+6.793952e-2_wp) 1237 748 zdm = (zdza*zh+zdzb)*zh+zdzc 1238 ! alphaPE749 ! ! alphaPE 1239 750 zddelta = ( 2._wp * (za*zdzc+zc*zdza) - zb*zdzb ) / zdelta 1240 751 zdAp = zAp * (zdzb/zb - zdza/za - zddelta) … … 1244 755 & zAp*zdCp/(1._wp+zsgn*zCp*zCp) ) / zh 1245 756 ! 1246 pab_pe(ji,jj,jk,jp_tem) = - ( zdn * (1._wp + zP ) + zn * zdP ) / rau0 * tmask(ji,jj,jk) 1247 ! 757 pab_pe(ji,jj,jk,jp_tem) = - ( zdn * (1._wp + zP ) + zn * zdP ) * r1_rau0 * tmask(ji,jj,jk) 1248 758 END DO 1249 759 END DO … … 1254 764 DO jj = 1, jpj 1255 765 DO ji = 1, jpi 1256 zt = pts(ji,jj,jk,jp_tem) - 10._wp ! potential temperature (t-T0) 1257 zh = fsdept(ji,jj,jk) ! depth in meters at t-point 1258 pab_pe(ji,jj,jk,jp_tem) = rn_alph0 * ( 1._wp + 0.5_wp*(rn_cabbel*zt + rn_thermo*zh) ) & 1259 & * tmask(ji,jj,jk) ! alphaPE 766 zt = pts(ji,jj,jk,jp_tem) - 10._wp ! potential temperature (t-T0) 767 zh = fsdept(ji,jj,jk) ! depth in meters at t-point 768 ! ! alphaPE 769 pab_pe(ji,jj,jk,jp_tem) = rn_alph0 * ( 1._wp + 0.5_wp*(rn_cabbel*zt + rn_thermo*zh) ) * tmask(ji,jj,jk) 770 ! ! Potential Energy 1260 771 pen(ji,jj,jk) = ( - rn_alph0 * ( 1._wp + rn_cabbel*zt + 0.5_wp*rn_thermo*zh ) *zt & & 1261 772 & + rn_beta0 * zs + 0.5_wp*rn_gamm0 * zh ) * tmask(ji,jj,jk) 1262 !1263 773 END DO 1264 774 END DO … … 1273 783 END SELECT 1274 784 ! 1275 CALL wrk_dealloc( jpi, jpj, jpk, zws ) 1276 ! 1277 IF( nn_timing == 1 ) CALL timing_stop('eos_pen') 785 IF( nn_timing == 1 ) CALL timing_stop('eos_pen') 1278 786 ! 1279 787 END SUBROUTINE eos_pen 1280 1281 1282 1283 SUBROUTINE eos_bn2_ab( pts, pn2 )1284 !!----------------------------------------------------------------------1285 !! *** ROUTINE eos_bn2_ab ***1286 !!1287 !! ** Purpose : Update alpha/beta then compute the local Brunt-Vaisala1288 !! frequency at the time-step of the input arguments1289 !!1290 !! ** Method : call eos_ab to obtain alpha and beta coefficients1291 !! then interpolate them vertically, and compute pn21292 !! Macro-tasked on horizontal slab (jk-loop)1293 !! N.B. N^2 is set to zero at the first level (JK=1) in inidtr1294 !! and is never used at this level.1295 !!1296 !! ** Action : - pn2 : the brunt-vaisala frequency1297 !!1298 !! References : Vallis, Atmospheric and Oceanic Fluid Dynamics, 20061299 !! Jackett and McDougall, J. Atmos. Ocean. Tech., 19951300 !! McDougall, J. Phys. Oceano., 19871301 !!----------------------------------------------------------------------1302 REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in ) :: pts ! 1 : potential temperature [Celcius]1303 ! ! 2 : salinity [psu]1304 REAL(wp), DIMENSION(jpi,jpj,jpk) , INTENT( out) :: pn2 ! Brunt-Vaisala frequency [s-1]1305 !!1306 INTEGER :: ji, jj, jk ! dummy loop indices1307 REAL(wp) :: zgde3w, za, zb, zr1 ! temporary scalars1308 #if defined key_zdfddm1309 REAL(wp) :: zds ! local scalars1310 #endif1311 REAL(wp), POINTER, DIMENSION(:,:,:,:) :: zab1312 !!----------------------------------------------------------------------1313 !1314 !1315 IF( nn_timing == 1 ) CALL timing_start('eos_bn2_ab')1316 !1317 ! pn2 : interior points only (2=< jk =< jpkm1 )1318 ! --------------------------1319 !1320 CALL wrk_alloc( jpi, jpj, jpk, jpts, zab )1321 !1322 CALL eos_ab( pts, zab )1323 !1324 DO jk = 2, jpkm11325 DO jj = 1, jpj1326 DO ji = 1, jpi1327 !1328 zgde3w = 1._wp / fse3w(ji,jj,jk)1329 zr1 = ( fsdepw(ji,jj,jk) - fsdept(ji,jj,jk) ) * zgde3w1330 !1331 za = ( zab(ji,jj,jk,jp_tem) + &1332 & ( zab(ji,jj,jk-1,jp_tem) - zab(ji,jj,jk,jp_tem) ) * zr1 ) &1333 & * tmask(ji,jj,jk)1334 zb = ( zab(ji,jj,jk,jp_sal) + &1335 & ( zab(ji,jj,jk-1,jp_sal) - zab(ji,jj,jk,jp_sal) ) * zr1 ) &1336 & * tmask(ji,jj,jk)1337 ! N21338 pn2(ji,jj,jk) = grav * zgde3w & ! N^21339 & * ( za * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) ) &1340 & - zb * ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) ) )1341 #if defined key_zdfddm1342 zds = zb * ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) ) ! Rrau = (alpha / beta) (dk[t] / dk[s])1343 IF ( ABS( zds) <= 1.e-20_wp ) zds = 1.e-20_wp1344 rrau(ji,jj,jk) = za * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) ) / zds1345 #endif1346 END DO1347 END DO1348 END DO1349 !1350 1351 IF(ln_ctl) CALL prt_ctl( tab3d_1=pn2, clinfo1=' bn2 : ', ovlap=1, kdim=jpk )1352 #if defined key_zdfddm1353 IF(ln_ctl) CALL prt_ctl( tab3d_1=rrau, clinfo1=' rrau : ', ovlap=1, kdim=jpk )1354 #endif1355 !1356 CALL wrk_dealloc( jpi, jpj, jpk, jpts, zab )1357 !1358 IF( nn_timing == 1 ) CALL timing_stop('eos_bn2_ab')1359 !1360 END SUBROUTINE eos_bn2_ab1361 1362 1363 SUBROUTINE eos_bn2( pts, pab, pn2 )1364 !!----------------------------------------------------------------------1365 !! *** ROUTINE eos_bn2 ***1366 !!1367 !! ** Purpose : Compute the local Brunt-Vaisala frequency at the1368 !! time-step of the input arguments1369 !!1370 !! ** Method : interpolate alpha/beta vertically, and compute pn21371 !! Macro-tasked on horizontal slab (jk-loop)1372 !! N.B. N^2 is set to zero at the first level (JK=1) in inidtr1373 !! and is never used at this level.1374 !!1375 !! ** Action : - pn2 : the brunt-vaisala frequency1376 !!1377 !! References : Vallis, Atmospheric and Oceanic Fluid Dynamics, 20061378 !! Jackett and McDougall, J. Atmos. Ocean. Tech., 19951379 !! McDougall, J. Phys. Oceano., 19871380 !!----------------------------------------------------------------------1381 REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in ) :: pts ! pot. temperature and salinity [Celcius,psu]1382 REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in ) :: pab ! partial derivative of in situ1383 ! ! density with respect to T and S [1/Celcius,1/psu]1384 REAL(wp), DIMENSION(jpi,jpj,jpk) , INTENT( out) :: pn2 ! Brunt-Vaisala frequency [1/s^2]1385 !!1386 INTEGER :: ji, jj, jk ! dummy loop indices1387 REAL(wp) :: zgde3w, za, zb, zr1 ! temporary scalars1388 #if defined key_zdfddm1389 REAL(wp) :: zds ! local scalars1390 #endif1391 !!----------------------------------------------------------------------1392 !1393 !1394 IF( nn_timing == 1 ) CALL timing_start('eos_bn2')1395 !1396 ! pn2 : interior points only (2=< jk =< jpkm1 )1397 ! --------------------------1398 !1399 DO jk = 2, jpkm11400 DO jj = 1, jpj1401 DO ji = 1, jpi1402 zgde3w = 1._wp / fse3w(ji,jj,jk)1403 zr1 = ( fsdepw(ji,jj,jk) - fsdept(ji,jj,jk) ) * zgde3w1404 !1405 za = ( pab(ji,jj,jk,jp_tem) + ( pab(ji,jj,jk-1,jp_tem) - pab(ji,jj,jk,jp_tem) ) * zr1 ) * tmask(ji,jj,jk)1406 zb = ( pab(ji,jj,jk,jp_sal) + ( pab(ji,jj,jk-1,jp_sal) - pab(ji,jj,jk,jp_sal) ) * zr1 ) * tmask(ji,jj,jk)1407 !1408 pn2(ji,jj,jk) = grav * zgde3w & ! N^21409 & * ( za * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) ) &1410 & - zb * ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) ) )1411 #if defined key_zdfddm1412 zds = zb * ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) ) ! Rrau = (alpha / beta) (dk[t] / dk[s])1413 IF ( ABS( zds) <= 1.e-20_wp ) zds = 1.e-20_wp1414 rrau(ji,jj,jk) = za * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) ) / zds1415 #endif1416 END DO1417 END DO1418 END DO1419 !1420 1421 IF(ln_ctl) CALL prt_ctl( tab3d_1=pn2, clinfo1=' bn2 : ', ovlap=1, kdim=jpk )1422 #if defined key_zdfddm1423 IF(ln_ctl) CALL prt_ctl( tab3d_1=rrau, clinfo1=' rrau : ', ovlap=1, kdim=jpk )1424 #endif1425 !1426 IF( nn_timing == 1 ) CALL timing_stop('eos_bn2')1427 !1428 END SUBROUTINE eos_bn21429 1430 1431 1432 FUNCTION tfreez( psal ) RESULT( ptf )1433 !!----------------------------------------------------------------------1434 !! *** ROUTINE tfreez ***1435 !!1436 !! ** Purpose : Compute the sea surface freezing temperature [Celcius]1437 !!1438 !! ** Method : UNESCO freezing point at the surface (pressure = 0???)1439 !! freezing point [Celcius]=(-.0575+1.710523e-3*sqrt(abs(s))-2.154996e-4*s)*s-7.53e-4*p1440 !! checkvalue: tf= -2.588567 Celsius for s=40.0psu, p=500. decibars1441 !!1442 !! Reference : UNESCO tech. papers in the marine science no. 28. 19781443 !!----------------------------------------------------------------------1444 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: psal ! salinity [psu]1445 ! Leave result array automatic rather than making explicitly allocated1446 REAL(wp), DIMENSION(jpi,jpj) :: ptf ! freezing temperature [Celcius]1447 !!----------------------------------------------------------------------1448 !1449 ptf(:,:) = ( - 0.0575_wp + 1.710523e-3_wp * SQRT( psal(:,:) ) &1450 & - 2.154996e-4_wp * psal(:,:) ) * psal(:,:)1451 !1452 END FUNCTION tfreez1453 788 1454 789 … … 1495 830 IF(lwp) WRITE(numout,*) ' use of simplified eos: rhd(T,S,z) = ' 1496 831 IF(lwp) WRITE(numout,*) ' -alph0*(1+cabbel*(T-T0)+thermo*z)*(T-T0) + beta0*(S-S0) + gamm0*z' 1497 IF( lk_zdfddm .AND. rn_beta0==0. ) & 1498 & CALL ctl_stop( ' double diffusive mixing parameterization requires that rn_beta<>0') 1499 offset = 1027._wp / rau0 - 1._wp 832 ! 833 offset = 1027._wp * r1_rau0 - 1._wp ! offset value so that rho(10,35,0)=1027kg/m3 834 ! 835 rab_b(:,:,:,jp_sal) = rn_beta0 * tmask(:,:,:) ! beta = constant : set one for all 836 rab_n(:,:,:,jp_sal) = rn_beta0 * tmask(:,:,:) 837 ! 838 IF(lwp) WRITE(numout,*) ' offset value so that rho(10,35,0)=1027kg/m3 : offset = ', offset 839 ! 840 IF( rn_beta0 == 0._wp ) THEN ! rho is not a function of salinty 841 IF( lk_ldfslp ) CALL ctl_stop( 'use of isoneutral mixing requires that rn_beta/=0') 842 IF( lk_zdfddm ) CALL ctl_stop( 'use of double diffusive mixing requires that rn_beta/=0') 843 ENDIF 1500 844 ! 1501 845 CASE DEFAULT !== ERROR in nn_eos ==! -
branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_cen2.F90
r3876 r3894 127 127 REAL(wp) :: zupsut, zcenut, zupst ! - - 128 128 REAL(wp) :: zupsvt, zcenvt, zcent, zice ! - - 129 REAL(wp), POINTER, DIMENSION(:,: ) :: ztfreez130 REAL(wp), POINTER, DIMENSION(:,:,:) :: zwz, zind 129 REAL(wp), POINTER, DIMENSION(:,:) :: zfzp ! 2D workspace 130 REAL(wp), POINTER, DIMENSION(:,:,:) :: zwz, zind ! 3D workspace 131 131 !!---------------------------------------------------------------------- 132 132 ! 133 133 IF( nn_timing == 1 ) CALL timing_start('tra_adv_cen2') 134 134 ! 135 CALL wrk_alloc( jpi, jpj, z tfreez)135 CALL wrk_alloc( jpi, jpj, zfzp ) 136 136 CALL wrk_alloc( jpi, jpj, jpk, zwz, zind ) 137 137 ! … … 168 168 !!gm not strickly exact : the freezing point should be computed at each ocean levels... 169 169 !!gm not a big deal since cen2 is no more used in global ice-ocean simulations 170 z tfreez(:,:) = tfreez( tsn(:,:,1,jp_sal) )170 zfzp(:,:) = eos_fzp( tsn(:,:,1,jp_sal) ) 171 171 DO jk = 1, jpk 172 172 DO jj = 1, jpj 173 173 DO ji = 1, jpi 174 174 ! ! below ice covered area (if tn < "freezing"+0.1 ) 175 IF( tsn(ji,jj,jk,jp_tem) <= z tfreez(ji,jj) + 0.1 ) THEN ; zice = 1.e0176 ELSE ; zice = 0.e0175 IF( tsn(ji,jj,jk,jp_tem) <= zfzp(ji,jj) + 0.1 ) THEN ; zice = 1._wp 176 ELSE ; zice = 0._wp 177 177 ENDIF 178 178 zind(ji,jj,jk) = MAX ( & … … 280 280 ENDIF 281 281 ! 282 CALL wrk_dealloc( jpi, jpj, z tfreez)282 CALL wrk_dealloc( jpi, jpj, zfzp ) 283 283 CALL wrk_dealloc( jpi, jpj, jpk, zwz, zind ) 284 284 ! -
branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_muscl.F90
r3876 r3894 30 30 USE lib_mpp ! distribued memory computing 31 31 USE lbclnk ! ocean lateral boundary condition (or mpp link) 32 !!gm USE trabbl ! tracers: bottom boundary layer33 !!gm USE eosbn2 ! equation of state34 32 35 33 IMPLICIT NONE 36 34 PRIVATE 37 35 38 PUBLIC tra_adv_muscl ! routine called by step.F90 39 40 LOGICAL :: l_trd ! flag to compute trends 41 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: upsmsk !: mixed upstream/centered scheme near some straits 42 ! ! and in closed seas (orca 2 and 4 configurations) 43 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: xind !: mixed upstream/centered index 36 PUBLIC tra_adv_muscl ! routine called by step.F90 37 38 LOGICAL :: l_trd ! flag to compute trends 39 40 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: upsmsk !: mixed upstream/centered scheme near some straits 41 ! ! and in closed seas (orca 2 and 4 configurations) 42 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: xind !: mixed upstream/centered index 43 44 44 !! * Substitutions 45 45 # include "domzgr_substitute.h90" … … 52 52 CONTAINS 53 53 54 SUBROUTINE tra_adv_muscl( kt, kit000, cdtype, p2dt, pun, pvn, pwn, &55 & ptb, pta, kjpt, ld_msc_ups )54 SUBROUTINE tra_adv_muscl( kt, kit000, cdtype, p2dt, pun, pvn, pwn, & 55 & ptb, pta, kjpt, ld_msc_ups ) 56 56 !!---------------------------------------------------------------------- 57 57 !! *** ROUTINE tra_adv_muscl *** … … 69 69 !! IPSL, Sept. 2000 (http://www.lodyc.jussieu.fr/opa) 70 70 !!---------------------------------------------------------------------- 71 USE oce 71 USE oce, ONLY: zwx => ua , zwy => va ! (ua,va) used as workspace 72 72 ! 73 73 INTEGER , INTENT(in ) :: kt ! ocean time-step index … … 80 80 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: ptb ! before tracer field 81 81 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: pta ! tracer trend 82 83 !84 INTEGER :: ji, jj, jk, jn ! dummy loop indices82 ! 83 INTEGER :: ji, jj, jk, jn ! dummy loop indices 84 INTEGER :: ierr ! local integer 85 85 REAL(wp) :: zu, z0u, zzwx, zw ! local scalars 86 86 REAL(wp) :: zv, z0v, zzwy, z0w ! - - 87 87 REAL(wp) :: ztra, zbtr, zdt, zalpha ! - - 88 REAL(wp), POINTER, DIMENSION(:,:,:) :: zslpx, zslpy 89 INTEGER :: ierr 88 REAL(wp), POINTER, DIMENSION(:,:,:) :: zslpx, zslpy ! workspace 90 89 !!---------------------------------------------------------------------- 91 90 ! … … 94 93 CALL wrk_alloc( jpi, jpj, jpk, zslpx, zslpy ) 95 94 ! 96 97 95 IF( kt == kit000 ) THEN 98 96 IF(lwp) WRITE(numout,*) … … 121 119 122 120 ! 123 ! Upstream / centeredscheme indicator121 ! Upstream / MUSCL scheme indicator 124 122 ! ------------------------------------ 125 xind(:,:,:) = 1._wp ! set equal to 1 where up-stream is not needed126 !127 123 IF( ld_msc_ups ) THEN 128 124 DO jk = 1, jpk 129 DO jj = 1, jpj 130 DO ji = 1, jpi 131 xind(ji,jj,jk) = 1 - MAX ( & 132 rnfmsk(ji,jj) * rnfmsk_z(jk), & ! near runoff mouths (& closed sea outflows) 133 upsmsk(ji,jj) ) * tmask(ji,jj,jk) ! some of some straits 134 END DO 135 END DO 125 xind(:,:,jk) = 1._wp & ! =>1 where up-stream is not needed 126 & - MAX ( rnfmsk(:,:) * rnfmsk_z(jk), & ! =>0 near runoff mouths (& closed sea outflows) 127 & upsmsk(:,:) ) * tmask(:,:,jk) ! =>0 near some straits 136 128 END DO 137 129 ENDIF -
branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/TRA/trabbl.F90
r3876 r3894 12 12 !! - ! 2010-06 (C. Ethe, G. Madec) merge TRA-TRC 13 13 !! - ! 2010-11 (G. Madec) add mbk. arrays associated to the deepest ocean level 14 !! - ! 2013-04 (F. Roquet, G. Madec) use of eosbn2 instead of local hard coded alpha and beta 14 15 !!---------------------------------------------------------------------- 15 16 #if defined key_trabbl || defined key_esopa … … 29 30 USE eosbn2 ! equation of state 30 31 USE trd_oce ! trends: ocean variables 31 USE trdtra ! trends manager: tracers 32 USE trdtra ! trends manager: tracers 33 ! 32 34 USE iom ! IOM server 33 35 USE in_out_manager ! I/O manager … … 36 38 USE wrk_nemo ! Memory Allocation 37 39 USE timing ! Timing 38 40 USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 39 41 40 42 IMPLICIT NONE … … 49 51 LOGICAL, PUBLIC, PARAMETER :: lk_trabbl = .TRUE. !: bottom boundary layer flag 50 52 51 ! !!* Namelist nambbl *52 INTEGER , PUBLIC :: nn_bbl_ldf = 0 !: =1 : diffusive bbl or not (=0)53 INTEGER , PUBLIC :: nn_bbl_adv = 0 !: =1/2 : advective bbl or not (=0)54 ! ! =1 : advective bbl using the bottom ocean velocity55 ! ! =2 : - - using utr_bbl proportional to grad(rho)56 REAL(wp), PUBLIC :: rn_ahtbbl = 1.e3_wp !: along slope bbl diffusive coefficient [m2/s]57 REAL(wp), PUBLIC :: rn_gambbl = 10.0_wp !: lateral coeff. for bottom boundary layer scheme [s]58 59 LOGICAL , PUBLIC :: l_bbl !: flag to compute bbl diffu. flux coefand transport53 ! !!* Namelist nambbl * 54 INTEGER , PUBLIC :: nn_bbl_ldf = 0 !: =1 : diffusive bbl or not (=0) 55 INTEGER , PUBLIC :: nn_bbl_adv = 0 !: =1/2 : advective bbl or not (=0) 56 ! ! =1 : advective bbl using the bottom ocean velocity 57 ! ! =2 : - - using utr_bbl proportional to grad(rho) 58 REAL(wp), PUBLIC :: rn_ahtbbl = 1.e3_wp !: along slope bbl diffusive coefficient [m2/s] 59 REAL(wp), PUBLIC :: rn_gambbl = 10.0_wp !: lateral coeff. for bottom boundary layer scheme [s] 60 61 LOGICAL , PUBLIC :: l_bbl !: flag to compute bbl diffu. flux coef. and transport 60 62 61 63 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC :: utr_bbl , vtr_bbl ! u- (v-) transport in the bottom boundary layer … … 105 107 !!---------------------------------------------------------------------- 106 108 INTEGER, INTENT( in ) :: kt ! ocean time-step 107 ! !109 ! 108 110 REAL(wp), POINTER, DIMENSION(:,:,:) :: ztrdt, ztrds 109 111 !!---------------------------------------------------------------------- … … 111 113 IF( nn_timing == 1 ) CALL timing_start( 'tra_bbl') 112 114 ! 113 IF( l_trdtra ) THEN !* Save ta and sa trends115 IF( l_trdtra ) THEN !* Save ta and sa trends 114 116 CALL wrk_alloc( jpi, jpj, jpk, ztrdt, ztrds ) 115 117 ztrdt(:,:,:) = tsa(:,:,:,jp_tem) … … 117 119 ENDIF 118 120 119 IF( l_bbl ) CALL bbl( kt, nit000, 'TRA' ) !* bbl coef. and transport (only if not already done in trcbbl)120 121 IF( nn_bbl_ldf == 1 ) THEN !* Diffusive bbl121 IF( l_bbl ) CALL bbl( kt, nit000, 'TRA' ) !* bbl coef. and transport (only if not already done in trcbbl) 122 123 IF( nn_bbl_ldf == 1 ) THEN !* Diffusive bbl 122 124 ! 123 125 CALL tra_bbl_dif( tsb, tsa, jpts ) 124 126 IF( ln_ctl ) & 125 127 CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' bbl_ldf - Ta: ', mask1=tmask, & 126 &tab3d_2=tsa(:,:,:,jp_sal), clinfo2= ' Sa: ', mask2=tmask, clinfo3='tra' )128 & tab3d_2=tsa(:,:,:,jp_sal), clinfo2= ' Sa: ', mask2=tmask, clinfo3='tra' ) 127 129 ! lateral boundary conditions ; just need for outputs 128 130 CALL lbc_lnk( ahu_bbl, 'U', 1. ) ; CALL lbc_lnk( ahv_bbl, 'V', 1. ) … … 132 134 END IF 133 135 134 IF( nn_bbl_adv /= 0 ) THEN !* Advective bbl136 IF( nn_bbl_adv /= 0 ) THEN !* Advective bbl 135 137 ! 136 138 CALL tra_bbl_adv( tsb, tsa, jpts ) 137 139 IF(ln_ctl) & 138 140 CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' bbl_adv - Ta: ', mask1=tmask, & 139 &tab3d_2=tsa(:,:,:,jp_sal), clinfo2= ' Sa: ', mask2=tmask, clinfo3='tra' )141 & tab3d_2=tsa(:,:,:,jp_sal), clinfo2= ' Sa: ', mask2=tmask, clinfo3='tra' ) 140 142 ! lateral boundary conditions ; just need for outputs 141 143 CALL lbc_lnk( utr_bbl, 'U', 1. ) ; CALL lbc_lnk( vtr_bbl, 'V', 1. ) … … 165 167 !! advection terms. 166 168 !! 167 !! ** Method : 168 !! * diffusive bbl (nn_bbl_ldf=1) : 169 !! ** Method : * diffusive bbl only (nn_bbl_ldf=1) : 169 170 !! When the product grad( rho) * grad(h) < 0 (where grad is an 170 171 !! along bottom slope gradient) an additional lateral 2nd order … … 196 197 DO jn = 1, kjpt ! tracer loop 197 198 ! ! =========== 198 # if defined key_vectopt_loop199 DO jj = 1, 1 ! vector opt. (forced unrolling)200 DO ji = 1, jpij201 #else202 199 DO jj = 1, jpj 203 200 DO ji = 1, jpi 204 #endif 205 ik = mbkt(ji,jj) ! bottom T-level index 206 zptb(ji,jj) = ptb(ji,jj,ik,jn) ! bottom before T and S 201 zptb(ji,jj) = ptb(ji,jj,mbkt(ji,jj),jn) ! bottom before T and S 207 202 END DO 208 203 END DO 209 ! ! Compute the trend 210 # if defined key_vectopt_loop 211 DO jj = 1, 1 ! vector opt. (forced unrolling) 212 DO ji = jpi+1, jpij-jpi-1 213 # else 214 DO jj = 2, jpjm1 204 ! 205 DO jj = 2, jpjm1 ! Compute the trend 215 206 DO ji = 2, jpim1 216 # endif 217 ik = mbkt(ji,jj) ! bottom T-level index 207 ik = mbkt(ji,jj) ! bottom T-level index 218 208 zbtr = r1_e1e2t(ji,jj) / fse3t(ji,jj,ik) 219 209 pta(ji,jj,ik,jn) = pta(ji,jj,ik,jn) & … … 264 254 DO jn = 1, kjpt ! tracer loop 265 255 ! ! =========== 266 # if defined key_vectopt_loop267 DO jj = 1, 1268 DO ji = 1, jpij-jpi-1 ! vector opt. (forced unrolling)269 # else270 256 DO jj = 1, jpjm1 271 257 DO ji = 1, jpim1 ! CAUTION start from i=1 to update i=2 when cyclic east-west 272 # endif273 258 IF( utr_bbl(ji,jj) /= 0.e0 ) THEN ! non-zero i-direction bbl advection 274 259 ! down-slope i/k-indices (deep) & up-slope i/k indices (shelf) … … 333 318 !! advection terms. 334 319 !! 335 !! ** Method : 336 !! * diffusive bbl (nn_bbl_ldf=1) : 320 !! ** Method : * diffusive bbl (nn_bbl_ldf=1) : 337 321 !! When the product grad( rho) * grad(h) < 0 (where grad is an 338 322 !! along bottom slope gradient) an additional lateral 2nd order … … 342 326 !! a downslope velocity of 20 cm/s if the condition for slope 343 327 !! convection is satified) 344 !! * advective bbl (nn_bbl_adv=1 or 2) :328 !! * advective bbl (nn_bbl_adv=1 or 2) : 345 329 !! nn_bbl_adv = 1 use of the ocean velocity as bbl velocity 346 330 !! nn_bbl_adv = 2 follow Campin and Goosse (1999) implentation … … 353 337 !! Campin, J.-M., and H. Goosse, 1999, Tellus, 412-430. 354 338 !!---------------------------------------------------------------------- 355 !356 339 INTEGER , INTENT(in ) :: kt ! ocean time-step index 357 INTEGER , INTENT(in ) :: kit000 340 INTEGER , INTENT(in ) :: kit000 ! first time step index 358 341 CHARACTER(len=3), INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator) 359 !! 360 INTEGER :: ji, jj ! dummy loop indices 361 INTEGER :: ik ! local integers 362 INTEGER :: iis , iid , ijs , ijd ! - - 363 INTEGER :: ikus, ikud, ikvs, ikvd ! - - 364 REAL(wp) :: zsign, zsigna, zgbbl ! local scalars 365 REAL(wp) :: zgdrho, zt, zs, zh ! - - 366 !! 367 REAL(wp) :: fsalbt, fsbeta, pft, pfs, pfh ! statement function 368 REAL(wp), POINTER, DIMENSION(:,:) :: zub, zvb, ztb, zsb, zdep 369 !!----------------------- zv_bbl----------------------------------------------- 370 ! ratio alpha/beta = fsalbt : ratio of thermal over saline expension coefficients 371 ! ================ pft : potential temperature in degrees celcius 372 ! pfs : salinity anomaly (s-35) in psu 373 ! pfh : depth in meters 374 ! nn_eos = 0 (Jackett and McDougall 1994 formulation) 375 fsalbt( pft, pfs, pfh ) = & ! alpha/beta 376 ( ( ( -0.255019e-07 * pft + 0.298357e-05 ) * pft & 377 - 0.203814e-03 ) * pft & 378 + 0.170907e-01 ) * pft & 379 + 0.665157e-01 & 380 +(-0.678662e-05 * pfs - 0.846960e-04 * pft + 0.378110e-02 ) * pfs & 381 + ( ( - 0.302285e-13 * pfh & 382 - 0.251520e-11 * pfs & 383 + 0.512857e-12 * pft * pft ) * pfh & 384 - 0.164759e-06 * pfs & 385 +( 0.791325e-08 * pft - 0.933746e-06 ) * pft & 386 + 0.380374e-04 ) * pfh 387 fsbeta( pft, pfs, pfh ) = & ! beta 388 ( ( -0.415613e-09 * pft + 0.555579e-07 ) * pft & 389 - 0.301985e-05 ) * pft & 390 + 0.785567e-03 & 391 + ( 0.515032e-08 * pfs & 392 + 0.788212e-08 * pft - 0.356603e-06 ) * pfs & 393 +( ( 0.121551e-17 * pfh & 394 - 0.602281e-15 * pfs & 395 - 0.175379e-14 * pft + 0.176621e-12 ) * pfh & 396 + 0.408195e-10 * pfs & 397 + ( - 0.213127e-11 * pft + 0.192867e-09 ) * pft & 398 - 0.121555e-07 ) * pfh 342 ! 343 INTEGER :: ji, jj ! dummy loop indices 344 INTEGER :: ik, ikip1, ikjp1 ! local integers 345 INTEGER :: iis, iid, ikus, ikud ! - - 346 INTEGER :: ijs, ijd, ikvs, ikvd ! - - 347 REAL(wp) :: za, zt, zsign , zgbbl ! local scalars 348 REAL(wp) :: zb, zs, zsigna, zgdrho, zh ! - - 349 REAL(wp), DIMENSION(jpi,jpj,jpts) :: zts, zab ! 3D workspace 350 REAL(wp), DIMENSION(jpi,jpj) :: zub, zvb ! 2D workspace 399 351 !!---------------------------------------------------------------------- 400 352 ! 401 353 IF( nn_timing == 1 ) CALL timing_start( 'bbl') 402 !403 CALL wrk_alloc( jpi, jpj, zub, zvb, ztb, zsb, zdep )404 354 ! 405 355 IF( kt == kit000 ) THEN … … 408 358 IF(lwp) WRITE(numout,*) '~~~~~~~~~~' 409 359 ENDIF 410 ! !* bottom temperature, salinity, velocity and depth 411 #if defined key_vectopt_loop 412 DO jj = 1, 1 ! vector opt. (forced unrolling) 413 DO ji = 1, jpij 414 #else 360 ! !* bottom variables (T, S, alpha, beta, depth, velocity) 415 361 DO jj = 1, jpj 416 362 DO ji = 1, jpi 417 #endif 418 ik = mbkt(ji,jj) ! bottom T-level index419 zt b (ji,jj) = tsb(ji,jj,ik,jp_tem) * tmask(ji,jj,1) ! bottom before T and S420 z sb (ji,jj) = tsb(ji,jj,ik,jp_sal) * tmask(ji,jj,1)421 z dep(ji,jj) = fsdept_0(ji,jj,ik) ! bottom T-level reference depth363 ik = mbkt(ji,jj) ! bottom T-level index 364 zts(ji,jj,jp_tem) = tsb (ji,jj,ik,jp_tem) ! bottom before T and S 365 zts(ji,jj,jp_sal) = tsb (ji,jj,ik,jp_sal) 366 zab(ji,jj,jp_tem) = rab_b(ji,jj,ik,jp_tem) ! bottom alpha and beta 367 zab(ji,jj,jp_sal) = rab_b(ji,jj,ik,jp_sal) 422 368 ! 423 zub(ji,jj) = un(ji,jj,mbku(ji,jj)) ! bottom velocity369 zub(ji,jj) = un(ji,jj,mbku(ji,jj)) ! bottom velocity 424 370 zvb(ji,jj) = vn(ji,jj,mbkv(ji,jj)) 425 371 END DO 426 372 END DO 427 373 ! 428 374 ! !-------------------! 429 375 IF( nn_bbl_ldf == 1 ) THEN ! diffusive bbl ! 430 376 ! !-------------------! 431 377 DO jj = 1, jpjm1 ! (criteria for non zero flux: grad(rho).grad(h) < 0 ) 432 DO ji = 1, jpim1 433 ! ! i-direction 434 zt = 0.5 * ( ztb (ji,jj) + ztb (ji+1,jj) ) ! T, S anomalie, and depth 435 zs = 0.5 * ( zsb (ji,jj) + zsb (ji+1,jj) ) - 35.0 436 zh = 0.5 * ( zdep(ji,jj) + zdep(ji+1,jj) ) 437 ! ! masked bbl i-gradient of density 438 zgdrho = ( fsalbt( zt, zs, zh ) * ( ztb(ji+1,jj) - ztb(ji,jj) ) & 439 & - ( zsb(ji+1,jj) - zsb(ji,jj) ) ) * umask(ji,jj,1) 378 DO ji = 1, fs_jpim1 ! vector opt. 379 ! ! i-direction 380 za = zab(ji+1,jj,jp_tem) + zab(ji,jj,jp_tem) ! 2*(alpha,beta) at u-point 381 zb = zab(ji+1,jj,jp_sal) + zab(ji,jj,jp_sal) 382 ! ! 2*masked bottom density gradient 383 zgdrho = ( za * ( zts(ji+1,jj,jp_tem) - zts(ji,jj,jp_tem) ) & 384 & - zb * ( zts(ji+1,jj,jp_sal) - zts(ji,jj,jp_sal) ) ) * umask(ji,jj,1) 440 385 ! 441 zsign = SIGN( 0.5, -zgdrho * REAL( mgrhu(ji,jj) ) ) ! sign of ( i-gradient * i-slope )442 ahu_bbl(ji,jj) = ( 0.5 - zsign ) * ahu_bbl_0(ji,jj) 386 zsign = SIGN( 0.5, -zgdrho * REAL( mgrhu(ji,jj) ) ) ! sign of ( i-gradient * i-slope ) 387 ahu_bbl(ji,jj) = ( 0.5 - zsign ) * ahu_bbl_0(ji,jj) ! masked diffusive flux coeff. 443 388 ! 444 ! ! j-direction 445 zt = 0.5 * ( ztb (ji,jj+1) + ztb (ji,jj) ) ! T, S anomalie, and depth 446 zs = 0.5 * ( zsb (ji,jj+1) + zsb (ji,jj) ) - 35.0 447 zh = 0.5 * ( zdep(ji,jj+1) + zdep(ji,jj) ) 448 ! ! masked bbl j-gradient of density 449 zgdrho = ( fsalbt( zt, zs, zh ) * ( ztb(ji,jj+1) - ztb(ji,jj) ) & 450 & - ( zsb(ji,jj+1) - zsb(ji,jj) ) ) * vmask(ji,jj,1) 389 ! ! j-direction 390 za = zab(ji,jj+1,jp_tem) + zab(ji,jj,jp_tem) ! 2*(alpha,beta) at v-point 391 zb = zab(ji,jj+1,jp_sal) + zab(ji,jj,jp_sal) 392 ! ! 2*masked bottom density gradient 393 zgdrho = ( za * ( zts(ji,jj+1,jp_tem) - zts(ji,jj,jp_tem) ) & 394 & - zb * ( zts(ji,jj+1,jp_sal) - zts(ji,jj,jp_sal) ) ) * vmask(ji,jj,1) 451 395 ! 452 zsign 396 zsign = SIGN( 0.5, -zgdrho * REAL( mgrhv(ji,jj) ) ) ! sign of ( j-gradient * j-slope ) 453 397 ahv_bbl(ji,jj) = ( 0.5 - zsign ) * ahv_bbl_0(ji,jj) 454 !455 398 END DO 456 399 END DO … … 466 409 DO jj = 1, jpjm1 ! criteria: grad(rho).grad(h)<0 and grad(rho).grad(h)<0 467 410 DO ji = 1, fs_jpim1 ! vector opt. 468 ! ! i-direction 469 zt = 0.5 * ( ztb (ji,jj) + ztb (ji+1,jj) ) ! T, S anomalie, and depth 470 zs = 0.5 * ( zsb (ji,jj) + zsb (ji+1,jj) ) - 35.0 471 zh = 0.5 * ( zdep(ji,jj) + zdep(ji+1,jj) ) 472 ! ! masked bbl i-gradient of density 473 zgdrho = ( fsalbt( zt, zs, zh ) * ( ztb(ji+1,jj) - ztb(ji,jj) ) & 474 & - ( zsb(ji+1,jj) - zsb(ji,jj) ) ) * umask(ji,jj,1) 475 ! 476 zsign = SIGN( 0.5, - zgdrho * REAL( mgrhu(ji,jj) ) ) ! sign of i-gradient * i-slope 477 zsigna= SIGN( 0.5, zub(ji,jj) * REAL( mgrhu(ji,jj) ) ) ! sign of u * i-slope 478 ! 479 ! ! bbl velocity 411 ! ! i-direction 412 za = zab(ji+1,jj,jp_tem) + zab(ji,jj,jp_tem) ! 2*(alpha,beta) at u-point 413 zb = zab(ji+1,jj,jp_sal) + zab(ji,jj,jp_sal) 414 ! ! 2*masked bottom density gradient 415 zgdrho = ( za * ( zts(ji+1,jj,jp_tem) - zts(ji,jj,jp_tem) ) & 416 - zb * ( zts(ji+1,jj,jp_sal) - zts(ji,jj,jp_sal) ) ) * umask(ji,jj,1) 417 ! 418 zsign = SIGN( 0.5, - zgdrho * REAL( mgrhu(ji,jj) ) ) ! sign of i-gradient * i-slope 419 zsigna= SIGN( 0.5, zub(ji,jj) * REAL( mgrhu(ji,jj) ) ) ! sign of u * i-slope 420 ! 421 ! ! bbl velocity 480 422 utr_bbl(ji,jj) = ( 0.5 + zsigna ) * ( 0.5 - zsign ) * e2u(ji,jj) * e3u_bbl_0(ji,jj) * zub(ji,jj) 481 423 ! 482 ! ! j-direction 483 zt = 0.5 * ( ztb (ji,jj+1) + ztb (ji,jj) ) ! T, S anomalie, and depth 484 zs = 0.5 * ( zsb (ji,jj+1) + zsb (ji,jj) ) - 35.0 485 zh = 0.5 * ( zdep(ji,jj+1) + zdep(ji,jj) ) 486 ! ! masked bbl j-gradient of density 487 zgdrho = ( fsalbt( zt, zs, zh ) * ( ztb(ji,jj+1) - ztb(ji,jj) ) & 488 & - ( zsb(ji,jj+1) - zsb(ji,jj) ) ) * vmask(ji,jj,1) 489 zsign = SIGN( 0.5, - zgdrho * REAL( mgrhv(ji,jj) ) ) ! sign of j-gradient * j-slope 490 zsigna= SIGN( 0.5, zvb(ji,jj) * REAL( mgrhv(ji,jj) ) ) ! sign of u * i-slope 491 ! 492 ! ! bbl velocity 424 ! ! j-direction 425 za = zab(ji,jj+1,jp_tem) + zab(ji,jj,jp_tem) ! 2*(alpha,beta) at v-point 426 zb = zab(ji,jj+1,jp_sal) + zab(ji,jj,jp_sal) 427 ! ! 2*masked bottom density gradient 428 zgdrho = ( za * ( zts(ji,jj+1,jp_tem) - zts(ji,jj,jp_tem) ) & 429 & - zb * ( zts(ji,jj+1,jp_sal) - zts(ji,jj,jp_sal) ) ) * vmask(ji,jj,1) 430 zsign = SIGN( 0.5, - zgdrho * REAL( mgrhv(ji,jj) ) ) ! sign of j-gradient * j-slope 431 zsigna= SIGN( 0.5, zvb(ji,jj) * REAL( mgrhv(ji,jj) ) ) ! sign of u * i-slope 432 ! 433 ! ! bbl transport 493 434 vtr_bbl(ji,jj) = ( 0.5 + zsigna ) * ( 0.5 - zsign ) * e1v(ji,jj) * e3v_bbl_0(ji,jj) * zvb(ji,jj) 494 435 END DO … … 499 440 DO jj = 1, jpjm1 ! criteria: rho_up > rho_down 500 441 DO ji = 1, fs_jpim1 ! vector opt. 501 ! ! i-direction442 ! ! i-direction 502 443 ! down-slope T-point i/k-index (deep) & up-slope T-point i/k-index (shelf) 503 444 iid = ji + MAX( 0, mgrhu(ji,jj) ) ; iis = ji + 1 - MAX( 0, mgrhu(ji,jj) ) 504 445 ikud = mbku_d(ji,jj) ; ikus = mbku(ji,jj) 505 446 ! 506 ! ! mid-depth density anomalie (up-slope minus down-slope) 507 zt = 0.5 * ( ztb (ji,jj) + ztb (ji+1,jj) ) ! mid slope depth of T, S, and depth 508 zs = 0.5 * ( zsb (ji,jj) + zsb (ji+1,jj) ) - 35.0 509 zh = 0.5 * ( zdep(ji,jj) + zdep(ji+1,jj) ) 510 zgdrho = fsbeta( zt, zs, zh ) & 511 & * ( fsalbt( zt, zs, zh ) * ( ztb(iid,jj) - ztb(iis,jj) ) & 512 & - ( zsb(iid,jj) - zsb(iis,jj) ) ) * umask(ji,jj,1) 513 zgdrho = MAX( 0.e0, zgdrho ) ! only if shelf is denser than deep 514 ! 515 ! ! bbl transport (down-slope direction) 447 za = zab(ji+1,jj,jp_tem) + zab(ji,jj,jp_tem) ! 2*(alpha,beta) at u-point 448 zb = zab(ji+1,jj,jp_sal) + zab(ji,jj,jp_sal) 449 ! ! masked bottom density gradient 450 zgdrho = 0.5 * ( za * ( zts(iid,jj,jp_tem) - zts(iis,jj,jp_tem) ) & 451 & - zb * ( zts(iid,jj,jp_sal) - zts(iis,jj,jp_sal) ) ) * umask(ji,jj,1) 452 zgdrho = MAX( 0.e0, zgdrho ) ! only if shelf is denser than deep 453 ! 454 ! ! bbl transport (down-slope direction) 516 455 utr_bbl(ji,jj) = e2u(ji,jj) * e3u_bbl_0(ji,jj) * zgbbl * zgdrho * REAL( mgrhu(ji,jj) ) 517 456 ! 518 ! ! j-direction457 ! ! j-direction 519 458 ! down-slope T-point j/k-index (deep) & of the up -slope T-point j/k-index (shelf) 520 459 ijd = jj + MAX( 0, mgrhv(ji,jj) ) ; ijs = jj + 1 - MAX( 0, mgrhv(ji,jj) ) 521 460 ikvd = mbkv_d(ji,jj) ; ikvs = mbkv(ji,jj) 522 461 ! 523 ! ! mid-depth density anomalie (up-slope minus down-slope) 524 zt = 0.5 * ( ztb (ji,jj) + ztb (ji,jj+1) ) ! mid slope depth of T, S, and depth 525 zs = 0.5 * ( zsb (ji,jj) + zsb (ji,jj+1) ) - 35.0 526 zh = 0.5 * ( zdep(ji,jj) + zdep(ji,jj+1) ) 527 zgdrho = fsbeta( zt, zs, zh ) & 528 & * ( fsalbt( zt, zs, zh ) * ( ztb(ji,ijd) - ztb(ji,ijs) ) & 529 & - ( zsb(ji,ijd) - zsb(ji,ijs) ) ) * vmask(ji,jj,1) 530 zgdrho = MAX( 0.e0, zgdrho ) ! only if shelf is denser than deep 531 ! 532 ! ! bbl transport (down-slope direction) 462 za = zab(ji,jj+1,jp_tem) + zab(ji,jj,jp_tem) ! 2*(alpha,beta) at v-point 463 zb = zab(ji,jj+1,jp_sal) + zab(ji,jj,jp_sal) 464 ! ! masked bottom density gradient 465 zgdrho = 0.5 * ( za * ( zts(ji,ijd,jp_tem) - zts(ji,ijs,jp_tem) ) & 466 & - zb * ( zts(ji,ijd,jp_sal) - zts(ji,ijs,jp_sal) ) ) * vmask(ji,jj,1) 467 zgdrho = MAX( 0.e0, zgdrho ) ! only if shelf is denser than deep 468 ! 469 ! ! bbl transport (down-slope direction) 533 470 vtr_bbl(ji,jj) = e1v(ji,jj) * e3v_bbl_0(ji,jj) * zgbbl * zgdrho * REAL( mgrhv(ji,jj) ) 534 471 END DO … … 537 474 ! 538 475 ENDIF 539 !540 CALL wrk_dealloc( jpi, jpj, zub, zvb, ztb, zsb, zdep )541 476 ! 542 477 IF( nn_timing == 1 ) CALL timing_stop( 'bbl') … … 606 541 607 542 !* sign of grad(H) at u- and v-points 608 mgrhu(jpi,:) = 0 . ; mgrhu(:,jpj) = 0. ; mgrhv(jpi,:) = 0. ; mgrhv(:,jpj) = 0.543 mgrhu(jpi,:) = 0 ; mgrhu(:,jpj) = 0 ; mgrhv(jpi,:) = 0 ; mgrhv(:,jpj) = 0 609 544 DO jj = 1, jpjm1 610 545 DO ji = 1, jpim1 -
branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/TRA/tranpc.F90
r3876 r3894 81 81 inpci = 0 82 82 83 CALL eos( tsa, rhd, zrhop ) ! Potential density83 CALL eos( tsa, fsdept(:,:,:), rhd, zrhop ) ! Potential density 84 84 85 85 IF( l_trdtra ) THEN !* Save ta and sa trends -
branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/TRA/zpshde.F90
r3294 r3894 74 74 !! Idem for di(s) and dj(s) 75 75 !! 76 !! For rho, we call eos _insitu_2d which will compute rd~(t~,s~) at77 !! the good depth zh from interpolated T and S for the different78 !! formulationof the equation of state (eos).76 !! For rho, we call eos which will compute rd~(t~,s~) at the right 77 !! depth zh from interpolated T and S for the different formulations 78 !! of the equation of state (eos). 79 79 !! Gradient formulation for rho : 80 !! di(rho) = rd~ - rd(i,j,k) orrd(i+1,j,k) - rd~80 !! di(rho) = rd~ - rd(i,j,k) or rd(i+1,j,k) - rd~ 81 81 !! 82 82 !! ** Action : - pgtu, pgtv: horizontal gradient of tracer at u- & v-points 83 83 !! - pgru, pgrv: horizontal gradient of rho (if present) at u- & v-points 84 84 !!---------------------------------------------------------------------- 85 !86 85 INTEGER , INTENT(in ) :: kt ! ocean time-step index 87 86 INTEGER , INTENT(in ) :: kjpt ! number of tracers 88 87 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: pta ! 4D tracers fields 89 88 REAL(wp), DIMENSION(jpi,jpj, kjpt), INTENT( out) :: pgtu, pgtv ! hor. grad. of ptra at u- & v-pts 90 REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ), OPTIONAL :: prd ! 3D density anomaly fields91 REAL(wp), DIMENSION(jpi,jpj ), INTENT( out), OPTIONAL :: pgru, pgrv ! hor. grad. of prd at u- & v-pts89 REAL(wp), DIMENSION(jpi,jpj,jpk) , INTENT(in ), OPTIONAL :: prd ! 3D density anomaly fields 90 REAL(wp), DIMENSION(jpi,jpj) , INTENT( out), OPTIONAL :: pgru, pgrv ! hor. grad. of prd at u- & v-pts 92 91 ! 93 92 INTEGER :: ji, jj, jn ! Dummy loop indices 94 93 INTEGER :: iku, ikv, ikum1, ikvm1 ! partial step level (ocean bottom level) at u- and v-points 95 94 REAL(wp) :: ze3wu, ze3wv, zmaxu, zmaxv ! temporary scalars 96 REAL(wp), POINTER, DIMENSION(:,: ) :: zri, zrj, zhi, zhj97 REAL(wp), POINTER, DIMENSION(:,:,:) :: zti, ztj ! interpolated value of tracer95 REAL(wp), DIMENSION(jpi,jpj,1) :: zri, zrj, zhi, zhj ! NB: 3rd dim=1 to use eos 96 REAL(wp), DIMENSION(jpi,jpj,1,kjpt) :: zti, ztj ! 98 97 !!---------------------------------------------------------------------- 99 98 ! 100 99 IF( nn_timing == 1 ) CALL timing_start( 'zps_hde') 101 !102 CALL wrk_alloc( jpi, jpj, zri, zrj, zhi, zhj )103 CALL wrk_alloc( jpi, jpj, kjpt, zti, ztj )104 100 ! 105 101 DO jn = 1, kjpt !== Interpolation of tracers at the last ocean level ==! … … 121 117 zmaxu = ze3wu / fse3w(ji+1,jj,iku) 122 118 ! interpolated values of tracers 123 zti(ji,jj, jn) = pta(ji+1,jj,iku,jn) + zmaxu * ( pta(ji+1,jj,ikum1,jn) - pta(ji+1,jj,iku,jn) )119 zti(ji,jj,1,jn) = pta(ji+1,jj,iku,jn) + zmaxu * ( pta(ji+1,jj,ikum1,jn) - pta(ji+1,jj,iku,jn) ) 124 120 ! gradient of tracers 125 pgtu(ji,jj,jn) = umask(ji,jj,1) * ( zti(ji,jj, jn) - pta(ji,jj,iku,jn) )121 pgtu(ji,jj,jn) = umask(ji,jj,1) * ( zti(ji,jj,1,jn) - pta(ji,jj,iku,jn) ) 126 122 ELSE ! case 2 127 123 zmaxu = -ze3wu / fse3w(ji,jj,iku) 128 124 ! interpolated values of tracers 129 zti(ji,jj, jn) = pta(ji,jj,iku,jn) + zmaxu * ( pta(ji,jj,ikum1,jn) - pta(ji,jj,iku,jn) )125 zti(ji,jj,1,jn) = pta(ji,jj,iku,jn) + zmaxu * ( pta(ji,jj,ikum1,jn) - pta(ji,jj,iku,jn) ) 130 126 ! gradient of tracers 131 pgtu(ji,jj,jn) = umask(ji,jj,1) * ( pta(ji+1,jj,iku,jn) - zti(ji,jj, jn) )127 pgtu(ji,jj,jn) = umask(ji,jj,1) * ( pta(ji+1,jj,iku,jn) - zti(ji,jj,1,jn) ) 132 128 ENDIF 133 129 ! … … 136 132 zmaxv = ze3wv / fse3w(ji,jj+1,ikv) 137 133 ! interpolated values of tracers 138 ztj(ji,jj, jn) = pta(ji,jj+1,ikv,jn) + zmaxv * ( pta(ji,jj+1,ikvm1,jn) - pta(ji,jj+1,ikv,jn) )134 ztj(ji,jj,1,jn) = pta(ji,jj+1,ikv,jn) + zmaxv * ( pta(ji,jj+1,ikvm1,jn) - pta(ji,jj+1,ikv,jn) ) 139 135 ! gradient of tracers 140 pgtv(ji,jj,jn) = vmask(ji,jj,1) * ( ztj(ji,jj, jn) - pta(ji,jj,ikv,jn) )136 pgtv(ji,jj,jn) = vmask(ji,jj,1) * ( ztj(ji,jj,1,jn) - pta(ji,jj,ikv,jn) ) 141 137 ELSE ! case 2 142 138 zmaxv = -ze3wv / fse3w(ji,jj,ikv) 143 139 ! interpolated values of tracers 144 ztj(ji,jj, jn) = pta(ji,jj,ikv,jn) + zmaxv * ( pta(ji,jj,ikvm1,jn) - pta(ji,jj,ikv,jn) )140 ztj(ji,jj,1,jn) = pta(ji,jj,ikv,jn) + zmaxv * ( pta(ji,jj,ikvm1,jn) - pta(ji,jj,ikv,jn) ) 145 141 ! gradient of tracers 146 pgtv(ji,jj,jn) = vmask(ji,jj,1) * ( pta(ji,jj+1,ikv,jn) - ztj(ji,jj, jn) )142 pgtv(ji,jj,jn) = vmask(ji,jj,1) * ( pta(ji,jj+1,ikv,jn) - ztj(ji,jj,1,jn) ) 147 143 ENDIF 148 144 # if ! defined key_vectopt_loop … … 167 163 ze3wu = fse3w(ji+1,jj ,iku) - fse3w(ji,jj,iku) 168 164 ze3wv = fse3w(ji ,jj+1,ikv) - fse3w(ji,jj,ikv) 169 IF( ze3wu >= 0._wp ) THEN ; zhi(ji,jj ) = fsdept(ji ,jj,iku) ! i-direction: case 1170 ELSE ; zhi(ji,jj ) = fsdept(ji+1,jj,iku) ! - - case 2171 ENDIF 172 IF( ze3wv >= 0._wp ) THEN ; zhj(ji,jj ) = fsdept(ji,jj ,ikv) ! j-direction: case 1173 ELSE ; zhj(ji,jj ) = fsdept(ji,jj+1,ikv) ! - - case 2165 IF( ze3wu >= 0._wp ) THEN ; zhi(ji,jj,1) = fsdept(ji ,jj,iku) ! i-direction: case 1 166 ELSE ; zhi(ji,jj,1) = fsdept(ji+1,jj,iku) ! - - case 2 167 ENDIF 168 IF( ze3wv >= 0._wp ) THEN ; zhj(ji,jj,1) = fsdept(ji,jj ,ikv) ! j-direction: case 1 169 ELSE ; zhj(ji,jj,1) = fsdept(ji,jj+1,ikv) ! - - case 2 174 170 ENDIF 175 171 # if ! defined key_vectopt_loop … … 195 191 ze3wu = fse3w(ji+1,jj ,iku) - fse3w(ji,jj,iku) 196 192 ze3wv = fse3w(ji ,jj+1,ikv) - fse3w(ji,jj,ikv) 197 IF( ze3wu >= 0._wp ) THEN ; pgru(ji,jj) = umask(ji,jj,1) * ( zri(ji ,jj ) - prd(ji,jj,iku) ) ! i: 1198 ELSE ; pgru(ji,jj) = umask(ji,jj,1) * ( prd(ji+1,jj,iku) - zri(ji,jj ) ) ! i: 2199 ENDIF 200 IF( ze3wv >= 0._wp ) THEN ; pgrv(ji,jj) = vmask(ji,jj,1) * ( zrj(ji,jj ) - prd(ji,jj,ikv) ) ! j: 1201 ELSE ; pgrv(ji,jj) = vmask(ji,jj,1) * ( prd(ji,jj+1,ikv) - zrj(ji,jj ) ) ! j: 2193 IF( ze3wu >= 0._wp ) THEN ; pgru(ji,jj) = umask(ji,jj,1) * ( zri(ji ,jj, 1 ) - prd(ji,jj,iku) ) ! i: 1 194 ELSE ; pgru(ji,jj) = umask(ji,jj,1) * ( prd(ji+1,jj,iku) - zri(ji,jj, 1 ) ) ! i: 2 195 ENDIF 196 IF( ze3wv >= 0._wp ) THEN ; pgrv(ji,jj) = vmask(ji,jj,1) * ( zrj(ji,jj , 1 ) - prd(ji,jj,ikv) ) ! j: 1 197 ELSE ; pgrv(ji,jj) = vmask(ji,jj,1) * ( prd(ji,jj+1,ikv) - zrj(ji,jj, 1 ) ) ! j: 2 202 198 ENDIF 203 199 # if ! defined key_vectopt_loop … … 209 205 END IF 210 206 ! 211 CALL wrk_dealloc( jpi, jpj, zri, zrj, zhi, zhj )212 CALL wrk_dealloc( jpi, jpj, kjpt, zti, ztj )213 !214 207 IF( nn_timing == 1 ) CALL timing_stop( 'zps_hde') 215 208 ! -
branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/TRD/trdglo.F90
r3876 r3894 212 212 zkepe(:,:,:) = 0._wp 213 213 214 CALL eos( tsn, rhd, rhop ) ! now potential and in situ densities214 CALL eos( tsn, fsdept(:,:,:), rhd, rhop ) ! now potential and in situ densities 215 215 216 216 zcof = 0.5_wp / rau0 ! Density flux at w-point -
branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfddm.F90
r3610 r3894 6 6 !! History : OPA ! 2000-08 (G. Madec) double diffusive mixing 7 7 !! NEMO 1.0 ! 2002-06 (G. Madec) F90: Free form and module 8 !! 3.3 ! 2010-10 (C. Ethe, G. Madec) reorganisation of initialisation phase 8 !! 3.3 ! 2010-10 (C. Ethe, G. Madec) reorganisation of initialisation phase 9 !! 3.6 ! 2013-04 (G. Madec, F. Roquet) zrau compute locally using interpolation of alpha & beta 9 10 !!---------------------------------------------------------------------- 10 11 #if defined key_zdfddm || defined key_esopa … … 18 19 USE dom_oce ! ocean space and time domain variables 19 20 USE zdf_oce ! ocean vertical physics variables 21 ! 20 22 USE in_out_manager ! I/O manager 21 23 USE lbclnk ! ocean lateral boundary conditions (or mpp link) … … 34 36 LOGICAL , PUBLIC, PARAMETER :: lk_zdfddm = .TRUE. !: double diffusive mixing flag 35 37 36 REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:) :: avs !: salinity vertical diffusivity coeff. at w-point 37 REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:) :: rrau !: heat/salt buoyancy flux ratio 38 REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:) :: avs !: salinity vertical diffusivity coeff. at w-point 38 39 39 40 ! !!* Namelist namzdf_ddm : double diffusive mixing * … … 42 43 43 44 !! * Substitutions 45 # include "domzgr_substitute.h90" 44 46 # include "vectopt_loop_substitute.h90" 45 47 !!---------------------------------------------------------------------- 46 !! NEMO/OPA 4.0 , NEMO Consortium (2011)48 !! NEMO/OPA 3.5 , NEMO Consortium (2013) 47 49 !! $Id$ 48 50 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) … … 54 56 !! *** ROUTINE zdf_ddm_alloc *** 55 57 !!---------------------------------------------------------------------- 56 ALLOCATE( avs(jpi,jpj,jpk) , rrau(jpi,jpj,jpk), STAT= zdf_ddm_alloc )58 ALLOCATE( avs(jpi,jpj,jpk) , STAT= zdf_ddm_alloc ) 57 59 ! 58 60 IF( lk_mpp ) CALL mpp_sum ( zdf_ddm_alloc ) … … 71 73 !! diffusive mixing (i.e. salt fingering and diffusive layering) 72 74 !! following Merryfield et al. (1999). The rate of double diffusive 73 !! mixing depend on the buoyancy ratio: Rrau=alpha/beta dk[T]/dk[S] 74 !! which is computed in rn2.F 75 !! mixing depend on the buoyancy ratio (R=alpha/beta dk[T]/dk[S]): 75 76 !! * salt fingering (Schmitt 1981): 76 !! for R rau > 1 and rn2 > 0 : zavfs = rn_avts / ( 1 + (Rrau/rn_hsbfr)^6 )77 !! for R rau> 1 and rn2 > 0 : zavfs = O78 !! otherwise : zavft = 0.7 zavs / R rau77 !! for R > 1 and rn2 > 0 : zavfs = rn_avts / ( 1 + (R/rn_hsbfr)^6 ) 78 !! for R > 1 and rn2 > 0 : zavfs = O 79 !! otherwise : zavft = 0.7 zavs / R 79 80 !! * diffusive layering (Federov 1988): 80 !! for 0< Rrau < 1 and rn2 > 0 : zavdt = 1.3635e-6 81 !! * exp( 4.6 exp(-0.54 (1/Rrau-1) ) ) 81 !! for 0< R < 1 and N^2 > 0 : zavdt = 1.3635e-6 * exp( 4.6 exp(-0.54 (1/R-1) ) ) 82 82 !! otherwise : zavdt = 0 83 !! for .5 < R rau < 1 and rn2 > 0 : zavds = zavdt (1.885 Rrau-0.85)84 !! for 0 < R rau <.5 and rn2 > 0 : zavds = zavdt 0.15 Rrau83 !! for .5 < R < 1 and N^2 > 0 : zavds = zavdt (1.885 R -0.85) 84 !! for 0 < R <.5 and N^2 > 0 : zavds = zavdt 0.15 R 85 85 !! otherwise : zavds = 0 86 86 !! * update the eddy diffusivity: … … 96 96 ! 97 97 INTEGER :: ji, jj , jk ! dummy loop indices 98 REAL(wp) :: zinr, zrr ! temporary scalars 99 REAL(wp) :: zavft, zavfs ! - - 100 REAL(wp) :: zavdt, zavds ! - - 101 REAL(wp), POINTER, DIMENSION(:,:) :: zmsks, zmskf, zmskd1, zmskd2, zmskd3 98 REAL(wp) :: zaw, zbw, zrw ! local scalars 99 REAL(wp) :: zdt, zds 100 REAL(wp) :: zinr, zrr ! - - 101 REAL(wp) :: zavft, zavfs ! - - 102 REAL(wp) :: zavdt, zavds ! - - 103 REAL(wp), POINTER, DIMENSION(:,:) :: zrau, zmsks, zmskf, zmskd1, zmskd2, zmskd3 102 104 !!---------------------------------------------------------------------- 103 105 ! 104 106 IF( nn_timing == 1 ) CALL timing_start('zdf_ddm') 105 107 ! 106 CALL wrk_alloc( jpi,jpj, z msks, zmskf, zmskd1, zmskd2, zmskd3 )108 CALL wrk_alloc( jpi,jpj, zrau, zmsks, zmskf, zmskd1, zmskd2, zmskd3 ) 107 109 108 110 ! ! =============== … … 111 113 ! Define the mask 112 114 ! --------------- 113 rrau(:,:,jk) = MAX( 1.e-20, rrau(:,:,jk) ) ! only retains positive value of rrau 115 DO jj = 1, jpj ! R=zrau = (alpha / beta) (dk[t] / dk[s]) 116 DO ji = 1, jpi 117 zrw = ( fsdepw(ji,jj,jk ) - fsdept(ji,jj,jk) ) & 118 & / ( fsdept(ji,jj,jk-1) - fsdept(ji,jj,jk) ) 119 ! 120 zaw = ( rab_n(ji,jj,jk,jp_tem) * (1. - zrw) + rab_n(ji,jj,jk-1,jp_tem) * zrw ) * tmask(ji,jj,jk) 121 zbw = ( rab_n(ji,jj,jk,jp_sal) * (1. - zrw) + rab_n(ji,jj,jk-1,jp_sal) * zrw ) * tmask(ji,jj,jk) 122 ! 123 zdt = zaw * ( tsn(ji,jj,jk-1,jp_tem) - tsn(ji,jj,jk,jp_tem) ) 124 zds = zbw * ( tsn(ji,jj,jk-1,jp_sal) - tsn(ji,jj,jk,jp_sal) ) 125 IF( ABS( zds) <= 1.e-20_wp ) zds = 1.e-20_wp 126 zrau(ji,jj) = MAX( 1.e-20, zdt / zds ) ! only retains positive value of zrau 127 END DO 128 END DO 114 129 115 130 DO jj = 1, jpj ! indicators: … … 119 134 ELSE ; zmsks(ji,jj) = 1._wp 120 135 ENDIF 121 ! salt fingering indicator: msksf=1 if rrau>1; 0 elsewhere122 IF( rrau(ji,jj,jk) <= 1.) THEN ; zmskf(ji,jj) = 0._wp136 ! salt fingering indicator: msksf=1 if R>1; 0 elsewhere 137 IF( zrau(ji,jj) <= 1. ) THEN ; zmskf(ji,jj) = 0._wp 123 138 ELSE ; zmskf(ji,jj) = 1._wp 124 139 ENDIF 125 140 ! diffusive layering indicators: 126 ! ! mskdl1=1 if 0< rrau<1; 0 elsewhere127 IF( rrau(ji,jj,jk) >= 1.) THEN ; zmskd1(ji,jj) = 0._wp141 ! ! mskdl1=1 if 0< R <1; 0 elsewhere 142 IF( zrau(ji,jj) >= 1. ) THEN ; zmskd1(ji,jj) = 0._wp 128 143 ELSE ; zmskd1(ji,jj) = 1._wp 129 144 ENDIF 130 ! ! mskdl2=1 if 0< rrau<0.5; 0 elsewhere131 IF( rrau(ji,jj,jk) >= 0.5) THEN ; zmskd2(ji,jj) = 0._wp145 ! ! mskdl2=1 if 0< R <0.5; 0 elsewhere 146 IF( zrau(ji,jj) >= 0.5 ) THEN ; zmskd2(ji,jj) = 0._wp 132 147 ELSE ; zmskd2(ji,jj) = 1._wp 133 148 ENDIF 134 ! mskdl3=1 if 0.5< rrau<1; 0 elsewhere135 IF( rrau(ji,jj,jk) <= 0.5 .OR. rrau(ji,jj,jk) >= 1. ) THEN ; zmskd3(ji,jj) = 0._wp136 ELSE 149 ! mskdl3=1 if 0.5< R <1; 0 elsewhere 150 IF( zrau(ji,jj) <= 0.5 .OR. zrau(ji,jj) >= 1. ) THEN ; zmskd3(ji,jj) = 0._wp 151 ELSE ; zmskd3(ji,jj) = 1._wp 137 152 ENDIF 138 153 END DO … … 149 164 !CDIR NOVERRCHK 150 165 DO ji = 1, jpi 151 zinr = 1. /rrau(ji,jj,jk)166 zinr = 1._wp / zrau(ji,jj) 152 167 ! salt fingering 153 zrr = rrau(ji,jj,jk)/rn_hsbfr168 zrr = zrau(ji,jj) / rn_hsbfr 154 169 zrr = zrr * zrr 155 170 zavfs = rn_avts / ( 1 + zrr*zrr*zrr ) * zmsks(ji,jj) * zmskf(ji,jj) … … 157 172 ! diffusive layering 158 173 zavdt = 1.3635e-6 * EXP( 4.6 * EXP( -0.54*(zinr-1.) ) ) * zmsks(ji,jj) * zmskd1(ji,jj) 159 zavds = zavdt * zmsks(ji,jj) * ( ( 1.85 * rrau(ji,jj,jk) - 0.85 ) * zmskd3(ji,jj) &160 & + 0.15 * rrau(ji,jj,jk) * zmskd2(ji,jj) )174 zavds = zavdt * zmsks(ji,jj) * ( ( 1.85 * zrau(ji,jj) - 0.85 ) * zmskd3(ji,jj) & 175 & + 0.15 * zrau(ji,jj) * zmskd2(ji,jj) ) 161 176 ! add to the eddy viscosity coef. previously computed 162 177 avs (ji,jj,jk) = avt(ji,jj,jk) + zavfs + zavds -
branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfkpp.F90
r3792 r3894 26 26 USE phycst ! physical constants 27 27 USE eosbn2 ! equation of state 28 USE zdfddm ! double diffusion mixing 28 USE zdfddm ! double diffusion mixing (avs array) 29 29 USE in_out_manager ! I/O manager 30 30 USE lib_mpp ! MPP library … … 32 32 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 33 33 USE prtctl ! Print control 34 USE trd mod_oce! ocean trends definition34 USE trd_oce ! ocean trends definition 35 35 USE trdtra ! tracers trends 36 36 USE timing ! Timing … … 246 246 REAL(wp) :: zdelta, zdelta2, zdzup, zdzdn, zdzh, zvath, zgat1, zdat1, zkm1m, zkm1t 247 247 #if defined key_zdfddm 248 REAL(wp) :: zrrau, zds, zavdds, zavddt,zinr ! double diffusion mixing 249 REAL(wp), POINTER, DIMENSION(:,:) :: zdifs 248 REAL(wp) :: zrw, zkm1s ! local scalars 249 REAL(wp) :: zrrau, zdt, zds, zavdds, zavddt, zinr ! double diffusion mixing 250 REAL(wp), POINTER, DIMENSION(:,:) :: zdifs 250 251 REAL(wp), POINTER, DIMENSION(:) :: za2s, za3s, zkmps 251 REAL(wp) :: zkm1s252 252 REAL(wp), POINTER, DIMENSION(:,:) :: zblcs 253 253 REAL(wp), POINTER, DIMENSION(:,:,:) :: zdiffus … … 332 332 avt (ji,jj,jk) = avt (ji,jj,jk) + rn_difri * zfri 333 333 ENDIF 334 ! 334 335 #if defined key_zdfddm 335 avs (ji,jj,jk) = avt (ji,jj,jk) 336 ! Double diffusion mixing ; NOT IN ROUTINE ZDFDDM.F90 337 ! ------------------------------------------------------------------ 338 ! only retains positive value of rrau 339 zrrau = MAX( rrau(ji,jj,jk), epsln ) 340 zds = tsn(ji,jj,jk-1,jp_sal) - tsn(ji,jj,jk,jp_sal) 341 IF( zrrau > 1. .AND. zds > 0.) THEN 342 ! 343 ! Salt fingering case. 344 !--------------------- 345 ! Compute interior diffusivity for double diffusive mixing of 346 ! salinity. Upper bound "zrrau" by "Rrho0"; (Rrho0=1.9, difcoefnuf=0.001). 347 ! After that set interior diffusivity for double diffusive mixing 348 ! of temperature 336 ! 337 ! Double diffusion mixing ; NOT in ROUTINE zdfddm.F90 338 ! ------------------------- 339 avs (ji,jj,jk) = avt (ji,jj,jk) 340 341 ! R=zrau = (alpha / beta) (dk[t] / dk[s]) 342 zrw = ( fsdepw(ji,jj,jk ) - fsdept(ji,jj,jk) ) & 343 & / ( fsdept(ji,jj,jk-1) - fsdept(ji,jj,jk) ) 344 ! 345 zaw = ( rab_n(ji,jj,jk,jp_tem) * (1. - zrw) + rab_n(ji,jj,jk-1,jp_tem) * zrw ) * tmask(ji,jj,jk) 346 zbw = ( rab_n(ji,jj,jk,jp_sal) * (1. - zrw) + rab_n(ji,jj,jk-1,jp_sal) * zrw ) * tmask(ji,jj,jk) 347 ! 348 zdt = zaw * ( tsn(ji,jj,jk-1,jp_tem) - tsn(ji,jj,jk,jp_tem) ) 349 zds = zbw * ( tsn(ji,jj,jk-1,jp_sal) - tsn(ji,jj,jk,jp_sal) ) 350 IF( ABS( zds) <= 1.e-20_wp ) zds = 1.e-20_wp 351 zrrau = MAX( epsln , zdt / zds ) ! only retains positive value of zrau 352 ! 353 IF( zrrau > 1. .AND. zds > 0.) THEN ! Salt fingering case. 354 ! !--------------------- 355 ! Compute interior diffusivity for double diffusive mixing of salinity. 356 ! Upper bound "zrrau" by "Rrho0"; (Rrho0=1.9, difcoefnuf=0.001). 357 ! After that set interior diffusivity for double diffusive mixing of temperature 349 358 zavdds = MIN( zrrau, Rrho0 ) 350 359 zavdds = ( zavdds - 1.0 ) / ( Rrho0 - 1.0 ) … … 353 362 zavdds = difssf * zavdds 354 363 zavddt = 0.7 * zavdds 355 ELSEIF( zrrau < 1. .AND. zrrau > 0. .AND. zds < 0.) THEN356 364 ! 357 ! Diffusive convection case. 358 !--------------------------- 359 ! Compute interior diffusivity for double diffusive mixing of 360 ! temperature (Marmorino and Caldwell, 1976); 365 ELSEIF( zrrau < 1. .AND. zrrau > 0. .AND. zds < 0.) THEN ! Diffusive convection case. 366 ! !--------------------------- 367 ! Compute interior diffusivity for double diffusive mixing of temperature (Marmorino and Caldwell, 1976); 361 368 ! Compute interior diffusivity for double diffusive mixing of salinity 362 369 zinr = 1. / zrrau 363 370 zavddt = 0.909 * EXP( 4.6 * EXP( -0.54* ( zinr - 1. ) ) ) 364 371 zavddt = difsdc * zavddt 365 IF( zrrau < 0.5) THEN 366 zavdds = zavddt * 0.15 * zrrau 367 ELSE 368 zavdds = zavddt * (1.85 * zrrau - 0.85 ) 372 IF( zrrau < 0.5) THEN ; zavdds = zavddt * 0.15 * zrrau 373 ELSE ; zavdds = zavddt * (1.85 * zrrau - 0.85 ) 369 374 ENDIF 370 375 ELSE … … 385 390 !--------------------------------------------------------------------- 386 391 DO jj = 2, jpjm1 387 DO ji = fs_2, fs_jpim1 388 IF( nn_eos < 1) THEN 389 zt = tsn(ji,jj,1,jp_tem) 390 zs = tsn(ji,jj,1,jp_sal) - 35.0 391 zh = fsdept(ji,jj,1) 392 ! potential volumic mass 393 zrhos = rhop(ji,jj,1) 394 zalbet = ( ( ( - 0.255019e-07 * zt + 0.298357e-05 ) * zt & ! ratio alpha/beta 395 & - 0.203814e-03 ) * zt & 396 & + 0.170907e-01 ) * zt & 397 & + 0.665157e-01 & 398 & + ( - 0.678662e-05 * zs & 399 & - 0.846960e-04 * zt + 0.378110e-02 ) * zs & 400 & + ( ( - 0.302285e-13 * zh & 401 & - 0.251520e-11 * zs & 402 & + 0.512857e-12 * zt * zt ) * zh & 403 & - 0.164759e-06 * zs & 404 & +( 0.791325e-08 * zt - 0.933746e-06 ) * zt & 405 & + 0.380374e-04 ) * zh 406 407 zbeta = ( ( -0.415613e-09 * zt + 0.555579e-07 ) * zt & ! beta 408 & - 0.301985e-05 ) * zt & 409 & + 0.785567e-03 & 410 & + ( 0.515032e-08 * zs & 411 & + 0.788212e-08 * zt - 0.356603e-06 ) * zs & 412 & +( ( 0.121551e-17 * zh & 413 & - 0.602281e-15 * zs & 414 & - 0.175379e-14 * zt + 0.176621e-12 ) * zh & 415 & + 0.408195e-10 * zs & 416 & + ( - 0.213127e-11 * zt + 0.192867e-09 ) * zt & 417 & - 0.121555e-07 ) * zh 418 419 zthermal = zbeta * zalbet / ( rcp * zrhos + epsln ) 420 zhalin = zbeta * tsn(ji,jj,1,jp_sal) * rcs 421 ELSE 422 zrhos = rhop(ji,jj,1) + rau0 * ( 1. - tmask(ji,jj,1) ) 423 zthermal = rn_alpha / ( rcp * zrhos + epsln ) 424 zhalin = rn_beta * tsn(ji,jj,1,jp_sal) * rcs 425 zbeta = rn_beta 426 ENDIF 392 DO ji = fs_2, fs_jpim1 393 zrhos = rau0 * ( 1._wp + rhd(ji,jj,1) ) * tmask(ji,jj,1) 394 zthermal = rab_n(ji,jj,1,jp_tem) / ( rcp * zrhos + epsln ) 395 zbeta = rab_n(ji,jj,1,jp_sal) 396 zhalin = zbeta * tsn(ji,jj,1,jp_sal) * rcs 397 ! 427 398 ! Radiative surface buoyancy force 428 399 zBosol(ji,jj) = grav * zthermal * qsr(ji,jj) … … 435 406 ws0(ji,jj) = - ( ( emp(ji,jj)-rnf(ji,jj) ) * tsn(ji,jj,1,jp_sal) & 436 407 & + sfx(ji,jj) ) * rcs * tmask(ji,jj,1) 437 END DO438 END DO408 END DO 409 END DO 439 410 440 411 zflageos = 0.5 + SIGN( 0.5, nn_eos - 1. ) … … 447 418 ! Friction velocity (zustar), at T-point : LMD94 eq. 2 448 419 zustar(ji,jj) = SQRT( taum(ji,jj) / ( zrhos + epsln ) ) 449 END DO450 END DO420 END DO 421 END DO 451 422 452 423 !CDIR NOVERRCHK … … 1270 1241 ztrds(:,:,:) = tsa(:,:,:,jp_sal) - ztrds(:,:,:) 1271 1242 !!bug gm jpttdzdf ==> jpttkpp 1272 CALL trd_tra( kt, 'TRA', jp_tem, jptra_ trd_zdf, ztrdt )1273 CALL trd_tra( kt, 'TRA', jp_sal, jptra_ trd_zdf, ztrds )1243 CALL trd_tra( kt, 'TRA', jp_tem, jptra_zdf, ztrdt ) 1244 CALL trd_tra( kt, 'TRA', jp_sal, jptra_zdf, ztrds ) 1274 1245 DEALLOCATE( ztrdt ) ; DEALLOCATE( ztrds ) 1275 1246 ENDIF -
branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftmx.F90
r3625 r3894 435 435 436 436 ztpc = 0.e0 437 zpc(:,:,:) = MAX( rn_n2min,rn2(:,:,:)) * zav_tide(:,:,:)437 zpc(:,:,:) = MAX( rn_n2min , rn2(:,:,:) ) * zav_tide(:,:,:) 438 438 DO jk= 2, jpkm1 439 439 DO jj = 1, jpj -
branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/oce.F90
r3893 r3894 29 29 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: rn2b , rn2 !: brunt-vaisala frequency**2 [s-2] 30 30 ! 31 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: rhd !: in situ density anomalierhd=(rho-rau0)/rau0 [no units]31 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: rhd !: density anomaly rhd=(rho-rau0)/rau0 [no units] 32 32 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: rhop !: potential volumic mass [kg/m3] 33 33 -
branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/step.F90
r3893 r3894 38 38 !! * Substitutions 39 39 # include "domzgr_substitute.h90" 40 # include "zdfddm_substitute.h90" 41 !!---------------------------------------------------------------------- 42 !! NEMO/OPA 3.3 , NEMO Consortium (2010) 40 !!---------------------------------------------------------------------- 41 !! NEMO/OPA 3.5 , NEMO Consortium (2010) 43 42 !! $Id$ 44 43 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) … … 104 103 ! Ocean physics update (ua, va, tsa used as workspace) 105 104 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 106 CALL eos( tsb, rab_b ) ! before thermal/haline expansion coef. 107 CALL bn2( tsb, rab_b, rn2b ) ! before Brunt-Vaisala frequency 108 CALL eos( tsn, rab_n ) ! now thermal/haline expansion coef. 109 CALL bn2( tsn, rab_n, rn2 ) ! now Brunt-Vaisala frequency 110 105 ! ! thermal/haline expansion coef. & Brunt-Vaisala frequency 106 CALL eos( tsb, rab_b, rn2b ) ! before 107 CALL eos( tsn, rab_n, rn2 ) ! now 111 108 ! 112 109 ! VERTICAL PHYSICS … … 142 139 ! 143 140 IF( lk_ldfslp ) THEN ! slope of lateral mixing 144 CALL eos( tsb, rhd )! before in situ density141 CALL eos( tsb, fsdept(:,:,:), rhd ) ! before in situ density 145 142 IF( ln_zps ) CALL zps_hde( kstp, jpts, tsb, gtsu, gtsv, & ! Partial steps: before horizontal gradient 146 143 & rhd, gru , grv ) ! of t, s, rd at the last ocean level … … 205 202 IF( ln_zdfnpc ) CALL tra_npc( kstp ) ! update after fields by non-penetrative convection 206 203 CALL tra_nxt( kstp ) ! tracer fields at next time step 207 CALL eos ( tsa, rhd, rhop )! Time-filtered in situ density for hpg computation204 CALL eos ( tsa, fsdept(:,:,:), rhd, rhop ) ! Time-filtered in situ density for hpg computation 208 205 IF( ln_zps ) CALL zps_hde( kstp, jpts, tsa, gtsu, gtsv, & ! zps: time filtered hor. derivative 209 206 & rhd, gru , grv ) ! of t, s, rd at the last ocean level 210 207 211 208 ELSE ! centered hpg (eos then time stepping) 212 CALL eos ( tsn, rhd, rhop )! now in situ density for hpg computation209 CALL eos ( tsn, fsdept(:,:,:), rhd, rhop ) ! now in situ density for hpg computation 213 210 IF( ln_zps ) CALL zps_hde( kstp, jpts, tsn, gtsu, gtsv, & ! zps: now hor. derivative 214 211 & rhd, gru , grv ) ! of t, s, rd at the last ocean level
Note: See TracChangeset
for help on using the changeset viewer.