Changeset 1559 for trunk/NEMO
- Timestamp:
- 2009-07-29T16:03:14+02:00 (15 years ago)
- Location:
- trunk/NEMO
- Files:
-
- 11 edited
-
LIM_SRC_2/limthd_2.F90 (modified) (6 diffs)
-
OPA_SRC/BDY/bdydyn.F90 (modified) (1 diff)
-
OPA_SRC/DIA/diaptr.F90 (modified) (16 diffs)
-
OPA_SRC/SBC/sbcana.F90 (modified) (9 diffs)
-
OPA_SRC/TRA/traadv_cen2.F90 (modified) (8 diffs)
-
OPA_SRC/TRA/traadv_qck.F90 (modified) (18 diffs)
-
OPA_SRC/ZDF/zdfini.F90 (modified) (5 diffs)
-
OPA_SRC/ZDF/zdfmxl.F90 (modified) (5 diffs)
-
OPA_SRC/daymod.F90 (modified) (6 diffs)
-
OPA_SRC/eosbn2.F90 (modified) (18 diffs)
-
OPA_SRC/lib_mpp.F90 (modified) (2 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMO/LIM_SRC_2/limthd_2.F90
r1482 r1559 4 4 !! LIM thermo ice model : ice thermodynamic 5 5 !!====================================================================== 6 !! History : 1.0 ! 7 !! 2.0 ! 8 !! 2.0 ! 9 !! - ! 6 !! History : 1.0 ! 2000-01 (LIM) 7 !! 2.0 ! 2002-07 (C. Ethe, G. Madec) F90 8 !! 2.0 ! 2003-08 (C. Ethe) add lim_thd_init 9 !! - ! 2008-2008 (A. Caubel, G. Madec, E. Maisonnave, S. Masson ) generic coupled interface 10 10 !!--------------------------------------------------------------------- 11 11 #if defined key_lim2 … … 47 47 # include "vectopt_loop_substitute.h90" 48 48 !!-------- ------------------------------------------------------------- 49 !! NEMO/LIM 2.0, UCL-LOCEAN-IPSL (2008)49 !! NEMO/LIM 3.2, UCL-LOCEAN-IPSL (2009) 50 50 !! $Id$ 51 51 !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) … … 82 82 REAL(wp) :: zinda ! switch for test. the val. of concen. 83 83 REAL(wp) :: zindb, zindg ! switches for test. the val of arg 84 REAL(wp) :: zfricp ! temporary scalar 84 85 REAL(wp) :: za , zh, zthsnice ! 85 86 REAL(wp) :: zfric_u ! friction velocity … … 87 88 REAL(wp) :: zfontn ! heat flux from snow thickness 88 89 REAL(wp) :: zfntlat, zpareff ! test. the val. of lead heat budget 89 REAL(wp) :: zfi ! temporary scalar 90 REAL(wp), DIMENSION(jpi,jpj) :: ztmp ! working array 90 REAL(wp), DIMENSION(jpi,jpj) :: ztmp ! 2D workspace 91 91 REAL(wp), DIMENSION(jpi,jpj) :: zqlbsbq ! link with lead energy budget qldif 92 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zmsk ! working array92 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zmsk ! 3D workspace 93 93 !!------------------------------------------------------------------- 94 94 … … 179 179 zindb = tms(ji,jj) * ( 1.0 - MAX( rzero , SIGN( rone , - zthsnice ) ) ) 180 180 pfrld(ji,jj) = frld(ji,jj) 181 zinda = 1.0 - MAX( rzero , SIGN( rone , - ( 1.0 - pfrld(ji,jj) ) ) ) 181 zfricp = 1.0 - frld(ji,jj) 182 zinda = 1.0 - MAX( rzero , SIGN( rone , - zfricp ) ) 182 183 183 184 ! solar irradiance transmission at the mixed layer bottom and used in the lead heat budget … … 192 193 ! partial computation of the lead energy budget (qldif) 193 194 #if defined key_coupled 194 zfi = 1.0 - pfrld(ji,jj)195 195 qldif(ji,jj) = tms(ji,jj) * rdt_ice & 196 & * ( ( qsr_tot(ji,jj) - qsr_ice(ji,jj,1) * zf i) * ( 1.0 - thcm(ji,jj) ) &197 & + ( qns_tot(ji,jj) - qns_ice(ji,jj,1) * zf i) &196 & * ( ( qsr_tot(ji,jj) - qsr_ice(ji,jj,1) * zfricp ) * ( 1.0 - thcm(ji,jj) ) & 197 & + ( qns_tot(ji,jj) - qns_ice(ji,jj,1) * zfricp ) & 198 198 & + frld(ji,jj) * ( fdtcn(ji,jj) + ( 1.0 - zindb ) * fsbbq(ji,jj) ) ) 199 199 #else -
trunk/NEMO/OPA_SRC/BDY/bdydyn.F90
r1502 r1559 183 183 WRITE(*,*) 'bdy_dyn_frs: You should not have seen this print! error?', kt 184 184 END SUBROUTINE bdy_dyn_frs 185 SUBROUTINE bdy_dyn_fla ! Empty routine 186 WRITE(*,*) 'bdy_dyn_fla: You should not have seen this print! error?' 185 SUBROUTINE bdy_dyn_fla( pssh ) ! Empty routine 186 REAL :: pssh(:,:) 187 WRITE(*,*) 'bdy_dyn_fla: You should not have seen this print! error?', pssh(1,1) 187 188 END SUBROUTINE bdy_dyn_fla 188 189 #endif -
trunk/NEMO/OPA_SRC/DIA/diaptr.F90
r1413 r1559 4 4 !! Ocean physics: Computes meridonal transports and zonal means 5 5 !!===================================================================== 6 !! History : 9.0 !03-09 (C. Talandier, G. Madec) Original code7 !! 9.0 !06-01 (A. Biastoch) Allow sub-basins computation8 !! 9.0 ! 03-09 (O. Marti) Add fields6 !! History : 1.0 ! 2003-09 (C. Talandier, G. Madec) Original code 7 !! 2.0 ! 2006-01 (A. Biastoch) Allow sub-basins computation 8 !! 3.2 ! 2003-03 (O. Marti, S. Flavoni) Add fields 9 9 !!---------------------------------------------------------------------- 10 10 … … 20 20 USE oce ! ocean dynamics and active tracers 21 21 USE dom_oce ! ocean space and time domain 22 USE daymod ! calandar 23 USE phycst ! physical constants 22 24 USE ldftra_oce ! ocean active tracers: lateral physics 23 USE lib_mpp24 USE in_out_manager25 25 USE dianam 26 USE phycst27 26 USE iom 28 27 USE ioipsl 29 USE daymod 28 USE in_out_manager 29 USE lib_mpp 30 30 31 31 IMPLICIT NONE … … 41 41 PUBLIC ptr_vjk ! call by tra_ldf & tra_adv routines 42 42 43 !!! ** init namelist (namptr) 44 LOGICAL , PUBLIC :: ln_diaptr = .FALSE. !: Poleward transport flag (T) or not (F) 45 LOGICAL , PUBLIC :: ln_subbas = .FALSE. !: Atlantic/Pacific/Indian basins calculation 46 LOGICAL , PUBLIC :: ln_diaznl = .FALSE. !: Add zonal means and meridional stream functions 47 LOGICAL , PUBLIC :: ln_ptrcomp = .FALSE. !: Add decomposition : overturning (and gyre, soon ...) 48 INTEGER , PUBLIC :: nf_ptr = 15 !: frequency of ptr computation 49 INTEGER , PUBLIC :: nf_ptr_wri = 15 !: frequency of ptr outputs 50 51 REAL(wp), DIMENSION(jpi,jpj), PUBLIC :: abasin, pbasin, ibasin, dbasin, sbasin !: Sub basin masks 52 53 REAL(wp), PUBLIC, DIMENSION(jpj) :: pht_adv, pst_adv !: heat and salt poleward transport: advection 54 REAL(wp), PUBLIC, DIMENSION(jpj) :: pht_ove_glo, pst_ove_glo, pht_ove_atl, pst_ove_atl, pht_ove_pac, pst_ove_pac, & 55 & pht_ove_ind, pst_ove_ind, pht_ove_ipc, pst_ove_ipc !: heat and salt poleward transport: overturning 56 REAL(wp), PUBLIC, DIMENSION(jpj) :: pht_ldf, pst_ldf !: heat and salt poleward transport: lateral diffusion 43 ! !!** namelist namptr ** 44 LOGICAL , PUBLIC :: ln_diaptr = .FALSE. !: Poleward transport flag (T) or not (F) 45 LOGICAL , PUBLIC :: ln_subbas = .FALSE. !: Atlantic/Pacific/Indian basins calculation 46 LOGICAL , PUBLIC :: ln_diaznl = .FALSE. !: Add zonal means and meridional stream functions 47 LOGICAL , PUBLIC :: ln_ptrcomp = .FALSE. !: Add decomposition : overturning (and gyre, soon ...) 48 INTEGER , PUBLIC :: nf_ptr = 15 !: frequency of ptr computation 49 INTEGER , PUBLIC :: nf_ptr_wri = 15 !: frequency of ptr outputs 50 51 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: abasin, pbasin, ibasin, dbasin, sbasin !: Sub basin masks 52 53 ! !!! poleward heat and salt transport 54 REAL(wp), PUBLIC, DIMENSION(jpj) :: pht_adv , pst_adv !: advection 55 REAL(wp), PUBLIC, DIMENSION(jpj) :: pht_ldf , pst_ldf !: lateral diffusion 56 REAL(wp), PUBLIC, DIMENSION(jpj) :: pht_ove_glo, pst_ove_glo !: global overturning 57 REAL(wp), PUBLIC, DIMENSION(jpj) :: pht_ove_atl, pst_ove_atl !: Atlantic overturning 58 REAL(wp), PUBLIC, DIMENSION(jpj) :: pht_ove_pac, pst_ove_pac !: Pacific overturning 59 REAL(wp), PUBLIC, DIMENSION(jpj) :: pht_ove_ind, pst_ove_ind !: Indian overturning 60 REAL(wp), PUBLIC, DIMENSION(jpj) :: pht_ove_ipc, pst_ove_ipc !: Indo-Pacific overturning 61 REAL(wp), PUBLIC, DIMENSION(jpj) :: ht_glo, ht_atl, ht_ind, ht_pac, ht_ipc !: heat 62 REAL(wp), PUBLIC, DIMENSION(jpj) :: st_glo, st_atl, st_ind, st_pac, st_ipc !: salt 63 64 INTEGER :: niter 65 INTEGER :: nidom_ptr 66 67 REAL(wp), DIMENSION(jpj,jpk) :: tn_jk_glo , sn_jk_glo ! global i-mean temperature and salinity 68 REAL(wp), DIMENSION(jpj,jpk) :: tn_jk_atl , sn_jk_atl ! Atlantic - - 69 REAL(wp), DIMENSION(jpj,jpk) :: tn_jk_pac , sn_jk_pac ! Pacific - - 70 REAL(wp), DIMENSION(jpj,jpk) :: tn_jk_ind , sn_jk_ind ! Indian - - 71 REAL(wp), DIMENSION(jpj,jpk) :: tn_jk_ipc , sn_jk_ipc ! Indo-Pacific - - 72 REAL(wp), DIMENSION(jpj,jpk) :: v_msf_glo ! global "meridional" Stream-Function 73 REAL(wp), DIMENSION(jpj,jpk) :: v_msf_atl ! Atlantic - - 74 REAL(wp), DIMENSION(jpj,jpk) :: v_msf_pac ! Pacific - - 75 REAL(wp), DIMENSION(jpj,jpk) :: v_msf_ind ! Indian - - 76 REAL(wp), DIMENSION(jpj,jpk) :: v_msf_ipc ! Indo-Pacific - - 77 REAL(wp), DIMENSION(jpj,jpk) :: surf_jk_glo, surf_jk_r_glo ! surface of global i-section and its inverse 78 REAL(wp), DIMENSION(jpj,jpk) :: surf_jk_atl, surf_jk_r_atl ! surface of Atlantic - - 79 REAL(wp), DIMENSION(jpj,jpk) :: surf_jk_pac, surf_jk_r_pac ! surface of Pacific - - 80 REAL(wp), DIMENSION(jpj,jpk) :: surf_jk_ind, surf_jk_r_ind ! surface of Indian - - 81 REAL(wp), DIMENSION(jpj,jpk) :: surf_jk_ipc, surf_jk_r_ipc ! surface of Indo-Pacific - - 57 82 #if defined key_diaeiv 58 REAL(wp), PUBLIC, DIMENSION(jpj) :: pht_eiv_glo, pst_eiv_glo, pht_eiv_atl, pst_eiv_atl, pht_eiv_pac, pst_eiv_pac, & 59 & pht_eiv_ind, pst_eiv_ind, pht_eiv_ipc, pst_eiv_ipc !: heat and salt poleward transport: bolus advection 60 #endif 61 REAL(wp), PUBLIC, DIMENSION(jpj) :: ht_glo,ht_atl,ht_ind,ht_pac,ht_ipc !: heat 62 REAL(wp), PUBLIC, DIMENSION(jpj) :: st_glo,st_atl,st_ind,st_pac,st_ipc !: salt 63 64 INTEGER :: niter 65 INTEGER :: nidom_ptr 66 67 REAL(wp), DIMENSION(jpj,jpk) :: tn_jk_glo, sn_jk_glo, & !: "zonal" mean temperature and salinity 68 & tn_jk_atl, sn_jk_atl, & 69 & tn_jk_pac, sn_jk_pac, & 70 & tn_jk_ind, sn_jk_ind, & 71 & tn_jk_ipc, sn_jk_ipc, & 72 & v_msf_glo , & !: "meridional" Stream-Function 73 & v_msf_atl , & 74 & v_msf_pac , & 75 & v_msf_ind , & 76 & v_msf_ipc , & 77 & surf_jk_glo , & !: Ocean "zonal" section surface 78 & surf_jk_atl , & 79 & surf_jk_pac , & 80 & surf_jk_ind , & 81 & surf_jk_ipc , & 82 & surf_jk_r_glo , & !: inverse of the ocean "zonal" section surface 83 & surf_jk_r_atl , & 84 & surf_jk_r_pac , & 85 & surf_jk_r_ind , & 86 & surf_jk_r_ipc 87 #if defined key_diaeiv 88 REAL(wp), DIMENSION(jpj,jpk) :: v_msf_eiv_glo, v_msf_eiv_atl, v_msf_eiv_pac, & 89 & v_msf_eiv_ind, v_msf_eiv_ipc !: bolus "meridional" Stream-Function 83 ! !!! eddy induced velocity (bolus) 84 REAL(wp), PUBLIC, DIMENSION(jpj) :: pht_eiv_glo, pst_eiv_glo !: global poleward heat and salt bolus advection 85 REAL(wp), PUBLIC, DIMENSION(jpj) :: pht_eiv_atl, pst_eiv_atl !: Atlantic - - 86 REAL(wp), PUBLIC, DIMENSION(jpj) :: pht_eiv_pac, pst_eiv_pac !: Pacific - - 87 REAL(wp), PUBLIC, DIMENSION(jpj) :: pht_eiv_ind, pst_eiv_ind !: Indian - - 88 REAL(wp), PUBLIC, DIMENSION(jpj) :: pht_eiv_ipc, pst_eiv_ipc !: Indo-Pacific - - 89 90 REAL(wp), DIMENSION(jpj,jpk) :: v_msf_eiv_glo ! global "meridional" bolus Stream-Function 91 REAL(wp), DIMENSION(jpj,jpk) :: v_msf_eiv_atl ! Atlantic - - 92 REAL(wp), DIMENSION(jpj,jpk) :: v_msf_eiv_pac ! Pacific - - 93 REAL(wp), DIMENSION(jpj,jpk) :: v_msf_eiv_ind ! Indian - - 94 REAL(wp), DIMENSION(jpj,jpk) :: v_msf_eiv_ipc ! Indo-Pacific - - 90 95 #endif 91 96 … … 94 99 # include "vectopt_loop_substitute.h90" 95 100 !!---------------------------------------------------------------------- 96 !! OPA 9.0 , LOCEAN-IPSL (2005)101 !! NEMO/OPA 3.2 , LOCEAN-IPSL (2009) 97 102 !! $Id$ 98 103 !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) … … 106 111 !! 107 112 !! ** Purpose : "zonal" and vertical sum computation of a "meridional" 108 !! flux array113 !! flux array 109 114 !! 110 115 !! ** Method : - i-k sum of pva using the interior 2D vmask (vmask_i). 111 !! pva is supposed to be a masked flux (i.e. * vmask*e1v*e3v)116 !! pva is supposed to be a masked flux (i.e. * vmask*e1v*e3v) 112 117 !! 113 118 !! ** Action : - p_fval: i-k-mean poleward flux of pva … … 182 187 !! ** Action : - p_fval: i-k-mean poleward flux of pva 183 188 !!---------------------------------------------------------------------- 184 REAL(wp) , INTENT(in), DIMENSION(jpi,jpj,jpk) :: pva! mask flux array at V-point185 REAL(wp) , INTENT(in), DIMENSION(jpi,jpj) , OPTIONAL :: bmask! Optional 2D basin mask189 REAL(wp) , INTENT(in), DIMENSION(jpi,jpj,jpk) :: pva ! mask flux array at V-point 190 REAL(wp) , INTENT(in), DIMENSION(jpi,jpj) , OPTIONAL :: bmask ! Optional 2D basin mask 186 191 !! 187 192 INTEGER :: ji, jj, jk ! dummy loop arguments 188 INTEGER , DIMENSION (1) :: ish189 INTEGER , DIMENSION (2) :: ish2190 REAL(wp), DIMENSION(jpj*jpk) :: zwork ! temporary vector for mpp_sum191 193 REAL(wp), DIMENSION(jpj,jpk) :: p_fval ! return function value 194 #if defined key_mpp_mpi 195 INTEGER, DIMENSION(1) :: ish 196 INTEGER, DIMENSION(2) :: ish2 197 REAL(wp), DIMENSION(jpj*jpk) :: zwork ! 1D workspace 198 #endif 192 199 !!-------------------------------------------------------------------- 193 ! 200 ! 194 201 p_fval(:,:) = 0.e0 195 202 ! 196 IF (PRESENT (bmask)) THEN203 IF( PRESENT( bmask ) ) THEN 197 204 DO jk = 1, jpkm1 198 205 DO jj = 2, jpjm1 206 !!gm here, use of tmask_i ==> no need of loop over nldi, nlei.... 199 207 DO ji = nldi, nlei ! No vector optimisation here. Better use a mask ? 200 208 p_fval(jj,jk) = p_fval(jj,jk) + pva(ji,jj,jk) * e1v(ji,jj) * fse3v(ji,jj,jk) & … … 215 223 ! 216 224 #if defined key_mpp_mpi 217 ish(1) = jpj*jpk ; ish2(1) = jpj ;ish2(2) = jpk218 zwork(:) = RESHAPE( p_fval, ish )225 ish(1) = jpj*jpk ; ish2(1) = jpj ; ish2(2) = jpk 226 zwork(:) = RESHAPE( p_fval, ish ) 219 227 CALL mpp_sum( zwork, jpj*jpk, ncomm_znl ) 220 p_fval(:,:) = RESHAPE( zwork, ish2 )228 p_fval(:,:) = RESHAPE( zwork, ish2 ) 221 229 #endif 222 230 ! 223 231 END FUNCTION ptr_vjk 232 224 233 225 234 FUNCTION ptr_tjk( pta, bmask ) RESULT ( p_fval ) … … 239 248 !! 240 249 INTEGER :: ji, jj, jk ! dummy loop arguments 241 INTEGER, DIMENSION (1) :: ish242 INTEGER, DIMENSION (2) :: ish2243 REAL(wp),DIMENSION(jpj*jpk) :: zwork ! temporary vector for mpp_sum244 250 REAL(wp),DIMENSION(jpj,jpk) :: p_fval ! return function value 251 #if defined key_mpp_mpi 252 INTEGER, DIMENSION(1) :: ish 253 INTEGER, DIMENSION(2) :: ish2 254 REAL(wp),DIMENSION(jpj*jpk) :: zwork ! 1D workspace 255 #endif 245 256 !!-------------------------------------------------------------------- 246 257 ! … … 285 296 INTEGER, INTENT(in) :: kt ! ocean time step index 286 297 !! 287 INTEGER :: jk, jj, ji 288 REAL(wp) :: zsverdrup , &! conversion from m3/s to Sverdrup289 & zpwatt, &! conversion from W to PW290 & zggram! conversion from g to Pg291 REAL(wp), DIMENSION(jpi,jpj,jpk) :: vt, vs298 INTEGER :: jk, jj, ji ! dummy loop 299 REAL(wp) :: zsverdrup ! conversion from m3/s to Sverdrup 300 REAL(wp) :: zpwatt ! conversion from W to PW 301 REAL(wp) :: zggram ! conversion from g to Pg 302 REAL(wp), DIMENSION(jpi,jpj,jpk) :: vt, vs ! 3D workspace 292 303 !!---------------------------------------------------------------------- 293 304 … … 325 336 v_msf_glo(:,:) = ptr_vjk( vn(:,:,:)+v_eiv(:,:,:) ) 326 337 IF( ln_subbas .AND. ln_diaznl ) THEN 327 v_msf_atl(:,:) = ptr_vjk( vn 328 v_msf_pac(:,:) = ptr_vjk( vn 329 v_msf_ind(:,:) = ptr_vjk( vn 330 v_msf_ipc(:,:) = ptr_vjk( vn 338 v_msf_atl(:,:) = ptr_vjk( vn(:,:,:)+v_eiv(:,:,:), abasin(:,:)*sbasin(:,:) ) 339 v_msf_pac(:,:) = ptr_vjk( vn(:,:,:)+v_eiv(:,:,:), pbasin(:,:)*sbasin(:,:) ) 340 v_msf_ind(:,:) = ptr_vjk( vn(:,:,:)+v_eiv(:,:,:), ibasin(:,:)*sbasin(:,:) ) 341 v_msf_ipc(:,:) = ptr_vjk( vn(:,:,:)+v_eiv(:,:,:), dbasin(:,:)*sbasin(:,:) ) 331 342 ENDIF 332 343 #else 333 344 v_msf_glo(:,:) = ptr_vjk( vn(:,:,:) ) 334 345 IF( ln_subbas .AND. ln_diaznl ) THEN 335 v_msf_atl(:,:) = ptr_vjk( vn 336 v_msf_pac(:,:) = ptr_vjk( vn 337 v_msf_ind(:,:) = ptr_vjk( vn 338 v_msf_ipc(:,:) = ptr_vjk( vn 346 v_msf_atl(:,:) = ptr_vjk( vn(:,:,:), abasin(:,:)*sbasin(:,:) ) 347 v_msf_pac(:,:) = ptr_vjk( vn(:,:,:), pbasin(:,:)*sbasin(:,:) ) 348 v_msf_ind(:,:) = ptr_vjk( vn(:,:,:), ibasin(:,:)*sbasin(:,:) ) 349 v_msf_ipc(:,:) = ptr_vjk( vn(:,:,:), dbasin(:,:)*sbasin(:,:) ) 339 350 ENDIF 340 351 #endif … … 474 485 ENDIF 475 486 ENDIF 476 477 ! outputs 478 CALL dia_ptr_wri( kt ) 479 487 ! 488 CALL dia_ptr_wri( kt ) ! outputs 489 ! 480 490 ENDIF 481 482 ! Close the file 483 IF( kt == nitend ) CALL histclo( numptr ) 491 ! 492 IF( kt == nitend ) CALL histclo( numptr ) ! Close the file 484 493 ! 485 494 END SUBROUTINE dia_ptr … … 492 501 !! ** Purpose : Initialization, namelist read 493 502 !!---------------------------------------------------------------------- 503 INTEGER :: inum ! temporary logical unit 504 #if defined key_mpp_mpi 505 INTEGER, DIMENSION(1) :: iglo, iloc, iabsf, iabsl, ihals, ihale, idid 506 #endif 507 !! 494 508 NAMELIST/namptr/ ln_diaptr, ln_diaznl, ln_subbas, ln_ptrcomp, nf_ptr, nf_ptr_wri 495 INTEGER :: inum ! temporary logical unit 496 #if defined key_mpp_mpi 497 INTEGER, DIMENSION (1) :: iglo, iloc, iabsf, iabsl, ihals, ihale, idid 498 #endif 499 !!---------------------------------------------------------------------- 500 501 ! Read Namelist namptr : poleward transport parameters 502 REWIND ( numnam ) 509 !!---------------------------------------------------------------------- 510 511 REWIND ( numnam ) ! Read Namelist namptr : poleward transport parameters 503 512 READ ( numnam, namptr ) 504 513 505 ! Control print 506 IF(lwp) THEN 514 IF(lwp) THEN ! Control print 507 515 WRITE(numout,*) 508 516 WRITE(numout,*) 'dia_ptr_init : poleward transport and msf initialization' 509 517 WRITE(numout,*) '~~~~~~~~~~~~' 510 WRITE(numout,*) ' 511 WRITE(numout,*) ' Switch for ptr diagnostic (T) or not (F) ln_diaptr= ', ln_diaptr512 WRITE(numout,*) ' Atla/Paci/Ind basins computation ln_subbas= ', ln_subbas513 WRITE(numout,*) ' Frequency of computation nf_ptr= ', nf_ptr514 WRITE(numout,*) ' Frequency of outputsnf_ptr_wri = ', nf_ptr_wri518 WRITE(numout,*) ' Namelist namptr : set ptr parameters' 519 WRITE(numout,*) ' Switch for ptr diagnostic (T) or not (F) ln_diaptr = ', ln_diaptr 520 WRITE(numout,*) ' Atl/Pac/Ind basins computation ln_subbas = ', ln_subbas 521 WRITE(numout,*) ' Frequency of computation nf_ptr = ', nf_ptr 522 WRITE(numout,*) ' Frequency of outputs nf_ptr_wri = ', nf_ptr_wri 515 523 ENDIF 516 524 517 ! 518 ! Define MPI communicator for zonal sum 519 ! 520 IF( lk_mpp ) THEN 521 CALL mpp_ini_znl 522 ENDIF 523 524 IF( ln_subbas ) THEN ! load sub-basin mask 525 IF( lk_mpp ) CALL mpp_ini_znl ! Define MPI communicator for zonal sum 526 527 IF( ln_subbas ) THEN ! load sub-basin mask 525 528 CALL iom_open( 'subbasins', inum ) 526 529 CALL iom_get( inum, jpdom_data, 'atlmsk', abasin ) ! Atlantic basin … … 533 536 ENDIF 534 537 538 !!gm CAUTION : this is only valid in fixed volume case ! 539 535 540 ! inverse of the ocean "zonal" v-point section 536 541 surf_jk_glo(:,:) = ptr_tjk( tmask(:,:,:) ) … … 580 585 nidom_ptr = FLIO_DOM_NONE 581 586 #endif 582 587 ! 583 588 END SUBROUTINE dia_ptr_init 584 589 … … 632 637 633 638 zdt = rdt 634 IF( nacc == 1 ) zdt = rdtmin639 IF( nacc == 1 ) zdt = rdtmin 635 640 636 641 ! Reference latitude … … 670 675 671 676 ! ! ======================= 672 ELSE ! OTHER configurations zjulian = zjulian - adatrj 673 ! set calendar origin to the beginning of the experiment 677 ELSE ! OTHER configurations 674 678 ! ! ======================= 675 679 zphi(:) = gphiv(1,:) ! assume lat/lon coordinate, select the first i-line -
trunk/NEMO/OPA_SRC/SBC/sbcana.F90
r1230 r1559 4 4 !! Ocean forcing: analytical momentum, heat and freshwater forcings 5 5 !!===================================================================== 6 !! History : 9.0 ! 06-06 (G. Madec) Original code 6 !! History : 3.0 ! 2006-06 (G. Madec) Original code 7 !! 3.2 ! 2009-07 (G. Madec) Style only 7 8 !!---------------------------------------------------------------------- 8 9 … … 26 27 PUBLIC sbc_gyre ! routine called in sbcmod module 27 28 28 ! ! * Namelist namsbc_ana29 ! !!* Namelist namsbc_ana * 29 30 INTEGER :: nn_tau000 = 1 ! nb of time-step during which the surface stress 30 !! increase from 0 to its nominal value31 ! ! increase from 0 to its nominal value 31 32 REAL(wp) :: rn_utau0 = 0.e0 ! constant wind stress value in i-direction 32 33 REAL(wp) :: rn_vtau0 = 0.e0 ! constant wind stress value in j-direction … … 39 40 # include "vectopt_loop_substitute.h90" 40 41 !!---------------------------------------------------------------------- 41 !! OPA 9.0 , LOCEAN-IPSL (2006)42 !! NEMO/OPA 3.2 , LOCEAN-IPSL (2009) 42 43 !! $Id$ 43 44 !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) … … 51 52 !! 52 53 !! ** Purpose : provide at each time-step the ocean surface boundary 53 !! condition, i.e. the momentum, heat and freshwater fluxes.54 !! condition, i.e. the momentum, heat and freshwater fluxes. 54 55 !! 55 56 !! ** Method : Constant and uniform surface forcing specified from 56 !! namsbc_ana namelist parameters. All the fluxes are time inde-57 !! pendant except the stresses which increase from zero during58 !! the first nn_tau000 time-step59 !! * C A U T I ON : never mask the surface stress field !57 !! namsbc_ana namelist parameters. All the fluxes are time 58 !! independant except the stresses which increase from zero 59 !! during the first nn_tau000 time-step 60 !! CAUTION : never mask the surface stress field ! 60 61 !! 61 62 !! ** Action : - set the ocean surface boundary condition, i.e. … … 64 65 INTEGER, INTENT(in) :: kt ! ocean time step 65 66 !! 66 INTEGER :: ji, jj ! dummy loop indices67 67 REAL(wp) :: zfacto ! local scalar 68 68 !! … … 100 100 ENDIF 101 101 102 ! Estimation of wind speed as a function of wind stress ( |tau|=rhoa*Cd*|U|^2 ) 103 CALL sbc_tau2wnd 102 CALL sbc_tau2wnd ! Estimation of wind speed as a function of wind stress ( |tau|=rhoa*Cd*|U|^2 ) 104 103 ! 105 104 END SUBROUTINE sbc_ana … … 110 109 !! *** ROUTINE sbc_ana *** 111 110 !! 112 !! ** Purpose : provide at each time-step the oceansurface boundary113 !! condition, i.e. the momentum, heat and freshwater fluxes.111 !! ** Purpose : provide at each time-step the GYRE surface boundary 112 !! condition, i.e. the momentum, heat and freshwater fluxes. 114 113 !! 115 114 !! ** Method : analytical seasonal cycle for GYRE configuration. 116 !! * C A U T I ON : never mask the surface stress field !115 !! CAUTION : never mask the surface stress field ! 117 116 !! 118 117 !! ** Action : - set the ocean surface boundary condition, i.e. … … 122 121 !!---------------------------------------------------------------------- 123 122 INTEGER, INTENT(in) :: kt ! ocean time step 124 123 !! 125 124 INTEGER :: ji, jj ! dummy loop indices 126 125 INTEGER :: zyear0 ! initial year … … 298 297 WRITE(numout,*)' raajj = ', raajj 299 298 ENDIF 300 299 ! 301 300 END SUBROUTINE sbc_gyre 302 301 -
trunk/NEMO/OPA_SRC/TRA/traadv_cen2.F90
r1537 r1559 4 4 !! Ocean active tracers: horizontal & vertical advective trend 5 5 !!====================================================================== 6 !! History : 8.2 ! 01-08 (G. Madec, E. Durand) trahad+trazad=traadv 7 !! 8.5 ! 02-06 (G. Madec) F90: Free form and module 8 !! 9.0 ! 04-08 (C. Talandier) New trends organization 9 !! " " ! 05-11 (V. Garnier) Surface pressure gradient organization 10 !! " " ! 06-04 (R. Benshila, G. Madec) Step reorganization 11 !! " " ! 06-07 (G. madec) add ups_orca_set routine 6 !! History : 8.2 ! 2001-08 (G. Madec, E. Durand) trahad+trazad=traadv 7 !! 1.0 ! 2002-06 (G. Madec) F90: Free form and module 8 !! 9.0 ! 2004-08 (C. Talandier) New trends organization 9 !! - ! 2005-11 (V. Garnier) Surface pressure gradient organization 10 !! 2.0 ! 2006-04 (R. Benshila, G. Madec) Step reorganization 11 !! - ! 2006-07 (G. madec) add ups_orca_set routine 12 !! 3.2 ! 2009-07 (G. Madec) add avmb, avtb in restart for cen2 advection 12 13 !!---------------------------------------------------------------------- 13 14 … … 53 54 # include "vectopt_loop_substitute.h90" 54 55 !!---------------------------------------------------------------------- 55 !! OPA 9.0 , LOCEAN-IPSL (2006)56 !! NEMO/OPA 3.2 , LOCEAN-IPSL (2009) 56 57 !! $Id$ 57 58 !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) … … 131 132 REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk) :: pwn ! ocean velocity w-component 132 133 !! 133 INTEGER :: ji, jj, jk 134 REAL(wp) :: z ta, zsa, zbtr, zhw, ze3tr, &! temporary scalars135 & zfp_ui, zfp_vj, zfp_w , zfui , & ! " "136 & zfm_ui, zfm_vj, zfm_w , zfvj , & ! " "137 & zcofi , zcofj , zcofk , & ! " "138 & zupsut, zupsus, zcenut, zcenus, & ! " "139 & zupsvt, zupsvs, zcenvt, zcenvs, & ! " "140 & zupst , zupss , zcent , zcens , & ! " "141 & z_hdivn_x, z_hdivn_y, z_hdivn142 REAL(wp) :: zice 134 INTEGER :: ji, jj, jk ! dummy loop indices 135 REAL(wp) :: zbtr, zhw, ze3tr ! temporary scalars 136 REAL(wp) :: zfp_ui, zfp_vj, zfp_w , zfui ! - - 137 REAL(wp) :: zfm_ui, zfm_vj, zfm_w , zfvj ! - - 138 REAL(wp) :: zcofi , zcofj , zcofk ! - - 139 REAL(wp) :: zupsut, zupsus, zcenut, zcenus ! - - 140 REAL(wp) :: zupsvt, zupsvs, zcenvt, zcenvs ! - - 141 REAL(wp) :: zupst , zupss , zcent , zcens ! - - 142 REAL(wp) :: z_hdivn_x, z_hdivn_y, z_hdivn ! - - 143 REAL(wp) :: zice ! - - 143 144 REAL(wp), DIMENSION(jpi,jpj) :: ztfreez ! 2D workspace 144 145 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zwz, ztrdt, zind ! 3D workspace … … 188 189 END DO 189 190 190 ! I. Horizontal advective fluxes 191 ! ------------------------------ 192 ! Second order centered tracer flux at u and v-points 193 ! ----------------------------------------------------- 194 ! ! =============== 195 DO jk = 1, jpkm1 ! Horizontal slab 196 ! ! =============== 191 ! I. Horizontal advection 192 ! ==================== 193 ! 194 DO jk = 1, jpkm1 195 ! ! Second order centered tracer flux at u- and v-points 197 196 DO jj = 1, jpjm1 198 197 DO ji = 1, fs_jpim1 ! vector opt. … … 201 200 zcofj = MAX( zind(ji,jj+1,jk), zind(ji,jj,jk) ) 202 201 ! volume fluxes * 1/2 203 #if defined key_zco204 zfui = 0.5 * e2u(ji,jj) * pun(ji,jj,jk)205 zfvj = 0.5 * e1v(ji,jj) * pvn(ji,jj,jk)206 #else207 202 zfui = 0.5 * e2u(ji,jj) * fse3u(ji,jj,jk) * pun(ji,jj,jk) 208 203 zfvj = 0.5 * e1v(ji,jj) * fse3v(ji,jj,jk) * pvn(ji,jj,jk) 209 #endif 204 ! 210 205 ! upstream scheme 211 206 zfp_ui = zfui + ABS( zfui ) … … 229 224 END DO 230 225 END DO 231 232 ! Tracer flux divergence at t-point added to the general trend 233 ! -------------------------------------------------------------- 226 ! ! Tracer flux divergence at t-point added to the general trend 234 227 DO jj = 2, jpjm1 235 228 DO ji = fs_2, fs_jpim1 ! vector opt. 236 #if defined key_zco237 zbtr = btr2(ji,jj)238 #else239 229 zbtr = btr2(ji,jj) / fse3t(ji,jj,jk) 240 #endif 241 ! horizontal advective trends 242 zta = - zbtr * ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk) & 243 & + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk) ) 244 zsa = - zbtr * ( zww(ji,jj,jk) - zww(ji-1,jj ,jk) & 245 & + zwz(ji,jj,jk) - zwz(ji ,jj-1,jk) ) 246 ! add it to the general tracer trends 247 ta(ji,jj,jk) = ta(ji,jj,jk) + zta 248 sa(ji,jj,jk) = sa(ji,jj,jk) + zsa 249 END DO 250 END DO 251 ! ! =============== 252 END DO ! End of slab 253 ! ! =============== 254 255 ! Save the horizontal advective trends for diagnostic 256 ! ----------------------------------------------------- 257 IF( l_trdtra ) THEN 258 ! T/S ZONAL advection trends 259 ztrdt(:,:,:) = 0.e0 ; ztrds(:,:,:) = 0.e0 230 ! 231 ta(ji,jj,jk) = ta(ji,jj,jk) - zbtr * ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk) & 232 & + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk) ) 233 sa(ji,jj,jk) = sa(ji,jj,jk) - zbtr * ( zww(ji,jj,jk) - zww(ji-1,jj ,jk) & 234 & + zwz(ji,jj,jk) - zwz(ji ,jj-1,jk) ) 235 END DO 236 END DO 237 END DO 238 239 240 IF( l_trdtra ) THEN ! Save the i- and j-advective trends for diagnostic (U.gradz(T) trends) 260 241 ! 261 242 DO jk = 1, jpkm1 … … 263 244 DO ji = fs_2, fs_jpim1 ! vector opt. 264 245 !-- Compute zonal divergence by splitting hdivn (see divcur.F90) 265 ! N.B. This computation is not valid along OBCs (if any) 266 #if defined key_zco 267 zbtr = btr2(ji,jj) 268 z_hdivn_x = ( e2u(ji ,jj) * pun(ji ,jj,jk) & 269 & - e2u(ji-1,jj) * pun(ji-1,jj,jk) ) * zbtr 270 #else 246 ! N.B. This computation is not valid with OBC, BDY, cla, eiv, advective bbl 271 247 zbtr = btr2(ji,jj) / fse3t(ji,jj,jk) 272 248 z_hdivn_x = ( e2u(ji ,jj) * fse3u(ji ,jj,jk) * pun(ji ,jj,jk) & 273 249 & - e2u(ji-1,jj) * fse3u(ji-1,jj,jk) * pun(ji-1,jj,jk) ) * zbtr 274 #endif 250 ! 275 251 ztrdt(ji,jj,jk) = - zbtr * ( zwx(ji,jj,jk) - zwx(ji-1,jj,jk) ) + tn(ji,jj,jk) * z_hdivn_x 276 252 ztrds(ji,jj,jk) = - zbtr * ( zww(ji,jj,jk) - zww(ji-1,jj,jk) ) + sn(ji,jj,jk) * z_hdivn_x … … 280 256 CALL trd_mod(ztrdt, ztrds, jptra_trd_xad, 'TRA', kt) 281 257 ! 282 ! T/S MERIDIONAL advection trends 258 DO jk = 1, jpkm1 ! T/S MERIDIONAL advection trends 259 DO jj = 2, jpjm1 260 DO ji = fs_2, fs_jpim1 ! vector opt. 261 zbtr = btr2(ji,jj) / fse3t(ji,jj,jk) 262 z_hdivn_y = ( e1v(ji, jj) * fse3v(ji,jj ,jk) * pvn(ji,jj ,jk) & 263 & - e1v(ji,jj-1) * fse3v(ji,jj-1,jk) * pvn(ji,jj-1,jk) ) * zbtr 264 ! 265 ztrdt(ji,jj,jk) = - zbtr * ( zwy(ji,jj,jk) - zwy(ji,jj-1,jk) ) + tn(ji,jj,jk) * z_hdivn_y 266 ztrds(ji,jj,jk) = - zbtr * ( zwz(ji,jj,jk) - zwz(ji,jj-1,jk) ) + sn(ji,jj,jk) * z_hdivn_y 267 END DO 268 END DO 269 END DO 270 CALL trd_mod(ztrdt, ztrds, jptra_trd_yad, 'TRA', kt) 271 ! 272 ztrdt(:,:,:) = ta(:,:,:) ; ztrds(:,:,:) = sa(:,:,:) ! Save the horizontal up-to-date ta/sa trends 273 ! 274 ENDIF 275 276 IF( ln_diaptr .AND. ( MOD( kt, nf_ptr ) == 0 ) ) THEN ! "zonal" mean advective heat and salt transport 277 pht_adv(:) = ptr_vj( zwy(:,:,:) ) 278 pst_adv(:) = ptr_vj( zwz(:,:,:) ) 279 ENDIF 280 281 IF(ln_ctl) CALL prt_ctl( tab3d_1=ta, clinfo1=' cen2 had - Ta: ', mask1=tmask, & 282 & tab3d_2=sa, clinfo2= ' Sa: ', mask2=tmask, clinfo3='tra' ) 283 284 285 ! II. Vertical advection 286 ! ================== 287 ! 288 zwx(:,:,jpk) = 0.e0 ; zwy(:,:,jpk) = 0.e0 ! Bottom value : flux set to zero 289 ! 290 IF( lk_vvl ) THEN ! Surface value : zero in variable volume 291 zwx(:,:, 1 ) = 0.e0 ; zwy(:,:, 1 ) = 0.e0 292 ELSE ! : linear free surface case 293 zwx(:,:, 1 ) = pwn(:,:,1) * tn(:,:,1) 294 zwy(:,:, 1 ) = pwn(:,:,1) * sn(:,:,1) 295 ENDIF 296 ! 297 DO jk = 2, jpk ! Second order centered tracer flux at w-point 298 DO jj = 2, jpjm1 299 DO ji = fs_2, fs_jpim1 ! vector opt. 300 zcofk = MAX( zind(ji,jj,jk-1), zind(ji,jj,jk) ) ! upstream indicator 301 zhw = 0.5 * pwn(ji,jj,jk) ! velocity * 1/2 302 ! 303 zfp_w = zhw + ABS( zhw ) ! upstream scheme 304 zfm_w = zhw - ABS( zhw ) 305 zupst = zfp_w * tb(ji,jj,jk) + zfm_w * tb(ji,jj,jk-1) 306 zupss = zfp_w * sb(ji,jj,jk) + zfm_w * sb(ji,jj,jk-1) 307 ! 308 zcent = zhw * ( tn(ji,jj,jk) + tn(ji,jj,jk-1) ) ! centered scheme 309 zcens = zhw * ( sn(ji,jj,jk) + sn(ji,jj,jk-1) ) 310 ! 311 zwx(ji,jj,jk) = zcofk * zupst + (1.-zcofk) * zcent ! mixed centered / upstream scheme 312 zwy(ji,jj,jk) = zcofk * zupss + (1.-zcofk) * zcens 313 END DO 314 END DO 315 END DO 316 ! 317 DO jk = 1, jpkm1 ! divergence of Tracer flux added to the general trend 318 DO jj = 2, jpjm1 319 DO ji = fs_2, fs_jpim1 ! vector opt. 320 ze3tr = 1. / fse3t(ji,jj,jk) 321 ta(ji,jj,jk) = ta(ji,jj,jk) - ze3tr * ( zwx(ji,jj,jk) - zwx(ji,jj,jk+1) ) 322 sa(ji,jj,jk) = sa(ji,jj,jk) - ze3tr * ( zwy(ji,jj,jk) - zwy(ji,jj,jk+1) ) 323 END DO 324 END DO 325 END DO 326 327 IF( l_trdtra ) THEN ! Save the vertical advective trends for diagnostic (W gradz(T) trends) 283 328 DO jk = 1, jpkm1 284 329 DO jj = 2, jpjm1 285 330 DO ji = fs_2, fs_jpim1 ! vector opt. 286 !-- Compute merid. divergence by splitting hdivn (see divcur.F90)287 ! N.B. This computation is not valid along OBCs (if any)288 #if defined key_zco289 zbtr = btr2(ji,jj)290 z_hdivn_y = ( e1v(ji,jj ) * pvn(ji,jj ,jk) &291 & - e1v(ji,jj-1) * pvn(ji,jj-1,jk) ) * zbtr292 #else293 zbtr = btr2(ji,jj) / fse3t(ji,jj,jk)294 z_hdivn_y = ( e1v(ji, jj) * fse3v(ji,jj ,jk) * pvn(ji,jj ,jk) &295 & - e1v(ji,jj-1) * fse3v(ji,jj-1,jk) * pvn(ji,jj-1,jk) ) * zbtr296 #endif297 ztrdt(ji,jj,jk) = - zbtr * ( zwy(ji,jj,jk) - zwy(ji,jj-1,jk) ) + tn(ji,jj,jk) * z_hdivn_y298 ztrds(ji,jj,jk) = - zbtr * ( zwz(ji,jj,jk) - zwz(ji,jj-1,jk) ) + sn(ji,jj,jk) * z_hdivn_y299 END DO300 END DO301 END DO302 CALL trd_mod(ztrdt, ztrds, jptra_trd_yad, 'TRA', kt)303 !304 ! Save the horizontal up-to-date ta/sa trends305 ztrdt(:,:,:) = ta(:,:,:)306 ztrds(:,:,:) = sa(:,:,:)307 ENDIF308 309 IF(ln_ctl) CALL prt_ctl( tab3d_1=ta, clinfo1=' cen2 had - Ta: ', mask1=tmask, &310 & tab3d_2=sa, clinfo2= ' Sa: ', mask2=tmask, clinfo3='tra' )311 312 ! "zonal" mean advective heat and salt transport313 ! ----------------------------------------------314 315 IF( ln_diaptr .AND. ( MOD( kt, nf_ptr ) == 0 ) ) THEN316 IF( lk_zco ) THEN317 DO jk = 1, jpkm1318 DO jj = 2, jpjm1319 DO ji = fs_2, fs_jpim1 ! vector opt.320 zwy(ji,jj,jk) = zwy(ji,jj,jk) * fse3v(ji,jj,jk)321 zwz(ji,jj,jk) = zwz(ji,jj,jk) * fse3v(ji,jj,jk)322 END DO323 END DO324 END DO325 ENDIF326 pht_adv(:) = ptr_vj( zwy(:,:,:) )327 pst_adv(:) = ptr_vj( zwz(:,:,:) )328 ENDIF329 330 ! II. Vertical advection331 ! ----------------------332 333 ! Bottom value : flux set to zero334 zwx(:,:,jpk) = 0.e0 ; zwy(:,:,jpk) = 0.e0335 336 ! Surface value337 IF( lk_vvl ) THEN338 ! variable volume: flux set to zero339 zwx(:,:, 1 ) = 0.e0 ; zwy(:,:, 1 ) = 0.e0340 ELSE341 ! free surface-constant volume342 zwx(:,:, 1 ) = pwn(:,:,1) * tn(:,:,1)343 zwy(:,:, 1 ) = pwn(:,:,1) * sn(:,:,1)344 ENDIF345 346 ! 1. Vertical advective fluxes347 ! ----------------------------348 ! Second order centered tracer flux at w-point349 DO jk = 2, jpk350 DO jj = 2, jpjm1351 DO ji = fs_2, fs_jpim1 ! vector opt.352 ! upstream indicator353 zcofk = MAX( zind(ji,jj,jk-1), zind(ji,jj,jk) )354 ! velocity * 1/2355 zhw = 0.5 * pwn(ji,jj,jk)356 ! upstream scheme357 zfp_w = zhw + ABS( zhw )358 zfm_w = zhw - ABS( zhw )359 zupst = zfp_w * tb(ji,jj,jk) + zfm_w * tb(ji,jj,jk-1)360 zupss = zfp_w * sb(ji,jj,jk) + zfm_w * sb(ji,jj,jk-1)361 ! centered scheme362 zcent = zhw * ( tn(ji,jj,jk) + tn(ji,jj,jk-1) )363 zcens = zhw * ( sn(ji,jj,jk) + sn(ji,jj,jk-1) )364 ! mixed centered / upstream scheme365 zwx(ji,jj,jk) = zcofk * zupst + (1.-zcofk) * zcent366 zwy(ji,jj,jk) = zcofk * zupss + (1.-zcofk) * zcens367 END DO368 END DO369 END DO370 371 ! 2. Tracer flux divergence at t-point added to the general trend372 ! -------------------------373 DO jk = 1, jpkm1374 DO jj = 2, jpjm1375 DO ji = fs_2, fs_jpim1 ! vector opt.376 ze3tr = 1. / fse3t(ji,jj,jk)377 ! vertical advective trends378 zta = - ze3tr * ( zwx(ji,jj,jk) - zwx(ji,jj,jk+1) )379 zsa = - ze3tr * ( zwy(ji,jj,jk) - zwy(ji,jj,jk+1) )380 ! add it to the general tracer trends381 ta(ji,jj,jk) = ta(ji,jj,jk) + zta382 sa(ji,jj,jk) = sa(ji,jj,jk) + zsa383 END DO384 END DO385 END DO386 387 ! 3. Save the vertical advective trends for diagnostic388 ! ----------------------------------------------------389 IF( l_trdtra ) THEN390 ! Recompute the vertical advection zta & zsa trends computed391 ! at the step 2. above in making the difference between the new392 ! trends and the previous one: ta()/sa - ztrdt()/ztrds() and substract393 ! the term tn()/sn()*hdivn() to recover the W gradz(T/S) trends394 395 DO jk = 1, jpkm1396 DO jj = 2, jpjm1397 DO ji = fs_2, fs_jpim1 ! vector opt.398 #if defined key_zco399 zbtr = btr2(ji,jj)400 z_hdivn_x = e2u(ji,jj)*pun(ji,jj,jk) - e2u(ji-1,jj)*pun(ji-1,jj,jk)401 z_hdivn_y = e1v(ji,jj)*pvn(ji,jj,jk) - e1v(ji,jj-1)*pvn(ji,jj-1,jk)402 #else403 331 zbtr = btr2(ji,jj) / fse3t(ji,jj,jk) 404 332 z_hdivn_x = e2u(ji,jj)*fse3u(ji,jj,jk)*pun(ji,jj,jk) - e2u(ji-1,jj)*fse3u(ji-1,jj,jk)*pun(ji-1,jj,jk) 405 333 z_hdivn_y = e1v(ji,jj)*fse3v(ji,jj,jk)*pvn(ji,jj,jk) - e1v(ji,jj-1)*fse3v(ji,jj-1,jk)*pvn(ji,jj-1,jk) 406 #endif 334 ! 407 335 z_hdivn = (z_hdivn_x + z_hdivn_y) * zbtr 408 336 ztrdt(ji,jj,jk) = ta(ji,jj,jk) - ztrdt(ji,jj,jk) - tn(ji,jj,jk) * z_hdivn -
trunk/NEMO/OPA_SRC/TRA/traadv_qck.F90
r1528 r1559 4 4 !! Ocean active tracers: horizontal & vertical advective trend 5 5 !!============================================================================== 6 !! History : 9.0 ! 2008-07 (G. Reffray) Original code6 !! History : 3.0 ! 2008-07 (G. Reffray) Original code 7 7 !!---------------------------------------------------------------------- 8 9 8 10 9 !!---------------------------------------------------------------------- 11 10 !! tra_adv_qck : update the tracer trend with the horizontal advection 12 11 !! trends using a 3rd order finite difference scheme 13 !! tra_adv_qck_ zon:14 !! tra_adv_qck_ mer:15 !! tra_adv_cen2_ ver: 2nd centered scheme for the vertical advection12 !! tra_adv_qck_i : 13 !! tra_adv_qck_j : 14 !! tra_adv_cen2_k : 2nd centered scheme for the vertical advection 16 15 !!---------------------------------------------------------------------- 17 !! * Modules used18 16 USE oce ! ocean dynamics and active tracers 19 17 USE dom_oce ! ocean space and time domain … … 31 29 PRIVATE 32 30 33 !! * Accessibility 34 PUBLIC tra_adv_qck ! routine called by step.F90 35 36 !! * Module variables 37 REAL(wp), DIMENSION(jpi,jpj), SAVE :: btr2 38 REAL(wp) , SAVE :: cst 39 !!---------------------------------------------------------------------- 31 PUBLIC tra_adv_qck ! routine called by step.F90 32 33 REAL(wp), DIMENSION(jpi,jpj) :: btr2 34 REAL(wp) :: r1_6 35 40 36 !! * Substitutions 41 37 # include "domzgr_substitute.h90" 42 38 # include "vectopt_loop_substitute.h90" 43 39 !!---------------------------------------------------------------------- 44 !! OPA 9.0 , LOCEAN-IPSL (2008)40 !! NEMO/OPA 3.2 , LOCEAN-IPSL (2009) 45 41 !! $Id$ 46 42 !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) … … 57 53 !! 58 54 !! ** Method : The advection is evaluated by a third order scheme 59 !! For a positive velocity u : 60 !! 61 !! 62 !! i-1 i i+1 i+2 63 !! 64 !! |--FU--|--FC--|--FD--|------| 65 !! u(i)>0 66 !! 67 !! For a negative velocity u : 68 !! 69 !! |------|--FD--|--FC--|--FU--| 70 !! u(i)<0 71 !! 72 !! FU is the second upwind point 73 !! FD is the first douwning point 74 !! FC is the central point (or the first upwind point) 75 !! 76 !! Flux(i) = u(i) * {0.5(FC+FD) -0.5C(i)(FD-FC) -((1-C(i))/6)(FU+FD-2FC)} 77 !! with C(i)=|u(i)|dx(i)/dt (Courant number) 55 !! For a positive velocity u : u(i)>0 56 !! |--FU--|--FC--|--FD--|------| 57 !! i-1 i i+1 i+2 58 !! 59 !! For a negative velocity u : u(i)<0 60 !! |------|--FD--|--FC--|--FU--| 61 !! i-1 i i+1 i+2 62 !! where FU is the second upwind point 63 !! FD is the first douwning point 64 !! FC is the central point (or the first upwind point) 65 !! 66 !! Flux(i) = u(i) * { 0.5(FC+FD) -0.5C(i)(FD-FC) -((1-C(i))/6)(FU+FD-2FC) } 67 !! with C(i)=|u(i)|dx(i)/dt (=Courant number) 78 68 !! 79 69 !! dt = 2*rdtra and the scalar values are tb and sb … … 81 71 !! On the vertical, the simple centered scheme used tn and sn 82 72 !! 83 !! The fluxes are bounded by the ULTIMATE limiter 84 !! toguarantee the monotonicity of the solution and to73 !! The fluxes are bounded by the ULTIMATE limiter to 74 !! guarantee the monotonicity of the solution and to 85 75 !! prevent the appearance of spurious numerical oscillations 86 76 !! … … 92 82 USE oce, ONLY : ztrdt => ua ! use ua as workspace 93 83 USE oce, ONLY : ztrds => va ! use va as workspace 94 !! * Arguments84 !! 95 85 INTEGER , INTENT(in) :: kt ! ocean time-step index 96 86 REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk) :: pun ! effective ocean velocity, u_component … … 109 99 IF(lwp) WRITE(numout,*) 110 100 btr2(:,:) = 1. / ( e1t(:,:) * e2t(:,:) ) 111 cst= 1. / 6.101 r1_6 = 1. / 6. 112 102 ENDIF 113 103 … … 119 109 !--------------------------------------------------------------------------- 120 110 121 CALL tra_adv_qck_ zon( kt,pun, tb, tn, ta, ztrdt, z2)122 CALL tra_adv_qck_ zon( kt,pun, sb, sn, sa, ztrds, z2)111 CALL tra_adv_qck_i( pun, tb, tn, ta, ztrdt, z2) 112 CALL tra_adv_qck_i( pun, sb, sn, sa, ztrds, z2) 123 113 124 114 IF( l_trdtra ) CALL trd_mod(ztrdt, ztrds, jptra_trd_xad, 'TRA', kt) 125 115 126 CALL tra_adv_qck_ mer( kt, pvn, tb, tn, ta, ztrdt, pht_adv, z2)127 CALL tra_adv_qck_ mer( kt, pvn, sb, sn, sa, ztrds, pst_adv, z2)116 CALL tra_adv_qck_j( kt, pvn, tb, tn, ta, ztrdt, pht_adv, z2) 117 CALL tra_adv_qck_j( kt, pvn, sb, sn, sa, ztrds, pst_adv, z2) 128 118 129 119 IF( l_trdtra ) THEN 130 120 CALL trd_mod(ztrdt, ztrds, jptra_trd_yad, 'TRA', kt) 131 121 ! 132 ! Save the horizontal up-to-date ta/sa trends 133 ztrdt(:,:,:) = ta(:,:,:) 122 ztrdt(:,:,:) = ta(:,:,:) ! Save the horizontal up-to-date ta/sa trends 134 123 ztrds(:,:,:) = sa(:,:,:) 135 124 END IF … … 141 130 !------------------------------------------------------------------------- 142 131 ! 143 CALL tra_adv_cen2_ ver( pwn, tn, ta )144 CALL tra_adv_cen2_ ver( pwn, sn, sa )132 CALL tra_adv_cen2_k( pwn, tn, ta ) 133 CALL tra_adv_cen2_k( pwn, sn, sa ) 145 134 ! 146 135 !Save the vertical advective trends for diagnostic … … 179 168 180 169 181 SUBROUTINE tra_adv_qck_zon ( kt, pun, tra, tran, traa, ztrdtra, z2 ) 182 !!---------------------------------------------------------------------- 183 !! 184 !!---------------------------------------------------------------------- 185 !! * Arguments 186 INTEGER, INTENT(in) :: kt ! ocean time-step index 170 SUBROUTINE tra_adv_qck_i ( pun, tra, tran, traa, ztrdtra, z2 ) 171 !!---------------------------------------------------------------------- 172 !! 173 !!---------------------------------------------------------------------- 187 174 REAL, INTENT(in) :: z2 188 175 REAL(wp), INTENT(in) , DIMENSION(jpi,jpj,jpk) :: pun, tra, tran ! horizontal effective velocity … … 196 183 REAL(wp), DIMENSION(jpi,jpj) :: zfu, zfc, zfd, zfho, zmskl, zsc_e 197 184 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zflux 198 !----------------------------------------------------------------------199 ! I. Part 1 : x-direction200 185 !---------------------------------------------------------------------- 201 186 … … 216 201 ! Mask at the T-points in the x-direction (mask=0 or mask=1) 217 202 zmask(ji,jj)=tmask(ji-1,jj,jk)+tmask(ji,jj,jk)+tmask(ji+1,jj,jk)-2 218 END DO219 END DO203 END DO 204 END DO 220 205 ! 221 206 !--- Lateral boundary conditions … … 244 229 zfd(ji,jj) = dir*tra (ji+1,jj,jk)+(1-dir)*tra (ji ,jj,jk) ! FD in the x-direction for T 245 230 zmskl(ji,jj) = dir*zmask(ji ,jj) +(1-dir)*zmask(ji+1,jj) 246 END DO247 END DO231 END DO 232 END DO 248 233 ! 249 234 !--- QUICKEST scheme … … 300 285 END IF 301 286 302 END SUBROUTINE tra_adv_qck_zon 303 304 305 SUBROUTINE tra_adv_qck_mer ( kt, pvn, tra, tran, traa, ztrdtra, trd_adv, z2 ) 306 !!---------------------------------------------------------------------- 307 !! 308 !!---------------------------------------------------------------------- 309 !! * Arguments 287 END SUBROUTINE tra_adv_qck_i 288 289 290 SUBROUTINE tra_adv_qck_j ( kt, pvn, tra, tran, traa, ztrdtra, trd_adv, z2 ) 291 !!---------------------------------------------------------------------- 292 !! 293 !!---------------------------------------------------------------------- 310 294 INTEGER, INTENT(in) :: kt ! ocean time-step index 311 295 REAL, INTENT(in) :: z2 … … 314 298 REAL(wp), INTENT(out) , DIMENSION(jpi,jpj,jpk) :: ztrdtra 315 299 REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: traa 316 ! 300 !! 317 301 INTEGER :: ji, jj, jk 318 302 REAL(wp) :: za, zbtr, dir, dx, dt ! temporary scalars … … 342 326 ! Mask at the T-points in the x-direction (mask=0 or mask=1) 343 327 zmask(ji,jj)=tmask(ji,jj-1,jk)+tmask(ji,jj,jk)+tmask(ji,jj+1,jk)-2 344 END DO345 END DO328 END DO 329 END DO 346 330 ! 347 331 !--- Lateral boundary conditions … … 370 354 zfd(ji,jj) = dir*tra (ji,jj+1,jk)+(1-dir)*tra (ji,jj ,jk) ! FD in the y-direction for T 371 355 zmskl(ji,jj) = dir*zmask(ji,jj )+(1-dir)*zmask(ji,jj+1) 372 END DO373 END DO356 END DO 357 END DO 374 358 ! 375 359 !--- QUICKEST scheme … … 440 424 ENDIF 441 425 442 END SUBROUTINE tra_adv_qck_mer 443 444 445 SUBROUTINE tra_adv_cen2_ver ( pwn, tra , traa ) 446 !!---------------------------------------------------------------------- 447 !! 448 !!---------------------------------------------------------------------- 449 !! * Arguments 450 REAL(wp), INTENT(in) , DIMENSION(jpi,jpj,jpk) :: pwn ! horizontal effective velocity 451 REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: tra, traa 452 ! 453 INTEGER ji, jj, jk 454 REAL(wp) :: za, ze3tr ! temporary scalars 455 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zflux 456 ! 457 ! Vertical advection 458 ! ------------------ 459 ! 460 ! 1. Vertical advective fluxes 461 ! ---------------------------- 462 ! 463 ! Bottom value : flux set to zero 464 zflux(:,:,jpk) = 0.e0 465 ! 466 ! Surface value 467 IF( lk_vvl ) THEN 468 ! variable volume : flux set to zero 469 zflux(:,:, 1 ) = 0.e0 470 ELSE 471 ! free surface-constant volume 472 zflux(:,:, 1 ) = pwn(:,:,1) * tra(:,:,1) 426 END SUBROUTINE tra_adv_qck_j 427 428 429 SUBROUTINE tra_adv_cen2_k ( pwn, ptn, pta ) 430 !!---------------------------------------------------------------------- 431 !! 432 !!---------------------------------------------------------------------- 433 REAL(wp), INTENT(in ), DIMENSION(jpi,jpj,jpk) :: pwn ! vertical effective velocity 434 REAL(wp), INTENT(in ), DIMENSION(jpi,jpj,jpk) :: ptn ! now tracer 435 REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: pta ! tracer general trend 436 !! 437 INTEGER :: ji, jj, jk ! dummy loop indices 438 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zflux ! 3D workspace 439 !!---------------------------------------------------------------------- 440 ! 441 ! !== Vertical advective fluxes ==! 442 zflux(:,:,jpk) = 0.e0 ! Bottom value : flux set to zero 443 ! 444 ! ! Surface value 445 IF( lk_vvl ) THEN ; zflux(:,:, 1 ) = 0.e0 ! Variable volume : flux set to zero 446 ELSE ; zflux(:,:, 1 ) = pwn(:,:,1) * ptn(:,:,1) ! Constant volume : advective flux through the surface 473 447 ENDIF 474 448 ! 475 ! Second order centered tracer flux at w-point 476 ! 477 DO jk = 2, jpkm1 449 DO jk = 2, jpkm1 ! Interior point: second order centered tracer flux at w-point 478 450 DO jj = 2, jpjm1 479 451 DO ji = fs_2, fs_jpim1 ! vector opt. 480 zflux(ji,jj,jk) = pwn(ji,jj,jk)*0.5*( tra(ji,jj,jk-1) + tra(ji,jj,jk) )452 zflux(ji,jj,jk) = 0.5 * pwn(ji,jj,jk) * ( ptn(ji,jj,jk-1) + ptn(ji,jj,jk) ) 481 453 END DO 482 454 END DO 483 455 END DO 484 456 ! 485 ! 2. Tracer flux divergence at t-point added to the general trend 486 ! --------------------------------------------------------------- 487 ! 488 DO jk = 1, jpkm1 457 DO jk = 1, jpkm1 !== Tracer flux divergence added to the general trend ==! 489 458 DO jj = 2, jpjm1 490 459 DO ji = fs_2, fs_jpim1 ! vector opt. 491 ze3tr = 1. / fse3t(ji,jj,jk) 492 ! vertical advective trends 493 za = - ze3tr * ( zflux(ji,jj,jk) - zflux(ji,jj,jk+1) ) 494 ! add it to the general tracer trends 495 traa(ji,jj,jk) = traa(ji,jj,jk) + za 460 pta(ji,jj,jk) = pta(ji,jj,jk) - ( zflux(ji,jj,jk) - zflux(ji,jj,jk+1) ) & 461 & / fse3t(ji,jj,jk) 496 462 END DO 497 463 END DO 498 464 END DO 499 465 ! 500 END SUBROUTINE tra_adv_cen2_ver 501 502 503 SUBROUTINE quickest(fu,fd,fc,fho,fc_cfl) 504 !!---------------------------------------------------------------------- 505 !! 506 !!---------------------------------------------------------------------- 507 !! * Arguments 466 END SUBROUTINE tra_adv_cen2_k 467 468 469 SUBROUTINE quickest( fu, fd, fc, fho, fc_cfl ) 470 !!---------------------------------------------------------------------- 471 !! 472 !!---------------------------------------------------------------------- 508 473 REAL(wp), INTENT(in) , DIMENSION(jpi,jpj) :: fu, fd, fc, fc_cfl 509 474 REAL(wp), INTENT(out) , DIMENSION(jpi,jpj) :: fho … … 513 478 zcoef1(:,:) = 0.5*( fc(:,:) + fd(:,:) ) 514 479 zcoef2(:,:) = 0.5*fc_cfl(:,:)*( fd(:,:) - fc(:,:) ) 515 zcoef3(:,:) = ( ( 1. - ( fc_cfl(:,:)*fc_cfl(:,:) ) )* cst)*zcurv(:,:)480 zcoef3(:,:) = ( ( 1. - ( fc_cfl(:,:)*fc_cfl(:,:) ) )*r1_6 )*zcurv(:,:) 516 481 fho (:,:) = zcoef1(:,:) - zcoef2(:,:) - zcoef3(:,:) ! phi_f QUICKEST 517 482 ! -
trunk/NEMO/OPA_SRC/ZDF/zdfini.F90
r1537 r1559 1 1 MODULE zdfini 2 2 !!====================================================================== 3 !! *** MODULE zdfini *** 4 !! Ocean physics : define vertical mixing variables 5 !!===================================================================== 3 !! *** MODULE zdfini *** 4 !! Ocean physics : read vertical mixing namelist and check consistancy 5 !!====================================================================== 6 !! History : 8.0 ! 1997-06 (G. Madec) Original code from inimix 7 !! 1.0 ! 2002-08 (G. Madec) F90 : free form 8 !! - ! 2005-06 (C. Ethe) KPP parameterization 9 !! - ! 2009-07 (G. Madec) add avmb, avtb in restart for cen2 advection 10 !!---------------------------------------------------------------------- 6 11 7 12 !!---------------------------------------------------------------------- 8 13 !! zdf_init : initialization, namelist read, and parameters control 9 14 !!---------------------------------------------------------------------- 10 !! * Modules used11 15 USE par_oce ! mesh and scale factors 12 16 USE ldftra_oce ! ocean active tracers: lateral physics … … 30 34 PRIVATE 31 35 32 !! * Routine accessibility33 PUBLIC zdf_init ! routine called by opa.F9036 PUBLIC zdf_init ! routine called by opa.F90 37 34 38 !!---------------------------------------------------------------------- 35 !! OPA 9.0 , LOCEAN-IPSL (2005)39 !! NEMO/OPA 3.2 , LOCEAN-IPSL (2009) 36 40 !! $Id$ 37 41 !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) … … 47 51 !! 48 52 !! ** Method : Read namelist namzdf, control logicals 49 !!50 !! History :51 !! ! 97-06 (G. Madec) Original code from inimix52 !! 8.5 ! 02-08 (G. Madec) F90 : free form53 !! 9.0 ! 05-06 (C. Ethe) KPP parameterization54 53 !!---------------------------------------------------------------------- 55 54 INTEGER :: ioptio ! temporary scalar … … 59 58 !!---------------------------------------------------------------------- 60 59 61 REWIND( numnam ) !Read nam_zdf namelist : vertical mixing parameters60 REWIND( numnam ) !* Read nam_zdf namelist : vertical mixing parameters 62 61 READ ( numnam, nam_zdf ) 63 62 64 IF(lwp) THEN !Parameter print63 IF(lwp) THEN !* Parameter print 65 64 WRITE(numout,*) 66 65 WRITE(numout,*) 'zdf_init: vertical physics' … … 81 80 ENDIF 82 81 83 ! Parameter & logicals controls 84 ! ----------------------------- 85 ! ... check of vertical mixing scheme on tracers 86 ! ==> will be done in trazdf module 87 88 ! ... check of mixing coefficient 82 ! !* Parameter & logical controls 83 ! ! ---------------------------- 84 ! 85 ! ! ... check of vertical mixing scheme on tracers 86 ! ==> will be done in trazdf module 87 ! 88 ! ! ... check of mixing coefficient 89 89 IF(lwp) WRITE(numout,*) 90 IF(lwp) WRITE(numout,*) ' 90 IF(lwp) WRITE(numout,*) ' vertical mixing option :' 91 91 ioptio = 0 92 92 IF( lk_zdfcst ) THEN 93 IF(lwp) WRITE(numout,*) ' 93 IF(lwp) WRITE(numout,*) ' constant eddy diffusion coefficients' 94 94 ioptio = ioptio+1 95 95 ENDIF 96 96 IF( lk_zdfric ) THEN 97 IF(lwp) WRITE(numout,*) ' 97 IF(lwp) WRITE(numout,*) ' Richardson dependent eddy coefficients' 98 98 ioptio = ioptio+1 99 99 ENDIF 100 100 IF( lk_zdftke_old ) THEN 101 IF(lwp) WRITE(numout,*) ' 101 IF(lwp) WRITE(numout,*) ' TKE dependent eddy coefficients' 102 102 ioptio = ioptio+1 103 103 ENDIF 104 104 IF( lk_zdftke ) THEN 105 IF(lwp) WRITE(numout,*) ' 105 IF(lwp) WRITE(numout,*) ' TKE dependent eddy coefficients' 106 106 ioptio = ioptio+1 107 107 ENDIF 108 108 IF( lk_zdfkpp ) THEN 109 IF(lwp) WRITE(numout,*) ' 109 IF(lwp) WRITE(numout,*) ' KPP dependent eddy coefficients' 110 110 ioptio = ioptio+1 111 111 ENDIF 112 112 IF( ioptio == 0 .OR. ioptio > 1 .AND. .NOT. lk_esopa ) & 113 CALL ctl_stop( ' one and only one vertical diffusion option has to be defined ' )114 115 ! ... Convection113 & CALL ctl_stop( ' one and only one vertical diffusion option has to be defined ' ) 114 ! 115 ! ! ... Convection 116 116 IF(lwp) WRITE(numout,*) 117 IF(lwp) WRITE(numout,*) ' 117 IF(lwp) WRITE(numout,*) ' convection :' 118 118 ioptio = 0 119 119 IF( ln_zdfnpc ) THEN 120 IF(lwp) WRITE(numout,*) ' 120 IF(lwp) WRITE(numout,*) ' use non penetrative convective scheme' 121 121 ioptio = ioptio+1 122 122 ENDIF 123 123 IF( ln_zdfevd ) THEN 124 IF(lwp) WRITE(numout,*) ' 124 IF(lwp) WRITE(numout,*) ' use enhanced vertical dif. scheme' 125 125 ioptio = ioptio+1 126 126 ENDIF 127 127 IF( lk_zdftke_old .OR. lk_zdftke ) THEN 128 IF(lwp) WRITE(numout,*) ' 128 IF(lwp) WRITE(numout,*) ' use the 1.5 turbulent closure' 129 129 ENDIF 130 130 IF( lk_zdfkpp ) THEN 131 IF(lwp) WRITE(numout,*) ' 131 IF(lwp) WRITE(numout,*) ' use the KPP closure scheme' 132 132 IF(lk_mpp) THEN 133 133 IF(lwp) WRITE(numout,cform_err) 134 IF(lwp) WRITE(numout,*) ' 134 IF(lwp) WRITE(numout,*) 'The KPP scheme is not ready to run in MPI' 135 135 ENDIF 136 136 ENDIF 137 IF ( ioptio > 1 .AND. .NOT. lk_esopa ) & 138 CALL ctl_stop( ' chose between ln_zdfnpc', ' and ln_zdfevd' ) 137 IF ( ioptio > 1 .AND. .NOT. lk_esopa ) CALL ctl_stop( ' chose between ln_zdfnpc and ln_zdfevd' ) 139 138 IF( ioptio == 0 .AND. .NOT.( lk_zdftke_old .OR. lk_zdftke .OR. lk_zdfkpp ) ) & 140 139 CALL ctl_stop( ' except for TKE sor KPP physics, a convection scheme is', & -
trunk/NEMO/OPA_SRC/ZDF/zdfmxl.F90
r1482 r1559 4 4 !! Ocean physics: mixed layer depth 5 5 !!====================================================================== 6 !! History : 7 !! 9.0 ! 03-08 (G. Madec) OpenMP/autotasking optimization6 !! History : 1.0 ! 2003-08 (G. Madec) original code 7 !! 3.2 ! 2009-07 (S. Masson, G. Madec) IOM + merge of DO-loop 8 8 !!---------------------------------------------------------------------- 9 9 !! zdf_mxl : Compute the turbocline and mixed layer depths. 10 10 !!---------------------------------------------------------------------- 11 !! * Modules used12 11 USE oce ! ocean dynamics and tracers variables 13 12 USE dom_oce ! ocean space and time domain variables … … 20 19 PRIVATE 21 20 22 !! * Routine accessibility 23 PUBLIC zdf_mxl ! called by step.F90 21 PUBLIC zdf_mxl ! called by step.F90 24 22 25 !! * Shared module variables 26 INTEGER, PUBLIC, DIMENSION(jpi,jpj) :: & !: 27 nmln !: number of level in the mixed layer 28 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: & !: 29 hmld , & !: mixing layer depth (turbocline) (m) 30 hmlp , & !: mixed layer depth (rho=rho0+zdcrit) (m) 31 hmlpt !: mixed layer depth at t-points (m) 23 INTEGER , PUBLIC, DIMENSION(jpi,jpj) :: nmln !: number of level in the mixed layer 24 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: hmld !: mixing layer depth (turbocline) [m] 25 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: hmlp !: mixed layer depth (rho=rho0+zdcrit) [m] 26 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: hmlpt !: mixed layer depth at t-points [m] 32 27 33 !! * module variables 34 REAL(wp) :: & 35 avt_c = 5.e-4_wp, & ! Kz criterion for the turbocline depth 36 rho_c = 0.01_wp ! density criterion for mixed layer depth 28 REAL(wp) :: avt_c = 5.e-4_wp ! Kz criterion for the turbocline depth 29 REAL(wp) :: rho_c = 0.01_wp ! density criterion for mixed layer depth 37 30 38 31 !! * Substitutions 39 32 # include "domzgr_substitute.h90" 40 33 !!---------------------------------------------------------------------- 41 !! OPA 9.0 , LOCEAN-IPSL (2005)34 !! NEMO/OPA 3.2 , LOCEAN-IPSL (2009) 42 35 !! $Id$ 43 !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt36 !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 44 37 !!---------------------------------------------------------------------- 45 38 … … 51 44 !! 52 45 !! ** Purpose : Compute the turbocline depth and the mixed layer depth 53 !! with density criteria.46 !! with density criteria. 54 47 !! 55 48 !! ** Method : The turbocline depth is the depth at which the vertical … … 58 51 !! value defined locally (avt_c here taken equal to 5 cm/s2) 59 52 !! 60 !! ** Action : 53 !! ** Action : nmln, hmld, hmlp, hmlpt 54 !!---------------------------------------------------------------------- 55 INTEGER, INTENT( in ) :: kt ! ocean time-step index 61 56 !! 62 !! History :63 !! ! 94-11 (M. Imbard) Original code64 !! 8.0 ! 96-01 (E. Guilyardi) sub mixed layer temp.65 !! 8.1 ! 97-07 (G. Madec) optimization66 !! 8.5 ! 02-06 (G. Madec) F90: Free form and module67 !!----------------------------------------------------------------------68 !! * Arguments69 INTEGER, INTENT( in ) :: kt ! ocean time-step index70 71 !! * Local declarations72 57 INTEGER :: ji, jj, jk ! dummy loop indices 73 INTEGER :: ik ! temporary integer 74 INTEGER, DIMENSION(jpi,jpj) :: & 75 imld ! temporary workspace 58 INTEGER :: iki, ikn ! temporary integer 59 INTEGER, DIMENSION(jpi,jpj) :: imld ! temporary workspace 76 60 !!---------------------------------------------------------------------- 77 61 … … 82 66 ENDIF 83 67 84 85 ! 1. Turbocline depth 86 ! ------------------- 87 ! last w-level at which avt<avt_c (starting from the bottom jk=jpk) 88 ! (since avt(.,.,jpk)=0, we have jpk=< imld =< 2 ) 89 DO jk = jpk, 2, -1 68 ! w-level of the mixing and mixed layers 69 nmln(:,:) = mbathy(:,:) ! Initialization to the number of w ocean point mbathy 70 imld(:,:) = mbathy(:,:) 71 DO jk = jpkm1, 2, -1 ! from the bottom to the surface 90 72 DO jj = 1, jpj 91 73 DO ji = 1, jpi 92 IF( avt(ji,jj,jk) < avt_c ) imld(ji,jj) = jk 74 IF( avt (ji,jj,jk) < avt_c ) imld(ji,jj) = jk ! Turbocline 75 IF( rhop(ji,jj,jk) > rhop(ji,jj,1) + rho_c ) nmln(ji,jj) = jk ! Mixed layer 93 76 END DO 94 77 END DO 95 78 END DO 96 97 ! Turbocline depth and sub-turbocline temperature 79 ! depth of the mixing and mixed layers 98 80 DO jj = 1, jpj 99 81 DO ji = 1, jpi 100 ik = imld(ji,jj) 101 hmld (ji,jj) = fsdepw(ji,jj,ik) * tmask(ji,jj,1) 82 iki = imld(ji,jj) 83 ikn = nmln(ji,jj) 84 hmld (ji,jj) = fsdepw(ji,jj,iki) * tmask(ji,jj,1) ! Turbocline depth 85 hmlp (ji,jj) = fsdepw(ji,jj,ikn) * tmask(ji,jj,1) ! Mixed layer depth 86 hmlpt(ji,jj) = fsdept(ji,jj,ikn-1) ! depth of the last T-point inside the mixed layer 102 87 END DO 103 88 END DO 104 89 CALL iom_put( "mldturb", hmld ) ! turbocline depth 105 106 !!gm idea 107 !! 108 !!gm DO jk = jpk, 2, -1 109 !!gm DO jj = 1, jpj 110 !!gm DO ji = 1, jpi 111 !!gm IF( avt(ji,jj,jk) < avt_c ) hmld(ji,jj) = fsdepw(ji,jj,jk) * tmask(ji,jj,1) 112 !!gm END DO 113 !!gm END DO 114 !!gm END DO 115 !!gm 116 117 ! 2. Mixed layer depth 118 ! -------------------- 119 ! Initialization to the number of w ocean point mbathy 120 nmln(:,:) = mbathy(:,:) 121 122 ! Last w-level at which rhop>=rho surf+rho_c (starting from jpk-1) 123 ! (rhop defined at t-point, thus jk-1 for w-level just above) 124 DO jk = jpkm1, 2, -1 125 DO jj = 1, jpj 126 DO ji = 1, jpi 127 IF( rhop(ji,jj,jk) > rhop(ji,jj,1) + rho_c ) nmln(ji,jj) = jk 128 END DO 129 END DO 130 END DO 131 132 ! Mixed layer depth 133 DO jj = 1, jpj 134 DO ji = 1, jpi 135 ik = nmln(ji,jj) 136 hmlp (ji,jj) = fsdepw(ji,jj,ik) * tmask(ji,jj,1) 137 hmlpt(ji,jj) = fsdept(ji,jj,ik-1) 138 END DO 139 END DO 140 CALL iom_put( "mld010", hmlp ) ! mixed layer depth 90 CALL iom_put( "mld010" , hmlp ) ! mixed layer depth 141 91 142 92 IF(ln_ctl) CALL prt_ctl( tab2d_1=REAL(nmln,wp), clinfo1=' nmln : ', tab2d_2=hmld, clinfo2=' hmld : ', ovlap=1 ) 143 93 ! 144 94 END SUBROUTINE zdf_mxl 145 95 -
trunk/NEMO/OPA_SRC/daymod.F90
r1359 r1559 4 4 !! Ocean : calendar 5 5 !!===================================================================== 6 !! History : !94-09 (M. Pontaud M. Imbard) Original code7 !! !97-03 (O. Marti)8 !! !97-05 (G. Madec)9 !! !97-08 (M. Imbard)10 !! 9.0 !03-09 (G. Madec) F90 + nyear, nmonth, nday11 !! !04-01 (A.M. Treguier) new calculation based on adatrj12 !! !06-08 (G. Madec) surface module major update6 !! History : OPA ! 1994-09 (M. Pontaud M. Imbard) Original code 7 !! ! 1997-03 (O. Marti) 8 !! ! 1997-05 (G. Madec) 9 !! ! 1997-08 (M. Imbard) 10 !! NEMO 1.0 ! 2003-09 (G. Madec) F90 + nyear, nmonth, nday 11 !! ! 2004-01 (A.M. Treguier) new calculation based on adatrj 12 !! ! 2006-08 (G. Madec) surface module major update 13 13 !!---------------------------------------------------------------------- 14 14 … … 30 30 USE in_out_manager ! I/O manager 31 31 USE iom ! 32 USE ioipsl, ONLY : ymds2ju 32 USE ioipsl, ONLY : ymds2ju ! for calendar 33 33 USE prtctl ! Print control 34 34 USE restart ! … … 37 37 PRIVATE 38 38 39 PUBLIC day ! called by step.F9040 PUBLIC day_init ! called by istate.F9039 PUBLIC day ! called by step.F90 40 PUBLIC day_init ! called by istate.F90 41 41 42 42 INTEGER , PUBLIC :: nyear !: current year … … 50 50 51 51 REAL(wp), PUBLIC :: fjulday !: julian day 52 REAL(wp), PUBLIC :: adatrj !: number of elapsed days since the begining of the run 53 ! !: it is the accumulated duration of previous runs 54 ! !: that may have been run with different time steps. 55 INTEGER , PUBLIC, DIMENSION(0:1) :: nyear_len !: length in days of the previous/current year 52 REAL(wp), PUBLIC :: adatrj !: number of elapsed days since the begining of the whole simulation 53 ! !: (cumulative duration of previous runs that may have used different time-step size) 54 INTEGER , PUBLIC, DIMENSION(0: 1) :: nyear_len !: length in days of the previous/current year 56 55 INTEGER , PUBLIC, DIMENSION(0:13) :: nmonth_len !: length in days of the months of the current year 57 REAL(wp), PUBLIC, DIMENSION(0:13) :: rmonth_half !: second since the beginning of the year and the halftof the months58 REAL(wp), PUBLIC, DIMENSION(0:13) :: rmonth_end !: second since the beginning of theyear and the end of the months59 REAL(wp), PUBLIC :: sec1jan000 !: second since Jan . 1st 00h of nit000 year and Jan. 1st 00h ofthe current year56 REAL(wp), PUBLIC, DIMENSION(0:13) :: rmonth_half !: second since Jan 1st 0h of the current year and the half of the months 57 REAL(wp), PUBLIC, DIMENSION(0:13) :: rmonth_end !: second since Jan 1st 0h of the current year and the end of the months 58 REAL(wp), PUBLIC :: sec1jan000 !: second since Jan 1st 0h of nit000 year and Jan 1st 0h the current year 60 59 61 60 ! this two variables are wrong DO NOT USE THEM !!! … … 65 64 & 31, 31, 30, 31, 30, 31 /) !: (365 days a year) 66 65 67 68 66 !!---------------------------------------------------------------------- 69 !! OPA 9.0 , LOCEAN-IPSL (2006)67 !! NEMO/OPA 3.2 , LOCEAN-IPSL (2009) 70 68 !! $Id$ 71 69 !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) … … 183 181 rmonth_end(jm) = rmonth_end(jm-1) + rday * REAL( nmonth_len(jm), wp ) 184 182 END DO 185 183 ! 186 184 END SUBROUTINE 187 185 -
trunk/NEMO/OPA_SRC/eosbn2.F90
r1493 r1559 5 5 !! - Brunt-Vaisala frequency 6 6 !!============================================================================== 7 !! History : !89-03 (O. Marti) Original code8 !! 6.0 ! 9 !! 6.0 ! 10 !! 7.0 ! 11 !! 8.1 ! 12 !! 8.1 ! 13 !! !99-02 (G. Madec, N. Grima) semi-implicit pressure gradient14 !! ! 01-09 (M. Ben Jelloul) bugfix onlinear eos15 !! 8.5 !02-10 (G. Madec) add eos_init16 !! 8.5 !02-11 (G. Madec, A. Bozec) partial step, eos_insitu_2d17 !! 9.0 !03-08 (G. Madec) F90, free form18 !! 9.0 !06-08 (G. Madec) add tfreez function7 !! History : OPA ! 1989-03 (O. Marti) Original code 8 !! 6.0 ! 1994-07 (G. Madec, M. Imbard) add bn2 9 !! 6.0 ! 1994-08 (G. Madec) Add Jackett & McDougall eos 10 !! 7.0 ! 1996-01 (G. Madec) statement function for e3 11 !! 8.1 ! 1997-07 (G. Madec) introduction of neos, OPA8.1 12 !! 8.1 ! 1997-07 (G. Madec) density instead of volumic mass 13 !! - ! 1999-02 (G. Madec, N. Grima) semi-implicit pressure gradient 14 !! 8.2 ! 2001-09 (M. Ben Jelloul) bugfix on linear eos 15 !! NEMO 1.0 ! 2002-10 (G. Madec) add eos_init 16 !! - ! 2002-11 (G. Madec, A. Bozec) partial step, eos_insitu_2d 17 !! - ! 2003-08 (G. Madec) F90, free form 18 !! 3.0 ! 2006-08 (G. Madec) add tfreez function 19 19 !!---------------------------------------------------------------------- 20 20 … … 51 51 PUBLIC tfreez ! called by sbcice_... modules 52 52 53 ! !* Namelist (nameos)53 ! !!* Namelist (nameos) * 54 54 INTEGER , PUBLIC :: neos = 0 !: = 0/1/2 type of eq. of state and Brunt-Vaisala frequ. 55 55 REAL(wp), PUBLIC :: ralpha = 2.0e-4 !: thermal expension coeff. (linear equation of state) 56 56 REAL(wp), PUBLIC :: rbeta = 7.7e-4 !: saline expension coeff. (linear equation of state) 57 REAL(wp), PUBLIC :: ralpbet !: alpha / beta ratio 57 58 58 59 !! * Substitutions … … 60 61 # include "vectopt_loop_substitute.h90" 61 62 !!---------------------------------------------------------------------- 62 !! OPA 9.0 , LOCEAN-IPSL (2006)63 !! NEMO/OPA 3.2 , LOCEAN-IPSL (2009) 63 64 !! $Id$ 64 65 !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) … … 107 108 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( out) :: prd ! in situ density 108 109 !! 109 INTEGER :: ji, jj, jk! dummy loop indices110 REAL(wp) :: &111 zt , zs , zh , zsr, & ! temporary scalars112 zr1, zr2, zr3, zr4, & ! " "113 zrhop, ze, zbw, zb, & ! " "114 zd , zc , zaw, za , & ! " "115 zb1, za1, zkw, zk0 ! " "110 INTEGER :: ji, jj, jk ! dummy loop indices 111 REAL(wp) :: zt , zs , zh , zsr ! temporary scalars 112 REAL(wp) :: zr1, zr2, zr3, zr4 ! - - 113 REAL(wp) :: zrhop, ze, zbw, zb ! - - 114 REAL(wp) :: zd , zc , zaw, za ! - - 115 REAL(wp) :: zb1, za1, zkw, zk0 ! - - 116 REAL(wp) :: zrau0r ! - - 116 117 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zws ! temporary workspace 117 118 !!---------------------------------------------------------------------- 118 119 119 SELECT CASE 120 ! 121 CASE ( 0 ) ! Jackett and McDougall (1994) formulation122 !120 SELECT CASE( neos ) 121 ! 122 CASE( 0 ) !== Jackett and McDougall (1994) formulation ==! 123 zrau0r = 1.e0 / rau0 123 124 !CDIR NOVERRCHK 124 125 zws(:,:,:) = SQRT( ABS( psal(:,:,:) ) ) 125 ! ! =============== 126 DO jk = 1, jpkm1 ! Horizontal slab 127 ! ! =============== 126 ! 127 DO jk = 1, jpkm1 128 128 DO jj = 1, jpj 129 129 DO ji = 1, jpi 130 zt = ptem(ji,jj,jk) 131 zs = psal(ji,jj,jk) 132 ! depth 133 zh = fsdept(ji,jj,jk) 134 ! square root salinity 135 zsr= zws(ji,jj,jk) 130 zt = ptem (ji,jj,jk) 131 zs = psal (ji,jj,jk) 132 zh = fsdept(ji,jj,jk) ! depth 133 zsr= zws (ji,jj,jk) ! square root salinity 134 ! 136 135 ! compute volumic mass pure water at atm pressure 137 136 zr1= ( ( ( ( 6.536332e-9*zt-1.120083e-6 )*zt+1.001685e-4)*zt & 138 -9.095290e-3 )*zt+6.793952e-2 )*zt+999.842594137 & -9.095290e-3 )*zt+6.793952e-2 )*zt+999.842594 139 138 ! seawater volumic mass atm pressure 140 zr2= ( ( ( 5.3875e-9*zt-8.2467e-7 ) *zt+7.6438e-5 ) *zt &141 -4.0899e-3 ) *zt+0.824493139 zr2= ( ( ( 5.3875e-9*zt-8.2467e-7 ) *zt+7.6438e-5 ) *zt & 140 & -4.0899e-3 ) *zt+0.824493 142 141 zr3= ( -1.6546e-6*zt+1.0227e-4 ) *zt-5.72466e-3 143 142 zr4= 4.8314e-4 144 143 ! 145 144 ! potential volumic mass (reference to the surface) 146 145 zrhop= ( zr4*zs + zr3*zsr + zr2 ) *zs + zr1 147 146 ! 148 147 ! add the compression terms 149 148 ze = ( -3.508914e-8*zt-1.248266e-8 ) *zt-2.595994e-6 150 149 zbw= ( 1.296821e-6*zt-5.782165e-9 ) *zt+1.045941e-4 151 150 zb = zbw + ze * zs 152 151 ! 153 152 zd = -2.042967e-2 154 153 zc = (-7.267926e-5*zt+2.598241e-3 ) *zt+0.1571896 155 154 zaw= ( ( 5.939910e-6*zt+2.512549e-3 ) *zt-0.1028859 ) *zt - 4.721788 156 155 za = ( zd*zsr + zc ) *zs + zaw 157 156 ! 158 157 zb1= (-0.1909078*zt+7.390729 ) *zt-55.87545 159 158 za1= ( ( 2.326469e-3*zt+1.553190)*zt-65.00517 ) *zt+1044.077 160 159 zkw= ( ( (-1.361629e-4*zt-1.852732e-2 ) *zt-30.41638 ) *zt + 2098.925 ) *zt+190925.6 161 160 zk0= ( zb1*zsr + za1 )*zs + zkw 162 161 ! 163 162 ! masked in situ density anomaly 164 163 prd(ji,jj,jk) = ( zrhop / ( 1.0 - zh / ( zk0 - zh * ( za - zh * zb ) ) ) & 165 - rau0 ) / rau0* tmask(ji,jj,jk)164 & - rau0 ) * zrau0r * tmask(ji,jj,jk) 166 165 END DO 167 166 END DO 168 ! ! =============== 169 END DO ! End of slab 170 ! ! =============== 171 ! 172 CASE ( 1 ) ! Linear formulation function of temperature only 173 ! 174 ! ! =============== 175 DO jk = 1, jpkm1 ! Horizontal slab 176 ! ! =============== 177 DO jj = 1, jpj 178 DO ji = 1, jpi 179 zt = ptem(ji,jj,jk) 180 zs = psal(ji,jj,jk) 181 ! ... density and potential volumic mass 182 prd(ji,jj,jk) = ( 0.0285 - ralpha * zt ) * tmask(ji,jj,jk) 183 END DO 184 END DO 185 ! ! =============== 186 END DO ! End of slab 187 ! ! =============== 188 ! 189 CASE ( 2 ) ! Linear formulation function of temperature and salinity 190 ! 191 ! ! =============== 192 DO jk = 1, jpkm1 ! Horizontal slab 193 ! ! =============== 194 DO jj = 1, jpj 195 DO ji = 1, jpi 196 zt = ptem(ji,jj,jk) 197 zs = psal(ji,jj,jk) 198 ! ... density and potential volumic mass 199 prd(ji,jj,jk) = ( rbeta * zs - ralpha * zt ) * tmask(ji,jj,jk) 200 END DO 201 END DO 202 ! ! =============== 203 END DO ! End of slab 204 ! ! =============== 167 END DO 168 ! 169 CASE( 1 ) !== Linear formulation function of temperature only ==! 170 DO jk = 1, jpkm1 171 prd(:,:,jk) = ( 0.0285 - ralpha * ptem(:,:,jk) ) * tmask(:,:,jk) 172 END DO 173 ! 174 CASE( 2 ) !== Linear formulation function of temperature and salinity ==! 175 DO jk = 1, jpkm1 176 prd(:,:,jk) = ( rbeta * psal(:,:,jk) - ralpha * ptem(:,:,jk) ) * tmask(:,:,jk) 177 END DO 205 178 ! 206 179 CASE DEFAULT 207 !208 180 WRITE(ctmp1,*) ' bad flag value for neos = ', neos 209 181 CALL ctl_stop( ctmp1 ) … … 211 183 END SELECT 212 184 ! 213 IF(ln_ctl) CALL prt_ctl( tab3d_1=prd, clinfo1=' eos : ', ovlap=1, kdim=jpk)185 IF(ln_ctl) CALL prt_ctl( tab3d_1=prd, clinfo1=' eos : ', ovlap=1, kdim=jpk ) 214 186 ! 215 187 END SUBROUTINE eos_insitu … … 267 239 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( out) :: prhop ! potential density (surface referenced) 268 240 269 INTEGER :: ji, jj, jk ! dummy loop indices 270 REAL(wp) :: & ! temporary scalars 271 zt, zs, zh, zsr, zr1, zr2, zr3, zr4, zrhop, ze, zbw, & 272 zb, zd, zc, zaw, za, zb1, za1, zkw, zk0 273 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zws 241 INTEGER :: ji, jj, jk ! dummy loop indices 242 REAL(wp) :: zt, zs, zh, zsr, zr1, zr2, zr3, zr4, zrhop, ze, zbw ! temporary scalars 243 REAL(wp) :: zb, zd, zc, zaw, za, zb1, za1, zkw, zk0, zrau0r ! - - 244 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zws ! 3D workspace 274 245 !!---------------------------------------------------------------------- 275 246 276 247 SELECT CASE ( neos ) 277 248 ! 278 CASE ( 0 ) ! Jackett and McDougall (1994) formulation279 !249 CASE( 0 ) !== Jackett and McDougall (1994) formulation ==! 250 zrau0r = 1.e0 / rau0 280 251 !CDIR NOVERRCHK 281 252 zws(:,:,:) = SQRT( ABS( psal(:,:,:) ) ) 282 ! ! =============== 283 DO jk = 1, jpkm1 ! Horizontal slab 284 ! ! =============== 253 ! 254 DO jk = 1, jpkm1 285 255 DO jj = 1, jpj 286 256 DO ji = 1, jpi 287 zt = ptem(ji,jj,jk) 288 zs = psal(ji,jj,jk) 289 ! depth 290 zh = fsdept(ji,jj,jk) 291 ! square root salinity 292 zsr= zws(ji,jj,jk) 257 zt = ptem (ji,jj,jk) 258 zs = psal (ji,jj,jk) 259 zh = fsdept(ji,jj,jk) ! depth 260 zsr= zws (ji,jj,jk) ! square root salinity 261 ! 293 262 ! compute volumic mass pure water at atm pressure 294 zr1= ( ( ( ( 6.536332e-9*zt-1.120083e-6 )*zt+1.001685e-4 )*zt &295 -9.095290e-3 )*zt+6.793952e-2 )*zt+999.842594263 zr1= ( ( ( ( 6.536332e-9*zt-1.120083e-6 )*zt+1.001685e-4 )*zt & 264 & -9.095290e-3 )*zt+6.793952e-2 )*zt+999.842594 296 265 ! seawater volumic mass atm pressure 297 266 zr2= ( ( ( 5.3875e-9*zt-8.2467e-7 ) *zt+7.6438e-5 ) *zt & 298 -4.0899e-3 ) *zt+0.824493267 & -4.0899e-3 ) *zt+0.824493 299 268 zr3= ( -1.6546e-6*zt+1.0227e-4 ) *zt-5.72466e-3 300 269 zr4= 4.8314e-4 301 270 ! 302 271 ! potential volumic mass (reference to the surface) 303 272 zrhop= ( zr4*zs + zr3*zsr + zr2 ) *zs + zr1 304 273 ! 305 274 ! save potential volumic mass 306 275 prhop(ji,jj,jk) = zrhop * tmask(ji,jj,jk) 307 276 ! 308 277 ! add the compression terms 309 278 ze = ( -3.508914e-8*zt-1.248266e-8 ) *zt-2.595994e-6 310 279 zbw= ( 1.296821e-6*zt-5.782165e-9 ) *zt+1.045941e-4 311 280 zb = zbw + ze * zs 312 281 ! 313 282 zd = -2.042967e-2 314 283 zc = (-7.267926e-5*zt+2.598241e-3 ) *zt+0.1571896 315 284 zaw= ( ( 5.939910e-6*zt+2.512549e-3 ) *zt-0.1028859 ) *zt - 4.721788 316 285 za = ( zd*zsr + zc ) *zs + zaw 317 286 ! 318 287 zb1= (-0.1909078*zt+7.390729 ) *zt-55.87545 319 288 za1= ( ( 2.326469e-3*zt+1.553190)*zt-65.00517 ) *zt+1044.077 320 289 zkw= ( ( (-1.361629e-4*zt-1.852732e-2 ) *zt-30.41638 ) *zt + 2098.925 ) *zt+190925.6 321 290 zk0= ( zb1*zsr + za1 )*zs + zkw 322 291 ! 323 292 ! masked in situ density anomaly 324 293 prd(ji,jj,jk) = ( zrhop / ( 1.0 - zh / ( zk0 - zh * ( za - zh * zb ) ) ) & 325 - rau0 ) / rau0* tmask(ji,jj,jk)294 & - rau0 ) * zrau0r * tmask(ji,jj,jk) 326 295 END DO 327 296 END DO 328 ! ! =============== 329 END DO ! End of slab 330 ! ! =============== 331 ! 332 CASE ( 1 ) ! Linear formulation function of temperature only 333 ! 334 ! ! =============== 335 DO jk = 1, jpkm1 ! Horizontal slab 336 ! ! =============== 337 DO jj = 1, jpj 338 DO ji = 1, jpi 339 zt = ptem(ji,jj,jk) 340 zs = psal(ji,jj,jk) 341 ! ... density and potential volumic mass 342 prd (ji,jj,jk) = ( 0.0285 - ralpha * zt ) * tmask(ji,jj,jk) 343 prhop(ji,jj,jk) = ( rau0 * prd(ji,jj,jk) + rau0 ) * tmask(ji,jj,jk) 344 END DO 345 END DO 346 ! ! =============== 347 END DO ! End of slab 348 ! ! =============== 349 ! 350 CASE ( 2 ) ! Linear formulation function of temperature and salinity 351 ! 352 ! ! =============== 353 DO jk = 1, jpkm1 ! Horizontal slab 354 ! ! =============== 355 DO jj = 1, jpj 356 DO ji = 1, jpi 357 zt = ptem(ji,jj,jk) 358 zs = psal(ji,jj,jk) 359 ! ... density and potential volumic mass 360 prd (ji,jj,jk) = ( rbeta * zs - ralpha * zt ) * tmask(ji,jj,jk) 361 prhop(ji,jj,jk) = ( rau0 * prd(ji,jj,jk) + rau0 ) * tmask(ji,jj,jk) 362 END DO 363 END DO 364 ! ! =============== 365 END DO ! End of slab 366 ! ! =============== 297 END DO 298 ! 299 CASE( 1 ) !== Linear formulation = F( temperature ) ==! 300 DO jk = 1, jpkm1 301 prd (:,:,jk) = ( 0.0285 - ralpha * ptem(:,:,jk) ) * tmask(:,:,jk) 302 prhop(:,:,jk) = ( 1.e0 + prd (:,:,jk) ) * rau0 * tmask(:,:,jk) 303 END DO 304 ! 305 CASE( 2 ) !== Linear formulation = F( temperature , salinity ) ==! 306 DO jk = 1, jpkm1 307 prd (:,:,jk) = ( rbeta * psal(:,:,jk) - ralpha * ptem(:,:,jk) ) * tmask(:,:,jk) 308 prhop(:,:,jk) = ( 1.e0 + prd (:,:,jk) ) * rau0 * tmask(:,:,jk) 309 END DO 367 310 ! 368 311 CASE DEFAULT 369 !370 312 WRITE(ctmp1,*) ' bad flag value for neos = ', neos 371 313 CALL ctl_stop( ctmp1 ) … … 419 361 REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: prd ! in situ density 420 362 !! 421 INTEGER :: ji, jj ! dummy loop indices 422 REAL(wp) :: & ! temporary scalars 423 zt, zs, zh, zsr, zr1, zr2, zr3, zr4, zrhop, ze, zbw, & 424 zb, zd, zc, zaw, za, zb1, za1, zkw, zk0, & 425 zmask 426 REAL(wp), DIMENSION(jpi,jpj) :: zws 363 INTEGER :: ji, jj ! dummy loop indices 364 REAL(wp) :: zt, zs, zh, zsr, zr1, zr2, zr3, zr4, zrhop, ze, zbw ! temporary scalars 365 REAL(wp) :: zb, zd, zc, zaw, za, zb1, za1, zkw, zk0, zmask ! - - 366 REAL(wp), DIMENSION(jpi,jpj) :: zws ! 2D workspace 427 367 !!---------------------------------------------------------------------- 428 368 429 369 prd(:,:) = 0.e0 430 370 431 SELECT CASE 432 ! 433 CASE ( 0 ) ! Jackett and McDougall (1994) formulation371 SELECT CASE( neos ) 372 ! 373 CASE( 0 ) !== Jackett and McDougall (1994) formulation ==! 434 374 ! 435 375 !CDIR NOVERRCHK … … 440 380 END DO 441 381 END DO 442 ! ! =============== 443 DO jj = 1, jpjm1 ! Horizontal slab 444 ! ! =============== 382 DO jj = 1, jpjm1 445 383 DO ji = 1, fs_jpim1 ! vector opt. 446 zmask = tmask(ji,jj,1) ! land/sea bottom mask = surf. mask 447 448 zt = ptem (ji,jj) ! interpolated T 449 zs = psal (ji,jj) ! interpolated S 450 zsr= zws(ji,jj) ! square root of interpolated S 451 zh = pdep(ji,jj) ! depth at the partial step level 452 384 zmask = tmask(ji,jj,1) ! land/sea bottom mask = surf. mask 385 zt = ptem (ji,jj) ! interpolated T 386 zs = psal (ji,jj) ! interpolated S 387 zsr = zws (ji,jj) ! square root of interpolated S 388 zh = pdep (ji,jj) ! depth at the partial step level 389 ! 453 390 ! compute volumic mass pure water at atm pressure 454 zr1 = ( ( ( ( 6.536332e-9*zt-1.120083e-6 )*zt+1.001685e-4 )*zt &455 -9.095290e-3 )*zt+6.793952e-2 )*zt+999.842594391 zr1 = ( ( ( ( 6.536332e-9*zt-1.120083e-6 )*zt+1.001685e-4 )*zt & 392 & -9.095290e-3 )*zt+6.793952e-2 )*zt+999.842594 456 393 ! seawater volumic mass atm pressure 457 zr2 = ( ( ( 5.3875e-9*zt-8.2467e-7 )*zt+7.6438e-5 ) *zt &458 -4.0899e-3 ) *zt+0.824493459 zr3 = ( -1.6546e-6*zt+1.0227e-4 ) *zt-5.72466e-3460 zr4 = 4.8314e-4461 394 zr2 = ( ( ( 5.3875e-9*zt-8.2467e-7 )*zt+7.6438e-5 ) *zt & 395 & -4.0899e-3 ) *zt+0.824493 396 zr3 = ( -1.6546e-6*zt+1.0227e-4 ) *zt-5.72466e-3 397 zr4 = 4.8314e-4 398 ! 462 399 ! potential volumic mass (reference to the surface) 463 400 zrhop= ( zr4*zs + zr3*zsr + zr2 ) *zs + zr1 464 401 ! 465 402 ! add the compression terms 466 403 ze = ( -3.508914e-8*zt-1.248266e-8 ) *zt-2.595994e-6 467 404 zbw= ( 1.296821e-6*zt-5.782165e-9 ) *zt+1.045941e-4 468 405 zb = zbw + ze * zs 469 406 ! 470 407 zd = -2.042967e-2 471 408 zc = (-7.267926e-5*zt+2.598241e-3 ) *zt+0.1571896 472 409 zaw= ( ( 5.939910e-6*zt+2.512549e-3 ) *zt-0.1028859 ) *zt -4.721788 473 410 za = ( zd*zsr + zc ) *zs + zaw 474 411 ! 475 412 zb1= (-0.1909078*zt+7.390729 ) *zt-55.87545 476 413 za1= ( ( 2.326469e-3*zt+1.553190)*zt-65.00517 ) *zt+1044.077 477 414 zkw= ( ( (-1.361629e-4*zt-1.852732e-2 ) *zt-30.41638 ) *zt & 478 +2098.925 ) *zt+190925.6415 & +2098.925 ) *zt+190925.6 479 416 zk0= ( zb1*zsr + za1 )*zs + zkw 480 417 ! 481 418 ! masked in situ density anomaly 482 prd(ji,jj) = ( zrhop / ( 1.0 - zh / ( zk0 - zh * ( za - zh * zb ) ) ) - rau0 ) & 483 & / rau0 * zmask 484 END DO 485 ! ! =============== 486 END DO ! End of slab 487 ! ! =============== 488 ! 489 CASE ( 1 ) ! Linear formulation function of temperature only 490 ! 491 ! ! =============== 492 DO jj = 1, jpjm1 ! Horizontal slab 493 ! ! =============== 419 prd(ji,jj) = ( zrhop / ( 1.0 - zh / ( zk0 - zh * ( za - zh * zb ) ) ) - rau0 ) / rau0 * zmask 420 END DO 421 END DO 422 ! 423 CASE( 1 ) !== Linear formulation = F( temperature ) ==! 424 DO jj = 1, jpjm1 494 425 DO ji = 1, fs_jpim1 ! vector opt. 495 426 prd(ji,jj) = ( 0.0285 - ralpha * ptem(ji,jj) ) * tmask(ji,jj,1) 496 427 END DO 497 ! ! =============== 498 END DO ! End of slab 499 ! ! =============== 500 ! 501 CASE ( 2 ) ! Linear formulation function of temperature and salinity 502 ! 503 ! ! =============== 504 DO jj = 1, jpjm1 ! Horizontal slab 505 ! ! =============== 428 END DO 429 ! 430 CASE( 2 ) !== Linear formulation = F( temperature , salinity ) ==! 431 DO jj = 1, jpjm1 506 432 DO ji = 1, fs_jpim1 ! vector opt. 507 433 prd(ji,jj) = ( rbeta * psal(ji,jj) - ralpha * ptem(ji,jj) ) * tmask(ji,jj,1) 508 434 END DO 509 ! ! =============== 510 END DO ! End of slab 511 ! ! =============== 435 END DO 512 436 ! 513 437 CASE DEFAULT 514 !515 438 WRITE(ctmp1,*) ' bad flag value for neos = ', neos 516 439 CALL ctl_stop( ctmp1 ) … … 556 479 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in ) :: psal ! salinity [psu] 557 480 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( out) :: pn2 ! Brunt-Vaisala frequency [s-1] 558 481 !! 559 482 INTEGER :: ji, jj, jk ! dummy loop indices 560 REAL(wp) :: zgde3w, zt, zs, zh, & ! temporary scalars 561 & zalbet, zbeta ! " " 483 REAL(wp) :: zgde3w, zt, zs, zh, zalbet, zbeta ! temporary scalars 562 484 #if defined key_zdfddm 563 REAL(wp) :: zds 485 REAL(wp) :: zds ! temporary scalars 564 486 #endif 565 487 !!---------------------------------------------------------------------- 566 567 ! pn2 : first and last levels568 ! ---------------------------569 ! bn^2=0. at jk=1 and jpk set in inidtr.F : no computation570 571 488 572 489 ! pn2 : interior points only (2=< jk =< jpkm1 ) 573 490 ! -------------------------- 574 575 SELECT CASE ( neos ) 576 ! 577 CASE ( 0 ) ! Jackett and McDougall (1994) formulation 578 ! 579 ! ! =============== 580 DO jk = 2, jpkm1 ! Horizontal slab 581 ! ! =============== 491 ! 492 SELECT CASE( neos ) 493 ! 494 CASE( 0 ) !== Jackett and McDougall (1994) formulation ==! 495 DO jk = 2, jpkm1 582 496 DO jj = 1, jpj 583 497 DO ji = 1, jpi … … 586 500 zs = 0.5 * ( psal(ji,jj,jk) + psal(ji,jj,jk-1) ) - 35.0 ! salinity anomaly (s-35) at w-point 587 501 zh = fsdepw(ji,jj,jk) ! depth in meters at w-point 588 502 ! 589 503 zalbet = ( ( ( - 0.255019e-07 * zt + 0.298357e-05 ) * zt & ! ratio alpha/beta 590 504 & - 0.203814e-03 ) * zt & … … 599 513 & +( 0.791325e-08 * zt - 0.933746e-06 ) * zt & 600 514 & + 0.380374e-04 ) * zh 601 515 ! 602 516 zbeta = ( ( -0.415613e-09 * zt + 0.555579e-07 ) * zt & ! beta 603 517 & - 0.301985e-05 ) * zt & … … 611 525 & + ( - 0.213127e-11 * zt + 0.192867e-09 ) * zt & 612 526 & - 0.121555e-07 ) * zh 613 527 ! 614 528 pn2(ji,jj,jk) = zgde3w * zbeta * tmask(ji,jj,jk) & ! N^2 615 529 & * ( zalbet * ( ptem(ji,jj,jk-1) - ptem(ji,jj,jk) ) & … … 623 537 END DO 624 538 END DO 625 ! ! =============== 626 END DO ! End of slab 627 ! ! =============== 628 ! 629 CASE ( 1 ) ! Linear formulation function of temperature only 630 ! 631 ! ! =============== 632 DO jk = 2, jpkm1 ! Horizontal slab 633 ! ! =============== 634 DO jj = 1, jpj 635 DO ji = 1, jpi 636 zgde3w = grav / fse3w(ji,jj,jk) * tmask(ji,jj,jk) 637 pn2(ji,jj,jk) = zgde3w * ralpha * ( ptem(ji,jj,jk-1) - ptem(ji,jj,jk) ) 638 END DO 639 END DO 640 ! ! =============== 641 END DO ! End of slab 642 ! ! =============== 643 ! 644 CASE ( 2 ) ! Linear formulation function of temperature and salinity 645 ! 646 ! ! =============== 647 DO jk = 2, jpkm1 ! Horizontal slab 648 ! ! =============== 649 DO jj = 1, jpj 650 DO ji = 1, jpi 651 zgde3w = grav / fse3w(ji,jj,jk) * tmask(ji,jj,jk) 652 pn2(ji,jj,jk) = zgde3w * ( ralpha * ( ptem(ji,jj,jk-1) - ptem(ji,jj,jk) ) & 653 & - rbeta * ( psal(ji,jj,jk-1) - psal(ji,jj,jk) ) ) 654 END DO 655 END DO 539 END DO 540 ! 541 CASE( 1 ) !== Linear formulation = F( temperature ) ==! 542 DO jk = 2, jpkm1 543 pn2(:,:,jk) = grav * ralpha * ( ptem(:,:,jk-1) - ptem(:,:,jk) ) / fse3w(:,:,jk) * tmask(:,:,jk) 544 END DO 545 ! 546 CASE( 2 ) !== Linear formulation = F( temperature , salinity ) ==! 547 DO jk = 2, jpkm1 548 pn2(:,:,jk) = grav * ( ralpha * ( ptem(:,:,jk-1) - ptem(:,:,jk) ) & 549 & - rbeta * ( psal(:,:,jk-1) - psal(:,:,jk) ) ) & 550 & / fse3w(:,:,jk) * tmask(:,:,jk) 551 END DO 656 552 #if defined key_zdfddm 657 ! ! Rrau = (alpha / beta) (dk[t] / dk[s]) 658 zalbet = ralpha / rbeta 553 DO jk = 2, jpkm1 ! Rrau = (alpha / beta) (dk[t] / dk[s]) 659 554 DO jj = 1, jpj 660 555 DO ji = 1, jpi 661 556 zds = ( psal(ji,jj,jk-1) - psal(ji,jj,jk) ) 662 557 IF ( ABS( zds ) <= 1.e-20 ) zds = 1.e-20 663 rrau(ji,jj,jk) = zalbet * ( ptem(ji,jj,jk-1) - ptem(ji,jj,jk) ) / zds558 rrau(ji,jj,jk) = ralpbet * ( ptem(ji,jj,jk-1) - ptem(ji,jj,jk) ) / zds 664 559 END DO 665 560 END DO 561 END DO 666 562 #endif 667 ! ! ===============668 END DO ! End of slab669 ! ! ===============670 563 ! 671 564 CASE DEFAULT 672 !673 565 WRITE(ctmp1,*) ' bad flag value for neos = ', neos 674 566 CALL ctl_stop( ctmp1 ) … … 676 568 END SELECT 677 569 678 IF(ln_ctl) CALL prt_ctl( tab3d_1=pn2, clinfo1=' bn2 : ', ovlap=1, kdim=jpk)570 IF(ln_ctl) CALL prt_ctl( tab3d_1=pn2, clinfo1=' bn2 : ', ovlap=1, kdim=jpk ) 679 571 #if defined key_zdfddm 680 IF(ln_ctl) CALL prt_ctl( tab3d_1=rrau, clinfo1=' rrau : ', ovlap=1, kdim=jpk)572 IF(ln_ctl) CALL prt_ctl( tab3d_1=rrau, clinfo1=' rrau : ', ovlap=1, kdim=jpk ) 681 573 #endif 682 574 ! … … 699 591 REAL(wp), DIMENSION(jpi,jpj) :: ptf ! freezing temperature [Celcius] 700 592 !!---------------------------------------------------------------------- 593 ! 701 594 ptf(:,:) = ( - 0.0575 + 1.710523e-3 * SQRT( psal(:,:) ) & 702 595 & - 2.154996e-4 * psal(:,:) ) * psal(:,:) 596 ! 703 597 END FUNCTION tfreez 704 598 … … 713 607 !!---------------------------------------------------------------------- 714 608 NAMELIST/nameos/ neos, ralpha, rbeta 715 609 !!---------------------------------------------------------------------- 610 ! 716 611 REWIND( numnam ) ! Read Namelist nameos : equation of state 717 612 READ ( numnam, nameos ) 718 719 ! Control print 720 IF(lwp) THEN 613 ! 614 IF(lwp) THEN ! Control print 721 615 WRITE(numout,*) 722 616 WRITE(numout,*) 'eos_init : equation of state' … … 727 621 WRITE(numout,*) ' saline exp. coef. (linear) rbeta = ', rbeta 728 622 ENDIF 729 730 SELECT CASE 731 732 CASE ( 0 ) ! Jackett and McDougall (1994) formulation623 ! 624 SELECT CASE( neos ) 625 ! 626 CASE( 0 ) !== Jackett and McDougall (1994) formulation ==! 733 627 IF(lwp) WRITE(numout,*) 734 628 IF(lwp) WRITE(numout,*) ' use of Jackett & McDougall (1994) equation of state and' 735 629 IF(lwp) WRITE(numout,*) ' McDougall (1987) Brunt-Vaisala frequency' 736 630 ! 737 CASE ( 1 ) ! Linear formulation function of temperature only631 CASE( 1 ) !== Linear formulation = F( temperature ) ==! 738 632 IF(lwp) WRITE(numout,*) 739 633 IF(lwp) WRITE(numout,*) ' use of linear eos rho(T) = rau0 * ( 1.0285 - ralpha * T )' … … 741 635 & ' that T and S are used as state variables' ) 742 636 ! 743 CASE ( 2 ) ! Linear formulation function of temperature and salinity 637 CASE( 2 ) !== Linear formulation = F( temperature , salinity ) ==! 638 ralpbet = ralpha / rbeta 744 639 IF(lwp) WRITE(numout,*) 745 640 IF(lwp) WRITE(numout,*) ' use of linear eos rho(T,S) = rau0 * ( rbeta * S - ralpha * T )' 746 641 ! 747 CASE DEFAULT ! E R R O R in neos642 CASE DEFAULT !== ERROR in neos ==! 748 643 WRITE(ctmp1,*) ' bad flag value for neos = ', neos 749 644 CALL ctl_stop( ctmp1 ) 645 ! 750 646 END SELECT 751 647 ! 752 648 END SUBROUTINE eos_init 753 649 -
trunk/NEMO/OPA_SRC/lib_mpp.F90
r1528 r1559 2291 2291 CONTAINS 2292 2292 2293 FUNCTION mynode( localComm) RESULT (function_value)2293 FUNCTION mynode( localComm ) RESULT (function_value) 2294 2294 INTEGER, OPTIONAL :: localComm 2295 function_value = 02295 IF( PRESENT( localComm ) .OR. .NOT.PRESENT( localComm ) ) function_value = 0 2296 2296 END FUNCTION mynode 2297 2297 … … 2446 2446 2447 2447 SUBROUTINE mpp_ini_znl 2448 INTEGER :: kcom2449 2448 WRITE(*,*) 'mpp_ini_znl: You should not have seen this print! error?' 2450 2449 END SUBROUTINE mpp_ini_znl
Note: See TracChangeset
for help on using the changeset viewer.
