Changeset 3876 for branches/2013/dev_r3858_CNRS3_Ediag
- Timestamp:
- 2013-04-19T08:48:21+02:00 (11 years ago)
- Location:
- branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM
- Files:
-
- 9 added
- 7 deleted
- 42 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/CONFIG/ORCA2_LIM/EXP00/namelist
r3795 r3876 848 848 !!====================================================================== 849 849 !! namnc4 netcdf4 chunking and compression settings ("key_netcdf4") 850 !! namtrd dynamics and/or tracer trends ("key_trddyn","key_trdtra","key_trdmld")850 !! namtrd dynamics and/or tracer trends 851 851 !! namflo float parameters ("key_float") 852 852 !! namptr Poleward Transport Diagnostics … … 866 866 / 867 867 !----------------------------------------------------------------------- 868 &namtrd ! diagnostics on dynamics and/or tracer trends ("key_trddyn" and/or "key_trdtra") 869 ! ! or mixed-layer trends or barotropic vorticity ("key_trdmld" or "key_trdvor") 870 !----------------------------------------------------------------------- 871 nn_trd = 365 ! time step frequency dynamics and tracers trends 868 &namtrd ! diagnostics on dynamics and/or tracer trends 869 ! ! or mixed-layer trends or barotropic vorticity 870 !----------------------------------------------------------------------- 871 ln_glo_trd = .FALSE. ! (T) global domain averaged diag for T, T^2, KE, and PE 872 ln_dyn_trd = .FALSE. ! (T) 3D momentum trend output 873 ln_dyn_mld = .FALSE. ! (T) 2D momentum trends averaged over the mixed layer 874 ln_vor_trd = .FALSE. ! (T) 2D barotropic vorticity trends 875 ln_KE_trd = .FALSE. ! (T) 3D Kinetic Energy trends 876 ln_PE_trd = .FALSE. ! (T) 3D Potential Energy trends 877 ln_tra_trd = .FALSE. ! (T) 3D tracer trend output 878 ln_tra_mld = .FALSE. ! (T) 2D tracer trends averaged over the mixed layer 879 nn_trd = 365 ! print frequency (ln_glo_trd=T) (unit=time step) 880 / 881 !----------------------------------------------------------------------- 882 &namtrd_mxl ! mixed layer averaged trends diagnosed on dynamics and/or tracer trends 883 !----------------------------------------------------------------------- 884 nn_trd = 365 ! print frequency (ln_glo_trd=T) (unit=time step) 885 rn_zho_c = 0.01 ! density criteria used to define the MLD (do not change unless changing the value used in zdfmxl) 872 886 nn_ctls = 0 ! control surface type in mixed-layer trends (0,1 or n<jpk) 873 887 rn_ucf = 1. ! unit conversion factor (=1 -> /seconds ; =86400. -> /day) -
branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/DYN/dynadv_cen2.F90
r3294 r3876 15 15 USE oce ! ocean dynamics and tracers 16 16 USE dom_oce ! ocean space and time domain 17 USE trd mod_oce ! ocean variables trends18 USE trd mod ! ocean dynamics trends17 USE trd_oce ! trends: ocean variables 18 USE trddyn ! trend manager: dynamics 19 19 USE in_out_manager ! I/O manager 20 20 USE lib_mpp ! MPP library … … 103 103 zfu_uw(:,:,:) = ua(:,:,:) - zfu_uw(:,:,:) 104 104 zfv_vw(:,:,:) = va(:,:,:) - zfv_vw(:,:,:) 105 CALL trd_ mod( zfu_uw, zfv_vw, jpdyn_trd_had, 'DYN', kt )105 CALL trd_dyn( zfu_uw, zfv_vw, jpdyn_keg, kt ) 106 106 zfu_t(:,:,:) = ua(:,:,:) 107 107 zfv_t(:,:,:) = va(:,:,:) … … 153 153 zfu_t(:,:,:) = ua(:,:,:) - zfu_t(:,:,:) 154 154 zfv_t(:,:,:) = va(:,:,:) - zfv_t(:,:,:) 155 CALL trd_ mod( zfu_t, zfv_t, jpdyn_trd_zad, 'DYN', kt )155 CALL trd_dyn( zfu_t, zfv_t, jpdyn_zad, kt ) 156 156 ENDIF 157 157 ! ! Control print -
branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/DYN/dynadv_ubs.F90
r3294 r3876 16 16 USE oce ! ocean dynamics and tracers 17 17 USE dom_oce ! ocean space and time domain 18 USE trd mod ! ocean dynamics trends19 USE trd mod_oce ! ocean variables trends18 USE trd_oce ! trends: ocean variables 19 USE trddyn ! trend manager: dynamics 20 20 USE in_out_manager ! I/O manager 21 21 USE prtctl ! Print control … … 196 196 zfu_uw(:,:,:) = ua(:,:,:) - zfu_uw(:,:,:) 197 197 zfv_vw(:,:,:) = va(:,:,:) - zfv_vw(:,:,:) 198 CALL trd_ mod( zfu_uw, zfv_vw, jpdyn_trd_had, 'DYN', kt )198 CALL trd_dyn( zfu_uw, zfv_vw, jpdyn_keg, kt ) 199 199 zfu_t(:,:,:) = ua(:,:,:) 200 200 zfv_t(:,:,:) = va(:,:,:) … … 245 245 zfu_t(:,:,:) = ua(:,:,:) - zfu_t(:,:,:) 246 246 zfv_t(:,:,:) = va(:,:,:) - zfv_t(:,:,:) 247 CALL trd_ mod( zfu_t, zfv_t, jpdyn_trd_zad, 'DYN', kt )247 CALL trd_dyn( zfu_t, zfv_t, jpdyn_zad, kt ) 248 248 ENDIF 249 249 ! ! Control print -
branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/DYN/dynbfr.F90
r3294 r3876 10 10 11 11 !!---------------------------------------------------------------------- 12 !! dyn_bfr : Update the momentum trend with the bottom friction contribution12 !! dyn_bfr : Update the momentum trend with the bottom friction contribution 13 13 !!---------------------------------------------------------------------- 14 USE oce 15 USE dom_oce 16 USE zdf_oce 17 USE zdfbfr 18 USE trd mod ! ocean active dynamics and tracers trends19 USE trd mod_oce ! ocean variables trends20 USE in_out_manager 21 USE prtctl 22 USE timing 23 USE wrk_nemo 14 USE oce ! ocean dynamics and tracers variables 15 USE dom_oce ! ocean space and time domain variables 16 USE zdf_oce ! ocean vertical physics variables 17 USE zdfbfr ! ocean bottom friction variables 18 USE trd_oce ! trends: ocean variables 19 USE trddyn ! trend manager: dynamics 20 USE in_out_manager ! I/O manager 21 USE prtctl ! Print control 22 USE timing ! Timing 23 USE wrk_nemo ! Memory Allocation 24 24 25 25 IMPLICIT NONE 26 26 PRIVATE 27 27 28 PUBLIC dyn_bfr 28 PUBLIC dyn_bfr ! routine called by step.F90 29 29 30 30 !! * Substitutions … … 57 57 IF( nn_timing == 1 ) CALL timing_start('dyn_bfr') 58 58 ! 59 !!gm issue: better to put the logical in step to control the call of zdf_bfr 60 !! ==> change the logical from ln_bfrimp to ln_bfr_exp !! 59 61 IF( .NOT.ln_bfrimp) THEN ! only for explicit bottom friction form 60 62 ! implicit bfr is implemented in dynzdf_imp 61 63 64 !!gm bug : time step is only rdt (not 2 rdt if euler start !) 62 65 zm1_2dt = - 1._wp / ( 2._wp * rdt ) 63 66 … … 89 92 ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 90 93 ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 91 CALL trd_ mod( ztrdu(:,:,:), ztrdv(:,:,:), jpdyn_trd_bfr, 'DYN', kt )94 CALL trd_dyn( ztrdu(:,:,:), ztrdv(:,:,:), jpdyn_bfr, kt ) 92 95 CALL wrk_dealloc( jpi,jpj,jpk, ztrdu, ztrdv ) 93 96 ENDIF -
branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90
r3764 r3876 31 31 USE dom_oce ! ocean space and time domain 32 32 USE phycst ! physical constants 33 USE trd mod ! ocean dynamics trends34 USE trd mod_oce ! ocean variables trends33 USE trd_oce ! trends: ocean variables 34 USE trddyn ! trend manager: dynamics 35 35 USE in_out_manager ! I/O manager 36 36 USE prtctl ! Print control 37 USE lbclnk ! lateral boundary condition 37 USE lbclnk ! lateral boundary condition 38 38 USE lib_mpp ! MPP library 39 39 USE wrk_nemo ! Memory Allocation … … 46 46 PUBLIC dyn_hpg_init ! routine called by opa module 47 47 48 ! !!* Namelist namdyn_hpg : hydrostatic pressure gradient 48 ! !!* Namelist namdyn_hpg : hydrostatic pressure gradient 49 49 LOGICAL , PUBLIC :: ln_hpg_zco = .TRUE. !: z-coordinate - full steps 50 50 LOGICAL , PUBLIC :: ln_hpg_zps = .FALSE. !: z-coordinate - partial steps (interpolation) … … 54 54 LOGICAL , PUBLIC :: ln_dynhpg_imp = .FALSE. !: semi-implicite hpg flag 55 55 56 INTEGER , PUBLIC :: nhpg = 0 ! = 0 to 7, type of pressure gradient scheme used ! (deduced from ln_hpg_... flags) (PUBLIC for TAM)56 INTEGER :: nhpg = 0 ! = 0 to 7, type of pressure gradient scheme used ! (deduced from ln_hpg_... flags) 57 57 58 58 !! * Substitutions … … 70 70 !! *** ROUTINE dyn_hpg *** 71 71 !! 72 !! ** Method : Call the hydrostatic pressure gradient routine 72 !! ** Method : Call the hydrostatic pressure gradient routine 73 73 !! using the scheme defined in the namelist 74 !! 74 !! 75 75 !! ** Action : - Update (ua,va) with the now hydrastatic pressure trend 76 !! - Save the trend(l_trddyn=T)76 !! - send trends to trd_dyn for futher diagnostics (l_trddyn=T) 77 77 !!---------------------------------------------------------------------- 78 78 INTEGER, INTENT(in) :: kt ! ocean time-step index … … 99 99 ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 100 100 ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 101 CALL trd_ mod( ztrdu, ztrdv, jpdyn_trd_hpg, 'DYN', kt )101 CALL trd_dyn( ztrdu, ztrdv, jpdyn_hpg, kt ) 102 102 CALL wrk_dealloc( jpi,jpj,jpk, ztrdu, ztrdv ) 103 103 ENDIF … … 427 427 END SUBROUTINE hpg_sco 428 428 429 429 430 SUBROUTINE hpg_djc( kt ) 430 431 !!--------------------------------------------------------------------- … … 664 665 !! 665 666 !! ** Action : - Update (ua,va) with the now hydrastatic pressure trend 666 !! - Save the trend (l_trddyn=T)667 !!668 667 !!---------------------------------------------------------------------- 669 668 INTEGER, PARAMETER :: polynomial_type = 1 ! 1: cubic spline, 2: linear … … 717 716 718 717 ! Transfer the depth of "T(:,:,:)" to vertical coordinate "zdept(:,:,:)" 719 DO jj = 1, jpj; DO ji = 1, jpi 720 zdept(ji,jj,1) = 0.5_wp * fse3w(ji,jj,1) - sshn(ji,jj) * znad 721 END DO ; END DO 722 723 DO jk = 2, jpk; DO jj = 1, jpj; DO ji = 1, jpi 724 zdept(ji,jj,jk) = zdept(ji,jj,jk-1) + fse3w(ji,jj,jk) 725 END DO ; END DO ; END DO 726 727 fsp(:,:,:) = zrhh(:,:,:) 718 DO jj = 1, jpj 719 DO ji = 1, jpi 720 zdept(ji,jj,1) = 0.5_wp * fse3w(ji,jj,1) - sshn(ji,jj) * znad 721 END DO 722 END DO 723 724 DO jk = 2, jpk 725 DO jj = 1, jpj 726 DO ji = 1, jpi 727 zdept(ji,jj,jk) = zdept(ji,jj,jk-1) + fse3w(ji,jj,jk) 728 END DO 729 END DO 730 END DO 731 732 fsp(:,:,:) = zrhh (:,:,:) 728 733 xsp(:,:,:) = zdept(:,:,:) 729 734 … … 926 931 END SUBROUTINE hpg_prj 927 932 933 928 934 SUBROUTINE cspline(fsp, xsp, asp, bsp, csp, dsp, polynomial_type) 929 935 !!---------------------------------------------------------------------- … … 934 940 !! ** Method : f(x) = asp + bsp*x + csp*x^2 + dsp*x^3 935 941 !! Reference: CJC Kruger, Constrained Cubic Spline Interpoltation 936 !!937 942 !!---------------------------------------------------------------------- 938 943 IMPLICIT NONE … … 942 947 INTEGER, INTENT(in) :: polynomial_type ! 1: cubic spline 943 948 ! 2: Linear 944 945 ! Local Variables 949 ! 946 950 INTEGER :: ji, jj, jk ! dummy loop indices 947 951 INTEGER :: jpi, jpj, jpkm1 … … 1033 1037 ENDIF 1034 1038 1035 1036 1039 END SUBROUTINE cspline 1037 1040 … … 1043 1046 !! ** Purpose : 1-d linear interpolation 1044 1047 !! 1045 !! ** Method : 1046 !! interpolation is straight forward 1048 !! ** Method : interpolation is straight forward 1047 1049 !! extrapolation is also permitted (no value limit) 1048 !!1049 1050 !!---------------------------------------------------------------------- 1050 1051 IMPLICIT NONE … … 1063 1064 END FUNCTION interp1 1064 1065 1066 1065 1067 FUNCTION interp2(x, a, b, c, d) RESULT(f) 1066 1068 !!---------------------------------------------------------------------- … … 1126 1128 END FUNCTION integ_spline 1127 1129 1128 1129 1130 !!====================================================================== 1130 1131 END MODULE dynhpg -
branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/DYN/dynkeg.F90
r3294 r3876 14 14 USE oce ! ocean dynamics and tracers 15 15 USE dom_oce ! ocean space and time domain 16 USE trd mod ! ocean dynamics trends17 USE trd mod_oce ! ocean variables trends16 USE trd_oce ! trends: ocean variables 17 USE trddyn ! trend manager: dynamics 18 18 USE in_out_manager ! I/O manager 19 19 USE lib_mpp ! MPP library … … 52 52 !! 53 53 !! ** Action : - Update the (ua, va) with the hor. ke gradient trend 54 !! - s ave this trends(l_trddyn=T) for post-processing54 !! - send this trends to trd_dyn (l_trddyn=T) for post-processing 55 55 !!---------------------------------------------------------------------- 56 56 INTEGER, INTENT( in ) :: kt ! ocean time-step index … … 131 131 ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 132 132 ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 133 CALL trd_ mod( ztrdu, ztrdv, jpdyn_trd_keg, 'DYN', kt )133 CALL trd_dyn( ztrdu, ztrdv, jpdyn_keg, kt ) 134 134 CALL wrk_dealloc( jpi,jpj,jpk, ztrdu, ztrdv ) 135 135 ENDIF -
branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf.F90
r3294 r3876 20 20 USE dynldf_iso ! lateral mixing (dyn_ldf_iso routine) 21 21 USE dynldf_lap ! lateral mixing (dyn_ldf_lap routine) 22 USE trd mod ! ocean dynamics and tracer trends23 USE trd mod_oce ! ocean variables trends22 USE trd_oce ! trends: ocean variables 23 USE trddyn ! trend manager: dynamics 24 24 USE prtctl ! Print control 25 25 USE in_out_manager ! I/O manager … … 54 54 !! ** Purpose : compute the lateral ocean dynamics physics. 55 55 !!---------------------------------------------------------------------- 56 !57 56 INTEGER, INTENT(in) :: kt ! ocean time-step index 58 57 ! … … 106 105 ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 107 106 ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 108 CALL trd_ mod( ztrdu, ztrdv, jpdyn_trd_ldf, 'DYN', kt )107 CALL trd_dyn( ztrdu, ztrdv, jpdyn_ldf, kt ) 109 108 CALL wrk_dealloc( jpi, jpj, jpk, ztrdu, ztrdv ) 110 109 ENDIF -
branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf_bilap.F90
r3634 r3876 19 19 USE dom_oce ! ocean space and time domain 20 20 USE ldfdyn_oce ! ocean dynamics: lateral physics 21 ! 21 22 USE in_out_manager ! I/O manager 22 USE trdmod ! ocean dynamics trends23 USE trdmod_oce ! ocean variables trends24 23 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 25 24 USE wrk_nemo ! Memory Allocation … … 70 69 !! Add this before trend to the general trend (ua,va): 71 70 !! (ua,va) = (ua,va) + (diffu,diffv) 72 !! 'key_trddyn' defined: the two components of the horizontal73 !! diffusion trend are saved.74 71 !! 75 72 !! ** Action : - Update (ua,va) with the before iso-level biharmonic -
branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf_bilapg.F90
r3634 r3876 20 20 USE ldfdyn_oce ! ocean dynamics lateral physics 21 21 USE zdf_oce ! ocean vertical physics 22 USE trdmod ! ocean dynamics trends23 USE trdmod_oce ! ocean variables trends24 22 USE ldfslp ! iso-neutral slopes available 23 ! 25 24 USE in_out_manager ! I/O manager 26 25 USE lib_mpp ! MPP library … … 81 80 !! -3- Add this trend to the general trend (ta,sa): 82 81 !! (ua,va) = (ua,va) + (zwk3,zwk4) 83 !! 'key_trddyn' defined: the trend is saved for diagnostics.84 82 !! 85 83 !! ** Action : - Update (ua,va) arrays with the before geopotential 86 84 !! biharmonic mixing trend. 87 !! - save the trend in (zwk3,zwk4) ('key_trddyn')88 85 !!---------------------------------------------------------------------- 89 86 INTEGER, INTENT( in ) :: kt ! ocean time-step index … … 176 173 !! pu and pv (all the components except 177 174 !! second order vertical derivative term) 178 !! 'key_trddyn' defined: the trend is saved for diagnostics.179 175 !!---------------------------------------------------------------------- 180 176 !! -
branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf_iso.F90
r3294 r3876 22 22 USE ldftra_oce ! ocean tracer lateral physics 23 23 USE zdf_oce ! ocean vertical physics 24 USE trdmod ! ocean dynamics trends25 USE trdmod_oce ! ocean variables trends26 24 USE ldfslp ! iso-neutral slopes 27 25 USE lbclnk ! ocean lateral boundary conditions (or mpp link) -
branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf_lap.F90
r3294 r3876 19 19 USE ldfdyn_oce ! ocean dynamics: lateral physics 20 20 USE zdf_oce ! ocean vertical physics 21 !!gm ?? USE ldfslp ! iso-neutral slopes 22 ! 21 23 USE in_out_manager ! I/O manager 22 USE trdmod ! ocean dynamics trends23 USE trdmod_oce ! ocean variables trends24 USE ldfslp ! iso-neutral slopes25 24 USE timing ! Timing 26 25 … … 57 56 !! Add this before trend to the general trend (ua,va): 58 57 !! (ua,va) = (ua,va) + (diffu,diffv) 59 !! 'key_trddyn' activated: the two components of the horizontal60 !! diffusion trend are saved.61 58 !! 62 59 !! ** Action : - Update (ua,va) with the before iso-level harmonic -
branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/DYN/dynnept.F90
r3723 r3876 5 5 !! recoded version of simplest case (u*, v* only) 6 6 !!====================================================================== 7 !! History : 1.0 ! 2007-06 (Michael Dunphy) Modular form: - new namelist parameters 8 !! - horizontal diffusion for Neptune 9 !! - vertical diffusion for gm in momentum eqns 10 !! - option to use Neptune in Coriolis eqn 7 !! History : 1.0 ! 2007-06 (Zeliang Wang, Michael Dunphy, BIO) Original code 8 !! Modular form: - new namelist parameters 9 !! - horizontal diffusion for Neptune 10 !! - vertical diffusion for gm in momentum eqns 11 !! - option to use Neptune in Coriolis eqn 11 12 !! 2011-08 (Jeff Blundell, NOCS) Simplified form for temporally invariant u*, v* 12 !! Horizontal and vertical diffusivity formulations removed 13 !! Dynamic allocation of storage added 14 !! Option of ramping Neptune vel. down 15 !! to zero added in shallow depths added 13 !! Horizontal and vertical diffusivity formulations removed 14 !! Dynamic allocation of storage added 15 !! Option of ramping Neptune vel. down to zero added in shallow depths added 16 !!---------------------------------------------------------------------- 17 16 18 !!---------------------------------------------------------------------- 17 19 !! dynnept_alloc : … … 30 32 USE phycst 31 33 USE lbclnk 32 USE wrk_nemo ! Memory Allocation34 USE wrk_nemo ! Memory Allocation 33 35 34 36 IMPLICIT NONE 35 37 PRIVATE 36 38 37 !! * Routine accessibility 38 PUBLIC dyn_nept_init ! routine called by nemogcm.F90 39 PUBLIC dyn_nept_cor ! routine called by step.F90 40 !! dynnept_alloc() is called only by dyn_nept_init, within this module 41 !! dyn_nept_div_cur_init is called only by dyn_nept_init, within this module 42 !! dyn_nept_vel is called only by dyn_nept_cor, within this module 43 44 !! * Shared module variables 45 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: zunep, zvnep ! Neptune u and v 46 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: zhdivnep ! hor. div for Neptune vel. 47 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: zmrotnep ! curl for Neptune vel. 48 49 50 !! * Namelist namdyn_nept variables 39 PUBLIC dyn_nept_init ! routine called by nemogcm.F90 40 PUBLIC dyn_nept_cor ! routine called by step.F90 41 42 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: zunep, zvnep ! Neptune u and v 43 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: zhdivnep ! hor. div for Neptune vel. 44 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: zmrotnep ! curl for Neptune vel. 45 46 ! !!* Namelist namdyn_nept variables 51 47 LOGICAL, PUBLIC :: ln_neptsimp = .FALSE. ! yes/no simplified neptune 52 48 … … 60 56 REAL(wp) :: rn_htrmax = 200.0 ! max. depth of transition range 61 57 62 !! * Module variables63 64 65 58 !! * Substitutions 66 59 # include "vectopt_loop_substitute.h90" 67 60 # include "domzgr_substitute.h90" 68 61 !!---------------------------------------------------------------------- 69 !! OPA 9.0 , implemented by Bedford Institute of Oceanography 70 !!---------------------------------------------------------------------- 71 62 !! NEMO/OPA 3.3 , NEMO Consortium (2011) 63 !! $Id: dynadv_cen2.F90 3316 2012-02-21 16:00:02Z gm $ 64 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 65 !!---------------------------------------------------------------------- 72 66 CONTAINS 73 67 … … 90 84 !! and compute the arrays zunep and zvnep 91 85 !! 92 !! ** Method : zunep = 93 !! zvnep = 94 !! 95 !! ** History : 1.0 ! 07-05 (Zeliang Wang) Original code for zunep, zvnep 96 !! 1.1 ! 07-06 (Michael Dunphy) namelist and initialisation 97 !! 2.0 ! 2011-07 (Jeff Blundell, NOCS) 98 !! ! Simplified form for temporally invariant u*, v* 99 !! ! Horizontal and vertical diffusivity formulations removed 100 !! ! Includes optional tapering-off in shallow depths 86 !! ** Method : Simplified form for temporally invariant u*, v* 87 !! Horizontal and vertical diffusivity formulations removed 88 !! Includes optional tapering-off in shallow depths 101 89 !!---------------------------------------------------------------------- 102 90 USE iom 103 ! !91 ! 104 92 INTEGER :: ji, jj, jk ! dummy loop indices 105 93 REAL(wp) :: unemin,unemax,vnemin,vnemax ! extrema of (u*, v*) fields … … 140 128 ENDIF 141 129 ! 142 IF( .NOT. ln_neptsimp ) RETURN130 IF( .NOT. ln_neptsimp ) RETURN 143 131 ! ! Dynamically allocate local work arrays 144 132 CALL wrk_alloc( jpi, jpj , ht, htn, tscale, tsp, hur_n, hvr_n, hu_n, hv_n ) … … 352 340 !! ** Action : - compute zhdivnep, the hor. divergence of (u*, v*) 353 341 !! - compute zmrotnep, the rel. vorticity of (u*, v*) 354 !! 355 !! History : OPA ! 1987-06 (P. Andrich, D. L Hostis) Original code 356 !! 4.0 ! 1991-11 (G. Madec) 357 !! 6.0 ! 1993-03 (M. Guyon) symetrical conditions 358 !! 7.0 ! 1996-01 (G. Madec) s-coordinates 359 !! 8.0 ! 1997-06 (G. Madec) lateral boundary cond., lbc 360 !! 8.1 ! 1997-08 (J.M. Molines) Open boundaries 361 !! 8.2 ! 2000-03 (G. Madec) no slip accurate 362 !! NEMO 1.0 ! 2002-09 (G. Madec, E. Durand) Free form, F90 363 !! - ! 2005-01 (J. Chanut) Unstructured open boundaries 364 !! - ! 2003-08 (G. Madec) merged of cur and div, free form, F90 365 !! - ! 2005-01 (J. Chanut, A. Sellar) unstructured open boundaries 366 !! 3.3 ! 2010-09 (D.Storkey and E.O'Dea) bug fixes for BDY module 367 !! ! 2011-06 (Jeff Blundell, NOCS) Adapt code from divcur.F90 368 !! ! to compute Neptune effect fields only 369 !!---------------------------------------------------------------------- 370 ! 342 !!---------------------------------------------------------------------- 371 343 INTEGER :: ji, jj, jk ! dummy loop indices 372 344 !!---------------------------------------------------------------------- … … 508 480 509 481 SUBROUTINE dyn_nept_smooth_vel( htold, htnew, ld_option ) 510 511 482 !!---------------------------------------------------------------------- 512 483 !! *** ROUTINE dyn_nept_smooth_vel *** … … 598 569 END SUBROUTINE dyn_nept_smooth_vel 599 570 571 !!============================================================================== 600 572 END MODULE dynnept -
branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90
r3764 r3876 17 17 !! 3.3 ! 2010-09 (D. Storkey, E.O'Dea) Bug fix for BDY module 18 18 !! 3.3 ! 2011-03 (P. Oddo) Bug fix for time-splitting+(BDY-OBC) and not VVL 19 !! 3.4 ! 2012-02 (G. Madec) add the diagnostic of the time filter trends 19 20 !!------------------------------------------------------------------------- 20 21 … … 29 30 USE dynadv ! dynamics: vector invariant versus flux form 30 31 USE domvvl ! variable volume 32 USE trd_oce ! trends: ocean variables 33 USE trddyn ! trend manager: dynamics 34 USE trdken ! trend manager: kinetic energy 31 35 USE obc_oce ! ocean open boundary conditions 32 36 USE obcdyn ! open boundary condition for momentum (obc_dyn routine) … … 37 41 USE bdydyn ! ocean open boundary conditions 38 42 USE bdyvol ! ocean open boundary condition (bdy_vol routines) 43 ! 39 44 USE in_out_manager ! I/O manager 45 USE iom ! I/O manager library 40 46 USE lbclnk ! lateral boundary condition (or mpp link) 41 47 USE lib_mpp ! MPP library 42 48 USE wrk_nemo ! Memory Allocation 43 49 USE prtctl ! Print control 50 USE timing ! Timing 44 51 #if defined key_agrif 45 52 USE agrif_opa_interp 46 53 #endif 47 USE timing ! Timing48 54 49 55 IMPLICIT NONE … … 81 87 !! at the local domain boundaries through lbc_lnk call, 82 88 !! at the one-way open boundaries (lk_obc=T), 83 !! at the AGRIF zoom 89 !! at the AGRIF zoom boundaries (lk_agrif=T) 84 90 !! 85 91 !! * Apply the time filter applied and swap of the dynamics … … 101 107 REAL(wp) :: z2dt ! temporary scalar 102 108 #endif 103 REAL(wp) :: zue3a, zue3n, zue3b, zuf, zec ! local scalars104 REAL(wp) :: zve3a, zve3n, zve3b, zvf 105 REAL(wp), POINTER, DIMENSION(:,:,:) :: ze3u_f, ze3v_f 109 REAL(wp) :: zue3a, zue3n, zue3b, zuf, zec ! local scalars 110 REAL(wp) :: zve3a, zve3n, zve3b, zvf, z1_2dt ! - - 111 REAL(wp), POINTER, DIMENSION(:,:,:) :: ze3u_f, ze3v_f, zua, zva 106 112 !!---------------------------------------------------------------------- 107 113 ! 108 IF( nn_timing == 1 ) CALL timing_start('dyn_nxt')109 ! 110 CALL wrk_alloc( jpi,jpj,jpk, ze3u_f, ze3v_f )114 IF( nn_timing == 1 ) CALL timing_start('dyn_nxt') 115 ! 116 CALL wrk_alloc( jpi,jpj,jpk, ze3u_f, ze3v_f, zua, zva ) 111 117 ! 112 118 IF( kt == nit000 ) THEN … … 156 162 # if defined key_obc 157 163 ! !* OBC open boundaries 158 IF( lk_obc ) CALL obc_dyn( kt )164 IF( lk_obc ) CALL obc_dyn( kt ) 159 165 ! 160 166 IF( .NOT. lk_dynspg_flt ) THEN … … 163 169 ! sshn_b (= after ssha_b) for time-splitting case (lk_dynspg_ts=T) 164 170 ! - Correct the barotropic velocities 165 IF( lk_obc ) CALL obc_dyn_bt( kt )171 IF( lk_obc ) CALL obc_dyn_bt( kt ) 166 172 ! 167 173 !!gm ERROR - potential BUG: sshn should not be modified at this stage !! ssh_nxt not alrady called … … 186 192 # endif 187 193 #endif 194 195 IF( l_trddyn ) THEN ! prepare the atf trend computation + some diagnostics 196 z1_2dt = 1._wp / (2. * rdt) ! Euler or leap-frog time step 197 IF( neuler == 0 .AND. kt == nit000 ) z1_2dt = 1._wp / rdt 198 ! 199 ! ! Kinetic energy and Conversion 200 IF( ln_KE_trd ) CALL trd_dyn( ua, va, jpdyn_ken, kt ) 201 ! 202 IF( ln_dyn_trd ) THEN ! 3D output: total momentum trends 203 zua(:,:,:) = ( ua(:,:,:) - ub(:,:,:) ) * z1_2dt 204 zva(:,:,:) = ( va(:,:,:) - vb(:,:,:) ) * z1_2dt 205 CALL iom_put( "utrd_tot", zua ) ! total momentum trends, except the asselin time filter 206 CALL iom_put( "vtrd_tot", zva ) 207 ENDIF 208 ! 209 zua(:,:,:) = un(:,:,:) ! save the now velocity before the asselin filter 210 zva(:,:,:) = vn(:,:,:) ! (caution: there will be a shift by 1 timestep in the 211 ! ! computation of the asselin filter trends) 212 ENDIF 188 213 189 214 ! Time filter and swap of dynamics arrays … … 274 299 ENDIF 275 300 301 IF( l_trddyn ) THEN ! 3D output: asselin filter trends on momentum 302 zua(:,:,:) = ( ub(:,:,:) - zua(:,:,:) ) * z1_2dt 303 zva(:,:,:) = ( vb(:,:,:) - zva(:,:,:) ) * z1_2dt 304 CALL trd_dyn( zua, zva, jpdyn_atf, kt ) 305 ENDIF 306 276 307 IF(ln_ctl) CALL prt_ctl( tab3d_1=un, clinfo1=' nxt - Un: ', mask1=umask, & 277 308 & tab3d_2=vn, clinfo2=' Vn: ' , mask2=vmask ) 278 309 ! 279 CALL wrk_dealloc( jpi,jpj,jpk, ze3u_f, ze3v_f )310 CALL wrk_dealloc( jpi,jpj,jpk, ze3u_f, ze3v_f, zua, zva ) 280 311 ! 281 312 IF( nn_timing == 1 ) CALL timing_stop('dyn_nxt') -
branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg.F90
r3625 r3876 22 22 USE dynspg_flt ! surface pressure gradient (dyn_spg_flt routine) 23 23 USE dynadv ! dynamics: vector invariant versus flux form 24 USE trd mod ! ocean dynamics trends25 USE trd mod_oce ! ocean variables trends24 USE trd_oce ! trends: ocean variables 25 USE trddyn ! trend manager: dynamics 26 26 USE prtctl ! Print control (prt_ctl routine) 27 27 USE in_out_manager ! I/O manager 28 28 USE lib_mpp ! MPP library 29 USE solver 30 USE wrk_nemo 31 USE timing 29 USE solver ! solver initialization 30 USE wrk_nemo ! Memory Allocation 31 USE timing ! Timing 32 32 33 33 … … 118 118 END DO 119 119 END DO 120 121 !!gm add here a call to dyn_trd for apr, the surf pressure trends ???? 122 120 123 ENDIF 121 124 … … 142 145 ! 143 146 CALL wrk_dealloc( jpi, jpj, zpice ) 147 148 !!gm add here a call to dyn_trd for ice pressure gradient, the surf pressure trends ???? 149 144 150 ENDIF 145 151 … … 170 176 CASE( 2 ) 171 177 z2dt = 2. * rdt 172 IF( neuler == 0 .AND. kt == nit000 ) z2dt = rdt178 IF( neuler == 0 .AND. kt == nit000 ) z2dt = rdt 173 179 ztrdu(:,:,:) = ( ua(:,:,:) - ub(:,:,:) ) / z2dt - ztrdu(:,:,:) 174 180 ztrdv(:,:,:) = ( va(:,:,:) - vb(:,:,:) ) / z2dt - ztrdv(:,:,:) 175 181 END SELECT 176 CALL trd_ mod( ztrdu, ztrdv, jpdyn_trd_spg, 'DYN', kt )182 CALL trd_dyn( ztrdu, ztrdv, jpdyn_spg, kt ) 177 183 ! 178 184 CALL wrk_dealloc( jpi, jpj, jpk, ztrdu, ztrdv ) -
branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_exp.F90
r3680 r3876 22 22 USE obc_par ! open boundary condition parameters 23 23 USE obcdta ! open boundary condition data (bdy_dta_bt routine) 24 ! 24 25 USE in_out_manager ! I/O manager 25 26 USE lib_mpp ! distributed memory computing library -
branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_flt.F90
r3765 r3876 13 13 !! - ! 2006-08 (J.Chanut, A.Sellar) Calls to BDY routines. 14 14 !! 3.2 ! 2009-03 (G. Madec, M. Leclair, R. Benshila) introduce sshwzv module 15 !! 3.5 ! 2012-05 (F. Roquet, G. Madec) add some trends diag 15 16 !!---------------------------------------------------------------------- 16 17 #if defined key_dynspg_flt || defined key_esopa … … 39 40 USE bdyvol ! ocean open boundary condition (bdy_vol routine) 40 41 USE cla ! cross land advection 42 USE trd_oce ! trends: ocean variables 43 USE trddyn ! trend manager: dynamics 44 ! 41 45 USE in_out_manager ! I/O manager 42 46 USE lib_mpp ! distributed memory computing library … … 46 50 USE iom 47 51 USE lib_fortran 52 USE timing ! Timing 48 53 #if defined key_agrif 49 54 USE agrif_opa_interp 50 55 #endif 51 USE timing ! Timing52 56 53 57 IMPLICIT NONE … … 102 106 !! ** Action : - Update (ua,va) with the surf. pressure gradient trend 103 107 !! 104 !! References : Roullet and Madec 1999, JGR.108 !! References : Roullet and Madec, JGR, 2000. 105 109 !!--------------------------------------------------------------------- 106 110 INTEGER, INTENT(in ) :: kt ! ocean time-step index 107 111 INTEGER, INTENT( out) :: kindic ! solver convergence flag (<0 if not converge) 108 ! !112 ! 109 113 INTEGER :: ji, jj, jk ! dummy loop indices 110 114 REAL(wp) :: z2dt, z2dtg, zgcb, zbtd, ztdgu, ztdgv ! local scalars 111 REAL(wp), POINTER, DIMENSION(:,:,:) :: zub, zvb 115 REAL(wp), POINTER, DIMENSION(:,:,:) :: ztrdu, ztrdv 116 REAL(wp), POINTER, DIMENSION(:,:) :: zpw 112 117 !!---------------------------------------------------------------------- 113 118 ! 114 119 IF( nn_timing == 1 ) CALL timing_start('dyn_spg_flt') 115 !116 CALL wrk_alloc( jpi,jpj,jpk, zub, zvb )117 120 ! 118 121 IF( kt == nit000 ) THEN … … 183 186 END DO 184 187 END DO 188 ! save the explicit SPG trends for further diagnostics 189 ! --------------------------------------------- 190 IF( l_trddyn ) THEN ! temporary save of ta and sa trends 191 CALL wrk_alloc( jpi, jpj, jpk, ztrdu, ztrdv ) 192 DO jk = 1, jpkm1 ! unweighted time stepping 193 DO jj = 2, jpjm1 194 DO ji = fs_2, fs_jpim1 ! vector opt. 195 ztrdu(ji,jj,jk) = spgu(ji,jj) * umask(ji,jj,jk) 196 ztrdv(ji,jj,jk) = spgv(ji,jj) * vmask(ji,jj,jk) 197 END DO 198 END DO 199 END DO 200 CALL trd_dyn( ztrdu, ztrdv, jpdyn_spgexp, kt ) 201 ENDIF 185 202 ! 186 203 ENDIF … … 208 225 END DO 209 226 210 ! vertical sum227 ! ! vertical sum 211 228 !CDIR NOLOOPCHG 212 229 IF( lk_vopt_loop ) THEN ! vector opt., forced unroll … … 217 234 END DO 218 235 END DO 219 ELSE ! No vector opt.236 ELSE ! No vector opt. 220 237 DO jk = 1, jpkm1 221 238 DO jj = 2, jpjm1 … … 228 245 ENDIF 229 246 230 ! transport: multiplied by the horizontal scale factor 231 DO jj = 2, jpjm1 247 DO jj = 2, jpjm1 ! transport: multiplied by the horizontal scale factor 232 248 DO ji = fs_2, fs_jpim1 ! vector opt. 233 249 spgu(ji,jj) = spgu(ji,jj) * e2u(ji,jj) … … 342 358 ENDIF 343 359 #endif 360 361 IF( l_trddyn ) THEN 362 ztrdu(:,:,:) = ua(:,:,:) ! save the after velocity before the filtered SPG 363 ztrdv(:,:,:) = va(:,:,:) 364 ! 365 CALL wrk_alloc( jpi, jpj, zpw ) 366 ! 367 zpw(:,:) = - z2dt * gcx(:,:) 368 CALL iom_put( "ssh_flt" , zpw ) ! output equivalent ssh modification due to implicit filter 369 ! 370 ! ! save surface pressure flux: -pw at z=0 371 zpw(:,:) = - rau0 * grav * sshn(:,:) * wn(:,:,1) * tmask(:,:,1) 372 CALL iom_put( "pw0_exp" , zpw ) 373 zpw(:,:) = wn(:,:,1) 374 CALL iom_put( "w0" , zpw ) 375 zpw(:,:) = rau0 * z2dtg * gcx(:,:) * wn(:,:,1) * tmask(:,:,1) 376 CALL iom_put( "pw0_flt" , zpw ) 377 ! 378 CALL wrk_dealloc( jpi, jpj, zpw ) 379 ! 380 ENDIF 381 344 382 ! Add the trends multiplied by z2dt to the after velocity 345 383 ! ------------------------------------------------------- … … 356 394 END DO 357 395 358 ! write filtered free surface arrays in restart file 359 ! -------------------------------------------------- 360 IF( lrst_oce ) CALL flt_rst( kt, 'WRITE' ) 361 ! 362 CALL wrk_dealloc( jpi,jpj,jpk, zub, zvb ) 396 IF( l_trddyn ) THEN ! save the explicit SPG trends for further diagnostics 397 ztrdu(:,:,:) = ( ua(:,:,:) - ztrdu(:,:,:) ) / z2dt 398 ztrdv(:,:,:) = ( va(:,:,:) - ztrdv(:,:,:) ) / z2dt 399 CALL trd_dyn( ztrdu, ztrdv, jpdyn_spgflt, kt ) 400 ! 401 CALL wrk_dealloc( jpi, jpj, jpk, ztrdu, ztrdv ) 402 ENDIF 403 404 IF( lrst_oce ) CALL flt_rst( kt, 'WRITE' ) ! write filtered free surface arrays in restart file 363 405 ! 364 406 IF( nn_timing == 1 ) CALL timing_stop('dyn_spg_flt') … … 373 415 !! ** Purpose : Read or write filtered free surface arrays in restart file 374 416 !!---------------------------------------------------------------------- 375 INTEGER , INTENT(in) :: kt 376 CHARACTER(len=*), INTENT(in) :: cdrw 417 INTEGER , INTENT(in) :: kt ! ocean time-step 418 CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag 377 419 !!---------------------------------------------------------------------- 378 420 ! -
branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90
r3680 r3876 400 400 ! !* Update the forcing (BDY and tides) 401 401 ! ! ------------------ 402 IF( lk_obc ) CALL obc_dta_bt ( kt, jn )403 IF( lk_bdy ) CALL bdy_dta ( kt, jit=jn, time_offset=+1 )404 IF ( ln_tide_pot .AND. lk_tide)CALL upd_tide( kt, jn )402 IF( lk_obc ) CALL obc_dta_bt ( kt, jn ) 403 IF( lk_bdy ) CALL bdy_dta ( kt, jit=jn, time_offset=+1 ) 404 IF( ln_tide_pot .AND. lk_tide ) CALL upd_tide( kt, jn ) 405 405 406 406 ! !* after ssh_e -
branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/DYN/dynvor.F90
r3802 r3876 15 15 !! 3.2 ! 2009-04 (R. Benshila) vvl: correction of een scheme 16 16 !! 3.3 ! 2010-10 (C. Ethe, G. Madec) reorganisation of initialisation phase 17 !! 3.5 ! 2012-02 (G. Madec) trend simplification: suppress jpdyn_trd_dat vorticity 17 18 !!---------------------------------------------------------------------- 18 19 … … 29 30 USE dommsk ! ocean mask 30 31 USE dynadv ! momentum advection (use ln_dynadv_vec value) 31 USE trd mod ! ocean dynamics trends32 USE trd mod_oce ! ocean variables trends32 USE trd_oce ! trends: ocean variables 33 USE trddyn ! trend manager: dynamics 33 34 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 34 35 USE prtctl ! Print control … … 73 74 !! ** Action : - Update (ua,va) with the now vorticity term trend 74 75 !! - save the trends in (ztrdu,ztrdv) in 2 parts (relative 75 !! and planetary vorticity trends) ( 'key_trddyn')76 !! and planetary vorticity trends) (l_trddyn=T) 76 77 !!---------------------------------------------------------------------- 77 78 INTEGER, INTENT( in ) :: kt ! ocean time-step index … … 108 109 ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 109 110 ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 110 CALL trd_ mod( ztrdu, ztrdv, jpdyn_trd_rvo, 'DYN', kt )111 CALL trd_dyn( ztrdu, ztrdv, jpdyn_rvo, kt ) 111 112 ztrdu(:,:,:) = ua(:,:,:) 112 113 ztrdv(:,:,:) = va(:,:,:) … … 114 115 ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 115 116 ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 116 CALL trd_mod( ztrdu, ztrdv, jpdyn_trd_pvo, 'DYN', kt ) 117 CALL trd_mod( ztrdu, ztrdv, jpdyn_trd_dat, 'DYN', kt ) 117 CALL trd_dyn( ztrdu, ztrdv, jpdyn_pvo, kt ) 118 118 ELSE 119 119 CALL vor_ene( kt, ntot, ua, va ) ! total vorticity … … 127 127 ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 128 128 ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 129 CALL trd_ mod( ztrdu, ztrdv, jpdyn_trd_rvo, 'DYN', kt )129 CALL trd_dyn( ztrdu, ztrdv, jpdyn_rvo, kt ) 130 130 ztrdu(:,:,:) = ua(:,:,:) 131 131 ztrdv(:,:,:) = va(:,:,:) … … 133 133 ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 134 134 ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 135 CALL trd_mod( ztrdu, ztrdv, jpdyn_trd_pvo, 'DYN', kt ) 136 CALL trd_mod( ztrdu, ztrdv, jpdyn_trd_dat, 'DYN', kt ) 135 CALL trd_dyn( ztrdu, ztrdv, jpdyn_pvo, kt ) 137 136 ELSE 138 137 CALL vor_ens( kt, ntot, ua, va ) ! total vorticity … … 146 145 ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 147 146 ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 148 CALL trd_ mod( ztrdu, ztrdv, jpdyn_trd_rvo, 'DYN', kt )147 CALL trd_dyn( ztrdu, ztrdv, jpdyn_rvo, kt ) 149 148 ztrdu(:,:,:) = ua(:,:,:) 150 149 ztrdv(:,:,:) = va(:,:,:) … … 152 151 ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 153 152 ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 154 CALL trd_mod( ztrdu, ztrdv, jpdyn_trd_pvo, 'DYN', kt ) 155 CALL trd_mod( ztrdu, ztrdv, jpdyn_trd_dat, 'DYN', kt ) 153 CALL trd_dyn( ztrdu, ztrdv, jpdyn_pvo, kt ) 156 154 ELSE 157 155 CALL vor_mix( kt ) ! total vorticity (mix=ens-ene) … … 165 163 ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 166 164 ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 167 CALL trd_ mod( ztrdu, ztrdv, jpdyn_trd_rvo, 'DYN', kt )165 CALL trd_dyn( ztrdu, ztrdv, jpdyn_rvo, kt ) 168 166 ztrdu(:,:,:) = ua(:,:,:) 169 167 ztrdv(:,:,:) = va(:,:,:) … … 171 169 ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 172 170 ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 173 CALL trd_mod( ztrdu, ztrdv, jpdyn_trd_pvo, 'DYN', kt ) 174 CALL trd_mod( ztrdu, ztrdv, jpdyn_trd_dat, 'DYN', kt ) 171 CALL trd_dyn( ztrdu, ztrdv, jpdyn_pvo, kt ) 175 172 ELSE 176 173 CALL vor_een( kt, ntot, ua, va ) ! total vorticity … … 211 208 !! 212 209 !! ** Action : - Update (ua,va) with the now vorticity term trend 213 !! - save the trends in (ztrdu,ztrdv) in 2 parts (relative214 !! and planetary vorticity trends) ('key_trddyn')215 210 !! 216 211 !! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689. … … 328 323 !! 329 324 !! ** Action : - Update (ua,va) arrays with the now vorticity term trend 330 !! - Save the trends in (ztrdu,ztrdv) in 2 parts (relative331 !! and planetary vorticity trends) ('key_trddyn')332 325 !! 333 326 !! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689. … … 444 437 !! 445 438 !! ** Action : - Update (ua,va) arrays with the now vorticity term trend 446 !! - Save the trends in (ztrdu,ztrdv) in 2 parts (relative447 !! and planetary vorticity trends) ('key_trddyn')448 439 !! 449 440 !! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689. … … 557 548 !! 558 549 !! ** Action : - Update (ua,va) with the now vorticity term trend 559 !! - save the trends in (ztrdu,ztrdv) in 2 parts (relative560 !! and planetary vorticity trends) ('key_trddyn')561 550 !! 562 551 !! References : Arakawa and Lamb 1980, Mon. Wea. Rev., 109, 18-36 -
branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/DYN/dynzad.F90
r3294 r3876 16 16 USE dom_oce ! ocean space and time domain 17 17 USE sbc_oce ! surface boundary condition: ocean 18 USE trdmod_oce ! ocean variables trends 19 USE trdmod ! ocean dynamics trends 18 USE trd_oce ! trends: ocean variables 19 USE trddyn ! trend manager: dynamics 20 ! 20 21 USE in_out_manager ! I/O manager 21 USE lib_mpp 22 USE lib_mpp ! MPP library 22 23 USE prtctl ! Print control 23 USE wrk_nemo 24 USE timing 24 USE wrk_nemo ! Memory Allocation 25 USE timing ! Timing 25 26 26 27 IMPLICIT NONE … … 53 54 !! 54 55 !! ** Action : - Update (ua,va) with the vert. momentum adv. trends 55 !! - S ave the trends in (ztrdu,ztrdv) ('key_trddyn')56 !! - Send the trends to trddyn for diagnostics (l_trddyn=T) 56 57 !!---------------------------------------------------------------------- 57 58 INTEGER, INTENT(in) :: kt ! ocean time-step inedx … … 118 119 ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 119 120 ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 120 CALL trd_ mod(ztrdu, ztrdv, jpdyn_trd_zad, 'DYN', kt)121 CALL trd_dyn( ztrdu, ztrdv, jpdyn_zad, kt ) 121 122 CALL wrk_dealloc( jpi, jpj, jpk, ztrdu, ztrdv ) 122 123 ENDIF -
branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf.F90
r3294 r3876 20 20 21 21 USE ldfdyn_oce ! ocean dynamics: lateral physics 22 USE trd mod ! ocean active dynamics and tracers trends23 USE trd mod_oce ! ocean variables trends22 USE trd_oce ! trends: ocean variables 23 USE trddyn ! trend manager: dynamics 24 24 USE in_out_manager ! I/O manager 25 25 USE lib_mpp ! MPP library … … 91 91 ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 92 92 ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 93 CALL trd_ mod( ztrdu, ztrdv, jpdyn_trd_zdf, 'DYN', kt )93 CALL trd_dyn( ztrdu, ztrdv, jpdyn_zdf, kt ) 94 94 ! 95 95 CALL wrk_dealloc( jpi, jpj, jpk, ztrdu, ztrdv ) -
branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf_imp.F90
r3625 r3876 94 94 IF( ln_bfrimp ) THEN 95 95 # if defined key_vectopt_loop 96 DO jj = 1, 197 DO ji = jpi+2, jpij-jpi-1 ! vector opt. (forced unrolling)96 DO jj = 1, 1 97 DO ji = jpi+2, jpij-jpi-1 ! vector opt. (forced unrolling) 98 98 # else 99 DO jj = 2, jpjm1100 DO ji = 2, jpim199 DO jj = 2, jpjm1 100 DO ji = 2, jpim1 101 101 # endif 102 ikbu = mbku(ji,jj) ! ocean bottom level at u- and v-points103 ikbv = mbkv(ji,jj) ! (deepest ocean u- and v-points)104 zavmu(ji,jj) = avmu(ji,jj,ikbu+1)105 zavmv(ji,jj) = avmv(ji,jj,ikbv+1)106 avmu(ji,jj,ikbu+1) = -bfrua(ji,jj) * fse3uw(ji,jj,ikbu+1)107 avmv(ji,jj,ikbv+1) = -bfrva(ji,jj) * fse3vw(ji,jj,ikbv+1)108 END DO109 END DO102 ikbu = mbku(ji,jj) ! ocean bottom level at u- and v-points 103 ikbv = mbkv(ji,jj) ! (deepest ocean u- and v-points) 104 zavmu(ji,jj) = avmu(ji,jj,ikbu+1) 105 zavmv(ji,jj) = avmv(ji,jj,ikbv+1) 106 avmu(ji,jj,ikbu+1) = -bfrua(ji,jj) * fse3uw(ji,jj,ikbu+1) 107 avmv(ji,jj,ikbv+1) = -bfrva(ji,jj) * fse3vw(ji,jj,ikbv+1) 108 END DO 109 END DO 110 110 ENDIF 111 111 … … 284 284 IF( ln_bfrimp ) THEN 285 285 # if defined key_vectopt_loop 286 DO jj = 1, 1287 DO ji = jpi+2, jpij-jpi-1 ! vector opt. (forced unrolling)286 DO jj = 1, 1 287 DO ji = jpi+2, jpij-jpi-1 ! vector opt. (forced unrolling) 288 288 # else 289 DO jj = 2, jpjm1290 DO ji = 2, jpim1289 DO jj = 2, jpjm1 290 DO ji = 2, jpim1 291 291 # endif 292 ikbu = mbku(ji,jj) ! ocean bottom level at u- and v-points293 ikbv = mbkv(ji,jj) ! (deepest ocean u- and v-points)294 avmu(ji,jj,ikbu+1) = zavmu(ji,jj)295 avmv(ji,jj,ikbv+1) = zavmv(ji,jj)296 END DO297 END DO292 ikbu = mbku(ji,jj) ! ocean bottom level at u- and v-points 293 ikbv = mbkv(ji,jj) ! (deepest ocean u- and v-points) 294 avmu(ji,jj,ikbu+1) = zavmu(ji,jj) 295 avmv(ji,jj,ikbv+1) = zavmv(ji,jj) 296 END DO 297 END DO 298 298 ENDIF 299 299 ! -
branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/TRA/eosbn2.F90
r3625 r3876 18 18 !! 3.3 ! 2010-05 (C. Ethe, G. Madec) merge TRC-TRA 19 19 !! - ! 2010-10 (G. Nurser, G. Madec) add eos_alpbet used in ldfslp 20 !! 3.5 ! 2012-03 (F. Roquet, G. Madec) add pen_ddt_dds used in trdpen 21 !! - ! 2012-05 (F. Roquet) add Vallis and original JM95 equation of state 22 !! - ! 2012-05 (F. Roquet) add eos_alpha_beta 20 23 !!---------------------------------------------------------------------- 21 24 … … 28 31 !! eos_bn2 : Compute the Brunt-Vaisala frequency 29 32 !! eos_alpbet : calculates the in situ thermal/haline expansion ratio 33 !! eos_expansion_ratio : calculates the partial derivative of in situ density with respect to T and S 30 34 !! tfreez : Compute the surface freezing temperature 31 35 !! eos_init : set eos parameters (namelist) … … 51 55 END INTERFACE 52 56 53 PUBLIC eos ! called by step, istate, tranpc and zpsgrd modules 54 PUBLIC eos_init ! called by istate module 55 PUBLIC bn2 ! called by step module 56 PUBLIC eos_alpbet ! called by ldfslp module 57 PUBLIC tfreez ! called by sbcice_... modules 57 PUBLIC eos ! called by step, istate, tranpc and zpsgrd modules 58 PUBLIC eos_init ! called by istate module 59 PUBLIC bn2 ! called by step module 60 PUBLIC eos_alpbet ! called by ldfslp module 61 PUBLIC eos_alpha_beta ! used for density diagnostics in dyntra 62 PUBLIC eos_pen ! used for pe diagnostics in trdpen 63 PUBLIC tfreez ! called by sbcice_... modules 58 64 59 65 ! !!* Namelist (nameos) * 60 INTEGER , PUBLIC :: nn_eos = 0 !: = 0/1/2 type of eq. of state and Brunt-Vaisala frequ.66 INTEGER , PUBLIC :: nn_eos = 0 !: = -1/0/1/2/3 type of eq. of state and associated Brunt-Vaisala frequ. 61 67 REAL(wp), PUBLIC :: rn_alpha = 2.0e-4_wp !: thermal expension coeff. (linear equation of state) 62 68 REAL(wp), PUBLIC :: rn_beta = 7.7e-4_wp !: saline expension coeff. (linear equation of state) 63 69 64 REAL(wp), PUBLIC :: ralpbet !: alpha / beta ratio 70 REAL(wp) :: vallis_ratio = 0 ! 1027/rau0 71 REAL(wp) :: vallis_diff = 0 ! (1027-rau0)/rau0 65 72 66 73 !! * Substitutions … … 82 89 !! defined through the namelist parameter nn_eos. 83 90 !! 84 !! ** Method : 3 cases: 85 !! nn_eos = 0 : Jackett and McDougall (1994) equation of state. 91 !! ** Method : 5 cases: 92 !! nn_eos = -1 : Jackett and McDougall (1995) equation of state. 93 !! Check value: rho = 1041.83267 kg/m**3 for p=3000 dbar, 94 !! t = 3 deg celcius, s=35.5 psu 95 !! nn_eos = 0 : standard NEMO equation of state, based on 96 !! Jackett and McDougall (1995) equation of state. 86 97 !! the in situ density is computed directly as a function of 87 98 !! potential temperature relative to the surface (the opa t … … 103 114 !! salinity 104 115 !! prd(t,s) = rn_beta * s - rn_alpha * tn - 1. 116 !! nn_eos = 3 : Vallis simplified equation of state (Vallis book, 2006) 117 !! prd(t,s,p) = ( rho(t,s,p) - rau0 ) / rau0 118 !! where the pressure p in decibars is approximated by the depth in meters. 105 119 !! Note that no boundary condition problem occurs in this routine 106 120 !! as pts are defined over the whole domain. … … 108 122 !! ** Action : compute prd , the in situ density (no units) 109 123 !! 110 !! References : Jackett and McDougall, J. Atmos. Ocean. Tech., 1994111 !! ----------------------------------------------------------------------112 !! 124 !! References : Vallis, Atmospheric and Oceanic Fluid Dynamics, 2006 125 !! Jackett and McDougall, J. Atmos. Ocean. Tech., 1995 126 !!---------------------------------------------------------------------- 113 127 REAL(wp), DIMENSION(:,:,:,:), INTENT(in ) :: pts ! 1 : potential temperature [Celcius] 114 128 ! ! 2 : salinity [psu] 115 129 REAL(wp), DIMENSION(:,:,:) , INTENT( out) :: prd ! in situ density [-] 116 ! !130 ! 117 131 INTEGER :: ji, jj, jk ! dummy loop indices 118 132 REAL(wp) :: zt , zs , zh , zsr ! local scalars … … 123 137 REAL(wp), POINTER, DIMENSION(:,:,:) :: zws 124 138 !!---------------------------------------------------------------------- 125 126 139 ! 127 140 IF( nn_timing == 1 ) CALL timing_start('eos') … … 131 144 SELECT CASE( nn_eos ) 132 145 ! 133 CASE( 0 ) !== Jackett and McDougall (1994) formulation ==! 146 CASE( -1 ) !== Jackett and McDougall (1995) formulation ==! 147 ! 134 148 !CDIR NOVERRCHK 135 149 zws(:,:,:) = SQRT( ABS( pts(:,:,:,jp_sal) ) ) … … 156 170 ! 157 171 ! add the compression terms 172 ze = ( 6.207323e-11_wp*zt+6.128773e-9_wp ) *zt-2.040237e-7_wp 173 zbw= ( 1.394680e-8_wp*zt-1.202016e-6_wp ) *zt+2.102898e-5_wp 174 zb = zbw + ze * zs 175 ! 176 zd = -1.480266e-4_wp 177 zc = (-2.059331e-7_wp*zt+1.847318e-4_wp ) *zt-6.704388e-3 178 zaw= ( ( -1.956415e-6_wp*zt+2.984642e-4_wp ) *zt-2.212276e-2_wp ) *zt -3.186519_wp 179 za = ( zd*zsr + zc ) *zs + zaw 180 ! 181 zb1= (-4.619924e-3_wp*zt+9.085835e-2_wp ) *zt+3.886640_wp 182 za1= ( ( -5.084188e-4_wp*zt+6.283263e-2_wp) *zt-3.101089_wp ) *zt+528.4855_wp 183 zkw= ( ( (-4.190253e-4_wp*zt+9.648704e-2_wp ) *zt-17.06103_wp ) *zt + 1444.304_wp ) *zt+196593.3_wp 184 zk0= ( zb1*zsr + za1 )*zs + zkw 185 ! 186 ! masked in situ density anomaly. Important: rau0=1035, should be 1027! 187 prd(ji,jj,jk) = ( zrhop / ( 1.0_wp - zh / ( zk0 - zh * ( za - zh * zb ) ) ) & 188 & - rau0 ) * r1_rau0 * tmask(ji,jj,jk) 189 END DO 190 END DO 191 END DO 192 ! 193 CASE( 0 ) !== modified Jackett and McDougall (1995) formulation ==! 194 ! 195 !CDIR NOVERRCHK 196 zws(:,:,:) = SQRT( ABS( pts(:,:,:,jp_sal) ) ) 197 ! 198 DO jk = 1, jpkm1 199 DO jj = 1, jpj 200 DO ji = 1, jpi 201 zt = pts (ji,jj,jk,jp_tem) 202 zs = pts (ji,jj,jk,jp_sal) 203 zh = fsdept(ji,jj,jk) ! depth 204 zsr= zws (ji,jj,jk) ! square root salinity 205 ! 206 ! compute volumic mass pure water at atm pressure 207 zr1= ( ( ( ( 6.536332e-9_wp *zt - 1.120083e-6_wp )*zt + 1.001685e-4_wp )*zt & 208 & -9.095290e-3_wp )*zt + 6.793952e-2_wp )*zt + 999.842594_wp 209 ! seawater volumic mass atm pressure 210 zr2= ( ( ( 5.3875e-9_wp*zt-8.2467e-7_wp ) *zt+7.6438e-5_wp ) *zt & 211 & -4.0899e-3_wp ) *zt+0.824493_wp 212 zr3= ( -1.6546e-6_wp*zt+1.0227e-4_wp ) *zt-5.72466e-3_wp 213 zr4= 4.8314e-4_wp 214 ! 215 ! potential volumic mass (reference to the surface) 216 zrhop= ( zr4*zs + zr3*zsr + zr2 ) *zs + zr1 217 ! 218 ! add the compression terms 158 219 ze = ( -3.508914e-8_wp*zt-1.248266e-8_wp ) *zt-2.595994e-6_wp 159 220 zbw= ( 1.296821e-6_wp*zt-5.782165e-9_wp ) *zt+1.045941e-4_wp … … 187 248 END DO 188 249 ! 250 CASE( 3 ) !== Vallis (2006) simplified EOS ==! 251 DO jk = 1, jpkm1 252 DO jj = 1, jpj 253 DO ji = 1, jpi 254 zt = pts (ji,jj,jk,jp_tem) - 10._wp 255 zs = pts (ji,jj,jk,jp_sal) - 35._wp 256 zh = fsdept(ji,jj,jk) ! depth 257 ! masked in situ density anomaly 258 prd(ji,jj,jk) = vallis_diff + vallis_ratio * ( & 259 & -(1.67e-4_wp*(1._wp+1.1e-4_wp*zh)+0.5e-5_wp*zt)*zt + 0.78e-3_wp*zs + 4.39e-6_wp*zh & 260 & ) * tmask(ji,jj,jk) 261 END DO 262 END DO 263 END DO 264 ! 189 265 END SELECT 190 266 ! … … 207 283 !! namelist parameter nn_eos. 208 284 !! 209 !! ** Method : 210 !! nn_eos = 0 : Jackett and McDougall (1994) equation of state. 285 !! ** Method : 5 cases: 286 !! nn_eos = -1 : Jackett and McDougall (1995) equation of state. 287 !! nn_eos = 0 : standard NEMO equation of state, based on 288 !! Jackett and McDougall (1995) equation of state. 211 289 !! the in situ density is computed directly as a function of 212 290 !! potential temperature relative to the surface (the opa t … … 235 313 !! = rn_beta * s - rn_alpha * tn - 1. 236 314 !! rhop(t,s) = rho(t,s) 315 !! nn_eos = 3 : Vallis simplified equation of state (Vallis book, 2006) 316 !! prd(t,s,p) = ( rho(t,s,p) - rau0 ) / rau0 317 !! rhop(t,s) = rho(t,s,0) 237 318 !! Note that no boundary condition problem occurs in this routine 238 319 !! as (tn,sn) or (ta,sa) are defined over the whole domain. … … 241 322 !! - prhop, the potential volumic mass (Kg/m3) 242 323 !! 243 !! References : Jackett and McDougall, J. Atmos. Ocean. Tech., 1994 324 !! References : Vallis, Atmospheric and Oceanic Fluid Dynamics, 2006 325 !! Jackett and McDougall, J. Atmos. Ocean. Tech., 1995 244 326 !! Brown and Campana, Mon. Weather Rev., 1978 245 327 !!---------------------------------------------------------------------- … … 262 344 SELECT CASE ( nn_eos ) 263 345 ! 264 CASE( 0 ) !== Jackett and McDougall (1994) formulation ==!346 CASE( -1 ) !== Jackett and McDougall (1995) formulation ==! 265 347 !CDIR NOVERRCHK 266 348 zws(:,:,:) = SQRT( ABS( pts(:,:,:,jp_sal) ) ) … … 290 372 ! 291 373 ! add the compression terms 374 ze = ( 6.207323e-11_wp*zt+6.128773e-9_wp ) *zt-2.040237e-7_wp 375 zbw= ( 1.394680e-8_wp*zt-1.202016e-6_wp ) *zt+2.102898e-5_wp 376 zb = zbw + ze * zs 377 ! 378 zd = -1.480266e-4_wp 379 zc = (-2.059331e-7_wp*zt+1.847318e-4_wp ) *zt-6.704388e-3 380 zaw= ( ( -1.956415e-6_wp*zt+2.984642e-4_wp ) *zt-2.212276e-2_wp ) *zt -3.186519_wp 381 za = ( zd*zsr + zc ) *zs + zaw 382 ! 383 zb1= (-4.619924e-3_wp*zt+9.085835e-2_wp ) *zt+3.886640_wp 384 za1= ( ( -5.084188e-4_wp*zt+6.283263e-2_wp) *zt-3.101089_wp ) *zt+528.4855_wp 385 zkw= ( ( (-4.190253e-4_wp*zt+9.648704e-2_wp ) *zt-17.06103_wp ) *zt + 1444.304_wp ) *zt+196593.3_wp 386 zk0= ( zb1*zsr + za1 )*zs + zkw 387 ! 388 ! masked in situ density anomaly 389 prd(ji,jj,jk) = ( zrhop / ( 1.0_wp - zh / ( zk0 - zh * ( za - zh * zb ) ) ) & 390 & - rau0 ) * r1_rau0 * tmask(ji,jj,jk) 391 END DO 392 END DO 393 END DO 394 ! 395 CASE( 0 ) !== modified Jackett and McDougall (1995) formulation ==! 396 !CDIR NOVERRCHK 397 zws(:,:,:) = SQRT( ABS( pts(:,:,:,jp_sal) ) ) 398 ! 399 DO jk = 1, jpkm1 400 DO jj = 1, jpj 401 DO ji = 1, jpi 402 zt = pts (ji,jj,jk,jp_tem) 403 zs = pts (ji,jj,jk,jp_sal) 404 zh = fsdept(ji,jj,jk) ! depth 405 zsr= zws (ji,jj,jk) ! square root salinity 406 ! 407 ! compute volumic mass pure water at atm pressure 408 zr1= ( ( ( ( 6.536332e-9_wp*zt-1.120083e-6_wp )*zt+1.001685e-4_wp )*zt & 409 & -9.095290e-3_wp )*zt+6.793952e-2_wp )*zt+999.842594_wp 410 ! seawater volumic mass atm pressure 411 zr2= ( ( ( 5.3875e-9_wp*zt-8.2467e-7_wp ) *zt+7.6438e-5_wp ) *zt & 412 & -4.0899e-3_wp ) *zt+0.824493_wp 413 zr3= ( -1.6546e-6_wp*zt+1.0227e-4_wp ) *zt-5.72466e-3_wp 414 zr4= 4.8314e-4_wp 415 ! 416 ! potential volumic mass (reference to the surface) 417 zrhop= ( zr4*zs + zr3*zsr + zr2 ) *zs + zr1 418 ! 419 ! save potential volumic mass 420 prhop(ji,jj,jk) = zrhop * tmask(ji,jj,jk) 421 ! 422 ! add the compression terms 292 423 ze = ( -3.508914e-8_wp*zt-1.248266e-8_wp ) *zt-2.595994e-6_wp 293 424 zbw= ( 1.296821e-6_wp*zt-5.782165e-9_wp ) *zt+1.045941e-4_wp … … 323 454 END DO 324 455 ! 456 CASE( 3 ) !== Vallis (2006) simplified EOS ==! 457 DO jk = 1, jpkm1 458 DO jj = 1, jpj 459 DO ji = 1, jpi 460 zt = pts (ji,jj,jk,jp_tem) - 10._wp 461 zs = pts (ji,jj,jk,jp_sal) - 35._wp 462 zh = fsdept(ji,jj,jk) ! depth 463 ! masked in situ density anomaly 464 prd(ji,jj,jk) = vallis_diff + vallis_ratio * ( & 465 & -(1.67e-4_wp*(1._wp+1.1e-4_wp*zh)+0.5e-5_wp*zt)*zt + 0.78e-3_wp*zs + 4.39e-6_wp*zh & 466 & ) * tmask(ji,jj,jk) 467 ! masked in situ density anomaly 468 prhop(ji,jj,jk) = ( 1.0_wp - ( 1.67e-4_wp + 0.5e-5_wp*zt ) *zt + 0.78e-3_wp *zs ) & 469 & * 1027._wp * tmask(ji,jj,jk) 470 ! 471 END DO 472 END DO 473 END DO 474 ! 325 475 END SELECT 326 476 ! … … 342 492 !! defined through the namelist parameter nn_eos. * 2D field case 343 493 !! 344 !! ** Method : 345 !! nn_eos = 0 : Jackett and McDougall (1994) equation of state. 494 !! ** Method : 5 cases: 495 !! nn_eos = -1 : Jackett and McDougall (1995) equation of state. 496 !! nn_eos = 0 : standard NEMO equation of state, based on 497 !! Jackett and McDougall (1995) equation of state. 346 498 !! the in situ density is computed directly as a function of 347 499 !! potential temperature relative to the surface (the opa t … … 363 515 !! salinity 364 516 !! prd(t,s) = rn_beta * s - rn_alpha * tn - 1. 517 !! nn_eos = 3 : Vallis simplified equation of state (Vallis book, 2006) 518 !! prd(t,s,p) = ( rho(t,s,p) - rau0 ) / rau0 519 !! where the pressure p in decibars is approximated by the depth in meters. 365 520 !! Note that no boundary condition problem occurs in this routine 366 521 !! as pts are defined over the whole domain. … … 368 523 !! ** Action : - prd , the in situ density (no units) 369 524 !! 370 !! References : Jackett and McDougall, J. Atmos. Ocean. Tech., 1994 525 !! References : Vallis, Atmospheric and Oceanic Fluid Dynamics, 2006 526 !! Jackett and McDougall, J. Atmos. Ocean. Tech., 1995 371 527 !!---------------------------------------------------------------------- 372 528 !! … … 391 547 SELECT CASE( nn_eos ) 392 548 ! 393 CASE( 0 ) !== Jackett and McDougall (1994) formulation ==!549 CASE( -1 ) !== Jackett and McDougall (1995) formulation ==! 394 550 ! 395 551 !CDIR NOVERRCHK … … 403 559 DO ji = 1, fs_jpim1 ! vector opt. 404 560 zmask = tmask(ji,jj,1) ! land/sea bottom mask = surf. mask 405 zt = pts (ji,jj,jp_tem) ! interpolated T 406 zs = pts (ji,jj,jp_sal) ! interpolated S 561 zt = pts (ji,jj,jp_tem) ! interpolated T 562 zs = pts (ji,jj,jp_sal) ! interpolated S 563 zsr = zws (ji,jj) ! square root of interpolated S 564 zh = pdep (ji,jj) ! depth at the partial step level 565 ! 566 ! compute volumic mass pure water at atm pressure 567 zr1 = ( ( ( ( 6.536332e-9_wp*zt-1.120083e-6_wp )*zt+1.001685e-4_wp )*zt & 568 & -9.095290e-3_wp )*zt+6.793952e-2_wp )*zt+999.842594_wp 569 ! seawater volumic mass atm pressure 570 zr2 = ( ( ( 5.3875e-9_wp*zt-8.2467e-7_wp )*zt+7.6438e-5_wp ) *zt & 571 & -4.0899e-3_wp ) *zt+0.824493_wp 572 zr3 = ( -1.6546e-6_wp*zt+1.0227e-4_wp ) *zt-5.72466e-3_wp 573 zr4 = 4.8314e-4_wp 574 ! 575 ! potential volumic mass (reference to the surface) 576 zrhop= ( zr4*zs + zr3*zsr + zr2 ) *zs + zr1 577 ! 578 ! add the compression terms 579 ze = ( 6.207323e-11_wp*zt+6.128773e-9_wp ) *zt-2.040237e-7_wp 580 zbw= ( 1.394680e-8_wp*zt-1.202016e-6_wp ) *zt+2.102898e-5_wp 581 zb = zbw + ze * zs 582 ! 583 zd = -1.480266e-4_wp 584 zc = (-2.059331e-7_wp*zt+1.847318e-4_wp ) *zt-6.704388e-3 585 zaw= ( ( -1.956415e-6_wp*zt+2.984642e-4_wp ) *zt-2.212276e-2_wp ) *zt -3.186519_wp 586 za = ( zd*zsr + zc ) *zs + zaw 587 ! 588 zb1= (-4.619924e-3_wp*zt+9.085835e-2_wp ) *zt+3.886640_wp 589 za1= ( ( -5.084188e-4_wp*zt+6.283263e-2_wp) *zt-3.101089_wp ) *zt+528.4855_wp 590 zkw= ( ( (-4.190253e-4_wp*zt+9.648704e-2_wp ) *zt-17.06103_wp ) *zt + 1444.304_wp ) *zt+196593.3_wp 591 zk0= ( zb1*zsr + za1 )*zs + zkw 592 ! 593 ! masked in situ density anomaly 594 prd(ji,jj) = ( zrhop / ( 1.0_wp - zh / ( zk0 - zh * ( za - zh * zb ) ) ) - rau0 ) / rau0 * zmask 595 END DO 596 END DO 597 ! 598 CASE( 0 ) !== modified Jackett and McDougall (1995) formulation ==! 599 ! 600 !CDIR NOVERRCHK 601 DO jj = 1, jpjm1 602 !CDIR NOVERRCHK 603 DO ji = 1, fs_jpim1 ! vector opt. 604 zws(ji,jj) = SQRT( ABS( pts(ji,jj,jp_sal) ) ) 605 END DO 606 END DO 607 DO jj = 1, jpjm1 608 DO ji = 1, fs_jpim1 ! vector opt. 609 zmask = tmask(ji,jj,1) ! land/sea bottom mask = surf. mask 610 zt = pts (ji,jj,jp_tem) ! interpolated T 611 zs = pts (ji,jj,jp_sal) ! interpolated S 407 612 zsr = zws (ji,jj) ! square root of interpolated S 408 613 zh = pdep (ji,jj) ! depth at the partial step level … … 455 660 END DO 456 661 ! 662 CASE( 3 ) !== Vallis (2006) simplified EOS ==! 663 DO jj = 1, jpjm1 664 DO ji = 1, fs_jpim1 ! vector opt. 665 zmask = tmask(ji,jj,1) ! land/sea bottom mask = surf. mask 666 zt = pts (ji,jj,jp_tem) - 10._wp ! interpolated T 667 zs = pts (ji,jj,jp_sal) - 35._wp ! interpolated S 668 zh = pdep (ji,jj) ! depth at the partial step level 669 ! masked in situ density anomaly 670 prd(ji,jj) = vallis_diff + vallis_ratio * ( & 671 & -(1.67e-4_wp*(1._wp+1.1e-4_wp*zh)+0.5e-5_wp*zt)*zt + 0.78e-3_wp*zs + 4.39e-6_wp*zh & 672 & ) * zmask 673 ! 674 END DO 675 END DO 676 ! 457 677 END SELECT 458 678 … … 473 693 !! step of the input arguments 474 694 !! 475 !! ** Method : 476 !! * nn_eos = 0 : UNESCO sea water properties 477 !! The brunt-vaisala frequency is computed using the polynomial 478 !! polynomial expression of McDougall (1987): 695 !! ** Method : 5 cases: 696 !! * nn_eos = -1 : The brunt-vaisala frequency is computed for 697 !! the Jackett and McDougall (1995) equation of state: 479 698 !! N^2 = grav * beta * ( alpha/beta*dk[ t ] - dk[ s ] )/e3w 699 !! If lk_zdfddm=T, the heat/salt buoyancy flux ratio Rrau is 700 !! computed and used in zdfddm module : 701 !! Rrau = alpha/beta * ( dk[ t ] / dk[ s ] ) 702 !! * nn_eos = 0 : The brunt-vaisala frequency is based on the 703 !! formulation of McDougall (1987). 480 704 !! If lk_zdfddm=T, the heat/salt buoyancy flux ratio Rrau is 481 705 !! computed and used in zdfddm module : … … 494 718 !! ** Action : - pn2 : the brunt-vaisala frequency 495 719 !! 496 !! References : McDougall, J. Phys. Oceanogr., 17, 1950-1964, 1987. 720 !! References : Vallis, Atmospheric and Oceanic Fluid Dynamics, 2006 721 !! Jackett and McDougall, J. Atmos. Ocean. Tech., 1995 722 !! McDougall, J. Phys. Oceano., 1987 497 723 !!---------------------------------------------------------------------- 498 724 REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in ) :: pts ! 1 : potential temperature [Celcius] … … 501 727 !! 502 728 INTEGER :: ji, jj, jk ! dummy loop indices 503 REAL(wp) :: zgde3w, zt, zs, zh, zalbet, zbeta ! local scalars 729 REAL(wp) :: zgde3w, zt, zs, zh, zsr, zr1, zr2, zr3, zr4, zrhop ! local scalars 730 REAL(wp) :: ze, zbw, zb, zd, zc, zaw, za, zb1, za1, zkw, zk0, zM ! - - 731 REAL(wp) :: zDenom, zrho, zdzk0, zdza, zdzb, zdrhop, zdM, zdDenom ! - - 732 REAL(wp) :: zalpbet, zalpha, zbeta ! - - 504 733 #if defined key_zdfddm 505 734 REAL(wp) :: zds ! local scalars 506 735 #endif 736 REAL(wp), POINTER, DIMENSION(:,:,:) :: zws 507 737 !!---------------------------------------------------------------------- 508 738 509 739 ! 510 740 IF( nn_timing == 1 ) CALL timing_start('bn2') 741 ! 742 CALL wrk_alloc( jpi, jpj, jpk, zws ) 511 743 ! 512 744 ! pn2 : interior points only (2=< jk =< jpkm1 ) … … 515 747 SELECT CASE( nn_eos ) 516 748 ! 517 CASE( 0 ) !== Jackett and McDougall (1994) formulation ==! 749 CASE( -1 ) !== Jackett and McDougall (1995) ==! 750 zws(:,:,:) = SQRT( ABS( pts(:,:,:,jp_sal) ) ) 751 DO jk = 2, jpkm1 752 DO jj = 1, jpj 753 DO ji = 1, jpi 754 zgde3w = grav / fse3w(ji,jj,jk) 755 zt = 0.5 * ( pts(ji,jj,jk,jp_tem) + pts(ji,jj,jk-1,jp_tem) ) ! potential temperature at w-pt 756 zs = 0.5 * ( pts(ji,jj,jk,jp_sal) + pts(ji,jj,jk-1,jp_sal) ) ! salinity at w-pt 757 zsr = 0.5 * ( zws(ji,jj,jk) + zws(ji,jj,jk-1) ) ! square root of S at w-pt 758 zh = fsdepw(ji,jj,jk) ! depth in meters at w-point 759 ! 760 ! compute volumic mass pure water at atm pressure 761 zr1= ( ( ( ( 6.536332e-9_wp *zt - 1.120083e-6_wp )*zt + 1.001685e-4_wp )*zt & 762 & -9.095290e-3_wp )*zt + 6.793952e-2_wp )*zt + 999.842594_wp 763 ! seawater volumic mass atm pressure 764 zr2= ( ( ( 5.3875e-9_wp*zt-8.2467e-7_wp ) *zt+7.6438e-5_wp ) *zt & 765 & -4.0899e-3_wp ) *zt+0.824493_wp 766 zr3= ( -1.6546e-6_wp*zt+1.0227e-4_wp ) *zt-5.72466e-3_wp 767 zr4= 4.8314e-4_wp 768 ! 769 ! potential volumic mass (reference to the surface) 770 zrhop= ( zr4*zs + zr3*zsr + zr2 ) *zs + zr1 771 ! 772 ! add the compression terms 773 ze = ( 6.207323e-11_wp*zt+6.128773e-9_wp ) *zt-2.040237e-7_wp 774 zbw= ( 1.394680e-8_wp*zt-1.202016e-6_wp ) *zt+2.102898e-5_wp 775 zb = zbw + ze * zs 776 ! 777 zd = -1.480266e-4_wp 778 zc = (-2.059331e-7_wp*zt+1.847318e-4_wp ) *zt-6.704388e-3 779 zaw= ( ( -1.956415e-6_wp*zt+2.984642e-4_wp ) *zt-2.212276e-2_wp ) *zt -3.186519_wp 780 za = ( zd*zsr + zc ) *zs + zaw 781 ! 782 zb1= (-4.619924e-3_wp*zt+9.085835e-2_wp ) *zt+3.886640_wp 783 za1= ( ( -5.084188e-4_wp*zt+6.283263e-2_wp) *zt-3.101089_wp ) *zt+528.4855_wp 784 zkw= ( ( (-4.190253e-4_wp*zt+9.648704e-2_wp ) *zt-17.06103_wp ) *zt + 1444.304_wp ) *zt+196593.3_wp 785 zk0= ( zb1*zsr + za1 )*zs + zkw 786 ! 787 zM= ( zb*zh - za )*zh + zk0 788 zDenom= 1._wp - zh / zM 789 ! compute in-situ density 790 zrho = zrhop/zDenom 791 ! alpha 792 zdzb = (2.593642e-6_wp*zt-5.782165e-9_wp) + (-7.017828e-8_wp*zt-1.248266e-8_wp) * zs 793 zdza = (-14.535852e-5_wp*zt+2.598241e-3)*zs+((17.81973e-6_wp*zt+5.025098e-3_wp)*zt-0.1028859_wp) 794 zdzk0 = ((-0.3818156_wp*zt+7.390729_wp)*zsr+((6.979407e-3_wp*zt+3.10638_wp)*zt-65.00517_wp))*zs & 795 & +(((-5.446516e-4_wp*zt-5.558196e-2_wp)*zt-60.83276_wp)*zt+2098.925_wp) 796 zdrhop = ((-3.3092e-6_wp*zt+1.0227e-4_wp)*zsr & 797 & +(((21.55e-9_wp*zt-24.7401e-7_wp)*zt+15.2876e-5_wp)*zt-4.0899e-3_wp))*zs & 798 & +((((32.68166e-9_wp*zt-4.480332e-6_wp)*zt+3.005055e-4_wp)*zt-18.19058e-3_wp)*zt+6.793952e-2_wp) 799 zdM = (zdzb*zh-zdza)*zh+zdzk0 800 zdDenom = zh * zdM / (zM*zM) 801 zalpha = - ( zdrhop - zrho*zdDenom ) / zDenom / rau0 802 ! beta 803 zdzb = ze 804 zdza = -3.0644505e-2_wp*zsr+zc 805 zdzk0 = 1.5_wp*zb1*zsr+za1 806 zdrhop = 9.6628e-4_wp*zs+1.5_wp*zr3*zsr+zr2 807 zdM = (zdzb*zh-zdza)*zh+zdzk0 808 zdDenom = zh * zdM / (zM*zM) 809 zbeta = ( zdrhop - zrho*zdDenom ) / zDenom / rau0 810 ! alpha/beta 811 zalpbet = zalpha / zbeta 812 ! N2 813 pn2(ji,jj,jk) = zgde3w * zbeta * tmask(ji,jj,jk) & ! N^2 814 & * ( zalpbet * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) ) & 815 & - ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) ) ) 816 #if defined key_zdfddm 817 ! !!bug **** caution a traiter zds=dk[S]= 0 !!!! 818 zds = ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) ) ! Rrau = (alpha / beta) (dk[t] / dk[s]) 819 IF ( ABS( zds) <= 1.e-20_wp ) zds = 1.e-20_wp 820 rrau(ji,jj,jk) = zalpbet * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) ) / zds 821 #endif 822 END DO 823 END DO 824 END DO 825 ! 826 CASE( 0 ) !== McDougall (1987) formulation ==! 518 827 DO jk = 2, jpkm1 519 828 DO jj = 1, jpj … … 524 833 zh = fsdepw(ji,jj,jk) ! depth in meters at w-point 525 834 ! 526 zal bet = ( ( ( - 0.255019e-07_wp * zt + 0.298357e-05_wp ) * zt& ! ratio alpha/beta835 zalpbet = ( ( ( - 0.255019e-07_wp * zt + 0.298357e-05_wp ) * zt & ! ratio alpha/beta 527 836 & - 0.203814e-03_wp ) * zt & 528 837 & + 0.170907e-01_wp ) * zt & … … 550 859 ! 551 860 pn2(ji,jj,jk) = zgde3w * zbeta * tmask(ji,jj,jk) & ! N^2 552 & * ( zal bet * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) ) &861 & * ( zalpbet * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) ) & 553 862 & - ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) ) ) 554 863 #if defined key_zdfddm … … 556 865 zds = ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) ) ! Rrau = (alpha / beta) (dk[t] / dk[s]) 557 866 IF ( ABS( zds) <= 1.e-20_wp ) zds = 1.e-20_wp 558 rrau(ji,jj,jk) = zal bet * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) ) / zds867 rrau(ji,jj,jk) = zalpbet * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) ) / zds 559 868 #endif 560 869 END DO … … 564 873 CASE( 1 ) !== Linear formulation = F( temperature ) ==! 565 874 DO jk = 2, jpkm1 566 pn2(:,:,jk) = grav * rn_alpha * ( pts(:,:,jk-1,jp_tem) - pts(:,:,jk,jp_tem) ) / fse3w(:,:,jk) * tmask(:,:,jk) 875 pn2(:,:,jk) = grav * rn_alpha * ( pts(:,:,jk-1,jp_tem) - pts(:,:,jk,jp_tem) ) & 876 & / fse3w(:,:,jk) * tmask(:,:,jk) 567 877 END DO 568 878 ! … … 579 889 zds = ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) ) 580 890 IF ( ABS( zds ) <= 1.e-20_wp ) zds = 1.e-20_wp 581 rrau(ji,jj,jk) = r alpbet* ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) ) / zds891 rrau(ji,jj,jk) = rn_alpha / rn_beta * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) ) / zds 582 892 END DO 583 893 END DO 584 894 END DO 585 895 #endif 896 ! 897 CASE( 3 ) !== Vallis (2006) simplified EOS ==! 898 DO jk = 2, jpkm1 899 DO jj = 1, jpj 900 DO ji = 1, jpi 901 zgde3w = grav / fse3w(ji,jj,jk) 902 zt = 0.5 * ( pts(ji,jj,jk,jp_tem) + pts(ji,jj,jk-1,jp_tem) ) - 10._wp ! potential temperature at w-pt 903 zs = 0.5 * ( pts(ji,jj,jk,jp_sal) + pts(ji,jj,jk-1,jp_sal) ) - 35._wp ! salinity anomaly (s-35) at w-pt 904 zh = fsdepw(ji,jj,jk) ! depth in meters at w-point 905 ! 906 zalpha = 1027._wp * ( 1.67e-4_wp * ( 1._wp + 1.1e-4_wp*zh ) + 1.e-5_wp*zt ) / rau0 ! alpha 907 zbeta = 1027._wp * 0.78e-3_wp / rau0 ! beta 908 zalpbet = zalpha / zbeta ! ratio alpha/beta 909 ! 910 pn2(ji,jj,jk) = zgde3w * zbeta * tmask(ji,jj,jk) & ! N^2 911 & * ( zalpbet * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) ) & 912 & - ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) ) ) 913 #if defined key_zdfddm 914 ! !!bug **** caution a traiter zds=dk[S]= 0 !!!! 915 zds = ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) ) ! Rrau = (alpha / beta) (dk[t] / dk[s]) 916 IF ( ABS( zds) <= 1.e-20_wp ) zds = 1.e-20_wp 917 rrau(ji,jj,jk) = zalpbet * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) ) / zds 918 #endif 919 END DO 920 END DO 921 END DO 922 ! 586 923 END SELECT 587 924 … … 591 928 #endif 592 929 ! 930 CALL wrk_dealloc( jpi, jpj, jpk, zws ) 931 ! 593 932 IF( nn_timing == 1 ) CALL timing_stop('bn2') 594 933 ! … … 603 942 !! 604 943 !! ** Method : calculates alpha / beta ratio at T-points 605 !! * nn_eos = 0 : UNESCO sea water properties 606 !! The alpha/beta ratio is returned as 3-D array palpbet using the polynomial 607 !! polynomial expression of McDougall (1987). 944 !! * nn_eos = -1 : alpha / beta ratio is computed for 945 !! the Jackett and McDougall (1995) equation of state: 946 !! Scalar beta0 is returned = 1. 947 !! * nn_eos = 0 : alpha / beta ratio using the formulation of McDougall (1987). 608 948 !! Scalar beta0 is returned = 1. 609 949 !! * nn_eos = 1 : linear equation of state (temperature only) … … 611 951 !! Scalar beta0 is returned = 0. 612 952 !! * nn_eos = 2 : linear equation of state (temperature & salinity) 613 !! The alpha/beta ratio is returned as ralpbet 953 !! The alpha/beta ratio is returned as palpbet 954 !! Scalar beta0 is returned = 1. 955 !! * nn_eos = 3 : Vallis equation of state (Vallis 2006) 614 956 !! Scalar beta0 is returned = 1. 615 957 !! … … 622 964 !! 623 965 INTEGER :: ji, jj, jk ! dummy loop indices 624 REAL(wp) :: zt, zs, zh ! local scalars 966 REAL(wp) :: zt, zs, zh, zsr, zr1, zr2, zr3, zr4, zrhop ! temporary scalars 967 REAL(wp) :: ze, zbw, zb, zd, zc, zaw, za, zb1, za1, zkw, zk0, zM ! - - 968 REAL(wp) :: zDenom, zrho, zdzk0, zdza, zdzb, zdrhop, zdM, zdDenom ! 969 REAL(wp) :: zalpha, zbeta ! local scalars 970 REAL(wp), POINTER, DIMENSION(:,:,:) :: zws 625 971 !!---------------------------------------------------------------------- 626 972 ! 627 973 IF( nn_timing == 1 ) CALL timing_start('eos_alpbet') 628 974 ! 975 CALL wrk_alloc( jpi, jpj, jpk, zws ) 976 ! 629 977 SELECT CASE ( nn_eos ) 630 978 ! 631 CASE ( 0 ) ! Jackett and McDougall (1994) formulation 979 CASE ( -1 ) ! Jackett and McDougall (1995) formulation 980 zws(:,:,:) = SQRT( ABS( pts(:,:,:,jp_sal) ) ) 981 DO jk = 1, jpkm1 982 DO jj = 1, jpj 983 DO ji = 1, jpi 984 zt = pts (ji,jj,jk,jp_tem) 985 zs = pts (ji,jj,jk,jp_sal) 986 zh = fsdept(ji,jj,jk) ! depth 987 zsr= zws (ji,jj,jk) ! square root salinity 988 ! 989 ! compute volumic mass pure water at atm pressure 990 zr1= ( ( ( ( 6.536332e-9_wp *zt - 1.120083e-6_wp )*zt + 1.001685e-4_wp )*zt & 991 & -9.095290e-3_wp )*zt + 6.793952e-2_wp )*zt + 999.842594_wp 992 ! seawater volumic mass atm pressure 993 zr2= ( ( ( 5.3875e-9_wp*zt-8.2467e-7_wp ) *zt+7.6438e-5_wp ) *zt & 994 & -4.0899e-3_wp ) *zt+0.824493_wp 995 zr3= ( -1.6546e-6_wp*zt+1.0227e-4_wp ) *zt-5.72466e-3_wp 996 zr4= 4.8314e-4_wp 997 ! 998 ! potential volumic mass (reference to the surface) 999 zrhop= ( zr4*zs + zr3*zsr + zr2 ) *zs + zr1 1000 ! 1001 ! add the compression terms 1002 ze = ( 6.207323e-11_wp*zt+6.128773e-9_wp ) *zt-2.040237e-7_wp 1003 zbw= ( 1.394680e-8_wp*zt-1.202016e-6_wp ) *zt+2.102898e-5_wp 1004 zb = zbw + ze * zs 1005 ! 1006 zd = -1.480266e-4_wp 1007 zc = (-2.059331e-7_wp*zt+1.847318e-4_wp ) *zt-6.704388e-3 1008 zaw= ( ( -1.956415e-6_wp*zt+2.984642e-4_wp ) *zt-2.212276e-2_wp ) *zt -3.186519_wp 1009 za = ( zd*zsr + zc ) *zs + zaw 1010 ! 1011 zb1= (-4.619924e-3_wp*zt+9.085835e-2_wp ) *zt+3.886640_wp 1012 za1= ( ( -5.084188e-4_wp*zt+6.283263e-2_wp) *zt-3.101089_wp ) *zt+528.4855_wp 1013 zkw= ( ( (-4.190253e-4_wp*zt+9.648704e-2_wp ) *zt-17.06103_wp ) *zt + 1444.304_wp ) *zt+196593.3_wp 1014 zk0= ( zb1*zsr + za1 )*zs + zkw 1015 ! 1016 zM= ( zb*zh - za )*zh + zk0 1017 zDenom= 1._wp - zh / zM 1018 ! compute in-situ density 1019 zrho = zrhop/zDenom 1020 ! alpha 1021 zdzb = (2.593642e-6_wp*zt-5.782165e-9_wp) + (-7.017828e-8_wp*zt-1.248266e-8_wp) * zs 1022 zdza = (-14.535852e-5_wp*zt+2.598241e-3)*zs+((17.81973e-6_wp*zt+5.025098e-3_wp)*zt-0.1028859_wp) 1023 zdzk0 = ((-0.3818156_wp*zt+7.390729_wp)*zsr+((6.979407e-3_wp*zt+3.10638_wp)*zt-65.00517_wp))*zs & 1024 & +(((-5.446516e-4_wp*zt-5.558196e-2_wp)*zt-60.83276_wp)*zt+2098.925_wp) 1025 zdrhop = ((-3.3092e-6_wp*zt+1.0227e-4_wp)*zsr & 1026 & +(((21.55e-9_wp*zt-24.7401e-7_wp)*zt+15.2876e-5_wp)*zt-4.0899e-3_wp))*zs & 1027 & +((((32.68166e-9_wp*zt-4.480332e-6_wp)*zt+3.005055e-4_wp)*zt-18.19058e-3_wp)*zt+6.793952e-2_wp) 1028 zdM = (zdzb*zh-zdza)*zh+zdzk0 1029 zdDenom = zh * zdM / (zM*zM) 1030 zalpha = - ( zdrhop - zrho*zdDenom ) / zDenom / rau0 1031 ! beta 1032 zdzb = ze 1033 zdza = -3.0644505e-2_wp*zsr+zc 1034 zdzk0 = 1.5_wp*zb1*zsr+za1 1035 zdrhop = 9.6628e-4_wp*zs+1.5_wp*zr3*zsr+zr2 1036 zdM = (zdzb*zh-zdza)*zh+zdzk0 1037 zdDenom = zh * zdM / (zM*zM) 1038 zbeta = ( zdrhop - zrho*zdDenom ) / zDenom / rau0 1039 ! alpha/beta 1040 palpbet(ji,jj,jk) = zalpha / zbeta * tmask(ji,jj,jk) 1041 END DO 1042 END DO 1043 END DO 1044 beta0 = 1._wp 1045 ! 1046 CASE ( 0 ) ! McDougall (1987) formulation 632 1047 DO jk = 1, jpk 633 1048 DO jj = 1, jpj … … 660 1075 ! 661 1076 CASE ( 2 ) !== Linear formulation = F( temperature , salinity ) ==! 662 palpbet(:,:,:) = ralpbet 1077 palpbet(:,:,:) = rn_alpha / rn_beta 1078 beta0 = 1._wp 1079 ! 1080 CASE( 3 ) !== Vallis (2006) simplified EOS ==! 1081 DO jk = 1, jpkm1 1082 DO jj = 1, jpj 1083 DO ji = 1, jpi 1084 zt = pts(ji,jj,jk,jp_tem) - 10._wp ! potential temperature 1085 zs = pts(ji,jj,jk,jp_sal) - 35._wp ! salinity anomaly (s-35) 1086 zh = fsdepw(ji,jj,jk) ! depth in meters at w-point 1087 ! 1088 zalpha = ( 1.67e-4_wp * ( 1._wp + 1.1e-4_wp*zh ) + 1.e-5_wp*zt ) ! alpha/vallis_ratio 1089 zbeta = 0.78e-3_wp ! beta/vallis_ratio 1090 ! 1091 palpbet(ji,jj,jk) = zalpha / zbeta * tmask(ji,jj,jk) 1092 END DO 1093 END DO 1094 END DO 663 1095 beta0 = 1._wp 664 1096 ! … … 670 1102 END SELECT 671 1103 ! 1104 CALL wrk_dealloc( jpi, jpj, jpk, zws ) 1105 ! 672 1106 IF( nn_timing == 1 ) CALL timing_stop('eos_alpbet') 673 1107 ! … … 675 1109 676 1110 1111 SUBROUTINE eos_alpha_beta( pts, palpha, pbeta ) 1112 !!---------------------------------------------------------------------- 1113 !! *** ROUTINE eos_alpha_beta *** 1114 !! 1115 !! ** Purpose : Calculates the in situ thermal/haline expansion terms at T-points 1116 !! 1117 !! ** Method : calculates alpha and beta at T-points 1118 !! * nn_eos = -1 : Jackett and McDougall (1995) equation of state 1119 !! * nn_eos = 0 : modified Jackett and McDougall (1995) equation of state 1120 !! * nn_eos = 1 : linear equation of state (temperature only) 1121 !! palpha = rn_alpha 1122 !! pbeta = 0 1123 !! * nn_eos = 2 : linear equation of state (temperature & salinity) 1124 !! palpha = rn_alpha 1125 !! pbeta = rn_beta 1126 !! * nn_eos = 3 : Vallis equation of state (Vallis 2006) 1127 !! alpha and beta ratios are returned as 3-D arrays. 1128 !! 1129 !! ** Action : - palpha : thermal expansion ratio at T-points 1130 !! : - pbeta : haline expansion ratio at T-points 1131 !!---------------------------------------------------------------------- 1132 REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in ) :: pts ! pot. temperature & salinity 1133 REAL(wp), DIMENSION(jpi,jpj,jpk) , INTENT( out) :: palpha ! alpha 1134 REAL(wp), DIMENSION(jpi,jpj,jpk) , INTENT( out) :: pbeta ! beta 1135 ! 1136 INTEGER :: ji, jj, jk ! dummy loop indices 1137 REAL(wp) :: zt, zs, zh, zsr, zr1, zr2, zr3, zr4, zrhop ! temporary scalars 1138 REAL(wp) :: ze, zbw, zb, zd, zc, zaw, za, zb1, za1, zkw, zk0, zM ! - - 1139 REAL(wp) :: zDenom, zrho, zdzk0, zdza, zdzb, zdrhop, zdM, zdDenom ! 1140 REAL(wp), POINTER, DIMENSION(:,:,:) :: zws 1141 !!---------------------------------------------------------------------- 1142 ! 1143 IF ( nn_timing == 1 ) CALL timing_start('eos_alpha_beta') 1144 ! 1145 CALL wrk_alloc( jpi, jpj, jpk, zws ) 1146 ! 1147 SELECT CASE ( nn_eos ) 1148 ! 1149 CASE ( -1 ) ! Jackett and McDougall (1995) formulation 1150 ! 1151 zws(:,:,:) = SQRT( ABS( pts(:,:,:,jp_sal) ) ) 1152 DO jk = 1, jpkm1 1153 DO jj = 1, jpj 1154 DO ji = 1, jpi 1155 zt = pts (ji,jj,jk,jp_tem) 1156 zs = pts (ji,jj,jk,jp_sal) 1157 zh = fsdept(ji,jj,jk) ! depth 1158 zsr= zws (ji,jj,jk) ! square root salinity 1159 ! 1160 ! compute volumic mass pure water at atm pressure 1161 zr1= ( ( ( ( 6.536332e-9_wp *zt - 1.120083e-6_wp )*zt + 1.001685e-4_wp )*zt & 1162 & -9.095290e-3_wp )*zt + 6.793952e-2_wp )*zt + 999.842594_wp 1163 ! seawater volumic mass atm pressure 1164 zr2= ( ( ( 5.3875e-9_wp*zt-8.2467e-7_wp ) *zt+7.6438e-5_wp ) *zt & 1165 & -4.0899e-3_wp ) *zt+0.824493_wp 1166 zr3= ( -1.6546e-6_wp*zt+1.0227e-4_wp ) *zt-5.72466e-3_wp 1167 zr4= 4.8314e-4_wp 1168 ! 1169 ! potential volumic mass (reference to the surface) 1170 zrhop= ( zr4*zs + zr3*zsr + zr2 ) *zs + zr1 1171 ! 1172 ! add the compression terms 1173 ze = ( 6.207323e-11_wp*zt+6.128773e-9_wp ) *zt-2.040237e-7_wp 1174 zbw= ( 1.394680e-8_wp*zt-1.202016e-6_wp ) *zt+2.102898e-5_wp 1175 zb = zbw + ze * zs 1176 ! 1177 zd = -1.480266e-4_wp 1178 zc = (-2.059331e-7_wp*zt+1.847318e-4_wp ) *zt-6.704388e-3 1179 zaw= ( ( -1.956415e-6_wp*zt+2.984642e-4_wp ) *zt-2.212276e-2_wp ) *zt -3.186519_wp 1180 za = ( zd*zsr + zc ) *zs + zaw 1181 ! 1182 zb1= (-4.619924e-3_wp*zt+9.085835e-2_wp ) *zt+3.886640_wp 1183 za1= ( ( -5.084188e-4_wp*zt+6.283263e-2_wp) *zt-3.101089_wp ) *zt+528.4855_wp 1184 zkw= ( ( (-4.190253e-4_wp*zt+9.648704e-2_wp ) *zt-17.06103_wp ) *zt + 1444.304_wp ) *zt+196593.3_wp 1185 zk0= ( zb1*zsr + za1 )*zs + zkw 1186 ! 1187 zM= ( zb*zh - za )*zh + zk0 1188 zDenom= 1._wp - zh / zM 1189 ! compute in-situ density 1190 zrho = zrhop/zDenom 1191 ! alpha 1192 zdzb = (2.593642e-6_wp*zt-5.782165e-9_wp) + (-7.017828e-8_wp*zt-1.248266e-8_wp) * zs 1193 zdza = (-14.535852e-5_wp*zt+2.598241e-3)*zs+((17.81973e-6_wp*zt+5.025098e-3_wp)*zt-0.1028859_wp) 1194 zdzk0 = ((-0.3818156_wp*zt+7.390729_wp)*zsr+((6.979407e-3_wp*zt+3.10638_wp)*zt-65.00517_wp))*zs & 1195 & +(((-5.446516e-4_wp*zt-5.558196e-2_wp)*zt-60.83276_wp)*zt+2098.925_wp) 1196 zdrhop = ((-3.3092e-6_wp*zt+1.0227e-4_wp)*zsr & 1197 & +(((21.55e-9_wp*zt-24.7401e-7_wp)*zt+15.2876e-5_wp)*zt-4.0899e-3_wp))*zs & 1198 & +((((32.68166e-9_wp*zt-4.480332e-6_wp)*zt+3.005055e-4_wp)*zt-18.19058e-3_wp)*zt+6.793952e-2_wp) 1199 zdM = (zdzb*zh-zdza)*zh+zdzk0 1200 zdDenom = zh * zdM / (zM*zM) 1201 palpha(ji,jj,jk) = - ( zdrhop - zrho*zdDenom ) / zDenom / rau0 * tmask(ji,jj,jk) 1202 ! beta 1203 zdzb = ze 1204 zdza = -3.0644505e-2_wp*zsr+zc 1205 zdzk0 = 1.5_wp*zb1*zsr+za1 1206 zdrhop = 9.6628e-4_wp*zs+1.5_wp*zr3*zsr+zr2 1207 zdM = (zdzb*zh-zdza)*zh+zdzk0 1208 zdDenom = zh * zdM / (zM*zM) 1209 pbeta(ji,jj,jk) = ( zdrhop - zrho*zdDenom ) / zDenom / rau0 * tmask(ji,jj,jk) 1210 END DO 1211 END DO 1212 END DO 1213 ! 1214 CASE ( 0 ) ! modified Jackett and McDougall (1995) formulation 1215 ! 1216 zws(:,:,:) = SQRT( ABS( pts(:,:,:,jp_sal) ) ) 1217 DO jk = 1, jpkm1 1218 DO jj = 1, jpj 1219 DO ji = 1, jpi 1220 zt = pts (ji,jj,jk,jp_tem) 1221 zs = pts (ji,jj,jk,jp_sal) 1222 zh = fsdept(ji,jj,jk) ! depth 1223 zsr= zws (ji,jj,jk) ! square root salinity 1224 ! 1225 ! compute volumic mass pure water at atm pressure 1226 zr1= ( ( ( ( 6.536332e-9_wp *zt - 1.120083e-6_wp )*zt + 1.001685e-4_wp )*zt & 1227 & -9.095290e-3_wp )*zt + 6.793952e-2_wp )*zt + 999.842594_wp 1228 ! seawater volumic mass atm pressure 1229 zr2= ( ( ( 5.3875e-9_wp*zt-8.2467e-7_wp ) *zt+7.6438e-5_wp ) *zt & 1230 & -4.0899e-3_wp ) *zt+0.824493_wp 1231 zr3= ( -1.6546e-6_wp*zt+1.0227e-4_wp ) *zt-5.72466e-3_wp 1232 zr4= 4.8314e-4_wp 1233 ! 1234 ! potential volumic mass (reference to the surface) 1235 zrhop= ( zr4*zs + zr3*zsr + zr2 ) *zs + zr1 1236 ! 1237 ! add the compression terms 1238 ze = ( -3.508914e-8_wp*zt-1.248266e-8_wp ) *zt-2.595994e-6_wp 1239 zbw= ( 1.296821e-6_wp*zt-5.782165e-9_wp ) *zt+1.045941e-4_wp 1240 zb = zbw + ze * zs 1241 ! 1242 zd = -2.042967e-2_wp 1243 zc = (-7.267926e-5_wp*zt+2.598241e-3_wp ) *zt+0.1571896_wp 1244 zaw= ( ( 5.939910e-6_wp*zt+2.512549e-3_wp ) *zt-0.1028859_wp ) *zt - 4.721788_wp 1245 za = ( zd*zsr + zc ) *zs + zaw 1246 ! 1247 zb1= (-0.1909078_wp*zt+7.390729_wp ) *zt-55.87545_wp 1248 za1= ( ( 2.326469e-3_wp*zt+1.553190_wp) *zt-65.00517_wp ) *zt+1044.077_wp 1249 zkw= ( ( (-1.361629e-4_wp*zt-1.852732e-2_wp ) *zt-30.41638_wp ) *zt + 2098.925_wp ) *zt+190925.6_wp 1250 zk0= ( zb1*zsr + za1 )*zs + zkw 1251 ! 1252 zM= ( zb*zh - za )*zh + zk0 1253 zDenom= 1._wp - zh / zM 1254 ! compute in-situ density 1255 zrho = zrhop/zDenom 1256 ! alpha 1257 zdzb = (2.593642e-6_wp*zt-5.782165e-9_wp) + (-7.017828e-8_wp*zt-1.248266e-8_wp) * zs 1258 zdza = (-14.535852e-5_wp*zt+2.598241e-3)*zs+((17.81973e-6_wp*zt+5.025098e-3_wp)*zt-0.1028859_wp) 1259 zdzk0 = ((-0.3818156_wp*zt+7.390729_wp)*zsr+((6.979407e-3_wp*zt+3.10638_wp)*zt-65.00517_wp))*zs & 1260 & +(((-5.446516e-4_wp*zt-5.558196e-2_wp)*zt-60.83276_wp)*zt+2098.925_wp) 1261 zdrhop = ((-3.3092e-6_wp*zt+1.0227e-4_wp)*zsr & 1262 & +(((21.55e-9_wp*zt-24.7401e-7_wp)*zt+15.2876e-5_wp)*zt-4.0899e-3_wp))*zs & 1263 & +((((32.68166e-9_wp*zt-4.480332e-6_wp)*zt+3.005055e-4_wp)*zt-18.19058e-3_wp)*zt+6.793952e-2_wp) 1264 zdM = (zdzb*zh-zdza)*zh+zdzk0 1265 zdDenom = zh * zdM / (zM*zM) 1266 palpha(ji,jj,jk) = - ( zdrhop - zrho*zdDenom ) / zDenom / rau0 * tmask(ji,jj,jk) 1267 ! beta 1268 zdzb = ze 1269 zdza = -3.0644505e-2_wp*zsr+zc 1270 zdzk0 = 1.5_wp*zb1*zsr+za1 1271 zdrhop = 9.6628e-4_wp*zs+1.5_wp*zr3*zsr+zr2 1272 zdM = (zdzb*zh-zdza)*zh+zdzk0 1273 zdDenom = zh * zdM / (zM*zM) 1274 pbeta(ji,jj,jk) = ( zdrhop - zrho*zdDenom ) / zDenom / rau0 * tmask(ji,jj,jk) 1275 END DO 1276 END DO 1277 END DO 1278 ! 1279 CASE ( 1 ) !== Linear formulation = F( temperature ) ==! 1280 palpha(:,:,:) = rn_alpha 1281 pbeta (:,:,:) = 0._wp 1282 ! 1283 CASE ( 2 ) !== Linear formulation = F( temperature , salinity ) ==! 1284 palpha(:,:,:) = rn_alpha 1285 pbeta (:,:,:) = rn_beta 1286 ! 1287 CASE( 3 ) !== Vallis (2006) simplified EOS ==! 1288 DO jk = 1, jpkm1 1289 DO jj = 1, jpj 1290 DO ji = 1, jpi 1291 zt = pts(ji,jj,jk,jp_tem) - 10._wp ! potential temperature (t-10) 1292 zh = fsdept(ji,jj,jk) ! depth in meters at w-point 1293 palpha(ji,jj,jk) = vallis_ratio * & 1294 & ( 1.67e-4_wp * ( 1._wp + 1.1e-4_wp*zh ) + 1.e-5_wp*zt ) * tmask(ji,jj,jk) ! alpha 1295 ! 1296 END DO 1297 END DO 1298 END DO 1299 pbeta (:,:,:) = vallis_ratio * 0.78e-3_wp * tmask(:,:,:) ! beta 1300 ! 1301 CASE DEFAULT 1302 IF(lwp) WRITE(numout,cform_err) 1303 IF(lwp) WRITE(numout,*) ' bad flag value for nn_eos = ', nn_eos 1304 nstop = nstop + 1 1305 ! 1306 END SELECT 1307 ! 1308 CALL wrk_dealloc( jpi, jpj, jpk, zws ) 1309 ! 1310 IF( nn_timing == 1 ) CALL timing_stop('eos_alpha_beta') 1311 ! 1312 END SUBROUTINE eos_alpha_beta 1313 1314 1315 1316 SUBROUTINE eos_pen( pts, palphaPE, pbetaPE, pPEanom ) 1317 !!---------------------------------------------------------------------- 1318 !! *** ROUTINE eos_pen *** 1319 !! 1320 !! ** Purpose : Calculates alpha_PE, beta_PE and PE at T-points 1321 !! 1322 !! ** Method : PE is defined analytically as the vertical 1323 !! primitive of EOS times -g integrated between 0 and z>0. 1324 !! PEanom is the PE anomaly: PEanom = ( PE - rau0 gz ) / rau0 gz 1325 !! = -1/z * /int_z^0 rd , where rd density anomaly 1326 !! dalphaPE and dbetaPE are partial derivatives of PE anomaly with respect to T and S: 1327 !! alphaPE = - 1/(rau0 gz) * dPE/dT = - dPEanom/dT 1328 !! betaPE = 1/(rau0 gz) * dPE/dS = - dPEanom/dS 1329 !! 1330 !! * nn_eos = -1 : Jackett and McDougall (1995) equation of state. 1331 !! * nn_eos = 0 : modified Jackett and McDougall (1995) equation of state. 1332 !! * nn_eos = 1 : linear equation of state (temperature only) 1333 !! palpha = rau0*g*rn_alpha*z 1334 !! pbeta = 0 1335 !! * nn_eos = 2 : linear equation of state (temperature & salinity) 1336 !! palpha = rau0*g*rn_alpha*z 1337 !! pbeta = -rau0*g*rn_beta*z 1338 !! * nn_eos = 3 : Vallis equation of state (Vallis 2006) 1339 !! 1340 !! ** Action : - pPEanom : PE anomaly given at T-points 1341 !! : - palphaPE : given at T-points 1342 !! : - pbetaPE : given at T-points 1343 !!---------------------------------------------------------------------- 1344 REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in ) :: pts ! pot. temperature & salinity 1345 REAL(wp), DIMENSION(jpi,jpj,jpk) , INTENT( out) :: palphaPE ! partial derivative of PE anom 1346 REAL(wp), DIMENSION(jpi,jpj,jpk) , INTENT( out) :: pbetaPE ! with respect to T and S, resp. 1347 REAL(wp), DIMENSION(jpi,jpj,jpk) , INTENT( out) :: pPEanom ! potential energy anomaly 1348 ! 1349 INTEGER :: ji, jj, jk ! dummy loop indices 1350 REAL(wp) :: zt, zs, zh, zsr, zr1, zr2, zr3, zr4, zrhop ! temporary scalars 1351 REAL(wp) :: ze, zbw, zb, zd, zc, zaw, za, zb1, za1, zkw, zk0, zM ! - - 1352 REAL(wp) :: zDenom, zrho, zdzk0, zdza, zdzb, zdrhop, zdM, zdDenom ! 1353 REAL(wp) :: zap1, zdelta, zsdelta, zAp, zBp, zlogBp, zCp, zatanCp ! 1354 REAL(wp) :: zddelta, zdAp, zdBp, zdlogBp, zdCp, zP, zdP, zEp ! 1355 REAL(wp), POINTER, DIMENSION(:,:,:) :: zws 1356 !!---------------------------------------------------------------------- 1357 ! 1358 IF ( nn_timing == 1 ) CALL timing_start('eos_pen') 1359 ! 1360 CALL wrk_alloc( jpi, jpj, jpk, zws ) 1361 ! 1362 SELECT CASE ( nn_eos ) 1363 ! 1364 CASE ( -1 ) ! Jackett and McDougall (1995) formulation 1365 ! 1366 zws(:,:,:) = SQRT( ABS( pts(:,:,:,jp_sal) ) ) 1367 DO jk = 1, jpkm1 1368 DO jj = 1, jpj 1369 DO ji = 1, jpi 1370 zt = pts (ji,jj,jk,jp_tem) 1371 zs = pts (ji,jj,jk,jp_sal) 1372 zh = fsdept(ji,jj,jk) ! depth 1373 zsr= zws (ji,jj,jk) ! square root salinity 1374 ! 1375 ! compute volumic mass pure water at atm pressure 1376 zr1= ( ( ( ( 6.536332e-9_wp *zt - 1.120083e-6_wp )*zt + 1.001685e-4_wp )*zt & 1377 & -9.095290e-3_wp )*zt + 6.793952e-2_wp )*zt + 999.842594_wp 1378 ! seawater volumic mass atm pressure 1379 zr2= ( ( ( 5.3875e-9_wp*zt-8.2467e-7_wp ) *zt+7.6438e-5_wp ) *zt & 1380 & -4.0899e-3_wp ) *zt+0.824493_wp 1381 zr3= ( -1.6546e-6_wp*zt+1.0227e-4_wp ) *zt-5.72466e-3_wp 1382 zr4= 4.8314e-4_wp 1383 ! 1384 ! potential volumic mass (reference to the surface) 1385 zrhop= ( zr4*zs + zr3*zsr + zr2 ) *zs + zr1 1386 ! 1387 ! add the compression terms 1388 ze = ( 6.207323e-11_wp*zt+6.128773e-9_wp ) *zt-2.040237e-7_wp 1389 zbw= ( 1.394680e-8_wp*zt-1.202016e-6_wp ) *zt+2.102898e-5_wp 1390 zb = zbw + ze * zs 1391 ! 1392 zd = -1.480266e-4_wp 1393 zc = (-2.059331e-7_wp*zt+1.847318e-4_wp ) *zt-6.704388e-3 1394 zaw= ( ( -1.956415e-6_wp*zt+2.984642e-4_wp ) *zt-2.212276e-2_wp ) *zt -3.186519_wp 1395 za = ( zd*zsr + zc ) *zs + zaw 1396 ! 1397 zb1= (-4.619924e-3_wp*zt+9.085835e-2_wp ) *zt+3.886640_wp 1398 za1= ( ( -5.084188e-4_wp*zt+6.283263e-2_wp) *zt-3.101089_wp ) *zt+528.4855_wp 1399 zkw= ( ( (-4.190253e-4_wp*zt+9.648704e-2_wp ) *zt-17.06103_wp ) *zt + 1444.304_wp ) *zt+196593.3_wp 1400 zk0= ( zb1*zsr + za1 )*zs + zkw 1401 ! 1402 zM= ( zb*zh - za )*zh + zk0 1403 zDenom= 1._wp - zh / zM 1404 ! compute in-situ density 1405 zrho = zrhop/zDenom 1406 ! 1407 zap1 = 1._wp + za 1408 zdelta = 4._wp*zb*zk0 - zap1*zap1 1409 zsdelta = sqrt( abs( zdelta ) ) 1410 zAp = zap1 / zb / zsdelta 1411 zBp = ( zM - zh ) / zk0 1412 zlogBp = log(abs(zBp)) 1413 zCp = zh * zsdelta / (2._wp*zk0-zh*zap1) 1414 zatanCp = atan(zCp) 1415 zP = zh + zAp * zatanCp + 0.5_wp*zlogBp/zb 1416 ! compute potential energy 1417 pPEanom(ji,jj,jk) = ( ( zrhop * zP ) / rau0 / zh - 1._wp ) * tmask(ji,jj,jk) !!wrong 1418 ! 1419 ! dPEdt 1420 zdzb = (2.593642e-6_wp*zt-5.782165e-9_wp) + (-7.017828e-8_wp*zt-1.248266e-8_wp) * zs 1421 zdza = (-14.535852e-5_wp*zt+2.598241e-3)*zs+((17.81973e-6_wp*zt+5.025098e-3_wp)*zt-0.1028859_wp) 1422 zdzk0 = ((-0.3818156_wp*zt+7.390729_wp)*zsr+((6.979407e-3_wp*zt+3.10638_wp)*zt-65.00517_wp))*zs & 1423 & +(((-5.446516e-4_wp*zt-5.558196e-2_wp)*zt-60.83276_wp)*zt+2098.925_wp) 1424 zdrhop = ((-3.3092e-6_wp*zt+1.0227e-4_wp)*zsr & 1425 & +(((21.55e-9_wp*zt-24.7401e-7_wp)*zt+15.2876e-5_wp)*zt-4.0899e-3_wp))*zs & 1426 & +((((32.68166e-9_wp*zt-4.480332e-6_wp)*zt+3.005055e-4_wp)*zt-18.19058e-3_wp)*zt+6.793952e-2_wp) 1427 zdM = (zdzb*zh-zdza)*zh+zdzk0 1428 zdDenom = zh * zdM / (zM*zM) 1429 ! 1430 zddelta = 4._wp * (zk0*zdzb+zb*zdzk0) - 2._wp*zap1*zdza 1431 zdAp = zAp * (zdza/zap1 - zdzb/zb - 0.5_wp*zddelta/zdelta) 1432 zdlogBp = zdM/(zM-zh) - zdzk0/zk0 1433 zdCp = zCp * (0.5_wp*zddelta/zdelta - (2._wp*zdzk0-zh*zdza)/(2._wp*zk0-zh*zap1)) 1434 zdP = 0.5_wp * ( zdlogBp - zlogBp*zdzb/zb ) / zb + zdAp*zatanCp + zAp*zdCp/(1._wp+zCp*zCp) 1435 ! 1436 palphaPE(ji,jj,jk) = - grav * ( zP*zdrhop + zrhop*zdP ) * tmask(ji,jj,jk) 1437 ! 1438 ! dPEds 1439 zdzb = ze 1440 zdza = -3.0644505e-2_wp*zsr+zc 1441 zdzk0 = 1.5_wp*zb1*zsr+za1 1442 zdrhop = 9.6628e-4_wp*zs+1.5_wp*zr3*zsr+zr2 1443 zdM = (zdzb*zh-zdza)*zh+zdzk0 1444 zdDenom = zh * zdM / (zM*zM) 1445 ! 1446 zddelta = 4._wp * (zk0*zdzb+zb*zdzk0) - 2._wp*zap1*zdza 1447 zdAp = zAp * (zdza/zap1 - zdzb/zb - 0.5_wp*zddelta/zdelta) 1448 zdlogBp = zdM/(zM-zh) - zdzk0/zk0 1449 zdCp = zCp * (0.5_wp*zddelta/zdelta - (2._wp*zdzk0-zh*zdza)/(2._wp*zk0-zh*zap1)) 1450 zdP = 0.5_wp * ( zdlogBp - zlogBp*zdzb/zb ) / zb + zdAp*zatanCp + zAp*zdCp/(1._wp+zCp*zCp) 1451 ! 1452 pbetaPE(ji,jj,jk) = - grav * ( zP*zdrhop + zrhop*zdP ) * tmask(ji,jj,jk) 1453 ! 1454 END DO 1455 END DO 1456 END DO 1457 ! 1458 CASE ( 0 ) ! modified Jackett and McDougall (1995) formulation 1459 ! 1460 zws(:,:,:) = SQRT( ABS( pts(:,:,:,jp_sal) ) ) 1461 DO jk = 1, jpkm1 1462 DO jj = 1, jpj 1463 DO ji = 1, jpi 1464 zt = pts (ji,jj,jk,jp_tem) 1465 zs = pts (ji,jj,jk,jp_sal) 1466 zh = fsdept(ji,jj,jk) ! depth 1467 zsr= zws (ji,jj,jk) ! square root salinity 1468 ! 1469 ! compute volumic mass pure water at atm pressure 1470 zr1= ( ( ( ( 6.536332e-9_wp *zt - 1.120083e-6_wp )*zt + 1.001685e-4_wp )*zt & 1471 & -9.095290e-3_wp )*zt + 6.793952e-2_wp )*zt + 999.842594_wp 1472 ! seawater volumic mass atm pressure 1473 zr2= ( ( ( 5.3875e-9_wp*zt-8.2467e-7_wp ) *zt+7.6438e-5_wp ) *zt & 1474 & -4.0899e-3_wp ) *zt+0.824493_wp 1475 zr3= ( -1.6546e-6_wp*zt+1.0227e-4_wp ) *zt-5.72466e-3_wp 1476 zr4= 4.8314e-4_wp 1477 ! 1478 ! potential volumic mass (reference to the surface) 1479 zrhop= ( zr4*zs + zr3*zsr + zr2 ) *zs + zr1 1480 ! 1481 ! add the compression terms 1482 ze = ( -3.508914e-8_wp*zt-1.248266e-8_wp ) *zt-2.595994e-6_wp 1483 zbw= ( 1.296821e-6_wp*zt-5.782165e-9_wp ) *zt+1.045941e-4_wp 1484 zb = zbw + ze * zs 1485 ! 1486 zd = -2.042967e-2_wp 1487 zc = (-7.267926e-5_wp*zt+2.598241e-3_wp ) *zt+0.1571896_wp 1488 zaw= ( ( 5.939910e-6_wp*zt+2.512549e-3_wp ) *zt-0.1028859_wp ) *zt - 4.721788_wp 1489 za = ( zd*zsr + zc ) *zs + zaw 1490 ! 1491 zb1= (-0.1909078_wp*zt+7.390729_wp ) *zt-55.87545_wp 1492 za1= ( ( 2.326469e-3_wp*zt+1.553190_wp) *zt-65.00517_wp ) *zt+1044.077_wp 1493 zkw= ( ( (-1.361629e-4_wp*zt-1.852732e-2_wp ) *zt-30.41638_wp ) *zt + 2098.925_wp ) *zt+190925.6_wp 1494 zk0= ( zb1*zsr + za1 )*zs + zkw 1495 ! 1496 zM= ( zb*zh - za )*zh + zk0 1497 zDenom= 1._wp - zh / zM 1498 ! compute in-situ density 1499 zrho = zrhop/zDenom 1500 ! 1501 zap1 = 1._wp + za 1502 zdelta = 4._wp*zb*zk0 - zap1*zap1 1503 zsdelta = sqrt( abs( zdelta ) ) 1504 zAp = zap1 / zb / zsdelta 1505 zBp = ( zM - zh ) / zk0 1506 zlogBp = log(abs(zBp)) 1507 zCp = zh * zsdelta / (2._wp*zk0-zh*zap1) 1508 zatanCp = atan(zCp) 1509 zP = zh + zAp * zatanCp + 0.5_wp*zlogBp/zb 1510 ! compute potential energy 1511 pPEanom(ji,jj,jk) = - grav * zrhop * zP * tmask(ji,jj,jk) 1512 ! 1513 ! dPEdt 1514 zdzb = (2.593642e-6_wp*zt-5.782165e-9_wp) + (-7.017828e-8_wp*zt-1.248266e-8_wp) * zs 1515 zdza = (-14.535852e-5_wp*zt+2.598241e-3)*zs+((17.81973e-6_wp*zt+5.025098e-3_wp)*zt-0.1028859_wp) 1516 zdzk0 = ((-0.3818156_wp*zt+7.390729_wp)*zsr+((6.979407e-3_wp*zt+3.10638_wp)*zt-65.00517_wp))*zs & 1517 & +(((-5.446516e-4_wp*zt-5.558196e-2_wp)*zt-60.83276_wp)*zt+2098.925_wp) 1518 zdrhop = ((-3.3092e-6_wp*zt+1.0227e-4_wp)*zsr & 1519 & +(((21.55e-9_wp*zt-24.7401e-7_wp)*zt+15.2876e-5_wp)*zt-4.0899e-3_wp))*zs & 1520 & +((((32.68166e-9_wp*zt-4.480332e-6_wp)*zt+3.005055e-4_wp)*zt-18.19058e-3_wp)*zt+6.793952e-2_wp) 1521 zdM = (zdzb*zh-zdza)*zh+zdzk0 1522 zdDenom = zh * zdM / (zM*zM) 1523 ! 1524 zddelta = 4._wp * (zk0*zdzb+zb*zdzk0) - 2._wp*zap1*zdza 1525 zdAp = zAp * (zdza/zap1 - zdzb/zb - 0.5_wp*zddelta/zdelta) 1526 zdlogBp = zdM/(zM-zh) - zdzk0/zk0 1527 zdCp = zCp * (0.5_wp*zddelta/zdelta - (2._wp*zdzk0-zh*zdza)/(2._wp*zk0-zh*zap1)) 1528 zdP = 0.5_wp * ( zdlogBp - zlogBp*zdzb/zb ) / zb + zdAp*zatanCp + zAp*zdCp/(1._wp+zCp*zCp) 1529 ! 1530 palphaPE(ji,jj,jk) = - grav * ( zP*zdrhop + zrhop*zdP ) * tmask(ji,jj,jk) 1531 ! 1532 ! dPEds 1533 zdzb = ze 1534 zdza = -3.0644505e-2_wp*zsr+zc 1535 zdzk0 = 1.5_wp*zb1*zsr+za1 1536 zdrhop = 9.6628e-4_wp*zs+1.5_wp*zr3*zsr+zr2 1537 zdM = (zdzb*zh-zdza)*zh+zdzk0 1538 zdDenom = zh * zdM / (zM*zM) 1539 ! 1540 zddelta = 4._wp * (zk0*zdzb+zb*zdzk0) - 2._wp*zap1*zdza 1541 zdAp = zAp * (zdza/zap1 - zdzb/zb - 0.5_wp*zddelta/zdelta) 1542 zdlogBp = zdM/(zM-zh) - zdzk0/zk0 1543 zdCp = zCp * (0.5_wp*zddelta/zdelta - (2._wp*zdzk0-zh*zdza)/(2._wp*zk0-zh*zap1)) 1544 zdP = 0.5_wp * ( zdlogBp - zlogBp*zdzb/zb ) / zb + zdAp*zatanCp + zAp*zdCp/(1._wp+zCp*zCp) 1545 ! 1546 pbetaPE(ji,jj,jk) = - grav * ( zP*zdrhop + zrhop*zdP ) * tmask(ji,jj,jk) 1547 ! 1548 END DO 1549 END DO 1550 END DO 1551 ! 1552 CASE ( 1 ) !== Linear formulation = F( temperature ) ==! 1553 palphaPE(:,:,:) = rn_alpha * tmask(:,:,:) 1554 pbetaPE (:,:,:) = 0._wp 1555 DO jk = 1, jpkm1 1556 pPEanom(:,:,jk) = ( 0.0285_wp - rn_alpha * pts(:,:,jk,jp_tem) ) * tmask(:,:,jk) 1557 END DO 1558 ! 1559 CASE ( 2 ) !== Linear formulation = F( temperature , salinity ) ==! 1560 palphaPE(:,:,:) = rn_alpha * tmask(:,:,:) 1561 pbetaPE (:,:,:) = rn_beta * tmask(:,:,:) 1562 DO jk = 1, jpkm1 1563 pPEanom(:,:,jk) = ( rn_beta * pts(:,:,jk,jp_sal) - rn_alpha * pts(:,:,jk,jp_tem) ) * tmask(:,:,jk) 1564 END DO 1565 ! 1566 CASE( 3 ) !== Vallis (2006) simplified EOS ==! 1567 DO jk = 1, jpkm1 1568 DO jj = 1, jpj 1569 DO ji = 1, jpi 1570 zt = pts(ji,jj,jk,jp_tem) - 10._wp ! potential temperature (t-10) 1571 zs = pts(ji,jj,jk,jp_sal) - 35._wp ! salinity anomaly (s-35) 1572 zh = fsdept(ji,jj,jk) ! depth in meters at w-point 1573 ! 1574 palphaPE(ji,jj,jk) = vallis_ratio * & 1575 & ( 1.67e-4_wp * ( 1._wp + 0.55e-4_wp*zh ) + 1.e-5_wp*zt ) * tmask(ji,jj,jk) ! alphaPE 1576 pPEanom(ji,jj,jk) = vallis_diff + vallis_ratio * & 1577 & ( - ( 1.67e-4_wp*(1._wp+0.55e-4_wp*zh) + 0.5e-5_wp*zt )*zt + 0.78e-3_wp*zs + 2.195e-6_wp*zh ) & 1578 & * tmask(ji,jj,jk) 1579 ! 1580 END DO 1581 END DO 1582 END DO 1583 pbetaPE (:,:,:) = vallis_ratio * 0.78e-3_wp * tmask(:,:,:) ! betaPE=beta 1584 ! 1585 CASE DEFAULT 1586 IF(lwp) WRITE(numout,cform_err) 1587 IF(lwp) WRITE(numout,*) ' bad flag value for nn_eos = ', nn_eos 1588 nstop = nstop + 1 1589 ! 1590 END SELECT 1591 ! 1592 CALL wrk_dealloc( jpi, jpj, jpk, zws ) 1593 ! 1594 IF( nn_timing == 1 ) CALL timing_stop('eos_pen') 1595 ! 1596 END SUBROUTINE eos_pen 1597 1598 1599 677 1600 FUNCTION tfreez( psal ) RESULT( ptf ) 678 1601 !!---------------------------------------------------------------------- 679 !! *** ROUTINE eos_init***1602 !! *** ROUTINE tfreez *** 680 1603 !! 681 1604 !! ** Purpose : Compute the sea surface freezing temperature [Celcius] … … 724 1647 SELECT CASE( nn_eos ) ! check option 725 1648 ! 726 CASE( 0 ) !== Jackett and McDougall (1994) formulation ==!1649 CASE( -1 ) !== Jackett and McDougall (1995) formulation ==! 727 1650 IF(lwp) WRITE(numout,*) 728 IF(lwp) WRITE(numout,*) ' use of Jackett & McDougall (1994) equation of state and' 729 IF(lwp) WRITE(numout,*) ' McDougall (1987) Brunt-Vaisala frequency' 1651 IF(lwp) WRITE(numout,*) ' use of Jackett & McDougall (1995) equation of state' 1652 ! 1653 CASE( 0 ) !== modified Jackett and McDougall (1995) formulation ==! 1654 IF(lwp) WRITE(numout,*) 1655 IF(lwp) WRITE(numout,*) ' use of modified Jackett & McDougall (1995) equation of state' 1656 IF(lwp) WRITE(numout,*) ' and McDougall (1987) Brunt-Vaisala frequency' 730 1657 ! 731 1658 CASE( 1 ) !== Linear formulation = F( temperature ) ==! … … 736 1663 ! 737 1664 CASE( 2 ) !== Linear formulation = F( temperature , salinity ) ==! 738 ralpbet = rn_alpha / rn_beta739 1665 IF(lwp) WRITE(numout,*) 740 1666 IF(lwp) WRITE(numout,*) ' use of linear eos rho(T,S) = rau0 * ( rn_beta * S - rn_alpha * T )' 1667 ! 1668 CASE( 3 ) !== Vallis (2006) formulation ==! 1669 IF(lwp) WRITE(numout,*) 1670 IF(lwp) WRITE(numout,*) ' use of Vallis (2006) equation of state' 1671 vallis_ratio = 1027._wp / rau0 1672 vallis_diff = (1027._wp-rau0) / rau0 741 1673 ! 742 1674 CASE DEFAULT !== ERROR in nn_eos ==! -
branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_cen2.F90
r3718 r3876 4 4 !! Ocean tracers: horizontal & vertical advective trend 5 5 !!====================================================================== 6 !! History : 8.2 ! 2001-08 (G. Madec, E. Durand)trahad+trazad=traadv7 !! 8 !! 9.0! 2004-08 (C. Talandier) New trends organization6 !! History : OPA ! 2001-08 (G. Madec, E. Durand) v8.2 trahad+trazad=traadv 7 !! NEMO 1.0 ! 2002-06 (G. Madec) F90: Free form and module 8 !! - ! 2004-08 (C. Talandier) New trends organization 9 9 !! - ! 2005-11 (V. Garnier) Surface pressure gradient organization 10 10 !! 2.0 ! 2006-04 (R. Benshila, G. Madec) Step reorganization … … 21 21 USE dom_oce ! ocean space and time domain 22 22 USE eosbn2 ! equation of state 23 USE trd mod_oce ! tracers trends24 USE trdtra ! tr acers trends23 USE trd_oce ! trends: ocean variables 24 USE trdtra ! trends manager: tracers 25 25 USE closea ! closed sea 26 26 USE sbcrnf ! river runoffs … … 37 37 PRIVATE 38 38 39 PUBLIC tra_adv_cen2 40 PUBLIC ups_orca_set 39 PUBLIC tra_adv_cen2 ! routine called by step.F90 40 PUBLIC ups_orca_set ! routine used by traadv_cen2_jki.F90 41 41 42 42 LOGICAL :: l_trd ! flag to compute trends … … 55 55 56 56 SUBROUTINE tra_adv_cen2( kt, kit000, cdtype, pun, pvn, pwn, & 57 & ptb, ptn, pta, kjpt )57 & ptb, ptn, pta, kjpt ) 58 58 !!---------------------------------------------------------------------- 59 59 !! *** ROUTINE tra_adv_cen2 *** … … 85 85 !! * Add this trend now to the general trend of tracer (ta,sa): 86 86 !! pta = pta + ztra 87 !! * trend diagnostic ( 'key_trdtra' defined): the trend is87 !! * trend diagnostic (l_trdtra=T or l_trctra=T): the trend is 88 88 !! saved for diagnostics. The trends saved is expressed as 89 !! Uh.gradh(T), i.e. 90 !! save trend = ztra + ptn divn 89 !! Uh.gradh(T), i.e. save trend = ztra + ptn divn 91 90 !! 92 91 !! Part II : vertical advection … … 104 103 !! Add this trend now to the general trend of tracer (ta,sa): 105 104 !! pta = pta + ztra 106 !! Trend diagnostic ( 'key_trdtra' defined): the trend is105 !! Trend diagnostic (l_trdtra=T or l_trctra=T): the trend is 107 106 !! saved for diagnostics. The trends saved is expressed as : 108 107 !! save trend = w.gradz(T) = ztra - ptn divn. … … 144 143 IF(lwp) WRITE(numout,*) 145 144 ! 146 IF 145 IF( .NOT. ALLOCATED( upsmsk ) ) THEN 147 146 ALLOCATE( upsmsk(jpi,jpj), STAT=ierr ) 148 147 IF( ierr /= 0 ) CALL ctl_stop('STOP', 'tra_adv_cen2: unable to allocate array') … … 260 259 END DO 261 260 262 ! ! trend diagnostics (contribution of upstream fluxes)261 ! ! trend diagnostics 263 262 IF( l_trd ) THEN 264 CALL trd_tra( kt, cdtype, jn, jptra_ trd_xad, zwx, pun, ptn(:,:,:,jn) )265 CALL trd_tra( kt, cdtype, jn, jptra_ trd_yad, zwy, pvn, ptn(:,:,:,jn) )266 CALL trd_tra( kt, cdtype, jn, jptra_ trd_zad, zwz, pwn, ptn(:,:,:,jn) )263 CALL trd_tra( kt, cdtype, jn, jptra_xad, zwx, pun, ptn(:,:,:,jn) ) 264 CALL trd_tra( kt, cdtype, jn, jptra_yad, zwy, pvn, ptn(:,:,:,jn) ) 265 CALL trd_tra( kt, cdtype, jn, jptra_zad, zwz, pwn, ptn(:,:,:,jn) ) 267 266 END IF 268 267 ! ! "Poleward" heat and salt transports (contribution of upstream fluxes) … … 272 271 ENDIF 273 272 ! 274 END DO273 END DO 275 274 276 275 ! --------------------------- required in restart file to ensure restartability) -
branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_muscl.F90
r3718 r3876 16 16 !!---------------------------------------------------------------------- 17 17 USE oce ! ocean dynamics and active tracers 18 USE trc_oce ! share passive tracers/Ocean variables 18 19 USE dom_oce ! ocean space and time domain 19 USE trdmod_oce ! tracers trends 20 USE trdtra ! tracers trends 21 USE in_out_manager ! I/O manager 20 USE trd_oce ! trends: ocean variables 21 USE trdtra ! tracers trends manager 22 22 USE dynspg_oce ! choice/control of key cpp for surface pressure gradient 23 USE trabbl ! tracers: bottom boundary layer 24 USE lib_mpp ! distribued memory computing 25 USE lbclnk ! ocean lateral boundary condition (or mpp link) 23 USE sbcrnf ! river runoffs 26 24 USE diaptr ! poleward transport diagnostics 27 USE trc_oce ! share passive tracers/Ocean variables25 ! 28 26 USE wrk_nemo ! Memory Allocation 29 27 USE timing ! Timing 30 28 USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 31 USE eosbn2 ! equation of state 32 USE sbcrnf ! river runoffs 29 USE in_out_manager ! I/O manager 30 USE lib_mpp ! distribued memory computing 31 USE lbclnk ! ocean lateral boundary condition (or mpp link) 32 !!gm USE trabbl ! tracers: bottom boundary layer 33 !!gm USE eosbn2 ! equation of state 33 34 34 35 IMPLICIT NONE … … 191 192 zalpha = 0.5 - z0u 192 193 zu = z0u - 0.5 * pun(ji,jj,jk) * zdt / ( e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) ) 193 zzwx = ptb(ji+1,jj,jk,jn) + xind(ji,jj,jk) * (zu * zslpx(ji+1,jj,jk))194 zzwy = ptb(ji ,jj,jk,jn) + xind(ji,jj,jk) * (zu * zslpx(ji ,jj,jk))194 zzwx = ptb(ji+1,jj,jk,jn) + xind(ji,jj,jk) * zu * zslpx(ji+1,jj,jk) 195 zzwy = ptb(ji ,jj,jk,jn) + xind(ji,jj,jk) * zu * zslpx(ji ,jj,jk) 195 196 zwx(ji,jj,jk) = pun(ji,jj,jk) * ( zalpha * zzwx + (1.-zalpha) * zzwy ) 196 197 ! … … 198 199 zalpha = 0.5 - z0v 199 200 zv = z0v - 0.5 * pvn(ji,jj,jk) * zdt / ( e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) ) 200 zzwx = ptb(ji,jj+1,jk,jn) + xind(ji,jj,jk) * (zv * zslpy(ji,jj+1,jk))201 zzwy = ptb(ji,jj ,jk,jn) + xind(ji,jj,jk) * (zv * zslpy(ji,jj ,jk))201 zzwx = ptb(ji,jj+1,jk,jn) + xind(ji,jj,jk) * zv * zslpy(ji,jj+1,jk) 202 zzwy = ptb(ji,jj ,jk,jn) + xind(ji,jj,jk) * zv * zslpy(ji,jj ,jk) 202 203 zwy(ji,jj,jk) = pvn(ji,jj,jk) * ( zalpha * zzwx + (1.-zalpha) * zzwy ) 203 204 END DO … … 222 223 ! ! trend diagnostics (contribution of upstream fluxes) 223 224 IF( l_trd ) THEN 224 CALL trd_tra( kt, cdtype, jn, jptra_ trd_xad, zwx, pun, ptb(:,:,:,jn) )225 CALL trd_tra( kt, cdtype, jn, jptra_ trd_yad, zwy, pvn, ptb(:,:,:,jn) )225 CALL trd_tra( kt, cdtype, jn, jptra_xad, zwx, pun, ptb(:,:,:,jn) ) 226 CALL trd_tra( kt, cdtype, jn, jptra_yad, zwy, pvn, ptb(:,:,:,jn) ) 226 227 END IF 227 228 ! ! "Poleward" heat and salt transports (contribution of upstream fluxes) … … 273 274 zalpha = 0.5 + z0w 274 275 zw = z0w - 0.5 * pwn(ji,jj,jk+1) * zdt * zbtr 275 zzwx = ptb(ji,jj,jk+1,jn) + xind(ji,jj,jk) * (zw * zslpx(ji,jj,jk+1))276 zzwy = ptb(ji,jj,jk ,jn) + xind(ji,jj,jk) * (zw * zslpx(ji,jj,jk ))276 zzwx = ptb(ji,jj,jk+1,jn) + xind(ji,jj,jk) * zw * zslpx(ji,jj,jk+1) 277 zzwy = ptb(ji,jj,jk ,jn) + xind(ji,jj,jk) * zw * zslpx(ji,jj,jk ) 277 278 zwx(ji,jj,jk+1) = pwn(ji,jj,jk+1) * ( zalpha * zzwx + (1.-zalpha) * zzwy ) 278 279 END DO … … 293 294 END DO 294 295 ! ! Save the vertical advective trends for diagnostic 295 IF( l_trd ) CALL trd_tra( kt, cdtype, jn, jptra_ trd_zad, zwx, pwn, ptb(:,:,:,jn) )296 IF( l_trd ) CALL trd_tra( kt, cdtype, jn, jptra_zad, zwx, pwn, ptb(:,:,:,jn) ) 296 297 ! 297 298 ENDDO -
branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_muscl2.F90
r3625 r3876 13 13 !!---------------------------------------------------------------------- 14 14 USE oce ! ocean dynamics and active tracers 15 USE trc_oce ! share passive tracers/Ocean variables 15 16 USE dom_oce ! ocean space and time domain 16 USE trd mod_oce ! tracers trends17 USE trdtra ! tr acers trends17 USE trd_oce ! trends: ocean variables 18 USE trdtra ! trends manager: tracers 18 19 USE in_out_manager ! I/O manager 19 20 USE dynspg_oce ! choice/control of key cpp for surface pressure gradient 20 USE trabbl ! tracers: bottom boundary layer 21 !!gm USE trabbl ! tracers: bottom boundary layer 22 USE diaptr ! poleward transport diagnostics 23 ! 21 24 USE lib_mpp ! distribued memory computing 22 25 USE lbclnk ! ocean lateral boundary condition (or mpp link) 23 USE diaptr ! poleward transport diagnostics24 USE trc_oce ! share passive tracers/Ocean variables25 26 USE wrk_nemo ! Memory Allocation 26 27 USE timing ! Timing … … 201 202 ! ! trend diagnostics (contribution of upstream fluxes) 202 203 IF( l_trd ) THEN 203 CALL trd_tra( kt, cdtype, jn, jptra_ trd_xad, zwx, pun, ptb(:,:,:,jn) )204 CALL trd_tra( kt, cdtype, jn, jptra_ trd_yad, zwy, pvn, ptb(:,:,:,jn) )204 CALL trd_tra( kt, cdtype, jn, jptra_xad, zwx, pun, ptb(:,:,:,jn) ) 205 CALL trd_tra( kt, cdtype, jn, jptra_yad, zwy, pvn, ptb(:,:,:,jn) ) 205 206 END IF 206 207 … … 284 285 END DO 285 286 ! ! trend diagnostics (contribution of upstream fluxes) 286 IF( l_trd ) CALL trd_tra( kt, cdtype, jn, jptra_ trd_zad, zwx, pwn, ptb(:,:,:,jn) )287 IF( l_trd ) CALL trd_tra( kt, cdtype, jn, jptra_zad, zwx, pwn, ptb(:,:,:,jn) ) 287 288 ! 288 289 END DO -
branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_qck.F90
r3625 r3876 17 17 USE oce ! ocean dynamics and active tracers 18 18 USE dom_oce ! ocean space and time domain 19 USE trdmod_oce ! ocean space and time domain 20 USE trdtra ! ocean tracers trends 21 USE trabbl ! advective term in the BBL 19 USE trc_oce ! share passive tracers/Ocean variables 20 USE trd_oce ! trends: ocean variables 21 USE trdtra ! trends manager: tracers 22 !!gm USE trabbl ! advective term in the BBL 23 USE dynspg_oce ! surface pressure gradient variables 24 USE diaptr ! poleward transport diagnostics 25 ! 22 26 USE lib_mpp ! distribued memory computing 23 27 USE lbclnk ! ocean lateral boundary condition (or mpp link) 24 USE dynspg_oce ! surface pressure gradient variables25 28 USE in_out_manager ! I/O manager 26 USE diaptr ! poleward transport diagnostics27 USE trc_oce ! share passive tracers/Ocean variables28 29 USE wrk_nemo ! Memory Allocation 29 30 USE timing ! Timing … … 233 234 END DO 234 235 ! ! trend diagnostics (contribution of upstream fluxes) 235 IF( l_trd ) CALL trd_tra( kt, cdtype, jn, jptra_ trd_xad, zwx, pun, ptn(:,:,:,jn) )236 IF( l_trd ) CALL trd_tra( kt, cdtype, jn, jptra_xad, zwx, pun, ptn(:,:,:,jn) ) 236 237 ! 237 238 END DO … … 359 360 END DO 360 361 ! ! trend diagnostics (contribution of upstream fluxes) 361 IF( l_trd ) CALL trd_tra( kt, cdtype, jn, jptra_ trd_yad, zwy, pvn, ptn(:,:,:,jn) )362 IF( l_trd ) CALL trd_tra( kt, cdtype, jn, jptra_yad, zwy, pvn, ptn(:,:,:,jn) ) 362 363 ! ! "Poleward" heat and salt transports (contribution of upstream fluxes) 363 364 IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nn_fptr ) == 0 ) ) THEN … … 422 423 END DO 423 424 ! ! Save the vertical advective trends for diagnostic 424 IF( l_trd ) CALL trd_tra( kt, cdtype, jn, jptra_ trd_zad, zwz, pwn, ptn(:,:,:,jn) )425 IF( l_trd ) CALL trd_tra( kt, cdtype, jn, jptra_zad, zwz, pwn, ptn(:,:,:,jn) ) 425 426 ! 426 427 END DO -
branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_tvd.F90
r3625 r3876 22 22 USE oce ! ocean dynamics and active tracers 23 23 USE dom_oce ! ocean space and time domain 24 USE trdmod_oce ! tracers trends 24 USE trc_oce ! share passive tracers/Ocean variables 25 USE trd_oce ! trends: ocean variables 25 26 USE trdtra ! tracers trends 26 USE in_out_manager ! I/O manager27 27 USE dynspg_oce ! choice/control of key cpp for surface pressure gradient 28 USE diaptr ! poleward transport diagnostics 29 ! 28 30 USE lib_mpp ! MPP library 29 31 USE lbclnk ! ocean lateral boundary condition (or mpp link) 30 USE diaptr ! poleward transport diagnostics 31 USE trc_oce ! share passive tracers/Ocean variables 32 USE in_out_manager ! I/O manager 32 33 USE wrk_nemo ! Memory Allocation 33 34 USE timing ! Timing … … 228 229 ztrdz(:,:,:) = ztrdz(:,:,:) + zwz(:,:,:) ! <<< Add to previously computed 229 230 230 CALL trd_tra( kt, cdtype, jn, jptra_ trd_xad, ztrdx, pun, ptn(:,:,:,jn) )231 CALL trd_tra( kt, cdtype, jn, jptra_ trd_yad, ztrdy, pvn, ptn(:,:,:,jn) )232 CALL trd_tra( kt, cdtype, jn, jptra_ trd_zad, ztrdz, pwn, ptn(:,:,:,jn) )231 CALL trd_tra( kt, cdtype, jn, jptra_xad, ztrdx, pun, ptn(:,:,:,jn) ) 232 CALL trd_tra( kt, cdtype, jn, jptra_yad, ztrdy, pvn, ptn(:,:,:,jn) ) 233 CALL trd_tra( kt, cdtype, jn, jptra_zad, ztrdz, pwn, ptn(:,:,:,jn) ) 233 234 END IF 234 235 ! ! "Poleward" heat and salt transports (contribution of upstream fluxes) -
branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_ubs.F90
r3787 r3876 14 14 USE oce ! ocean dynamics and active tracers 15 15 USE dom_oce ! ocean space and time domain 16 USE trdmod_oce ! ocean space and time domain 17 USE trdtra 18 USE lib_mpp 16 USE trc_oce ! share passive tracers/Ocean variables 17 USE trd_oce ! trends: ocean variables 18 USE trdtra ! trends manager: tracers 19 USE dynspg_oce ! choice/control of key cpp for surface pressure gradient 20 USE diaptr ! poleward transport diagnostics 21 ! 22 USE lib_mpp ! I/O library 19 23 USE lbclnk ! ocean lateral boundary condition (or mpp link) 20 24 USE in_out_manager ! I/O manager 21 USE diaptr ! poleward transport diagnostics22 USE dynspg_oce ! choice/control of key cpp for surface pressure gradient23 USE trc_oce ! share passive tracers/Ocean variables24 25 USE wrk_nemo ! Memory Allocation 25 26 USE timing ! Timing … … 51 52 !! and add it to the general trend of passive tracer equations. 52 53 !! 53 !! ** Method : The upstream biased 3rd order scheme (UBS) is based on an54 !! ** Method : The upstream biased scheme (UBS) is based on a 3rd order 54 55 !! upstream-biased parabolic interpolation (Shchepetkin and McWilliams 2005) 55 56 !! It is only used in the horizontal direction. … … 182 183 ! ! trend diagnostics (contribution of upstream fluxes) 183 184 IF( l_trd ) THEN 184 CALL trd_tra( kt, cdtype, jn, jptra_ trd_xad, zwx, pun, ptn(:,:,:,jn) )185 CALL trd_tra( kt, cdtype, jn, jptra_ trd_yad, zwy, pvn, ptn(:,:,:,jn) )185 CALL trd_tra( kt, cdtype, jn, jptra_xad, zwx, pun, ptn(:,:,:,jn) ) 186 CALL trd_tra( kt, cdtype, jn, jptra_yad, zwy, pvn, ptn(:,:,:,jn) ) 186 187 END IF 187 188 ! ! "Poleward" heat and salt transports (contribution of upstream fluxes) … … 265 266 END DO 266 267 END DO 267 CALL trd_tra( kt, cdtype, jn, jptra_ trd_zad, zltv )268 CALL trd_tra( kt, cdtype, jn, jptra_zad, zltv ) 268 269 ENDIF 269 270 ! -
branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/TRA/trabbc.F90
r3625 r3876 18 18 USE dom_oce ! domain: ocean 19 19 USE phycst ! physical constants 20 USE trd mod_oce ! trends: ocean variables21 USE trdtra ! trends : activetracers20 USE trd_oce ! trends: ocean variables 21 USE trdtra ! trends manager: tracers 22 22 USE in_out_manager ! I/O manager 23 23 USE prtctl ! Print control … … 99 99 IF( l_trdtra ) THEN ! Save the geothermal heat flux trend for diagnostics 100 100 ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:) 101 CALL trd_tra( kt, 'TRA', jp_tem, jptra_ trd_bbc, ztrdt )101 CALL trd_tra( kt, 'TRA', jp_tem, jptra_bbc, ztrdt ) 102 102 CALL wrk_dealloc( jpi, jpj, jpk, ztrdt ) 103 103 ENDIF -
branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/TRA/trabbl.F90
r3764 r3876 28 28 USE phycst ! physical constant 29 29 USE eosbn2 ! equation of state 30 USE trd mod_oce! trends: ocean variables31 USE trdtra ! trends : active tracers32 USE iom ! IOM server 30 USE trd_oce ! trends: ocean variables 31 USE trdtra ! trends manager: tracers 32 USE iom ! IOM server 33 33 USE in_out_manager ! I/O manager 34 34 USE lbclnk ! ocean lateral boundary conditions … … 148 148 ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:) 149 149 ztrds(:,:,:) = tsa(:,:,:,jp_sal) - ztrds(:,:,:) 150 CALL trd_tra( kt, 'TRA', jp_tem, jptra_ trd_bbl, ztrdt )151 CALL trd_tra( kt, 'TRA', jp_sal, jptra_ trd_bbl, ztrds )152 CALL wrk_dealloc( jpi, jpj, jpk, ztrdt, ztrds ) 150 CALL trd_tra( kt, 'TRA', jp_tem, jptra_bbl, ztrdt ) 151 CALL trd_tra( kt, 'TRA', jp_sal, jptra_bbl, ztrds ) 152 CALL wrk_dealloc( jpi, jpj, jpk, ztrdt, ztrds ) 153 153 ENDIF 154 154 ! … … 180 180 !! Campin, J.-M., and H. Goosse, 1999, Tellus, 412-430. 181 181 !!---------------------------------------------------------------------- 182 !183 182 INTEGER , INTENT(in ) :: kjpt ! number of tracers 184 183 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: ptb ! before and now tracer fields … … 399 398 - 0.121555e-07 ) * pfh 400 399 !!---------------------------------------------------------------------- 401 402 400 ! 403 401 IF( nn_timing == 1 ) CALL timing_start( 'bbl') … … 405 403 CALL wrk_alloc( jpi, jpj, zub, zvb, ztb, zsb, zdep ) 406 404 ! 407 408 405 IF( kt == kit000 ) THEN 409 406 IF(lwp) WRITE(numout,*) … … 411 408 IF(lwp) WRITE(numout,*) '~~~~~~~~~~' 412 409 ENDIF 413 414 410 ! !* bottom temperature, salinity, velocity and depth 415 411 #if defined key_vectopt_loop -
branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/TRA/tradmp.F90
r3294 r3876 27 27 USE oce ! ocean: variables 28 28 USE dom_oce ! ocean: domain variables 29 USE trd mod_oce ! ocean: trendvariables30 USE trdtra ! active tracers: trends29 USE trd_oce ! trends: ocean variables 30 USE trdtra ! trends manager: tracers 31 31 USE zdf_oce ! ocean: vertical physics 32 32 USE phycst ! physical constants … … 47 47 PUBLIC dtacof_zoom ! routine called by in both tradmp.F90 and trcdmp.F90 48 48 49 ! !!* Namelist namtra_dmp : T & S newtonian damping *49 ! !!* Namelist namtra_dmp : T & S newtonian damping * 50 50 LOGICAL, PUBLIC :: ln_tradmp = .TRUE. !: internal damping flag 51 51 INTEGER :: nn_hdmp = -1 ! = 0/-1/'latitude' for damping over T and S … … 111 111 ! 112 112 CALL wrk_alloc( jpi, jpj, jpk, jpts, zts_dta ) 113 ! 113 114 ! !== input T-S data at kt ==! 114 115 CALL dta_tsd( kt, zts_dta ) ! read and interpolates T-S data at kt … … 171 172 ! 172 173 IF( l_trdtra ) THEN ! trend diagnostic 173 CALL trd_tra( kt, 'TRA', jp_tem, jptra_ trd_dmp, ttrdmp )174 CALL trd_tra( kt, 'TRA', jp_sal, jptra_ trd_dmp, strdmp )174 CALL trd_tra( kt, 'TRA', jp_tem, jptra_dmp, ttrdmp ) 175 CALL trd_tra( kt, 'TRA', jp_sal, jptra_dmp, strdmp ) 175 176 ENDIF 176 177 ! ! Control print -
branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/TRA/traldf.F90
r3294 r3876 23 23 USE traldf_iso_grif ! lateral mixing (tra_ldf_iso_grif routine) 24 24 USE traldf_lap ! lateral mixing (tra_ldf_lap routine) 25 USE trd mod_oce ! ocean space and time domain26 USE trdtra ! ocean active tracers trends25 USE trd_oce ! trends: ocean variables 26 USE trdtra ! trends manager: tracers 27 27 USE prtctl ! Print control 28 28 USE in_out_manager ! I/O manager … … 35 35 PRIVATE 36 36 37 PUBLIC tra_ldf 38 PUBLIC tra_ldf_init 37 PUBLIC tra_ldf ! called by step.F90 38 PUBLIC tra_ldf_init ! called by opa.F90 39 39 ! 40 40 INTEGER :: nldf = 0 ! type of lateral diffusion used defined from ln_traldf_... namlist logicals) … … 112 112 ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:) 113 113 ztrds(:,:,:) = tsa(:,:,:,jp_sal) - ztrds(:,:,:) 114 CALL trd_tra( kt, 'TRA', jp_tem, jptra_ trd_ldf, ztrdt )115 CALL trd_tra( kt, 'TRA', jp_sal, jptra_ trd_ldf, ztrds )114 CALL trd_tra( kt, 'TRA', jp_tem, jptra_ldf, ztrdt ) 115 CALL trd_tra( kt, 'TRA', jp_sal, jptra_ldf, ztrds ) 116 116 CALL wrk_dealloc( jpi, jpj, jpk, ztrdt, ztrds ) 117 117 ENDIF -
branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/TRA/tranpc.F90
r3294 r3876 17 17 USE dom_oce ! ocean space and time domain 18 18 USE zdf_oce ! ocean vertical physics 19 USE trd mod_oce ! ocean active tracer trends20 USE trdtra ! ocean active tracer trends19 USE trd_oce ! trends: ocean variables 20 USE trdtra ! trends manager: tracers 21 21 USE eosbn2 ! equation of state (eos routine) 22 22 USE lbclnk ! lateral boundary conditions (or mpp link) … … 54 54 !! 55 55 !! ** Action : - (ta,sa) after the application od the npc scheme 56 !! - save the associated trends ( ttrd,strd) ('key_trdtra')56 !! - save the associated trends (l_trdtra=T) 57 57 !! 58 58 !! References : Madec, et al., 1991, JPO, 21, 9, 1349-1371. … … 196 196 ! ! =============== 197 197 ! 198 199 ! Lateral boundary conditions on ( ta, sa ) ( Unchanged sign) 200 CALL lbc_lnk( tsa(:,:,:,jp_tem), 'T', 1. ) ; CALL lbc_lnk( tsa(:,:,:,jp_sal), 'T', 1. ) 201 198 202 IF( l_trdtra ) THEN ! save the Non penetrative mixing trends for diagnostic 199 203 ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:) 200 204 ztrds(:,:,:) = tsa(:,:,:,jp_sal) - ztrds(:,:,:) 201 CALL trd_tra( kt, 'TRA', jp_tem, jptra_ trd_npc, ztrdt )202 CALL trd_tra( kt, 'TRA', jp_sal, jptra_ trd_npc, ztrds )205 CALL trd_tra( kt, 'TRA', jp_tem, jptra_npc, ztrdt ) 206 CALL trd_tra( kt, 'TRA', jp_sal, jptra_npc, ztrds ) 203 207 CALL wrk_dealloc( jpi, jpj, jpk, ztrdt, ztrds ) 204 208 ENDIF 205 209 206 ! Lateral boundary conditions on ( ta, sa ) ( Unchanged sign) 207 ! ------------------------------============ 208 CALL lbc_lnk( tsa(:,:,:,jp_tem), 'T', 1. ) ; CALL lbc_lnk( tsa(:,:,:,jp_sal), 'T', 1. ) 209 210 211 ! 2. non penetrative convective scheme statistics 212 ! ----------------------------------------------- 210 ! non penetrative convective scheme statistics 213 211 IF( nn_npcp /= 0 .AND. MOD( kt, nn_npcp ) == 0 ) THEN 214 212 IF(lwp) WRITE(numout,*)' kt=',kt, ' number of statically instable', & -
branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/TRA/tranxt.F90
r3294 r3876 27 27 USE dom_oce ! ocean space and time domain variables 28 28 USE sbc_oce ! surface boundary condition: ocean 29 USE zdf_oce ! ???29 USE zdf_oce ! ocean vertical mixing 30 30 USE domvvl ! variable volume 31 31 USE dynspg_oce ! surface pressure gradient variables 32 32 USE dynhpg ! hydrostatic pressure gradient 33 USE trdmod_oce ! ocean space and time domain variables 34 USE trdtra ! ocean active tracers trends 35 USE phycst 36 USE obc_oce 33 USE trd_oce ! trends: ocean variables 34 USE trdtra ! trends manager: tracers 35 USE traqsr ! penetrative solar radiation (needed for nksr) 36 USE phycst ! physical constant 37 USE ldftra_oce ! lateral physics on tracers 38 USE bdy_oce ! BDY open boundary condition variables 39 USE bdytra ! BDY open boundary condition (bdy_tra routine) 40 USE obc_oce ! open boundary condition variables 37 41 USE obctra ! open boundary condition (obc_tra routine) 38 USE bdy_oce 39 USE bdytra ! open boundary condition (bdy_tra routine) 42 ! 40 43 USE in_out_manager ! I/O manager 41 44 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 42 45 USE prtctl ! Print control 43 USE traqsr ! penetrative solar radiation (needed for nksr) 46 USE wrk_nemo ! Memory allocation 47 USE timing ! Timing 44 48 #if defined key_agrif 45 49 USE agrif_opa_update 46 50 USE agrif_opa_interp 47 51 #endif 48 USE wrk_nemo ! Memory allocation49 USE timing ! Timing50 52 51 53 IMPLICIT NONE … … 132 134 ztrdt(:,:,:) = tsn(:,:,:,jp_tem) 133 135 ztrds(:,:,:) = tsn(:,:,:,jp_sal) 136 IF( ln_traldf_iso ) THEN ! diagnose the "pure" Kz diffusive trend 137 CALL trd_tra( kt, 'TRA', jp_tem, jptra_zdfp, ztrdt ) 138 CALL trd_tra( kt, 'TRA', jp_sal, jptra_zdfp, ztrds ) 139 ENDIF 134 140 ENDIF 135 141 … … 159 165 ztrds(:,:,jk) = ( tsb(:,:,jk,jp_sal) - ztrds(:,:,jk) ) * zfact 160 166 END DO 161 CALL trd_tra( kt, 'TRA', jp_tem, jptra_ trd_atf, ztrdt )162 CALL trd_tra( kt, 'TRA', jp_sal, jptra_ trd_atf, ztrds )167 CALL trd_tra( kt, 'TRA', jp_tem, jptra_atf, ztrdt ) 168 CALL trd_tra( kt, 'TRA', jp_sal, jptra_atf, ztrds ) 163 169 CALL wrk_dealloc( jpi, jpj, jpk, ztrdt, ztrds ) 164 170 END IF … … 168 174 & tab3d_2=tsn(:,:,:,jp_sal), clinfo2= ' Sn: ', mask2=tmask ) 169 175 ! 170 ! 171 IF( nn_timing == 1 ) CALL timing_stop('tra_nxt') 176 IF( nn_timing == 1 ) CALL timing_stop('tra_nxt') 172 177 ! 173 178 END SUBROUTINE tra_nxt -
branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/TRA/traqsr.F90
r3680 r3876 13 13 14 14 !!---------------------------------------------------------------------- 15 !! tra_qsr : trend due to the solar radiation penetration16 !! tra_qsr_init : solar radiation penetration initialization15 !! tra_qsr : trend due to the solar radiation penetration 16 !! tra_qsr_init : solar radiation penetration initialization 17 17 !!---------------------------------------------------------------------- 18 USE oce ! ocean dynamics and active tracers 19 USE dom_oce ! ocean space and time domain 20 USE sbc_oce ! surface boundary condition: ocean 21 USE trc_oce ! share SMS/Ocean variables 22 USE trdmod_oce ! ocean variables trends 23 USE trdtra ! ocean active tracers trends 24 USE in_out_manager ! I/O manager 25 USE phycst ! physical constants 26 USE prtctl ! Print control 27 USE iom ! I/O manager 28 USE fldread ! read input fields 29 USE lib_mpp ! MPP library 18 USE oce ! ocean dynamics and active tracers 19 USE dom_oce ! ocean space and time domain 20 USE sbc_oce ! surface boundary condition: ocean 21 USE trc_oce ! share SMS/Ocean variables 22 USE trd_oce ! trends: ocean variables 23 USE trdtra ! trends manager: tracers 24 USE phycst ! physical constants 25 ! 26 USE in_out_manager ! I/O manager 27 USE prtctl ! Print control 28 USE iom ! I/O manager 29 USE fldread ! read input fields 30 USE lib_mpp ! MPP library 30 31 USE wrk_nemo ! Memory Allocation 31 32 USE timing ! Timing 32 33 33 34 34 IMPLICIT NONE … … 48 48 REAL(wp), PUBLIC :: rn_si1 = 23.0_wp !: deepest depth of extinction (water type I) (2 bands) 49 49 50 ! Module variables 51 REAL(wp) :: xsi0r !: inverse of rn_si0 52 REAL(wp) :: xsi1r !: inverse of rn_si1 50 INTEGER , PUBLIC :: nksr !: levels below which the light cannot penetrate ( depth larger than 391 m) 51 52 REAL(wp) :: xsi0r, xsi1r ! inverse of rn_si0 and rn_si1, resp. 53 REAL(wp), DIMENSION(3,61) :: rkrgb ! tabulated attenuation coefficients for RGB absorption 53 54 TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_chl ! structure of input Chl (file informations, fields read) 54 INTEGER, PUBLIC :: nksr ! levels below which the light cannot penetrate ( depth larger than 391 m)55 REAL(wp), DIMENSION(3,61) :: rkrgb !: tabulated attenuation coefficients for RGB absorption56 55 57 56 !! * Substitutions … … 87 86 !! 88 87 !! ** Action : - update ta with the penetrative solar radiation trend 89 !! - s ave the trend in ttrd ('key_trdtra')88 !! - send the trend to trdtra (l_trdtra=T) 90 89 !! 91 90 !! Reference : Jerlov, N. G., 1968 Optical Oceanography, Elsevier, 194pp. 92 91 !! Lengaigne et al. 2007, Clim. Dyn., V28, 5, 503-516. 93 92 !!---------------------------------------------------------------------- 94 !95 93 INTEGER, INTENT(in) :: kt ! ocean time-step 96 94 ! … … 116 114 ENDIF 117 115 118 IF( l_trdtra ) THEN ! Save t a and sa trends116 IF( l_trdtra ) THEN ! Save temperature trend 119 117 CALL wrk_alloc( jpi, jpj, jpk, ztrdt ) 120 118 ztrdt(:,:,:) = tsa(:,:,:,jp_tem) … … 141 139 ! Compute now qsr tracer content field 142 140 ! ************************************ 143 144 141 ! ! ============================================== ! 145 142 IF( lk_qsr_bio .AND. ln_qsr_bio ) THEN ! bio-model fluxes : all vertical coordinates ! … … 168 165 IF( nn_chldta == 1 .OR. lk_vvl ) THEN !* Variable Chlorophyll or ocean volume 169 166 ! 170 IF( nn_chldta == 1 ) THEN ! *Variable Chlorophyll167 IF( nn_chldta == 1 ) THEN !- Variable Chlorophyll 171 168 ! 172 169 CALL fld_read( kt, 1, sf_chl ) ! Read Chl data and provides it at the current time step … … 184 181 END DO 185 182 END DO 186 ELSE !Variable ocean volume but constant chrlorophyll187 zchl = 0.05 ! constant chlorophyll183 ELSE !- Variable ocean volume but constant chrlorophyll 184 zchl = 0.05 ! constant chlorophyll 188 185 irgb = NINT( 41 + 20.*LOG10( zchl ) + 1.e-15 ) 189 zekb(:,:) = rkrgb(1,irgb) ! Separation in R-G-B depending of the chlorophyll186 zekb(:,:) = rkrgb(1,irgb) ! Separation in R-G-B depending of the chlorophyll 190 187 zekg(:,:) = rkrgb(2,irgb) 191 188 zekr(:,:) = rkrgb(3,irgb) 192 189 ENDIF 193 190 ! 194 zcoef = ( 1. - rn_abs ) / 3.e0 !equi-partition in R-G-B195 ze0(:,:,1) = rn_abs 196 ze1(:,:,1) = zcoef * qsr(:,:)197 ze2(:,:,1) = zcoef * qsr(:,:)198 ze3(:,:,1) = zcoef * qsr(:,:)199 zea(:,:,1) = qsr(:,:)191 zcoef = ( 1. - rn_abs ) / 3.e0 !- equi-partition in R-G-B 192 ze0(:,:,1) = rn_abs * qsr(:,:) 193 ze1(:,:,1) = zcoef * qsr(:,:) 194 ze2(:,:,1) = zcoef * qsr(:,:) 195 ze3(:,:,1) = zcoef * qsr(:,:) 196 zea(:,:,1) = qsr(:,:) 200 197 ! 201 198 DO jk = 2, nksr+1 … … 283 280 IF( l_trdtra ) THEN ! qsr tracers trends saved for diagnostics 284 281 ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:) 285 CALL trd_tra( kt, 'TRA', jp_tem, jptra_ trd_qsr, ztrdt )282 CALL trd_tra( kt, 'TRA', jp_tem, jptra_qsr, ztrdt ) 286 283 CALL wrk_dealloc( jpi, jpj, jpk, ztrdt ) 287 284 ENDIF -
branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/TRA/trasbc.F90
r3764 r3876 18 18 USE dom_oce ! ocean space domain variables 19 19 USE phycst ! physical constant 20 USE sbcmod ! ln_rnf 21 USE sbcrnf ! River runoff 20 22 USE traqsr ! solar radiation penetration 21 USE trdmod_oce ! ocean trends 22 USE trdtra ! ocean trends 23 USE trd_oce ! trends: ocean variables 24 USE trdtra ! trends manager: tracers 25 ! 23 26 USE in_out_manager ! I/O manager 24 27 USE prtctl ! Print control 25 USE sbcrnf ! River runoff 26 USE sbcmod ! ln_rnf 27 USE iom 28 USE iom ! I/O library 28 29 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 29 30 USE wrk_nemo ! Memory Allocation … … 91 92 !! where emp, the surface freshwater budget (evaporation minus 92 93 !! precipitation minus runoff) given in kg/m2/s is divided 93 !! by rau0 = 10 20kg/m3 (density of sea water) to obtain m/s.94 !! by rau0 = 1035 kg/m3 (density of sea water) to obtain m/s. 94 95 !! Note: even though Fwe does not appear explicitly for 95 96 !! temperature in this routine, the heat carried by the water … … 107 108 !! ** Action : - Update the 1st level of (ta,sa) with the trend associated 108 109 !! with the tracer surface boundary condition 109 !! - s ave the trend it in ttrd ('key_trdtra')110 !! - send trends to trdtra module (l_trdtra=T) 110 111 !!---------------------------------------------------------------------- 111 112 INTEGER, INTENT(in) :: kt ! ocean time-step index … … 124 125 ENDIF 125 126 126 IF( l_trdtra ) 127 IF( l_trdtra ) THEN !* Save ta and sa trends 127 128 CALL wrk_alloc( jpi, jpj, jpk, ztrdt, ztrds ) 128 129 ztrdt(:,:,:) = tsa(:,:,:,jp_tem) … … 137 138 138 139 !---------------------------------------- 139 ! EMP, EMPSand QNS effects140 ! EMP, SFX and QNS effects 140 141 !---------------------------------------- 141 142 ! Set before sbc tracer content fields … … 146 147 & iom_varid( numror, 'sbc_hc_b', ldstop = .FALSE. ) > 0 ) THEN 147 148 IF(lwp) WRITE(numout,*) ' nit000-1 surface tracer content forcing fields red in the restart file' 148 zfact = 0.5 e0149 zfact = 0.5_wp 149 150 CALL iom_get( numror, jpdom_autoglo, 'sbc_hc_b', sbc_tsc_b(:,:,jp_tem) ) ! before heat content sbc trend 150 151 CALL iom_get( numror, jpdom_autoglo, 'sbc_sc_b', sbc_tsc_b(:,:,jp_sal) ) ! before salt content sbc trend 151 152 ELSE ! No restart or restart not found: Euler forward time stepping 152 zfact = 1. e0153 sbc_tsc_b(:,:,:) = 0. e0153 zfact = 1._wp 154 sbc_tsc_b(:,:,:) = 0._wp 154 155 ENDIF 155 156 ELSE ! Swap of forcing fields 156 157 ! ! ---------------------- 157 zfact = 0.5 e0158 zfact = 0.5_wp 158 159 sbc_tsc_b(:,:,:) = sbc_tsc(:,:,:) 159 160 ENDIF … … 226 227 ENDIF 227 228 228 IF( l_trdtra ) THEN ! save the horizontal diffusivetrends for further diagnostics229 IF( l_trdtra ) THEN ! save the trends for further diagnostics 229 230 ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:) 230 231 ztrds(:,:,:) = tsa(:,:,:,jp_sal) - ztrds(:,:,:) 231 CALL trd_tra( kt, 'TRA', jp_tem, jptra_ trd_nsr, ztrdt )232 CALL trd_tra( kt, 'TRA', jp_sal, jptra_ trd_nsr, ztrds )232 CALL trd_tra( kt, 'TRA', jp_tem, jptra_nsr, ztrdt ) 233 CALL trd_tra( kt, 'TRA', jp_sal, jptra_nsr, ztrds ) 233 234 CALL wrk_dealloc( jpi, jpj, jpk, ztrdt, ztrds ) 234 235 ENDIF -
branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/TRA/trazdf.F90
r3294 r3876 19 19 USE sbc_oce ! surface boundary condition: ocean 20 20 USE dynspg_oce 21 22 21 USE trazdf_exp ! vertical diffusion: explicit (tra_zdf_exp routine) 23 22 USE trazdf_imp ! vertical diffusion: implicit (tra_zdf_imp routine) 24 25 23 USE ldftra_oce ! ocean active tracers: lateral physics 26 USE trdmod_oce ! ocean active tracers: lateral physics 27 USE trdtra ! ocean tracers trends 24 USE trd_oce ! trends: ocean variables 25 USE trdtra ! trends manager: tracers 26 ! 28 27 USE in_out_manager ! I/O manager 29 28 USE prtctl ! Print control … … 96 95 ztrds(:,:,jk) = ( ( tsa(:,:,jk,jp_sal) - tsb(:,:,jk,jp_sal) ) / r2dtra(jk) ) - ztrds(:,:,jk) 97 96 END DO 98 CALL trd_tra( kt, 'TRA', jp_tem, jptra_ trd_zdf, ztrdt )99 CALL trd_tra( kt, 'TRA', jp_sal, jptra_ trd_zdf, ztrds )97 CALL trd_tra( kt, 'TRA', jp_tem, jptra_zdf, ztrdt ) 98 CALL trd_tra( kt, 'TRA', jp_sal, jptra_zdf, ztrds ) 100 99 CALL wrk_dealloc( jpi, jpj, jpk, ztrdt, ztrds ) 101 100 ENDIF -
branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/TRD/trdtra.F90
r3632 r3876 2 2 !!====================================================================== 3 3 !! *** MODULE trdtra *** 4 !! Ocean diagnostics: ocean tracers trends 4 !! Ocean diagnostics: ocean tracers trends pre-processing 5 5 !!===================================================================== 6 !! History : 1.0 ! 2004-08 (C. Talandier) Original code 7 !! 2.0 ! 2005-04 (C. Deltel) Add Asselin trend in the ML budget 8 !! 3.3 ! 2010-06 (C. Ethe) merge TRA-TRC 9 !!---------------------------------------------------------------------- 10 #if defined key_trdtra || defined key_trdtrc || defined key_trdmld || defined key_trdmld_trc 11 !!---------------------------------------------------------------------- 12 !! trd_tra : Call the trend to be computed 13 !!---------------------------------------------------------------------- 14 USE dom_oce ! ocean domain 15 USE trdmod_oce ! ocean active mixed layer tracers trends 16 USE trdmod ! ocean active mixed layer tracers trends 17 USE trdmod_trc ! ocean passive mixed layer tracers trends 18 USE in_out_manager ! I/O manager 19 USE lib_mpp ! MPP library 20 USE wrk_nemo ! Memory allocation 21 6 !! History : 3.3 ! 2010-06 (C. Ethe) creation for the TRA/TRC merge 7 !! 3.5 ! 2012-02 (G. Madec) update the comments 8 !!---------------------------------------------------------------------- 9 10 !!---------------------------------------------------------------------- 11 !! trd_tra : pre-process the tracer trends 12 !! trd_tra_adv : transform a div(U.T) trend into a U.grad(T) trend 13 !! trd_tra_mng : tracer trend manager: dispatch to the diagnostic modules 14 !! trd_tra_iom : output 3D tracer trends using IOM 15 !!---------------------------------------------------------------------- 16 USE oce ! ocean dynamics and tracers variables 17 USE dom_oce ! ocean domain 18 USE sbc_oce ! surface boundary condition: ocean 19 USE zdf_oce ! ocean vertical physics 20 USE trd_oce ! trends: ocean variables 21 USE trdmod_trc ! ocean passive mixed layer tracers trends 22 USE trdglo ! trends: global domain averaged 23 USE trdpen ! trends: Potential ENergy 24 USE trdmxl ! ocean active mixed layer tracers trends 25 USE ldftra_oce ! ocean active tracers lateral physics 26 USE zdfddm ! vertical physics: double diffusion 27 USE phycst ! physical constants 28 USE in_out_manager ! I/O manager 29 USE iom ! I/O manager library 30 USE lib_mpp ! MPP library 31 USE wrk_nemo ! Memory allocation 22 32 23 33 IMPLICIT NONE 24 34 PRIVATE 25 35 26 PUBLIC trd_tra ! called by all traXX modules 27 28 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: trdtx, trdty, trdt !: 36 PUBLIC trd_tra ! called by all tra_... modules 37 38 REAL(wp) :: r2dt ! time-step, = 2 rdttra except at nit000 (=rdttra) if neuler=0 39 40 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: trdtx, trdty, trdt ! use to store the temperature trends 29 41 30 42 !! * Substitutions 31 43 # include "domzgr_substitute.h90" 44 # include "zdfddm_substitute.h90" 32 45 # include "vectopt_loop_substitute.h90" 33 46 !!---------------------------------------------------------------------- 34 !! NEMO/OPA 4.0 , NEMO Consortium (2011)47 !! NEMO/OPA 3.3 , NEMO Consortium (2010) 35 48 !! $Id$ 36 49 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) … … 39 52 40 53 INTEGER FUNCTION trd_tra_alloc() 41 !!--------------------------------------------------------------------- -------54 !!--------------------------------------------------------------------- 42 55 !! *** FUNCTION trd_tra_alloc *** 43 !!--------------------------------------------------------------------- -------56 !!--------------------------------------------------------------------- 44 57 ALLOCATE( trdtx(jpi,jpj,jpk) , trdty(jpi,jpj,jpk) , trdt(jpi,jpj,jpk) , STAT= trd_tra_alloc ) 45 58 ! … … 53 66 !! *** ROUTINE trd_tra *** 54 67 !! 55 !! ** Purpose : Dispatch all trends computation, e.g. vorticity, mld or 56 !! integral constraints 68 !! ** Purpose : pre-process tracer trends 57 69 !! 58 !! ** Method /usage : For the mixed-layer trend, the control surface can be either59 !! a mixed layer depth (time varying) or a fixed surface (jk level or bowl).60 !! Choose control surface with nn_ctls in namelist NAMTRD :61 !! nn_ctls = 0 : use mixed layer with density criterion62 !! nn_ctls = 1 : read index from file 'ctlsurf_idx'63 !! nn_ctls > 1 : use fixed level surface jk = nn_ctls64 !!---------------------------------------------------------------------- 65 !66 INTEGER , INTENT(in) :: kt ! time step67 CHARACTER(len=3) , INTENT(in) :: ctype ! tracers trends type 'TRA'/'TRC'68 INTEGER , INTENT(in) :: ktra ! tracerindex69 INTEGER , INTENT(in) :: ktrd ! tracer trend index70 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in) :: ptrd ! tracer trend or flux71 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in), OPTIONAL :: pun ! velocity72 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in), OPTIONAL :: ptra ! Tracer variablea73 !74 REAL(wp), POINTER, DIMENSION(:,:,:) :: ztrds75 !!---------------------------------------------------------------------- 76 70 !! ** Method : - mask the trend 71 !! - advection (ptra present) converte the incoming flux (U.T) 72 !! into trend (U.T => -U.grat(T)=div(U.T)-T.div(U)) through a 73 !! call to trd_tra_adv 74 !! - 'TRA' case : regroup T & S trends 75 !! - send the trends to trd_tra_mng (trdmod_trc) for further processing 76 !!---------------------------------------------------------------------- 77 INTEGER , INTENT(in) :: kt ! time step 78 CHARACTER(len=3) , INTENT(in) :: ctype ! tracers trends type 'TRA'/'TRC' 79 INTEGER , INTENT(in) :: ktra ! tracer index 80 INTEGER , INTENT(in) :: ktrd ! tracer trend index 81 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in) :: ptrd ! tracer trend or flux 82 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in), OPTIONAL :: pun ! now velocity 83 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in), OPTIONAL :: ptra ! now tracer variable 84 ! 85 INTEGER :: jk ! loop indices 86 REAL(wp), POINTER, DIMENSION(:,:,:) :: zwt, zws, ztrdt, ztrds ! 3D workspace 87 !!---------------------------------------------------------------------- 88 ! 77 89 CALL wrk_alloc( jpi, jpj, jpk, ztrds ) 78 79 IF( .NOT. ALLOCATED( trdtx ) ) THEN 90 ! 91 IF( .NOT. ALLOCATED( trdtx ) ) THEN ! allocate trdtra arrays 80 92 IF( trd_tra_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'trd_tra : unable to allocate arrays' ) 81 93 ENDIF 82 83 ! Control of optional arguments 84 IF( ctype == 'TRA' .AND. ktra == jp_tem ) THEN 85 IF( PRESENT( ptra ) ) THEN 86 SELECT CASE( ktrd ) ! shift depending on the direction 87 CASE( jptra_trd_xad ) ; CALL trd_tra_adv( ptrd, pun, ptra, 'X', trdtx ) 88 CASE( jptra_trd_yad ) ; CALL trd_tra_adv( ptrd, pun, ptra, 'Y', trdty ) 89 CASE( jptra_trd_zad ) ; CALL trd_tra_adv( ptrd, pun, ptra, 'Z', trdt ) 90 END SELECT 91 ELSE 92 trdt(:,:,:) = ptrd(:,:,:) 93 IF( ktrd == jptra_trd_bbc .OR. ktrd == jptra_trd_qsr ) THEN 94 ztrds(:,:,:) = 0. 95 CALL trd_mod( trdt, ztrds, ktrd, ctype, kt ) 96 END IF 97 END IF 98 END IF 99 100 IF( ctype == 'TRA' .AND. ktra == jp_sal ) THEN 101 IF( PRESENT( ptra ) ) THEN 102 SELECT CASE( ktrd ) ! shift depending on the direction 103 CASE( jptra_trd_xad ) 104 CALL trd_tra_adv( ptrd, pun, ptra, 'X', ztrds ) 105 CALL trd_mod( trdtx, ztrds, ktrd, ctype, kt ) 106 CASE( jptra_trd_yad ) 107 CALL trd_tra_adv( ptrd, pun, ptra, 'Y', ztrds ) 108 CALL trd_mod( trdty, ztrds, ktrd, ctype, kt ) 109 CASE( jptra_trd_zad ) 110 CALL trd_tra_adv( ptrd, pun, ptra, 'Z', ztrds ) 111 CALL trd_mod( trdt , ztrds, ktrd, ctype, kt ) 112 END SELECT 113 ELSE 114 ztrds(:,:,:) = ptrd(:,:,:) 115 CALL trd_mod( trdt, ztrds, ktrd, ctype, kt ) 116 END IF 117 END IF 118 119 IF( ctype == 'TRC' ) THEN 120 ! 121 IF( PRESENT( ptra ) ) THEN 122 SELECT CASE( ktrd ) ! shift depending on the direction 123 CASE( jptra_trd_xad ) 124 CALL trd_tra_adv( ptrd, pun, ptra, 'X', ztrds ) 125 CALL trd_mod_trc( ztrds, ktra, ktrd, kt ) 126 CASE( jptra_trd_yad ) 127 CALL trd_tra_adv( ptrd, pun, ptra, 'Y', ztrds ) 128 CALL trd_mod_trc( ztrds, ktra, ktrd, kt ) 129 CASE( jptra_trd_zad ) 130 CALL trd_tra_adv( ptrd, pun, ptra, 'Z', ztrds ) 131 CALL trd_mod_trc( ztrds, ktra, ktrd, kt ) 132 END SELECT 133 ELSE 134 ztrds(:,:,:) = ptrd(:,:,:) 135 CALL trd_mod_trc( ztrds, ktra, ktrd, kt ) 136 END IF 94 95 IF( ctype == 'TRA' .AND. ktra == jp_tem ) THEN !== Temperature trend ==! 96 ! 97 SELECT CASE( ktrd ) 98 ! ! advection: transform the advective flux into a trend 99 CASE( jptra_xad ) ; CALL trd_tra_adv( ptrd, pun, ptra, 'X', trdtx ) 100 CASE( jptra_yad ) ; CALL trd_tra_adv( ptrd, pun, ptra, 'Y', trdty ) 101 CASE( jptra_zad ) ; CALL trd_tra_adv( ptrd, pun, ptra, 'Z', trdt ) 102 CASE( jptra_bbc, & ! qsr, bbc: on temperature only, send to trd_tra_mng 103 & jptra_qsr ) ; trdt(:,:,:) = ptrd(:,:,:) * tmask(:,:,:) 104 ztrds(:,:,:) = 0._wp 105 CALL trd_tra_mng( trdt, ztrds, ktrd, kt ) 106 CASE DEFAULT ! other trends: masked trends 107 trdt(:,:,:) = ptrd(:,:,:) * tmask(:,:,:) ! mask & store 108 END SELECT 109 ! 110 ENDIF 111 112 IF( ctype == 'TRA' .AND. ktra == jp_sal ) THEN !== Salinity trends ==! 113 ! 114 SELECT CASE( ktrd ) 115 ! ! advection: transform the advective flux into a trend 116 ! ! and send T & S trends to trd_tra_mng 117 CASE( jptra_xad ) ; CALL trd_tra_adv( ptrd , pun , ptra, 'X' , ztrds ) 118 CALL trd_tra_mng( trdtx, ztrds, ktrd, kt ) 119 CASE( jptra_yad ) ; CALL trd_tra_adv( ptrd , pun , ptra, 'Y' , ztrds ) 120 CALL trd_tra_mng( trdty, ztrds, ktrd, kt ) 121 CASE( jptra_zad ) ; CALL trd_tra_adv( ptrd , pun , ptra, 'Z' , ztrds ) 122 CALL trd_tra_mng( trdt , ztrds, ktrd, kt ) 123 CASE( jptra_zdfp ) ! diagnose the "PURE" Kz trend (here: just before the swap) 124 ! ! iso-neutral diffusion case otherwise jptra_zdf is "PURE" 125 CALL wrk_alloc( jpi, jpj, jpk, zwt, zws, ztrdt ) 126 ! 127 zwt(:,:, 1 ) = 0._wp ; zws(:,:, 1 ) = 0._wp ! vertical diffusive fluxes 128 zwt(:,:,jpk) = 0._wp ; zws(:,:,jpk) = 0._wp 129 DO jk = 2, jpk 130 zwt(:,:,jk) = avt(:,:,jk) * ( tsa(:,:,jk-1,jp_tem) - tsa(:,:,jk,jp_tem) ) / fse3w(:,:,jk) * tmask(:,:,jk) 131 zws(:,:,jk) = fsavs(:,:,jk) * ( tsa(:,:,jk-1,jp_sal) - tsa(:,:,jk,jp_sal) ) / fse3w(:,:,jk) * tmask(:,:,jk) 132 END DO 133 ! 134 ztrdt(:,:,jpk) = 0._wp ; ztrds(:,:,jpk) = 0._wp 135 DO jk = 1, jpkm1 136 ztrdt(:,:,jk) = ( zwt(:,:,jk) - zwt(:,:,jk+1) ) / fse3t(:,:,jk) 137 ztrds(:,:,jk) = ( zws(:,:,jk) - zws(:,:,jk+1) ) / fse3t(:,:,jk) 138 END DO 139 CALL trd_tra_mng( ztrdt, ztrds, jptra_zdfp, kt ) 140 ! 141 CALL wrk_dealloc( jpi, jpj, jpk, zwt, zws, ztrdt ) 142 ! 143 CASE DEFAULT ! other trends: mask and send T & S trends to trd_tra_mng 144 ztrds(:,:,:) = ptrd(:,:,:) * tmask(:,:,:) 145 CALL trd_tra_mng( trdt, ztrds, ktrd, kt ) 146 END SELECT 147 ENDIF 148 149 IF( ctype == 'TRC' ) THEN !== passive tracer trend ==! 150 ! 151 SELECT CASE( ktrd ) 152 ! ! advection: transform the advective flux into a masked trend 153 CASE( jptra_xad ) ; CALL trd_tra_adv( ptrd , pun , ptra, 'X', ztrds ) 154 CASE( jptra_yad ) ; CALL trd_tra_adv( ptrd , pun , ptra, 'Y', ztrds ) 155 CASE( jptra_zad ) ; CALL trd_tra_adv( ptrd , pun , ptra, 'Z', ztrds ) 156 CASE DEFAULT ! other trends: just masked 157 ztrds(:,:,:) = ptrd(:,:,:) * tmask(:,:,:) 158 END SELECT 159 ! ! send trend to trd_mod_trc 160 CALL trd_mod_trc( ztrds, ktra, ktrd, kt ) 137 161 ! 138 162 ENDIF … … 147 171 !! *** ROUTINE trd_tra_adv *** 148 172 !! 149 !! ** Purpose : transformed the i-, j- or k-advective flux into thes 150 !! i-, j- or k-advective trends, resp. 151 !! ** Method : i-advective trends = -un. di-1[T] = -( di-1[fi] - tn di-1[un] ) 152 !! k-advective trends = -un. di-1[T] = -( dj-1[fi] - tn dj-1[un] ) 153 !! k-advective trends = -un. di+1[T] = -( dk+1[fi] - tn dk+1[un] ) 154 !!---------------------------------------------------------------------- 155 REAL(wp) , INTENT(in ), DIMENSION(jpi,jpj,jpk) :: pf ! advective flux in one direction 156 REAL(wp) , INTENT(in ), DIMENSION(jpi,jpj,jpk) :: pun ! now velocity in one direction 157 REAL(wp) , INTENT(in ), DIMENSION(jpi,jpj,jpk) :: ptn ! now or before tracer 158 CHARACTER(len=1), INTENT(in ) :: cdir ! X/Y/Z direction 159 REAL(wp) , INTENT(out), DIMENSION(jpi,jpj,jpk) :: ptrd ! advective trend in one direction 173 !! ** Purpose : transformed a advective flux into a masked advective trends 174 !! 175 !! ** Method : use the following transformation: -div(U.T) = - U grad(T) + T.div(U) 176 !! i-advective trends = -un. di-1[T] = -( di-1[fi] - tn di-1[un] ) 177 !! j-advective trends = -un. di-1[T] = -( dj-1[fi] - tn dj-1[un] ) 178 !! k-advective trends = -un. di+1[T] = -( dk+1[fi] - tn dk+1[un] ) 179 !! where fi is the incoming advective flux. 180 !!---------------------------------------------------------------------- 181 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in ) :: pf ! advective flux in one direction 182 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in ) :: pun ! now velocity in one direction 183 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in ) :: ptn ! now or before tracer 184 CHARACTER(len=1) , INTENT(in ) :: cdir ! X/Y/Z direction 185 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( out) :: ptrd ! advective trend in one direction 160 186 ! 161 187 INTEGER :: ji, jj, jk ! dummy loop indices 162 INTEGER :: ii, ij, ik ! index shift function of the direction 163 REAL(wp) :: zbtr ! local scalar 164 !!---------------------------------------------------------------------- 165 166 SELECT CASE( cdir ) ! shift depending on the direction 167 CASE( 'X' ) ; ii = 1 ; ij = 0 ; ik = 0 ! i-advective trend 168 CASE( 'Y' ) ; ii = 0 ; ij = 1 ; ik = 0 ! j-advective trend 169 CASE( 'Z' ) ; ii = 0 ; ij = 0 ; ik =-1 ! k-advective trend 188 INTEGER :: ii, ij, ik ! index shift as function of the direction 189 !!---------------------------------------------------------------------- 190 ! 191 SELECT CASE( cdir ) ! shift depending on the direction 192 CASE( 'X' ) ; ii = 1 ; ij = 0 ; ik = 0 ! i-trend 193 CASE( 'Y' ) ; ii = 0 ; ij = 1 ; ik = 0 ! j-trend 194 CASE( 'Z' ) ; ii = 0 ; ij = 0 ; ik =-1 ! k-trend 170 195 END SELECT 171 172 ! ! set to zero uncomputed values 173 ptrd(jpi,:,:) = 0.e0 ; ptrd(1,:,:) = 0.e0 174 ptrd(:,jpj,:) = 0.e0 ; ptrd(:,1,:) = 0.e0 175 ptrd(:,:,jpk) = 0.e0 176 ! 177 ! 178 DO jk = 1, jpkm1 196 ! 197 ! ! set to zero uncomputed values 198 ptrd(jpi,:,:) = 0._wp ; ptrd(1,:,:) = 0._wp 199 ptrd(:,jpj,:) = 0._wp ; ptrd(:,1,:) = 0._wp 200 ptrd(:,:,jpk) = 0._wp 201 ! 202 DO jk = 1, jpkm1 ! advective trend 179 203 DO jj = 2, jpjm1 180 204 DO ji = fs_2, fs_jpim1 ! vector opt. 181 zbtr = 1.e0/ ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )182 ptrd(ji,jj,jk) = - zbtr * ( pf (ji,jj,jk) - pf (ji-ii,jj-ij,jk-ik)&183 & - ( pun(ji,jj,jk) - pun(ji-ii,jj-ij,jk-ik) ) * ptn(ji,jj,jk))205 ptrd(ji,jj,jk) = - ( pf (ji,jj,jk) - pf (ji-ii,jj-ij,jk-ik) & 206 & - ( pun(ji,jj,jk) - pun(ji-ii,jj-ij,jk-ik) ) * ptn(ji,jj,jk) ) & 207 & / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) * tmask(ji,jj,jk) 184 208 END DO 185 209 END DO … … 188 212 END SUBROUTINE trd_tra_adv 189 213 190 # else 191 !!---------------------------------------------------------------------- 192 !! Default case : Dummy module No trend diagnostics 193 !!---------------------------------------------------------------------- 194 USE par_oce ! ocean variables trends 195 CONTAINS 196 SUBROUTINE trd_tra( kt, ctype, ktra, ktrd, ptrd, pu, ptra ) 197 !!---------------------------------------------------------------------- 198 INTEGER , INTENT(in) :: kt ! time step 199 CHARACTER(len=3) , INTENT(in) :: ctype ! tracers trends type 'TRA'/'TRC' 200 INTEGER , INTENT(in) :: ktra ! tracer index 201 INTEGER , INTENT(in) :: ktrd ! tracer trend index 202 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in) :: ptrd ! tracer trend 203 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in), OPTIONAL :: pu ! velocity 204 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in), OPTIONAL :: ptra ! Tracer variable 205 WRITE(*,*) 'trd_3d: You should not have seen this print! error ?', ptrd(1,1,1), ptra(1,1,1), pu(1,1,1), & 206 & ktrd, ktra, ctype, kt 207 END SUBROUTINE trd_tra 208 # endif 214 215 SUBROUTINE trd_tra_mng( ptrdx, ptrdy, ktrd, kt ) 216 !!--------------------------------------------------------------------- 217 !! *** ROUTINE trd_tra_mng *** 218 !! 219 !! ** Purpose : Dispatch all tracer trends computation, e.g. 3D output, 220 !! integral constraints, potential energy, and/or 221 !! mixed layer budget. 222 !!---------------------------------------------------------------------- 223 REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: ptrdx ! Temperature or U trend 224 REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: ptrdy ! Salinity or V trend 225 INTEGER , INTENT(in ) :: ktrd ! tracer trend index 226 INTEGER , INTENT(in ) :: kt ! time step 227 !!---------------------------------------------------------------------- 228 229 IF( neuler == 0 .AND. kt == nit000 ) THEN ; r2dt = rdt ! = rdtra (restart with Euler time stepping) 230 ELSEIF( kt <= nit000 + 1) THEN ; r2dt = 2. * rdt ! = 2 rdttra (leapfrog) 231 ENDIF 232 233 ! ! 3D output of tracers trends using IOM interface 234 IF( ln_tra_trd ) CALL trd_tra_iom ( ptrdx, ptrdy, ktrd, kt ) 235 236 ! ! Integral Constraints Properties for tracers trends !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 237 IF( ln_glo_trd ) CALL trd_glo( ptrdx, ptrdy, ktrd, 'TRA', kt ) 238 239 ! ! Potential ENergy trends 240 IF( ln_PE_trd ) CALL trd_pen( ptrdx, ptrdy, ktrd, kt, r2dt ) 241 242 ! ! Mixed layer trends for active tracers 243 IF( ln_tra_mxl ) THEN 244 !----------------------------------------------------------------------------------------------- 245 ! W.A.R.N.I.N.G : 246 ! jptra_ldf : called by traldf.F90 247 ! at this stage we store: 248 ! - the lateral geopotential diffusion (here, lateral = horizontal) 249 ! - and the iso-neutral diffusion if activated 250 ! jptra_zdf : called by trazdf.F90 251 ! * in case of iso-neutral diffusion we store the vertical diffusion component in the 252 ! lateral trend including the K_z contrib, which will be removed later (see trd_mxl) 253 !----------------------------------------------------------------------------------------------- 254 255 SELECT CASE ( ktrd ) 256 CASE ( jptra_xad ) ; CALL trd_mxl_zint( ptrdx, ptrdy, jpmxl_xad, '3D' ) ! zonal advection 257 CASE ( jptra_yad ) ; CALL trd_mxl_zint( ptrdx, ptrdy, jpmxl_yad, '3D' ) ! merid. advection 258 CASE ( jptra_zad ) ; CALL trd_mxl_zint( ptrdx, ptrdy, jpmxl_zad, '3D' ) ! vertical advection 259 CASE ( jptra_ldf ) ; CALL trd_mxl_zint( ptrdx, ptrdy, jpmxl_ldf, '3D' ) ! lateral diffusion 260 CASE ( jptra_bbl ) ; CALL trd_mxl_zint( ptrdx, ptrdy, jpmxl_bbl, '3D' ) ! bottom boundary layer 261 CASE ( jptra_zdf ) 262 IF( ln_traldf_iso ) THEN ; CALL trd_mxl_zint( ptrdx, ptrdy, jpmxl_ldf, '3D' ) ! lateral diffusion (K_z) 263 ELSE ; CALL trd_mxl_zint( ptrdx, ptrdy, jpmxl_zdf, '3D' ) ! vertical diffusion (K_z) 264 ENDIF 265 CASE ( jptra_dmp ) ; CALL trd_mxl_zint( ptrdx, ptrdy, jpmxl_dmp, '3D' ) ! internal 3D restoring (tradmp) 266 CASE ( jptra_qsr ) ; CALL trd_mxl_zint( ptrdx, ptrdy, jpmxl_for, '3D' ) ! air-sea : penetrative sol radiat 267 CASE ( jptra_nsr ) ; ptrdx(:,:,2:jpk) = 0._wp ; ptrdy(:,:,2:jpk) = 0._wp 268 CALL trd_mxl_zint( ptrdx, ptrdy, jpmxl_for, '2D' ) ! air-sea : non penetr sol radiation 269 CASE ( jptra_bbc ) ; CALL trd_mxl_zint( ptrdx, ptrdy, jpmxl_bbc, '3D' ) ! bottom bound cond (geoth flux) 270 CASE ( jptra_npc ) ; CALL trd_mxl_zint( ptrdx, ptrdy, jpmxl_npc, '3D' ) ! non penetr convect adjustment 271 CASE ( jptra_atf ) ; CALL trd_mxl_zint( ptrdx, ptrdy, jpmxl_atf, '3D' ) ! asselin time filter (last trend) 272 ! 273 CALL trd_mxl( kt, r2dt ) ! trends: Mixed-layer (output) 274 END SELECT 275 ! 276 ENDIF 277 ! 278 END SUBROUTINE trd_tra_mng 279 280 281 SUBROUTINE trd_tra_iom( ptrdx, ptrdy, ktrd, kt ) 282 !!--------------------------------------------------------------------- 283 !! *** ROUTINE trd_tra_iom *** 284 !! 285 !! ** Purpose : output 3D tracer trends using IOM 286 !!---------------------------------------------------------------------- 287 REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: ptrdx ! Temperature or U trend 288 REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: ptrdy ! Salinity or V trend 289 INTEGER , INTENT(in ) :: ktrd ! tracer trend index 290 INTEGER , INTENT(in ) :: kt ! time step 291 !! 292 INTEGER :: ji, jj, jk ! dummy loop indices 293 INTEGER :: ikbu, ikbv ! local integers 294 REAL(wp), POINTER, DIMENSION(:,:) :: z2dx, z2dy ! 2D workspace 295 !!---------------------------------------------------------------------- 296 ! 297 !!gm Rq: mask the trends already masked in trd_tra, but lbc_lnk should probably be added 298 ! 299 SELECT CASE( ktrd ) 300 CASE( jptra_xad ) ; CALL iom_put( "ttrd_xad" , ptrdx ) ! x- horizontal advection 301 CALL iom_put( "strd_xad" , ptrdy ) 302 CASE( jptra_yad ) ; CALL iom_put( "ttrd_yad" , ptrdx ) ! y- horizontal advection 303 CALL iom_put( "strd_yad" , ptrdy ) 304 CASE( jptra_zad ) ; CALL iom_put( "ttrd_zad" , ptrdx ) ! z- vertical advection 305 CALL iom_put( "strd_zad" , ptrdy ) 306 IF( .NOT. lk_vvl ) THEN ! cst volume : adv flux through z=0 surface 307 CALL wrk_alloc( jpi, jpj, z2dx, z2dy ) 308 z2dx(:,:) = wn(:,:,1) * tsn(:,:,1,jp_tem) / fse3t(:,:,1) 309 z2dy(:,:) = wn(:,:,1) * tsn(:,:,1,jp_sal) / fse3t(:,:,1) 310 CALL iom_put( "ttrd_sad", z2dx ) 311 CALL iom_put( "strd_sad", z2dy ) 312 CALL wrk_dealloc( jpi, jpj, z2dx, z2dy ) 313 ENDIF 314 CASE( jptra_ldf ) ; CALL iom_put( "ttrd_ldf" , ptrdx ) ! lateral diffusion 315 CALL iom_put( "strd_ldf" , ptrdy ) 316 CASE( jptra_zdf ) ; CALL iom_put( "ttrd_zdf" , ptrdx ) ! vertical diffusion (including Kz contribution) 317 CALL iom_put( "strd_zdf" , ptrdy ) 318 CASE( jptra_zdfp ) ; CALL iom_put( "ttrd_zdfp", ptrdx ) ! PURE vertical diffusion (no isoneutral contribution) 319 CALL iom_put( "strd_zdfp", ptrdy ) 320 CASE( jptra_dmp ) ; CALL iom_put( "ttrd_dmp" , ptrdx ) ! internal restoring (damping) 321 CALL iom_put( "strd_dmp" , ptrdy ) 322 CASE( jptra_bbl ) ; CALL iom_put( "ttrd_bbl" , ptrdx ) ! bottom boundary layer 323 CALL iom_put( "strd_bbl" , ptrdy ) 324 CASE( jptra_npc ) ; CALL iom_put( "ttrd_npc" , ptrdx ) ! static instability mixing 325 CALL iom_put( "strd_npc" , ptrdy ) 326 CASE( jptra_nsr ) ; CALL iom_put( "ttrd_qns" , ptrdx ) ! surface forcing + runoff (ln_rnf=T) 327 CALL iom_put( "strd_cdt" , ptrdy ) 328 CASE( jptra_qsr ) ; CALL iom_put( "ttrd_qsr" , ptrdx ) ! penetrative solar radiat. (only on temperature) 329 CASE( jptra_bbc ) ; CALL iom_put( "ttrd_bbc" , ptrdx ) ! geothermal heating (only on temperature) 330 CASE( jptra_atf ) ; CALL iom_put( "ttrd_atf" , ptrdx ) ! asselin time Filter 331 CALL iom_put( "strd_atf" , ptrdy ) 332 END SELECT 333 ! 334 END SUBROUTINE trd_tra_iom 335 209 336 !!====================================================================== 210 337 END MODULE trdtra -
branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/TRD/trdvor.F90
r3294 r3876 4 4 !! Ocean diagnostics: momentum trends 5 5 !!===================================================================== 6 !! History : 1.0 ! 04-2006 (L. Brunier, A-M. Treguier) Original code 7 !! 2.0 ! 04-2008 (C. Talandier) New trends organization 6 !! History : 1.0 ! 2006-01 (L. Brunier, A-M. Treguier) Original code 7 !! 2.0 ! 2008-04 (C. Talandier) New trends organization 8 !! 3.5 ! 2012-02 (G. Madec) regroup beta.V computation with pvo trend 8 9 !!---------------------------------------------------------------------- 9 #if defined key_trdvor || defined key_esopa 10 !!---------------------------------------------------------------------- 11 !! 'key_trdvor' : momentum trend diagnostics 10 12 11 !!---------------------------------------------------------------------- 13 12 !! trd_vor : momentum trends averaged over the depth … … 17 16 USE oce ! ocean dynamics and tracers variables 18 17 USE dom_oce ! ocean space and time domain variables 19 USE trd mod_oce ! ocean variables trends18 USE trd_oce ! trends: ocean variables 20 19 USE zdf_oce ! ocean vertical physics 21 USE in_out_manager ! I/O manager20 USE sbc_oce ! surface boundary condition: ocean 22 21 USE phycst ! Define parameters for the routines 23 22 USE ldfdyn_oce ! ocean active tracers: lateral physics 24 23 USE dianam ! build the name of file (routine) 25 24 USE zdfmxl ! mixed layer depth 25 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 26 USE in_out_manager ! I/O manager 26 27 USE ioipsl ! NetCDF library 27 USE lbclnk ! ocean lateral boundary conditions (or mpp link)28 28 USE lib_mpp ! MPP library 29 29 USE wrk_nemo ! Memory allocation 30 31 30 32 31 IMPLICIT NONE … … 37 36 END INTERFACE 38 37 39 PUBLIC trd_vor ! routine called by step.F90 40 PUBLIC trd_vor_zint ! routine called by dynamics routines 38 PUBLIC trd_vor ! routine called by trddyn.F90 41 39 PUBLIC trd_vor_init ! routine called by opa.F90 42 40 PUBLIC trd_vor_alloc ! routine called by nemogcm.F90 … … 80 78 IF( trd_vor_alloc /= 0 ) CALL ctl_warn('trd_vor_alloc: failed to allocate arrays') 81 79 END FUNCTION trd_vor_alloc 80 81 82 SUBROUTINE trd_vor( putrd, pvtrd, ktrd, kt ) 83 !!---------------------------------------------------------------------- 84 !! *** ROUTINE trd_vor *** 85 !! 86 !! ** Purpose : computation of cumulated trends over analysis period 87 !! and make outputs (NetCDF or DIMG format) 88 !!---------------------------------------------------------------------- 89 REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: putrd, pvtrd ! U and V trends 90 INTEGER , INTENT(in ) :: ktrd ! trend index 91 INTEGER , INTENT(in ) :: kt ! time step 92 ! 93 INTEGER :: ji, jj ! dummy loop indices 94 REAL(wp), POINTER, DIMENSION(:,:) :: ztswu, ztswv ! 2D workspace 95 !!---------------------------------------------------------------------- 96 97 CALL wrk_alloc( jpi, jpj, ztswu, ztswv ) 98 99 SELECT CASE( ktrd ) 100 CASE( jpdyn_hpg ) ; CALL trd_vor_zint( putrd, pvtrd, jpvor_prg ) ! Hydrostatique Pressure Gradient 101 CASE( jpdyn_keg ) ; CALL trd_vor_zint( putrd, pvtrd, jpvor_keg ) ! KE Gradient 102 CASE( jpdyn_rvo ) ; CALL trd_vor_zint( putrd, pvtrd, jpvor_rvo ) ! Relative Vorticity 103 CASE( jpdyn_pvo ) ; CALL trd_vor_zint( putrd, pvtrd, jpvor_pvo ) ! Planetary Vorticity Term 104 CASE( jpdyn_ldf ) ; CALL trd_vor_zint( putrd, pvtrd, jpvor_ldf ) ! Horizontal Diffusion 105 CASE( jpdyn_zad ) ; CALL trd_vor_zint( putrd, pvtrd, jpvor_zad ) ! Vertical Advection 106 CASE( jpdyn_spg ) ; CALL trd_vor_zint( putrd, pvtrd, jpvor_spg ) ! Surface Pressure Grad. 107 CASE( jpdyn_zdf ) ! Vertical Diffusion 108 ztswu(:,:) = 0.e0 ; ztswv(:,:) = 0.e0 109 DO jj = 2, jpjm1 ! wind stress trends 110 DO ji = fs_2, fs_jpim1 ! vector opt. 111 ztswu(ji,jj) = 0.5 * ( utau_b(ji,jj) + utau(ji,jj) ) / ( fse3u(ji,jj,1) * rau0 ) 112 ztswv(ji,jj) = 0.5 * ( vtau_b(ji,jj) + vtau(ji,jj) ) / ( fse3v(ji,jj,1) * rau0 ) 113 END DO 114 END DO 115 ! 116 CALL trd_vor_zint( putrd, pvtrd, jpvor_zdf ) ! zdf trend including surf./bot. stresses 117 CALL trd_vor_zint( ztswu, ztswv, jpvor_swf ) ! surface wind stress 118 CASE( jpdyn_bfr ) 119 CALL trd_vor_zint( putrd, pvtrd, jpvor_bfr ) ! Bottom stress 120 ! 121 CASE( jpdyn_atf ) ! last trends: perform the output of 2D vorticity trends 122 CALL trd_vor_iom( kt ) 123 END SELECT 124 ! 125 CALL wrk_dealloc( jpi, jpj, ztswu, ztswv ) 126 ! 127 END SUBROUTINE trd_vor 82 128 83 129 … … 109 155 !! trends output in netCDF format using ioipsl 110 156 !!---------------------------------------------------------------------- 111 !112 157 INTEGER , INTENT(in ) :: ktrd ! ocean trend index 113 158 REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: putrdvor ! u vorticity trend … … 131 176 ! ===================================== 132 177 133 SELECT CASE (ktrd)134 ! 135 CASE (jpvor_bfr) ! bottom friction178 SELECT CASE( ktrd ) 179 ! 180 CASE( jpvor_bfr ) ! bottom friction 136 181 DO jj = 2, jpjm1 137 182 DO ji = fs_2, fs_jpim1 … … 143 188 END DO 144 189 ! 145 CASE (jpvor_swf) ! wind stress190 CASE( jpvor_swf ) ! wind stress 146 191 zudpvor(:,:) = putrdvor(:,:) * fse3u(:,:,1) * e1u(:,:) * umask(:,:,1) 147 192 zvdpvor(:,:) = pvtrdvor(:,:) * fse3v(:,:,1) * e2v(:,:) * vmask(:,:,1) … … 154 199 155 200 ! Curl 156 DO ji =1,jpim1157 DO jj =1,jpjm1201 DO ji = 1, jpim1 202 DO jj = 1, jpjm1 158 203 vortrd(ji,jj,ktrd) = ( zvdpvor(ji+1,jj) - zvdpvor(ji,jj) & 159 204 & - ( zudpvor(ji,jj+1) - zudpvor(ji,jj) ) ) / ( e1f(ji,jj) * e2f(ji,jj) ) … … 229 274 END DO 230 275 231 ! Save Beta.V term to avoid average before Curl232 ! Beta.V : intergration, noaverage233 IF( ktrd == jpvor_ bev) THEN276 ! Planetary vorticity: 2nd computation (Beta.V term) store the vertical sum 277 ! as Beta.V term need intergration, not average 278 IF( ktrd == jpvor_pvo ) THEN 234 279 zubet(:,:) = zudpvor(:,:) 235 280 zvbet(:,:) = zvdpvor(:,:) 236 ENDIF 237 238 ! Average except for Beta.V 281 DO ji = 1, jpim1 282 DO jj = 1, jpjm1 283 vortrd(ji,jj,jpvor_bev) = ( zvbet(ji+1,jj) - zvbet(ji,jj) & 284 & - ( zubet(ji,jj+1) - zubet(ji,jj) ) ) / ( e1f(ji,jj) * e2f(ji,jj) ) 285 END DO 286 END DO 287 ! Average of the Curl and Surface mask 288 vortrd(:,:,jpvor_bev) = vortrd(:,:,jpvor_bev) * hur(:,:) * fmask(:,:,1) 289 ENDIF 290 ! 291 ! Average 239 292 zudpvor(:,:) = zudpvor(:,:) * hur(:,:) 240 293 zvdpvor(:,:) = zvdpvor(:,:) * hvr(:,:) 241 294 ! 242 295 ! Curl 243 296 DO ji=1,jpim1 … … 247 300 END DO 248 301 END DO 249 250 302 ! Surface mask 251 303 vortrd(:,:,ktrd) = vortrd(:,:,ktrd) * fmask(:,:,1) 252 253 ! Special treatement for the Beta.V term254 ! Compute the Curl of the Beta.V term which is not averaged255 IF( ktrd == jpvor_bev ) THEN256 DO ji=1,jpim1257 DO jj=1,jpjm1258 vortrd(ji,jj,jpvor_bev) = ( zvbet(ji+1,jj) - zvbet(ji,jj) &259 & - ( zubet(ji,jj+1) - zubet(ji,jj) ) ) / ( e1f(ji,jj) * e2f(ji,jj) )260 END DO261 END DO262 263 ! Average on the Curl264 vortrd(:,:,jpvor_bev) = vortrd(:,:,jpvor_bev) * hur(:,:)265 266 ! Surface mask267 vortrd(:,:,jpvor_bev) = vortrd(:,:,jpvor_bev) * fmask(:,:,1)268 ENDIF269 304 270 305 IF( ndebug /= 0 ) THEN … … 278 313 279 314 280 SUBROUTINE trd_vor ( kt )315 SUBROUTINE trd_vor_iom( kt ) 281 316 !!---------------------------------------------------------------------- 282 317 !! *** ROUTINE trd_vor *** … … 285 320 !! and make outputs (NetCDF or DIMG format) 286 321 !!---------------------------------------------------------------------- 287 ! 288 INTEGER, INTENT(in) :: kt ! ocean time-step index 322 INTEGER , INTENT(in ) :: kt ! time step 289 323 ! 290 324 INTEGER :: ji, jj, jk, jl ! dummy loop indices … … 305 339 306 340 IF( kt > nit000 ) vor_avrb(:,:) = vor_avr(:,:) 307 308 IF( ndebug /= 0 ) THEN309 WRITE(numout,*) ' debuging trd_vor: I.1 done '310 CALL FLUSH(numout)311 ENDIF312 341 313 342 ! I.2 vertically integrated vorticity … … 330 359 331 360 ! Curl 332 DO ji =1,jpim1333 DO jj =1,jpjm1361 DO ji = 1, jpim1 362 DO jj = 1, jpjm1 334 363 vor_avr(ji,jj) = ( ( zvn(ji+1,jj) - zvn(ji,jj) ) & 335 364 & - ( zun(ji,jj+1) - zun(ji,jj) ) ) / ( e1f(ji,jj) * e2f(ji,jj) ) * fmask(ji,jj,1) … … 337 366 END DO 338 367 339 IF( ndebug /= 0 ) THEN340 WRITE(numout,*) ' debuging trd_vor: I.2 done'341 CALL FLUSH(numout)342 ENDIF343 344 368 ! ================================= 345 369 ! II. Cumulated trends … … 351 375 vor_avrbb(:,:) = vor_avrb(:,:) 352 376 vor_avrbn(:,:) = vor_avr (:,:) 353 ENDIF354 355 IF( ndebug /= 0 ) THEN356 WRITE(numout,*) ' debuging trd_vor: I1.1 done'357 CALL FLUSH(numout)358 377 ENDIF 359 378 … … 371 390 ENDIF 372 391 373 IF( ndebug /= 0 ) THEN374 WRITE(numout,*) ' debuging trd_vor: II.2 done'375 CALL FLUSH(numout)376 ENDIF377 378 392 ! ============================================= 379 393 ! III. Output in netCDF + residual computation … … 391 405 vor_avrtot(:,:) = ( vor_avr(:,:) - vor_avrbn(:,:) + vor_avrb(:,:) - vor_avrbb(:,:) ) * zmean 392 406 393 IF( ndebug /= 0 ) THEN394 WRITE(numout,*) ' zmean = ',zmean395 WRITE(numout,*) ' debuging trd_vor: III.1 done'396 CALL FLUSH(numout)397 ENDIF398 407 399 408 ! III.2 compute residual … … 406 415 CALL lbc_lnk( vor_avrres, 'F', 1. ) 407 416 408 IF( ndebug /= 0 ) THEN409 WRITE(numout,*) ' debuging trd_vor: III.2 done'410 CALL FLUSH(numout)411 ENDIF412 417 413 418 ! III.3 time evolution array swap … … 415 420 vor_avrbb(:,:) = vor_avrb(:,:) 416 421 vor_avrbn(:,:) = vor_avr (:,:) 417 418 IF( ndebug /= 0 ) THEN419 WRITE(numout,*) ' debuging trd_vor: III.3 done'420 CALL FLUSH(numout)421 ENDIF422 422 ! 423 423 nmoydpvor = 0 … … 463 463 CALL wrk_dealloc( jpi, jpj, zun, zvn ) 464 464 ! 465 END SUBROUTINE trd_vor 465 END SUBROUTINE trd_vor_iom 466 466 467 467 … … 587 587 END SUBROUTINE trd_vor_init 588 588 589 #else590 !!----------------------------------------------------------------------591 !! Default option : Empty module592 !!----------------------------------------------------------------------593 INTERFACE trd_vor_zint594 MODULE PROCEDURE trd_vor_zint_2d, trd_vor_zint_3d595 END INTERFACE596 CONTAINS597 SUBROUTINE trd_vor( kt ) ! Empty routine598 WRITE(*,*) 'trd_vor: You should not have seen this print! error?', kt599 END SUBROUTINE trd_vor600 SUBROUTINE trd_vor_zint_2d( putrdvor, pvtrdvor, ktrd )601 REAL, DIMENSION(:,:), INTENT( inout ) :: putrdvor, pvtrdvor602 INTEGER, INTENT( in ) :: ktrd ! ocean trend index603 WRITE(*,*) 'trd_vor_zint_2d: You should not have seen this print! error?', putrdvor(1,1), pvtrdvor(1,1), ktrd604 END SUBROUTINE trd_vor_zint_2d605 SUBROUTINE trd_vor_zint_3d( putrdvor, pvtrdvor, ktrd )606 REAL, DIMENSION(:,:,:), INTENT( inout ) :: putrdvor, pvtrdvor607 INTEGER, INTENT( in ) :: ktrd ! ocean trend index608 WRITE(*,*) 'trd_vor_zint_3d: You should not have seen this print! error?', putrdvor(1,1,1), pvtrdvor(1,1,1), ktrd609 END SUBROUTINE trd_vor_zint_3d610 SUBROUTINE trd_vor_init ! Empty routine611 WRITE(*,*) 'trd_vor_init: You should not have seen this print! error?'612 END SUBROUTINE trd_vor_init613 #endif614 589 !!====================================================================== 615 590 END MODULE trdvor -
branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/TRD/trdvor_oce.F90
r2715 r3876 4 4 !! Ocean trends : set vorticity trend variables 5 5 !!====================================================================== 6 !! History : 9.0 ! 04-2006 (L. Brunier, A-M. Treguier) Original code6 !! History : 1.0 ! 04-2006 (L. Brunier, A-M. Treguier) Original code 7 7 !!---------------------------------------------------------------------- 8 9 !!---------------------------------------------------------------------- 8 10 9 USE par_oce ! ocean parameters 11 10 … … 13 12 PRIVATE 14 13 15 #if defined key_trdvor16 LOGICAL, PUBLIC, PARAMETER :: lk_trdvor = .TRUE. !: momentum trend flag17 #else18 LOGICAL, PUBLIC, PARAMETER :: lk_trdvor = .FALSE. !: momentum trend flag19 #endif20 14 ! !!* vorticity trends index 21 15 INTEGER, PUBLIC, PARAMETER :: jpltot_vor = 11 !: Number of vorticity trend terms -
branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfmxl.F90
r3294 r3876 6 6 !! History : 1.0 ! 2003-08 (G. Madec) original code 7 7 !! 3.2 ! 2009-07 (S. Masson, G. Madec) IOM + merge of DO-loop 8 !! 3.5 ! 2012-03 (G. Madec) make public the density criteria for trdmxl 8 9 !!---------------------------------------------------------------------- 9 10 !! zdf_mxl : Compute the turbocline and mixed layer depths. … … 29 30 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hmlp !: mixed layer depth (rho=rho0+zdcrit) [m] 30 31 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hmlpt !: mixed layer depth at t-points [m] 32 33 REAL(wp), PUBLIC :: rho_c = 0.01_wp !: density criterion for mixed layer depth 34 REAL(wp) :: avt_c = 5.e-4_wp ! Kz criterion for the turbocline depth 31 35 32 36 !! * Substitutions … … 63 67 !! ** Method : The mixed layer depth is the shallowest W depth with 64 68 !! the density of the corresponding T point (just bellow) bellow a 65 !! given value defined locally as rho(10m) + zrho_c69 !! given value defined locally as rho(10m) + rho_c 66 70 !! The turbocline depth is the depth at which the vertical 67 71 !! eddy diffusivity coefficient (resulting from the vertical physics 68 72 !! alone, not the isopycnal part, see trazdf.F) fall below a given 69 !! value defined locally (avt_c here taken equal to 5 cm/s2 )73 !! value defined locally (avt_c here taken equal to 5 cm/s2 by default) 70 74 !! 71 75 !! ** Action : nmln, hmld, hmlp, hmlpt … … 76 80 INTEGER :: iikn, iiki ! temporary integer within a do loop 77 81 INTEGER, POINTER, DIMENSION(:,:) :: imld ! temporary workspace 78 REAL(wp) :: zrho_c = 0.01_wp ! density criterion for mixed layer depth79 REAL(wp) :: zavt_c = 5.e-4_wp ! Kz criterion for the turbocline depth80 82 !!---------------------------------------------------------------------- 81 83 ! … … 98 100 DO jj = 1, jpj 99 101 DO ji = 1, jpi 100 IF( rhop(ji,jj,jk) > rhop(ji,jj,nla10) + zrho_c ) nmln(ji,jj) = jk ! Mixed layer101 IF( avt (ji,jj,jk) < zavt_c ) imld(ji,jj) = jk ! Turbocline102 IF( rhop(ji,jj,jk) > rhop(ji,jj,nla10) + rho_c ) nmln(ji,jj) = jk ! Mixed layer 103 IF( avt (ji,jj,jk) < avt_c ) imld(ji,jj) = jk ! Turbocline 102 104 END DO 103 105 END DO -
branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90
r3769 r3876 59 59 USE zdfini ! vertical physics setting (zdf_init routine) 60 60 USE phycst ! physical constant (par_cst routine) 61 USE trd mod ! momentum/tracers trends (trd_mod_init routine)61 USE trdini ! dyn/tra trends initialization (trd_init routine) 62 62 USE asminc ! assimilation increments 63 63 USE asmbkg ! writing out state trajectory … … 395 395 IF( lk_diadct ) CALL dia_dct_init ! Sections tranports 396 396 CALL dia_hsb_init ! heat content, salt content and volume budgets 397 CALL trd_mod_init ! Mixed-layer/Vorticity/Integral constraintstrends397 CALL trd_init ! momentum/tracer trends 398 398 IF( lk_diaobs ) THEN ! Observation & model comparison 399 399 CALL dia_obs_init ! Initialize observational data … … 574 574 !! ** Method : 575 575 !!---------------------------------------------------------------------- 576 INTEGER, INTENT(in) :: num_pes! The number of MPI processes we have576 INTEGER, INTENT(in) :: num_pes ! The number of MPI processes we have 577 577 ! 578 578 INTEGER, PARAMETER :: nfactmax = 20 … … 583 583 INTEGER, DIMENSION(nfactmax) :: ifact ! Array of factors 584 584 !!---------------------------------------------------------------------- 585 585 ! 586 586 ierr = 0 587 587 ! 588 588 CALL factorise( ifact, nfactmax, nfact, num_pes, ierr ) 589 589 ! 590 590 IF( nfact <= 1 ) THEN 591 591 WRITE (numout, *) 'WARNING: factorisation of number of PEs failed' … … 629 629 INTEGER, PARAMETER :: ntest = 14 630 630 INTEGER :: ilfax(ntest) 631 631 ! 632 632 ! lfax contains the set of allowed factors. 633 633 data (ilfax(jl),jl=1,ntest) / 16384, 8192, 4096, 2048, 1024, 512, 256, & … … 680 680 681 681 #if defined key_mpp_mpi 682 682 683 SUBROUTINE nemo_northcomms 683 684 !!====================================================================== 684 685 !! *** ROUTINE nemo_northcomms *** 685 !! nemo_northcomms 686 !! nemo_northcomms : Setup for north fold exchanges with explicit peer to peer messaging 686 687 !!===================================================================== 687 688 !!---------------------------------------------------------------------- … … 691 692 !! 1.0 ! 2011-10 (A. C. Coward, NOCS & J. Donners, PRACE) 692 693 !!---------------------------------------------------------------------- 693 694 694 INTEGER :: ji, jj, jk, ij, jtyp ! dummy loop indices 695 695 INTEGER :: ijpj ! number of rows involved in north-fold exchange 696 696 INTEGER :: northcomms_alloc ! allocate return status 697 REAL(wp), ALLOCATABLE, DIMENSION ( :,: ) :: znnbrs ! workspace 698 LOGICAL, ALLOCATABLE, DIMENSION ( : ) :: lrankset ! workspace 697 REAL(wp), ALLOCATABLE, DIMENSION (:,:) :: znnbrs ! workspace 698 LOGICAL, ALLOCATABLE, DIMENSION (:) :: lrankset ! workspace 699 !!---------------------------------------------------------------------- 699 700 700 701 IF(lwp) WRITE(numout,*) … … 702 703 IF(lwp) WRITE(numout,*) '~~~~~~~~~~' 703 704 704 !!----------------------------------------------------------------------705 705 ALLOCATE( znnbrs(jpi,jpj), stat = northcomms_alloc ) 706 706 ALLOCATE( lrankset(jpnij), stat = northcomms_alloc ) … … 724 724 ! Exchange and store ranks on northern rows 725 725 726 DO jtyp = 1, 4726 DO jtyp = 1, 4 727 727 728 728 lrankset = .FALSE. 729 znnbrs = narea729 znnbrs = narea 730 730 SELECT CASE (jtyp) 731 731 CASE(1) … … 739 739 END SELECT 740 740 741 IF ( njmppt(narea) .EQ.MAXVAL( njmppt ) ) THEN741 IF ( njmppt(narea) == MAXVAL( njmppt ) ) THEN 742 742 DO jj = nlcj-ijpj+1, nlcj 743 743 ij = jj - nlcj + ijpj … … 817 817 818 818 END SUBROUTINE nemo_northcomms 819 819 820 #else 820 821 SUBROUTINE nemo_northcomms ! Dummy routine … … 822 823 END SUBROUTINE nemo_northcomms 823 824 #endif 825 824 826 !!====================================================================== 825 827 END MODULE nemogcm
Note: See TracChangeset
for help on using the changeset viewer.