Changeset 508 for trunk/NEMO/OPA_SRC
- Timestamp:
- 2006-10-03T17:58:55+02:00 (18 years ago)
- Location:
- trunk/NEMO/OPA_SRC
- Files:
-
- 15 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMO/OPA_SRC/DIA/diaptr.F90
r460 r508 5 5 !! (please no more than 2 lines) 6 6 !!===================================================================== 7 !! History : 9.0 ! 03-09 (C. Talandir, G. Madec) Original code 8 !! 9.0 ! 06-01 (A. Biastoch) Allow sub-basins computation 9 !!---------------------------------------------------------------------- 10 7 11 !!---------------------------------------------------------------------- 8 12 !! dia_ptr : Poleward Transport Diagnostics module … … 14 18 !! : flux array; Generic interface: ptr_vj_3d, ptr_vj_2d 15 19 !!---------------------------------------------------------------------- 16 !! History :17 !! 9.0 ! 03-09 (C. Talandir, G. Madec) Original code18 !! 9.0 ! 06-01 (A. Biastoch) Allow sub-basins computation19 !!----------------------------------------------------------------------20 !! * Modules used21 20 USE oce ! ocean dynamics and active tracers 22 21 USE dom_oce ! ocean space and time domain … … 26 25 USE dianam 27 26 USE phycst 28 USE ioipsl ! NetCDF IPSL library 27 USE iom 28 USE ioipsl 29 29 USE daymod 30 30 … … 36 36 END INTERFACE 37 37 38 !! * Routine accessibility 39 PUBLIC dia_ptr_init ! call in opa module 40 PUBLIC dia_ptr ! call in step module 41 PUBLIC ptr_vj ! call by tra_ldf & tra_adv routines 42 PUBLIC ptr_vjk ! call by tra_ldf & tra_adv routines 43 44 !! * Share Module variables 45 LOGICAL, PUBLIC :: & !!! ** init namelist (namptr) ** 46 ln_diaptr = .FALSE., & !: Poleward transport flag (T) or not (F) 47 ln_subbas = .FALSE. !: Atlantic/Pacific/Indian basins calculation 48 INTEGER, PUBLIC :: & !!: ** ptr namelist (namptr) ** 49 nf_ptr = 15 !: frequency of ptr computation 50 REAL(wp), PUBLIC, DIMENSION(jpj) :: & !!: poleward transport 51 pht_adv, pst_adv, & !: heat and salt: advection 52 pht_ove, pst_ove, & !: heat and salt: overturning 53 pht_ldf, pst_ldf, & !: heat and salt: lateral diffusion 54 #if defined key_diaeiv 55 pht_eiv, pst_eiv, & !: heat and salt: bolus advection 56 #endif 57 ht_atl,ht_ind,ht_pac, & !: heat 58 st_atl,st_ind,st_pac !: salt 59 REAL(wp),DIMENSION(jpi,jpj) :: & 60 abasin,pbasin,ibasin !: return function value 38 PUBLIC dia_ptr_init ! call in opa module 39 PUBLIC dia_ptr ! call in step module 40 PUBLIC ptr_vj ! call by tra_ldf & tra_adv routines 41 PUBLIC ptr_vjk ! call by tra_ldf & tra_adv routines 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 INTEGER , PUBLIC :: nf_ptr = 15 !: frequency of ptr computation 47 48 REAL(wp), PUBLIC, DIMENSION(jpj) :: pht_adv, pst_adv !: heat and salt poleward transport: advection 49 REAL(wp), PUBLIC, DIMENSION(jpj) :: pht_ove, pst_ove !: heat and salt poleward transport: overturning 50 REAL(wp), PUBLIC, DIMENSION(jpj) :: pht_ldf, pst_ldf !: heat and salt poleward transport: lateral diffusion 51 #if defined key_diaeiv 52 REAL(wp), PUBLIC, DIMENSION(jpj) :: pht_eiv, pst_eiv !: heat and salt poleward transport: bolus advection 53 #endif 54 REAL(wp), PUBLIC, DIMENSION(jpj) :: ht_atl,ht_ind,ht_pac !: heat 55 REAL(wp), PUBLIC, DIMENSION(jpj) :: st_atl,st_ind,st_pac !: salt 56 61 57 62 58 63 !! Module variables 64 REAL(wp), DIMENSION(jpj,jpk) :: & 65 tn_jk , sn_jk , & !: "zonal" mean temperature and salinity 66 v_msf_atl , & !: "meridional" Stream-Function 67 v_msf_glo , & !: "meridional" Stream-Function 68 v_msf_ipc , & !: "meridional" Stream-Function 69 #if defined key_diaeiv 70 v_msf_eiv , & !: bolus "meridional" Stream-Function 71 #endif 72 surf_jk_r !: inverse of the ocean "zonal" section surface 59 REAL(wp), DIMENSION(jpj,jpk) :: tn_jk , sn_jk , & !: "zonal" mean temperature and salinity 60 & v_msf_atl , & !: "meridional" Stream-Function 61 & v_msf_glo , & !: "meridional" Stream-Function 62 & v_msf_ipc , & !: "meridional" Stream-Function 63 & surf_jk_r !: inverse of the ocean "zonal" section surface 64 #if defined key_diaeiv 65 REAL(wp), DIMENSION(jpj,jpk) :: v_msf_eiv !: bolus "meridional" Stream-Function 66 #endif 67 REAL(wp), DIMENSION(jpi,jpj) :: abasin, pbasin, ibasin !: return function value 73 68 74 69 !! * Substitutions … … 78 73 !! OPA 9.0 , LOCEAN-IPSL (2005) 79 74 !! $Header$ 80 !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt75 !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 81 76 !!---------------------------------------------------------------------- 82 77 … … 94 89 !! 95 90 !! ** Action : - p_fval: i-k-mean poleward flux of pva 96 !! 97 !!---------------------------------------------------------------------- 98 !! * arguments 99 REAL(wp) , INTENT(in), DIMENSION(jpi,jpj,jpk) :: & 100 pva ! mask flux array at V-point 101 102 !! * local declarations 103 INTEGER :: ji, jj, jk ! dummy loop arguments 104 INTEGER :: ijpj ! ??? 105 REAL(wp),DIMENSION(jpj) :: & 106 p_fval ! function value 91 !!---------------------------------------------------------------------- 92 REAL(wp) , INTENT(in), DIMENSION(jpi,jpj,jpk) :: pva ! mask flux array at V-point 93 !! 94 INTEGER :: ji, jj, jk ! dummy loop arguments 95 INTEGER :: ijpj ! ??? 96 REAL(wp), DIMENSION(jpj) :: p_fval ! function value 107 97 !!-------------------------------------------------------------------- 108 98 ! 109 99 ijpj = jpj 110 100 p_fval(:) = 0.e0 … … 116 106 END DO 117 107 END DO 118 108 ! 119 109 IF( lk_mpp ) CALL mpp_sum( p_fval, ijpj ) !!bug I presume 120 110 ! 121 111 END FUNCTION ptr_vj_3d 122 123 112 124 113 … … 134 123 !! 135 124 !! ** Action : - p_fval: i-k-mean poleward flux of pva 136 !! 137 !!---------------------------------------------------------------------- 138 !! * arguments 139 REAL(wp) , INTENT(in), DIMENSION(jpi,jpj) :: & 140 pva ! mask flux array at V-point 141 142 !! * local declarations 143 INTEGER :: ji,jj ! dummy loop arguments 144 INTEGER :: ijpj ! ??? 145 REAL(wp),DIMENSION(jpj) :: & 146 p_fval ! function value 125 !!---------------------------------------------------------------------- 126 REAL(wp) , INTENT(in), DIMENSION(jpi,jpj) :: pva ! mask flux array at V-point 127 !! 128 INTEGER :: ji,jj ! dummy loop arguments 129 INTEGER :: ijpj ! ??? 130 REAL(wp), DIMENSION(jpj) :: p_fval ! function value 147 131 !!-------------------------------------------------------------------- 148 132 ! 149 133 ijpj = jpj 150 134 p_fval(:) = 0.e0 … … 154 138 END DO 155 139 END DO 156 140 ! 157 141 IF( lk_mpp ) CALL mpp_sum( p_fval, ijpj ) !!bug I presume 158 159 END FUNCTION ptr_vj_2d 160 142 ! 143 END FUNCTION ptr_vj_2d 161 144 162 145 … … 171 154 !! 172 155 !! ** Action : - p_fval: i-k-mean poleward flux of pva 173 !! 174 !!---------------------------------------------------------------------- 175 !! * arguments 176 REAL(wp) , INTENT(in), DIMENSION(jpi,jpj,jpk) :: & 177 pva ! mask flux array at V-point 178 179 !! * local declarations 180 INTEGER :: ji, jj, jk ! dummy loop arguments 181 INTEGER, DIMENSION (1) :: ish 182 INTEGER, DIMENSION (2) :: ish2 183 REAL(wp),DIMENSION(jpj*jpk) :: & 184 zwork ! temporary vector for mpp_sum 185 REAL(wp),DIMENSION(jpj,jpk) :: & 186 p_fval ! return function value 156 !!---------------------------------------------------------------------- 157 REAL(wp) , INTENT(in), DIMENSION(jpi,jpj,jpk) :: pva ! mask flux array at V-point 158 !! 159 INTEGER :: ji, jj, jk ! dummy loop arguments 160 INTEGER , DIMENSION (1) :: ish 161 INTEGER , DIMENSION (2) :: ish2 162 REAL(wp), DIMENSION(jpj*jpk) :: zwork ! temporary vector for mpp_sum 163 REAL(wp), DIMENSION(jpj,jpk) :: p_fval ! return function value 187 164 !!-------------------------------------------------------------------- 188 165 ! 189 166 p_fval(:,:) = 0.e0 190 167 ! 191 168 DO jk = 1, jpkm1 192 169 DO jj = 2, jpjm1 … … 197 174 END DO 198 175 END DO 199 176 ! 200 177 IF(lk_mpp) THEN 201 178 ish(1) = jpj*jpk ; ish2(1) = jpj ; ish2(2) = jpk … … 204 181 p_fval(:,:)= RESHAPE( zwork, ish2 ) 205 182 END IF 206 183 ! 207 184 END FUNCTION ptr_vjk 185 208 186 209 187 FUNCTION ptr_vtjk( pva ) RESULT ( p_fval ) … … 218 196 !! 219 197 !! ** Action : - p_fval: i-k-mean poleward flux of pva 220 !! 221 !!---------------------------------------------------------------------- 222 !! * arguments 223 REAL(wp) , INTENT(in), DIMENSION(jpi,jpj,jpk) :: & 224 pva ! mask flux array at V-point 225 226 !! * local declarations 227 INTEGER :: ji, jj, jk ! dummy loop arguments 228 INTEGER, DIMENSION (1) :: ish 229 INTEGER, DIMENSION (2) :: ish2 230 REAL(wp),DIMENSION(jpj*jpk) :: & 231 zwork ! temporary vector for mpp_sum 232 REAL(wp),DIMENSION(jpj,jpk) :: & 233 p_fval ! return function value 198 !!---------------------------------------------------------------------- 199 REAL(wp) , INTENT(in), DIMENSION(jpi,jpj,jpk) :: pva ! mask flux array at V-point 200 !! 201 INTEGER :: ji, jj, jk ! dummy loop arguments 202 INTEGER, DIMENSION (1) :: ish 203 INTEGER, DIMENSION (2) :: ish2 204 REAL(wp),DIMENSION(jpj*jpk) :: zwork ! temporary vector for mpp_sum 205 REAL(wp),DIMENSION(jpj,jpk) :: p_fval ! return function value 234 206 !!-------------------------------------------------------------------- 235 207 ! 236 208 p_fval(:,:) = 0.e0 237 209 DO jk = 1, jpkm1 … … 251 223 p_fval(:,:)= RESHAPE(zwork,ish2) 252 224 END IF 253 225 ! 254 226 END FUNCTION ptr_vtjk 255 227 … … 259 231 !! *** ROUTINE dia_ptr *** 260 232 !!---------------------------------------------------------------------- 261 !! * Moudules used262 USE ioipsl263 264 !! * Argument265 233 INTEGER, INTENT(in) :: kt ! ocean time step index 266 267 !! * Local variables 268 INTEGER :: jk,jj,ji ! dummy loop 269 REAL(wp) :: & 270 zsverdrup, & ! conversion from m3/s to Sverdrup 271 zpwatt, & ! conversion from W to PW 272 zggram ! conversion from g to Pg 234 !! 235 INTEGER :: jk, jj, ji ! dummy loop 236 REAL(wp) :: zsverdrup, & ! conversion from m3/s to Sverdrup 237 & zpwatt, & ! conversion from W to PW 238 & zggram ! conversion from g to Pg 273 239 274 240 REAL(wp), DIMENSION(jpi,jpj,jpk) :: & … … 277 243 vs_atl, vs_pac, vs_ind, & 278 244 zv_eiv 279 CHARACTER (len=32) :: & 280 clnam = 'subbasins.nc' 281 INTEGER :: itime,inum,ipi,ipj,ipk ! temporary integer 282 INTEGER, DIMENSION (1) :: istep 283 REAL(wp) :: zdate0,zsecond,zdt ! temporary scalars 284 REAL(wp), DIMENSION(jpidta,jpjdta) :: & 285 zlamt, zphit, zdta ! temporary workspace (NetCDF read) 286 REAL(wp), DIMENSION(jpk) :: & 287 zdept ! temporary workspace (NetCDF read) 245 INTEGER :: inum ! temporary logical unit 288 246 !!---------------------------------------------------------------------- 289 247 … … 293 251 zpwatt = 1.e-15 294 252 zggram = 1.e-6 295 ipi = jpidta296 ipj = jpjdta297 ipk = 1298 itime = 1299 zsecond = 0.e0300 zdate0 = 0.e0301 253 302 254 # if defined key_diaeiv … … 315 267 IF( ln_subbas ) THEN ! Basins computation 316 268 317 IF( kt == nit000 ) THEN ! load basin mask 318 itime = 1 319 ipi = jpidta 320 ipj = jpjdta 321 ipk = 1 322 zdt = 0.e0 323 istep = 0 324 clnam = 'subbasins.nc' 325 326 CALL flinopen(clnam,1,jpidta,1,jpjdta,.FALSE.,ipi,ipj, & 327 & ipk,zlamt,zphit,zdept,itime,istep,zdate0,zdt,inum) 328 329 ! get basins: 330 abasin (:,:) = 0.e0 331 pbasin (:,:) = 0.e0 332 ibasin (:,:) = 0.e0 333 334 ! Atlantic basin 335 CALL flinget(inum,'atlmsk',jpidta,jpjdta,1,itime,1, & 336 & 0,1,jpidta,1,jpjdta,zdta(:,:)) 337 DO jj = 1, nlcj ! interior values 338 DO ji = 1, nlci 339 abasin (ji,jj) = zdta( mig(ji), mjg(jj) ) 340 END DO 341 END DO 342 343 ! Pacific basin 344 CALL flinget(inum,'pacmsk',jpidta,jpjdta,1,itime,1, & 345 & 0,1,jpidta,1,jpjdta,zdta(:,:)) 346 DO jj = 1, nlcj ! interior values 347 DO ji = 1, nlci 348 pbasin (ji,jj) = zdta( mig(ji), mjg(jj) ) 349 END DO 350 END DO 351 352 ! Indian basin 353 CALL flinget(inum,'indmsk',jpidta,jpjdta,1,itime,1, & 354 & 0,1,jpidta,1,jpjdta,zdta(:,:)) 355 DO jj = 1, nlcj ! interior values 356 DO ji = 1, nlci 357 ibasin (ji,jj) = zdta( mig(ji), mjg(jj) ) 358 END DO 359 END DO 360 361 CALL flinclo(inum) 362 269 IF( kt == nit000 ) THEN ! load sub-basin mask 270 CALL iom_open( 'subbasins', inum ) 271 CALL iom_get( inum, jpdom_data, 'atlmsk', abasin ) ! Atlantic basin 272 CALL iom_get( inum, jpdom_data, 'pacmsk', pbasin ) ! Pacific basin 273 CALL iom_get( inum, jpdom_data, 'indmsk', ibasin ) ! Indian basin 274 CALL iom_close( inum ) 363 275 ENDIF 364 276 … … 396 308 #endif 397 309 IF( ln_subbas ) THEN 398 v_msf_atl(:,:) = ptr_vjk( v_atl (:,:,:) )399 v_msf_ipc(:,:) = ptr_vjk( v_ipc (:,:,:) )400 ht_atl(:) = SUM( ptr_vjk( vt_atl(:,:,:)),2 )401 ht_pac(:) = SUM( ptr_vjk( vt_pac(:,:,:)),2 )402 ht_ind(:) = SUM( ptr_vjk( vt_ind(:,:,:)),2 )403 st_atl(:) = SUM( ptr_vjk( vs_atl(:,:,:)),2 )404 st_pac(:) = SUM( ptr_vjk( vs_pac(:,:,:)),2 )405 st_ind(:) = SUM( ptr_vjk( vs_ind(:,:,:)),2 )310 v_msf_atl(:,:) = ptr_vjk( v_atl (:,:,:) ) 311 v_msf_ipc(:,:) = ptr_vjk( v_ipc (:,:,:) ) 312 ht_atl(:) = SUM( ptr_vjk( vt_atl(:,:,:)), 2 ) 313 ht_pac(:) = SUM( ptr_vjk( vt_pac(:,:,:)), 2 ) 314 ht_ind(:) = SUM( ptr_vjk( vt_ind(:,:,:)), 2 ) 315 st_atl(:) = SUM( ptr_vjk( vs_atl(:,:,:)), 2 ) 316 st_pac(:) = SUM( ptr_vjk( vs_pac(:,:,:)), 2 ) 317 st_ind(:) = SUM( ptr_vjk( vs_ind(:,:,:)), 2 ) 406 318 ENDIF 407 319 … … 466 378 ! Close the file 467 379 IF( kt == nitend ) CALL histclo( numptr ) 468 380 ! 469 381 END SUBROUTINE dia_ptr 470 382 … … 475 387 !! 476 388 !! ** Purpose : Initialization, namelist read 477 !!478 389 !!---------------------------------------------------------------------- 479 390 REAL(wp), DIMENSION(jpi,jpj,jpk) :: z_1 ! temporary workspace … … 485 396 REWIND ( numnam ) 486 397 READ ( numnam, namptr ) 487 488 398 489 399 ! Control print … … 513 423 !! 514 424 !! ** Method : NetCDF file 515 !! 516 !!---------------------------------------------------------------------- 517 !! * Arguments 425 !!---------------------------------------------------------------------- 518 426 INTEGER, INTENT(in) :: kt ! ocean time-step index 519 520 !! * Save variables 427 !! 521 428 INTEGER, SAVE :: nhoridz, ndepidzt, ndepidzw, ndex(1) 522 429 523 !! * Local variables 524 CHARACTER (len=40) :: & 525 clhstnam, clop ! temporary names 526 INTEGER :: iline, it, ji ! 527 REAL(wp) :: & 528 zsto, zout, zdt, zmax, & ! temporary scalars 529 zjulian 430 CHARACTER (len=40) :: clhstnam, clop ! temporary names 431 INTEGER :: iline, it, ji ! 432 REAL(wp) :: zsto, zout, zdt, zmax, zjulian ! temporary scalars 530 433 REAL(wp), DIMENSION(jpj) :: zphi, zfoo 531 434 !!---------------------------------------------------------------------- … … 720 623 721 624 ENDIF 722 625 ! 723 626 END SUBROUTINE dia_ptr_wri 724 627 -
trunk/NEMO/OPA_SRC/DYN/dynspg_flt.F90
r474 r508 4 4 !! Ocean dynamics: surface pressure gradient trend 5 5 !!====================================================================== 6 #if ( defined key_dynspg_flt && ! defined key_mpp_omp ) || defined key_esopa 6 !! History 8.0 ! 98-05 (G. Roullet) free surface 7 !! ! 98-10 (G. Madec, M. Imbard) release 8.2 8 !! 8.5 ! 02-08 (G. Madec) F90: Free form and module 9 !! " " ! 02-11 (C. Talandier, A-M Treguier) Open boundaries 10 !! 9.0 ! 04-08 (C. Talandier) New trends organization 11 !! " " ! 05-11 (V. Garnier) Surface pressure gradient organization 12 !! " " ! 06-07 (S. Masson) distributed restart using iom 13 !!---------------------------------------------------------------------- 14 #if ( defined key_dynspg_flt && ! defined key_mpp_omp ) || defined key_esopa 7 15 !!---------------------------------------------------------------------- 8 16 !! 'key_dynspg_flt' filtered free surface 9 17 !! NOT 'key_mpp_omp' k-j-i loop (vector opt.) 18 !!---------------------------------------------------------------------- 10 19 !!---------------------------------------------------------------------- 11 20 !! dyn_spg_flt : update the momentum trend with the surface pressure 12 21 !! gradient in the filtered free surface case with 13 22 !! vector optimization 14 !! ----------------------------------------------------------------------15 !! * Modules used23 !! flt_rst : read/write the time-splitting restart fields in the ocean restart file 24 !!---------------------------------------------------------------------- 16 25 USE oce ! ocean dynamics and tracers 17 26 USE dom_oce ! ocean space and time domain … … 21 30 USE flxrnf ! ocean runoffs 22 31 USE sol_oce ! ocean elliptic solver 32 USE solver ! solver initialization 23 33 USE solpcg ! preconditionned conjugate gradient solver 24 34 USE solsor ! Successive Over-relaxation solver … … 32 42 USE cla_dynspg ! cross land advection 33 43 USE prtctl ! Print control 34 USE in_out_manager ! I/O manager35 44 USE solmat ! matrix construction for elliptic solvers 36 45 USE agrif_opa_interp 46 USE in_out_manager ! I/O manager 47 USE iom 48 USE restart ! only for lrst_oce 37 49 38 50 IMPLICIT NONE 39 51 PRIVATE 40 52 41 !! * Accessibility42 53 PUBLIC dyn_spg_flt ! routine called by step.F90 43 54 … … 48 59 !! OPA 9.0 , LOCEAN-IPSL (2005) 49 60 !! $Header$ 50 !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt61 !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 51 62 !!---------------------------------------------------------------------- 52 63 … … 96 107 !! ** Action : - Update (ua,va) with the surf. pressure gradient trend 97 108 !! 98 !! References : 99 !! Roullet and Madec 1999, JGR. 100 !! 101 !! History : 102 !! ! 98-05 (G. Roullet) Original code 103 !! ! 98-10 (G. Madec, M. Imbard) release 8.2 104 !! 8.5 ! 02-08 (G. Madec) F90: Free form and module 105 !! ! 02-11 (C. Talandier, A-M Treguier) Open boundaries 106 !! 9.0 ! 04-08 (C. Talandier) New trends organization 107 !! " ! 05-11 (V. Garnier) Surface pressure gradient organization 109 !! References : Roullet and Madec 1999, JGR. 108 110 !!--------------------------------------------------------------------- 109 !! * Arguments110 111 INTEGER, INTENT( in ) :: kt ! ocean time-step index 111 INTEGER, INTENT( out ) :: kindic ! solver convergence flag 112 ! if the solver doesn t converge 113 ! the flag is < 0 114 !! * Local declarations 115 INTEGER :: ji, jj, jk ! dummy loop indices 116 REAL(wp) :: & 117 z2dt, z2dtg, zraur, znugdt, & ! temporary scalars 118 znurau, zssha, zgcb, zbtd, & ! " " 119 ztdgu, ztdgv ! " " 112 INTEGER, INTENT( out ) :: kindic ! solver convergence flag (<0 if not converge) 113 !! 114 INTEGER :: ji, jj, jk ! dummy loop indices 115 REAL(wp) :: z2dt, z2dtg, zraur, znugdt, & ! temporary scalars 116 & znurau, zssha, zgcb, zbtd, & ! " " 117 & ztdgu, ztdgv ! " " 120 118 !!---------------------------------------------------------------------- 121 119 ! 122 120 IF( kt == nit000 ) THEN 123 121 IF(lwp) WRITE(numout,*) … … 128 126 spgu(:,:) = 0.e0 ! surface pressure gradient (i-direction) 129 127 spgv(:,:) = 0.e0 ! surface pressure gradient (j-direction) 128 CALL solver_init( nit000 ) ! Elliptic solver initialisation 129 130 ! read filtered free surface arrays in restart file 131 CALL flt_rst( nit000, 'READ' ) ! read or initialize the following fields: 132 ! ! gcx, gcxb, sshb, sshn 130 133 ENDIF 131 134 … … 168 171 169 172 #if defined key_obc 170 ! Update velocities on each open boundary with the radiation algorithm 171 CALL obc_dyn( kt ) 172 ! Correction of the barotropic componant velocity to control the volume of the system 173 CALL obc_vol( kt ) 173 CALL obc_dyn( kt ) ! Update velocities on each open boundary with the radiation algorithm 174 CALL obc_vol( kt ) ! Correction of the barotropic componant velocity to control the volume of the system 174 175 #endif 175 176 #if defined key_agrif 176 ! Update velocities on each coarse/fine interfaces 177 178 CALL Agrif_dyn( kt ) 177 CALL Agrif_dyn( kt ) ! Update velocities on each coarse/fine interfaces 179 178 #endif 180 179 #if defined key_orca_r2 … … 243 242 IF( .NOT. AGRIF_ROOT() ) THEN 244 243 ! add contribution of gradient of after barotropic transport divergence 245 IF( (nbondi == -1) .OR. (nbondi == 2) ) gcb(3,:) = gcb(3,:)&246 & -znugdt * z2dt*laplacu(2,:)*gcdprc(3,:)*hu(2,:)*e2u(2,:)247 IF( (nbondi == 1) .OR. (nbondi == 2) ) gcb(nlci-2,:) = gcb(nlci-2,:)&248 & +znugdt * z2dt*laplacu(nlci-2,:)*gcdprc(nlci-2,:)*hu(nlci-2,:)*e2u(nlci-2,:)249 IF( (nbondj == -1) .OR. (nbondj == 2) ) gcb(:,3) = gcb(:,3)&250 & -znugdt * z2dt*laplacv(:,2)*gcdprc(:,3)*hv(:,2)*e1v(:,2)251 IF( (nbondj == 1) .OR. (nbondj == 2) ) gcb(:,nlcj-2) = gcb(:,nlcj-2)&252 & +znugdt * z2dt*laplacv(:,nlcj-2)*gcdprc(:,nlcj-2)*hv(:,nlcj-2)*e1v(:,nlcj-2)244 IF( nbondi == -1 .OR. nbondi == 2 ) gcb(3 ,:) = & 245 & gcb(3 ,:) - znugdt * z2dt * laplacu(2 ,:) * gcdprc(3 ,:) * hu(2 ,:) * e2u(2 ,:) 246 IF( nbondi == 1 .OR. nbondi == 2 ) gcb(nlci-2,:) = & 247 & gcb(nlci-2,:) + znugdt * z2dt * laplacu(nlci-2,:) * gcdprc(nlci-2,:) * hu(nlci-2,:) * e2u(nlci-2,:) 248 IF( nbondj == -1 .OR. nbondj == 2 ) gcb(: ,3) = & 249 & gcb(:,3 ) - znugdt * z2dt * laplacv(:,2 ) * gcdprc(:,3 ) * hv(:,2 ) * e1v(:,2 ) 250 IF( nbondj == 1 .OR. nbondj == 2 ) gcb(:,nlcj-2) = & 251 & gcb(:,nlcj-2) + znugdt * z2dt * laplacv(:,nlcj-2) * gcdprc(:,nlcj-2) * hv(:,nlcj-2) * e1v(:,nlcj-2) 253 252 ENDIF 254 253 #endif … … 263 262 epsr = eps * eps * rnorme 264 263 ncut = 0 265 ! if rnorme is 0, the solution is 0, the solver is n't called264 ! if rnorme is 0, the solution is 0, the solver is not called 266 265 IF( rnorme == 0.e0 ) THEN 267 266 gcx(:,:) = 0.e0 … … 313 312 IF( .NOT. Agrif_Root() ) THEN 314 313 ! caution : grad D (fine) = grad D (coarse) at coarse/fine interface 315 IF( (nbondi == -1) .OR. (nbondi == 2) ) spgu(2,:) = znugdt * z2dt * laplacu(2,:) * umask(2,:,1)316 IF( (nbondi == 1) .OR. (nbondi == 2)) spgu(nlci-2,:) = znugdt * z2dt * laplacu(nlci-2,:) * umask(nlci-2,:,1)317 IF( (nbondj == -1) .OR. (nbondj == 2) ) spgv(:,2) = znugdt * z2dt * laplacv(:,2) * vmask(:,2,1)318 IF( (nbondj == 1) .OR. (nbondj == 2)) spgv(:,nlcj-2) = znugdt * z2dt * laplacv(:,nlcj-2) * vmask(:,nlcj-2,1)314 IF( nbondi == -1 .OR. nbondi == 2 ) spgu(2 ,:) = znugdt * z2dt * laplacu(2 ,:) * umask(2 ,:,1) 315 IF( nbondi == 1 .OR. nbondi == 2 ) spgu(nlci-2,:) = znugdt * z2dt * laplacu(nlci-2,:) * umask(nlci-2,:,1) 316 IF( nbondj == -1 .OR. nbondj == 2 ) spgv(:,2 ) = znugdt * z2dt * laplacv(:,2 ) * vmask(: ,2,1) 317 IF( nbondj == 1 .OR. nbondj == 2 ) spgv(:,nlcj-2) = znugdt * z2dt * laplacv(:,nlcj-2) * vmask(:,nlcj-2,1) 319 318 ENDIF 320 319 #endif … … 323 322 ! ( c a u t i o n : (ua,va) here are the after velocity not the 324 323 ! trend, the leap-frog time stepping will not 325 ! be done in dynnxt.F routine)324 ! be done in dynnxt.F90 routine) 326 325 DO jk = 1, jpkm1 327 326 DO jj = 2, jpjm1 … … 332 331 END DO 333 332 END DO 334 335 333 336 334 ! Sea surface elevation time stepping … … 358 356 ENDIF 359 357 360 ! ! print sum trends (used for debugging) 361 IF(ln_ctl) CALL prt_ctl( tab2d_1=sshn, clinfo1=' spg - ssh: ', mask1=tmask ) 362 358 ! write filtered free surface arrays in restart file 359 ! -------------------------------------------------- 360 IF( lrst_oce ) CALL flt_rst( kt, 'WRITE' ) 361 362 ! print sum trends (used for debugging) 363 IF(ln_ctl) CALL prt_ctl( tab2d_1=sshn, clinfo1=' spg - ssh: ', mask1=tmask ) 364 ! 363 365 END SUBROUTINE dyn_spg_flt 366 367 368 SUBROUTINE flt_rst( kt, cdrw ) 369 !!--------------------------------------------------------------------- 370 !! *** ROUTINE ts_rst *** 371 !! 372 !! ** Purpose : Read or write filtered free surface arrays in restart file 373 !!---------------------------------------------------------------------- 374 INTEGER , INTENT(in) :: kt ! ocean time-step 375 CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag 376 !!---------------------------------------------------------------------- 377 378 IF( TRIM(cdrw) == 'READ' ) THEN 379 IF( iom_varid( numror, 'gcx' ) > 0 ) THEN 380 ! Caution : extra-hallow 381 ! gcx and gcxb are defined as: DIMENSION(1-jpr2di:jpi+jpr2di,1-jpr2dj:jpj+jpr2dj) 382 CALL iom_get( numror, jpdom_local, 'gcx' , gcx (1:jpi,1:jpj) ) 383 CALL iom_get( numror, jpdom_local, 'gcxb', gcxb(1:jpi,1:jpj) ) 384 CALL iom_get( numror, jpdom_local, 'sshb', sshb(:,:) ) 385 CALL iom_get( numror, jpdom_local, 'sshn', sshn(:,:) ) 386 IF( neuler == 0 ) THEN 387 sshb(:,:) = sshn(:,:) 388 gcxb(:,:) = gcx (:,:) 389 ENDIF 390 ELSE 391 gcx (:,:) = 0.e0 392 gcxb(:,:) = 0.e0 393 sshb(:,:) = 0.e0 394 sshn(:,:) = 0.e0 395 ENDIF 396 ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN 397 ! Caution : extra-hallow 398 ! gcx and gcxb are defined as: DIMENSION(1-jpr2di:jpi+jpr2di,1-jpr2dj:jpj+jpr2dj) 399 CALL iom_rstput( kt, nitrst, numrow, 'gcx' , gcx( 1:jpi,1:jpj) ) 400 CALL iom_rstput( kt, nitrst, numrow, 'gcxb', gcxb(1:jpi,1:jpj) ) 401 CALL iom_rstput( kt, nitrst, numrow, 'sshb', sshb(:,:) ) 402 CALL iom_rstput( kt, nitrst, numrow, 'sshn', sshn(:,:) ) 403 ENDIF 404 ! 405 END SUBROUTINE flt_rst 364 406 365 407 #else -
trunk/NEMO/OPA_SRC/DYN/dynspg_flt_jki.F90
r503 r508 25 25 USE solfet ! FETI solver 26 26 USE solsor_e ! Successive Over-relaxation solver with MPP optimization 27 USE solver 27 28 USE obc_oce ! Lateral open boundary condition 28 29 USE obcdyn ! ocean open boundary condition (obc_dyn routines) … … 35 36 USE solmat ! matrix construction for elliptic solvers 36 37 USE agrif_opa_interp 38 USE restart ! only for lrst_oce 39 USE iom 37 40 38 41 IMPLICIT NONE … … 112 115 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~~~ (free surface constant volume, autotasking case)' 113 116 114 ! set to zero free surface specific arrays 115 spgu(:,:) = 0.e0 ! surface pressure gradient (i-direction) 116 spgv(:,:) = 0.e0 ! surface pressure gradient (j-direction) 117 ! set to zero free surface specific arrays 118 spgu(:,:) = 0.e0 ! surface pressure gradient (i-direction) 119 120 spgv(:,:) = 0.e0 ! surface pressure gradient (j-direction) 121 122 CALL solver_init( nit000 ) ! Elliptic solver initialisation 123 124 ! read filtered free surface arrays in restart file 125 CALL flt_rst( nit000, 'READ' ) ! read or initialize the following fields: 126 ! ! gcx, gcxb, sshb, sshn 117 127 ENDIF 118 128 … … 354 364 CALL lbc_lnk( sshn, 'T', 1. ) 355 365 366 ! write filtered free surface arrays in restart file 367 ! -------------------------------------------------- 368 IF( lrst_oce ) CALL flt_rst( kt, 'WRITE' ) 369 356 370 IF(ln_ctl) THEN ! print sum trends (used for debugging) 357 371 CALL prt_ctl( tab3d_1=ua , clinfo1=' spg - Ua : ', mask1=umask, & … … 362 376 END SUBROUTINE dyn_spg_flt_jki 363 377 378 SUBROUTINE flt_rst( kt, cdrw ) 379 !!--------------------------------------------------------------------- 380 !! *** ROUTINE ts_rst *** 381 !! 382 !! ** Purpose : Read or write filtered free surface arrays in restart file 383 !!---------------------------------------------------------------------- 384 INTEGER , INTENT(in) :: kt ! ocean time-step 385 CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag 386 !!---------------------------------------------------------------------- 387 388 IF( TRIM(cdrw) == 'READ' ) THEN 389 IF( iom_varid( numror, 'gcx' ) > 0 ) THEN 390 ! Caution : extra-hallow 391 ! gcx and gcxb are defined as: DIMENSION(1-jpr2di:jpi+jpr2di,1-jpr2dj:jpj+jpr2dj) 392 CALL iom_get( numror, jpdom_local, 'gcx' , gcx (1:jpi,1:jpj) ) 393 CALL iom_get( numror, jpdom_local, 'gcxb', gcxb(1:jpi,1:jpj) ) 394 CALL iom_get( numror, jpdom_local, 'sshb', sshb(:,:) ) 395 CALL iom_get( numror, jpdom_local, 'sshn', sshn(:,:) ) 396 IF( neuler == 0 ) THEN 397 sshb(:,:) = sshn(:,:) 398 gcxb(:,:) = gcx (:,:) 399 ENDIF 400 ELSE 401 gcx (:,:) = 0.e0 402 gcxb(:,:) = 0.e0 403 sshb(:,:) = 0.e0 404 sshn(:,:) = 0.e0 405 ENDIF 406 ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN 407 ! Caution : extra-hallow 408 ! gcx and gcxb are defined as: DIMENSION(1-jpr2di:jpi+jpr2di,1-jpr2dj:jpj+jpr2dj) 409 CALL iom_rstput( kt, nitrst, numrow, 'gcx' , gcx( 1:jpi,1:jpj) ) 410 CALL iom_rstput( kt, nitrst, numrow, 'gcxb', gcxb(1:jpi,1:jpj) ) 411 CALL iom_rstput( kt, nitrst, numrow, 'sshb', sshb(:,:) ) 412 CALL iom_rstput( kt, nitrst, numrow, 'sshn', sshn(:,:) ) 413 ENDIF 414 ! 415 END SUBROUTINE flt_rst 416 364 417 #else 365 418 !!---------------------------------------------------------------------- -
trunk/NEMO/OPA_SRC/DYN/dynspg_rl.F90
r474 r508 4 4 !! Ocean dynamics: surface pressure gradient trend 5 5 !!====================================================================== 6 !! History : 7.0 ! 96-05 (G. Madec, M. Imbard, M. Guyon) rewitting in 1 7 !! routine, without pointers, and with the same matrix 8 !! for sor and pcg, mpp exchange, and symmetric conditions 9 !! " " ! 96-07 (A. Weaver) Euler forward step 10 !! " " ! 96-11 (A. Weaver) correction to preconditioning: 11 !! 8.0 ! 98-02 (M. Guyon) FETI method 12 !! " " ! 97-09 (J.-M. Molines) Open boundaries 13 !! 8.5 ! 02-08 (G. Madec) F90: Free form and module 14 !! ! 02-11 (C. Talandier, A-M Treguier) Open boundaries 15 !! 9.0 ! 04-08 (C. Talandier) New trends organization 16 !! " " ! 05-11 (V. Garnier) Surface pressure gradient organization 17 !! " " ! 06-07 (S. Masson) distributed restart using iom 18 !!--------------------------------------------------------------------- 6 19 #if defined key_dynspg_rl || defined key_esopa 7 20 !!---------------------------------------------------------------------- … … 10 23 !! dyn_spg_rl : update the momentum trend with the surface pressure 11 24 !! for the rigid-lid case. 12 !! ----------------------------------------------------------------------13 !! * Modules used25 !! rl_rst : read/write the rigid-lid restart fields in the ocean restart file 26 !!---------------------------------------------------------------------- 14 27 USE oce ! ocean dynamics and tracers 15 28 USE dom_oce ! ocean space and time domain … … 19 32 USE zdf_oce ! ocean vertical physics 20 33 USE sol_oce ! ocean elliptic solver 34 USE solver ! solver initialization 21 35 USE solpcg ! preconditionned conjugate gradient solver 22 36 USE solsor ! Successive Over-relaxation solver … … 28 42 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 29 43 USE in_out_manager ! I/O manager 44 USE iom 45 USE restart ! only for lrst_oce 30 46 31 47 IMPLICIT NONE 32 48 PRIVATE 33 49 34 !! * Accessibility35 50 PUBLIC dyn_spg_rl ! called by step.F90 36 51 … … 42 57 !! OPA 9.0 , LOCEAN-IPSL (2005) 43 58 !! $Header$ 44 !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt59 !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 45 60 !!---------------------------------------------------------------------- 46 61 … … 76 91 !! ** Action : - Update (ua,va) with the surf. pressure gradient trend 77 92 !! 78 !! References : 79 !! Madec et al. 1988, ocean modelling, issue 78, 1-6. 80 !! 81 !! History : 82 !! ! 96-05 (G. Madec, M. Imbard, M. Guyon) rewitting in 1 83 !! routine, without pointers, and with the same matrix 84 !! for sor and pcg, mpp exchange, and symmetric conditions 85 !! ! 96-07 (A. Weaver) Euler forward step 86 !! ! 96-11 (A. Weaver) correction to preconditioning: 87 !! ! 98-02 (M. Guyon) FETI method 88 !! ! 98-05 (G. Roullet) free surface 89 !! ! 97-09 (J.-M. Molines) Open boundaries 90 !! 8.5 ! 02-08 (G. Madec) F90: Free form and module 91 !! ! 02-11 (C. Talandier, A-M Treguier) Open boundaries 92 !! 9.0 ! 04-08 (C. Talandier) New trends organization 93 !! " ! 05-11 (V. Garnier) Surface pressure gradient organization 93 !! References : Madec et al. 1988, ocean modelling, issue 78, 1-6. 94 94 !!--------------------------------------------------------------------- 95 !! * Arguments96 95 INTEGER, INTENT( in ) :: kt ! ocean time-step index 97 INTEGER, INTENT( out ) :: kindic ! solver flag, take a negative value 98 ! ! when the solver doesnot converge 99 !! * Local declarations 96 INTEGER, INTENT( out ) :: kindic ! solver flag (<0 when the solver does not converge) 100 97 INTEGER :: ji, jj, jk ! dummy loop indices 101 98 REAL(wp) :: zbsfa, zgcx, z2dt … … 114 111 115 112 ! set to zero rigid-lid specific arrays 116 spgu(:,:) = 0.e0 ! surface pressure gradient (i-direction) 117 spgv(:,:) = 0.e0 ! surface pressure gradient (j-direction) 118 ENDIF 119 120 ! 0. Initializations: 121 ! ------------------- 122 # if defined key_obc 123 ! space index on boundary arrays 124 ib = 1 125 ibm = 2 126 ibm2 = 3 127 ! time index on boundary arrays 128 it = 1 129 itm = 2 130 itm2 = 3 131 # endif 132 113 spgu(:,:) = 0.e0 ! surface pressure gradient (i-direction) 114 spgv(:,:) = 0.e0 ! surface pressure gradient (j-direction) 115 116 CALL solver_init( nit000 ) ! Elliptic solver initialisation 117 118 ! read rigid lid arrays in restart file 119 CALL rl_rst( nit000, 'READ' ) ! read or initialize the following fields: 120 ! ! gcx, gcxb, bsfb, bsfn, bsfd 121 ENDIF 122 123 ! Vertically averaged momentum trend 124 ! ------------------------------------ 133 125 ! ! =============== 134 126 DO jj = 2, jpjm1 ! Vertical slab 135 127 ! ! =============== 136 137 ! 1. Vertically averaged momentum trend 138 ! ------------------------------------- 139 ! initialization to zero 140 spgu(:,jj) = 0. 128 129 spgu(:,jj) = 0. ! initialization to zero 141 130 spgv(:,jj) = 0. 142 ! vertical sum 143 DO jk = 1, jpk 131 DO jk = 1, jpk ! vertical sum 144 132 DO ji = 2, jpim1 145 133 spgu(ji,jj) = spgu(ji,jj) + ua(ji,jj,jk) * fse3u(ji,jj,jk) * umask(ji,jj,jk) … … 147 135 END DO 148 136 END DO 149 ! divide by the depth 150 spgu(:,jj) = spgu(:,jj) * hur(:,jj) 137 spgu(:,jj) = spgu(:,jj) * hur(:,jj) ! divide by the depth 151 138 spgv(:,jj) = spgv(:,jj) * hvr(:,jj) 152 139 … … 155 142 ! ! =============== 156 143 157 !,,,,,,,,,,,,,,,,,,,,,,,,,,,,,synchro,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,158 159 144 ! Boundary conditions on (spgu,spgv) 160 145 CALL lbc_lnk( spgu, 'U', -1. ) 161 146 CALL lbc_lnk( spgv, 'V', -1. ) 162 147 163 !,,,,,,,,,,,,,,,,,,,,,,,,,,,,,synchro,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,, 164 165 ! 2. Barotropic streamfunction trend (bsfd) 148 ! Barotropic streamfunction trend (bsfd) 166 149 ! ---------------------------------- 167 168 ! Curl of the vertically averaged velocity 169 DO jj = 2, jpjm1 150 DO jj = 2, jpjm1 ! Curl of the vertically averaged velocity 170 151 DO ji = 2, jpim1 171 152 gcb(ji,jj) = -gcdprc(ji,jj) & 172 173 153 & *( ( e2v(ji+1,jj )*spgv(ji+1,jj ) - e2v(ji,jj)*spgv(ji,jj) ) & 154 & -( e1u(ji ,jj+1)*spgu(ji ,jj+1) - e1u(ji,jj)*spgu(ji,jj) ) ) 174 155 END DO 175 156 END DO 176 157 177 158 # if defined key_obc 178 ! Open boundary contribution 179 DO jj = 2, jpjm1 159 DO jj = 2, jpjm1 ! Open boundary contribution 180 160 DO ji = 2, jpim1 181 161 gcb(ji,jj) = gcb(ji,jj) - gcdprc(ji,jj) * gcbob(ji,jj) … … 198 178 ! applied the lateral boundary conditions 199 179 IF( nsolv == 4) CALL lbc_lnk_e( gcb, c_solver_pt, 1. ) 200 201 !,,,,,,,,,,,,,,,,,,,,,,,,,,,,,synchro,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,202 180 203 181 ! Relative precision (computation on one processor) … … 234 212 ENDIF 235 213 236 237 214 ! bsf trend update 238 215 ! ---------------- 239 240 216 bsfd(1:nlci,1:nlcj) = gcx(1:nlci,1:nlcj) 241 242 217 243 218 ! update bsf trend with islands trend 244 219 ! ----------------------------------- 245 246 220 IF( lk_isl ) CALL isl_dyn_spg ! update bsfd 247 248 221 249 222 # if defined key_obc … … 337 310 ! 4. Barotropic stream function and array swap 338 311 ! -------------------------------------------- 339 340 312 ! Leap-frog time scheme, time filter and array swap 341 313 IF( neuler == 0 .AND. kt == nit000 ) THEN … … 362 334 363 335 # if defined key_obc 336 ib = 1 ! space index on boundary arrays 337 ibm = 2 338 ibm2 = 3 339 it = 1 ! time index on boundary arrays 340 itm = 2 341 itm2 = 3 342 364 343 ! Swap of boundary arrays 365 344 IF( lp_obc_east ) THEN … … 499 478 ENDIF 500 479 # endif 501 !502 !,,,,,,,,,,,,,,,,,,,,,,,,,,,,,synchro,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,503 480 504 481 ! add the surface pressure trend to the general trend 505 482 ! ----------------------------------------------------- 506 507 483 DO jj = 2, jpjm1 508 509 484 ! update the surface pressure gradient with the barotropic trend 510 485 DO ji = 2, jpim1 … … 519 494 END DO 520 495 END DO 521 522 END DO 496 END DO 497 498 ! write rigid lid arrays in restart file 499 ! -------------------------------------- 500 IF( lrst_oce ) CALL rl_rst( kt, 'WRITE' ) 523 501 524 502 END SUBROUTINE dyn_spg_rl 503 504 505 SUBROUTINE rl_rst( kt, cdrw ) 506 !!--------------------------------------------------------------------- 507 !! *** ROUTINE rl_rst *** 508 !! 509 !! ** Purpose : Read or write rigid-lid arrays in ocean restart file 510 !!---------------------------------------------------------------------- 511 INTEGER , INTENT(in) :: kt ! ocean time-step 512 CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag 513 !!---------------------------------------------------------------------- 514 ! 515 IF( TRIM(cdrw) == 'READ' ) THEN 516 IF( iom_varid( numror, 'gcx' ) > 0 ) THEN 517 ! Caution : extra-hallow 518 ! gcx and gcxb are defined as: DIMENSION(1-jpr2di:jpi+jpr2di,1-jpr2dj:jpj+jpr2dj) 519 CALL iom_get( numror, jpdom_local, 'gcx' , gcx (1:jpi,1:jpj) ) 520 CALL iom_get( numror, jpdom_local, 'gcxb', gcxb(1:jpi,1:jpj) ) 521 CALL iom_get( numror, jpdom_local, 'bsfb', bsfb(:,:) ) 522 CALL iom_get( numror, jpdom_local, 'bsfn', bsfn(:,:) ) 523 CALL iom_get( numror, jpdom_local, 'bsfd', bsfd(:,:) ) 524 IF( neuler == 0 ) THEN 525 gcxb(:,:) = gcx (:,:) 526 bsfb(:,:) = bsfn(:,:) 527 ENDIF 528 ELSE 529 gcx (:,:) = 0.e0 530 gcxb(:,:) = 0.e0 531 bsfb(:,:) = 0.e0 532 bsfn(:,:) = 0.e0 533 bsfd(:,:) = 0.e0 534 ENDIF 535 ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN 536 ! Caution : extra-hallow, gcx and gcxb are defined as: DIMENSION(1-jpr2di:jpi+jpr2di,1-jpr2dj:jpj+jpr2dj) 537 CALL iom_rstput( kt, nitrst, numrow, 'gcx' , gcx (1:jpi,1:jpj) ) 538 CALL iom_rstput( kt, nitrst, numrow, 'gcxb', gcxb(1:jpi,1:jpj) ) 539 CALL iom_rstput( kt, nitrst, numrow, 'bsfb', bsfb(:,:) ) 540 CALL iom_rstput( kt, nitrst, numrow, 'bsfn', bsfn(:,:) ) 541 CALL iom_rstput( kt, nitrst, numrow, 'bsfd', bsfd(:,:) ) 542 ENDIF 543 ! 544 END SUBROUTINE rl_rst 525 545 526 546 #else -
trunk/NEMO/OPA_SRC/DYN/dynspg_ts.F90
r455 r508 4 4 !! Ocean dynamics: surface pressure gradient trend 5 5 !!====================================================================== 6 !! History : 9.0 ! 04-12 (L. Bessieres, G. Madec) Original code 7 !! " " ! 05-11 (V. Garnier, G. Madec) optimization 8 !! 9.0 ! 06-08 (S. Masson) distributed restart using iom 9 !!--------------------------------------------------------------------- 6 10 #if ( defined key_dynspg_ts && ! defined key_mpp_omp ) || defined key_esopa 7 11 !!---------------------------------------------------------------------- … … 9 13 !! NOT 'key_mpp_omp' k-j-i loop (vector opt.) 10 14 !!---------------------------------------------------------------------- 15 !!---------------------------------------------------------------------- 11 16 !! dyn_spg_ts : compute surface pressure gradient trend using a time- 12 17 !! splitting scheme and add to the general trend 18 !! ts_rst : read/write the time-splitting restart fields in the ocean restart file 13 19 !!---------------------------------------------------------------------- 14 20 !! * Modules used … … 27 33 USE dynspg_oce ! surface pressure gradient variables 28 34 USE in_out_manager ! I/O manager 35 USE iom 36 USE restart ! only for lrst_oce 29 37 30 38 IMPLICIT NONE 31 39 PRIVATE 32 40 33 !! * Accessibility34 41 PUBLIC dyn_spg_ts ! routine called by step.F90 42 43 REAL(wp), DIMENSION(jpi,jpj) :: ftnw, ftne, & ! triad of coriolis parameter 44 & ftsw, ftse ! (only used with een vorticity scheme) 45 35 46 36 47 !! * Substitutions … … 74 85 !! ** Action : - Update (ua,va) with the surf. pressure gradient trend 75 86 !! 76 !! References : 77 !! Griffies et al., (2003): A technical guide to MOM4. NOAA/GFDL 78 !! 79 !! History : 80 !! 9.0 ! 04-12 (L. Bessieres, G. Madec) Original code 81 !! ! 05-11 (V. Garnier, G. Madec) optimization 87 !! References : Griffies et al., (2003): A technical guide to MOM4. NOAA/GFDL 82 88 !!--------------------------------------------------------------------- 83 !! * Arguments84 89 INTEGER, INTENT( in ) :: kt ! ocean time-step index 85 90 … … 97 102 zsshb_e, zub_e, zvb_e, & ! " " 98 103 zun_e, zvn_e ! " " 99 REAL(wp), DIMENSION(jpi,jpj),SAVE :: &100 ztnw, ztne, ztsw, ztse101 104 !!---------------------------------------------------------------------- 102 105 … … 109 112 110 113 IF( kt == nit000 ) THEN 111 114 ! 112 115 IF(lwp) WRITE(numout,*) 113 116 IF(lwp) WRITE(numout,*) 'dyn_spg_ts : surface pressure gradient trend' … … 115 118 IF(lwp) WRITE(numout,*) ' Number of sub cycle in 1 time-step (2 rdt) : icycle = ', FLOOR( 2*rdt/rdtbt ) 116 119 117 IF( .NOT. ln_rstart ) THEN 118 ! initialize barotropic specific arrays 119 sshb_b(:,:) = sshb(:,:) 120 sshn_b(:,:) = sshn(:,:) 121 un_b(:,:) = 0.e0 122 vn_b(:,:) = 0.e0 123 ! vertical sum 124 IF( lk_vopt_loop ) THEN ! vector opt., forced unroll 125 DO jk = 1, jpkm1 126 DO ji = 1, jpij 127 un_b(ji,1) = un_b(ji,1) + fse3u(ji,1,jk) * un(ji,1,jk) 128 vn_b(ji,1) = vn_b(ji,1) + fse3v(ji,1,jk) * vn(ji,1,jk) 129 END DO 130 END DO 131 ELSE ! No vector opt. 132 DO jk = 1, jpkm1 133 un_b(:,:) = un_b(:,:) + fse3u(:,:,jk) * un(:,:,jk) 134 vn_b(:,:) = vn_b(:,:) + fse3v(:,:,jk) * vn(:,:,jk) 135 END DO 136 ENDIF 137 ENDIF 120 CALL ts_rst( nit000, 'READ' ) ! read or initialize the following fields: 121 ! ! sshb, sshn, sshb_b, sshn_b, un_b, vn_b 122 138 123 ssha_e(:,:) = sshn(:,:) 139 124 ua_e(:,:) = un_b(:,:) … … 141 126 142 127 IF( ln_dynvor_een ) THEN 143 ztne(1,:) = 0.e0 ; ztnw(1,:) = 0.e0 ; ztse(1,:) = 0.e0 ; ztsw(1,:) = 0.e0128 ftne(1,:) = 0.e0 ; ftnw(1,:) = 0.e0 ; ftse(1,:) = 0.e0 ; ftsw(1,:) = 0.e0 144 129 DO jj = 2, jpj 145 130 DO ji = fs_2, jpi ! vector opt. 146 ztne(ji,jj) = ( ff(ji-1,jj ) + ff(ji ,jj ) + ff(ji ,jj-1) ) / 3.147 ztnw(ji,jj) = ( ff(ji-1,jj-1) + ff(ji-1,jj ) + ff(ji ,jj ) ) / 3.148 ztse(ji,jj) = ( ff(ji ,jj ) + ff(ji ,jj-1) + ff(ji-1,jj-1) ) / 3.149 ztsw(ji,jj) = ( ff(ji ,jj-1) + ff(ji-1,jj-1) + ff(ji-1,jj ) ) / 3.131 ftne(ji,jj) = ( ff(ji-1,jj ) + ff(ji ,jj ) + ff(ji ,jj-1) ) / 3. 132 ftnw(ji,jj) = ( ff(ji-1,jj-1) + ff(ji-1,jj ) + ff(ji ,jj ) ) / 3. 133 ftse(ji,jj) = ( ff(ji ,jj ) + ff(ji ,jj-1) + ff(ji-1,jj-1) ) / 3. 134 ftsw(ji,jj) = ( ff(ji ,jj-1) + ff(ji-1,jj-1) + ff(ji-1,jj ) ) / 3. 150 135 END DO 151 136 END DO 152 137 ENDIF 153 138 ! 154 139 ENDIF 155 140 156 141 ! Local constant initialization 157 142 ! -------------------------------- … … 216 201 END DO 217 202 END DO 218 203 ! 219 204 ELSEIF ( ln_dynvor_ens ) THEN ! enstrophy conserving scheme 220 205 DO jj = 2, jpjm1 … … 228 213 END DO 229 214 END DO 230 215 ! 231 216 ELSEIF ( ln_dynvor_een ) THEN ! enstrophy and energy conserving scheme 232 217 zfac25 = 0.25 … … 241 226 END DO 242 227 END DO 243 228 ! 244 229 ENDIF 245 230 … … 300 285 DO jit = 1, icycle ! sub-time-step loop ! 301 286 ! ! ==================== ! 302 303 287 z2dt_e = 2. * rdtbt 304 288 IF ( jit == 1 ) z2dt_e = rdtbt … … 360 344 END DO 361 345 END DO 362 346 ! 363 347 ELSEIF ( ln_dynvor_ens ) THEN ! enstrophy conserving scheme 364 348 DO jj = 2, jpjm1 … … 379 363 END DO 380 364 END DO 381 365 ! 382 366 ELSEIF ( ln_dynvor_een ) THEN ! energy and enstrophy conserving scheme 383 367 zfac25 = 0.25 … … 397 381 END DO 398 382 END DO 383 ! 399 384 ENDIF 400 385 … … 504 489 END DO 505 490 506 IF(ln_ctl) THEN ! print sum trends (used for debugging) 507 CALL prt_ctl(tab2d_1=sshn, clinfo1=' ssh : ', mask1=tmask) 491 ! write filtered free surface arrays in restart file 492 ! -------------------------------------------------- 493 IF( lrst_oce ) CALL ts_rst( kt, 'WRITE' ) 494 495 ! print sum trends (used for debugging) 496 IF( ln_ctl ) CALL prt_ctl( tab2d_1=sshn, clinfo1=' ssh : ', mask1=tmask ) 497 ! 498 END SUBROUTINE dyn_spg_ts 499 500 501 SUBROUTINE ts_rst( kt, cdrw ) 502 !!--------------------------------------------------------------------- 503 !! *** ROUTINE ts_rst *** 504 !! 505 !! ** Purpose : Read or write time-splitting arrays in restart file 506 !!---------------------------------------------------------------------- 507 INTEGER , INTENT(in) :: kt ! ocean time-step 508 CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag 509 ! 510 INTEGER :: ji, jk ! dummy loop indices 511 !!---------------------------------------------------------------------- 512 ! 513 IF( TRIM(cdrw) == 'READ' ) THEN 514 IF( iom_varid( numror, 'sshn' ) > 0 ) THEN 515 CALL iom_get( numror, jpdom_local, 'sshb' , sshb(:,:) ) 516 CALL iom_get( numror, jpdom_local, 'sshn' , sshn(:,:) ) 517 IF( neuler == 0 ) sshb(:,:) = sshn(:,:) 518 ELSE 519 sshb(:,:) = 0.e0 520 sshn(:,:) = 0.e0 521 ENDIF 522 IF( iom_varid( numror, 'sshn_b' ) > 0 ) THEN 523 CALL iom_get( numror, jpdom_local, 'sshb_b', sshb_b(:,:) ) ! free surface issued 524 CALL iom_get( numror, jpdom_local, 'sshn_b', sshn_b(:,:) ) ! from time-splitting loop 525 CALL iom_get( numror, jpdom_local, 'un_b' , un_b (:,:) ) ! horizontal transports issued 526 CALL iom_get( numror, jpdom_local, 'vn_b' , vn_b (:,:) ) ! from barotropic loop 527 IF( neuler == 0 ) sshb_b(:,:) = sshn_b(:,:) 528 ELSE 529 sshb_b(:,:) = sshb(:,:) 530 sshn_b(:,:) = sshn(:,:) 531 un_b (:,:) = 0.e0 532 vn_b (:,:) = 0.e0 533 ! vertical sum 534 IF( lk_vopt_loop ) THEN ! vector opt., forced unroll 535 DO jk = 1, jpkm1 536 DO ji = 1, jpij 537 un_b(ji,1) = un_b(ji,1) + fse3u(ji,1,jk) * un(ji,1,jk) 538 vn_b(ji,1) = vn_b(ji,1) + fse3v(ji,1,jk) * vn(ji,1,jk) 539 END DO 540 END DO 541 ELSE ! No vector opt. 542 DO jk = 1, jpkm1 543 un_b(:,:) = un_b(:,:) + fse3u(:,:,jk) * un(:,:,jk) 544 vn_b(:,:) = vn_b(:,:) + fse3v(:,:,jk) * vn(:,:,jk) 545 END DO 546 ENDIF 547 ENDIF 548 ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN 549 CALL iom_rstput( kt, nitrst, numrow, 'sshb' , sshb (:,:) ) 550 CALL iom_rstput( kt, nitrst, numrow, 'sshn' , sshn (:,:) ) 551 CALL iom_rstput( kt, nitrst, numrow, 'sshb_b', sshb_b(:,:) ) ! free surface issued 552 CALL iom_rstput( kt, nitrst, numrow, 'sshn_b', sshn_b(:,:) ) ! from barotropic loop 553 CALL iom_rstput( kt, nitrst, numrow, 'un_b' , un_b (:,:) ) ! horizontal transports issued 554 CALL iom_rstput( kt, nitrst, numrow, 'vn_b' , vn_b (:,:) ) ! from barotropic loop 508 555 ENDIF 509 510 END SUBROUTINE dyn_spg_ts 556 ! 557 END SUBROUTINE ts_rst 558 511 559 #else 512 560 !!---------------------------------------------------------------------- -
trunk/NEMO/OPA_SRC/SOL/solmat.F90
r413 r508 4 4 !! solver : construction of the matrix 5 5 !!====================================================================== 6 !! History : 1.0 ! 88-04 (G. Madec) Original code 7 !! ! 93-03 (M. Guyon) symetrical conditions 8 !! ! 93-06 (M. Guyon) suppress pointers 9 !! ! 96-05 (G. Madec) merge sor and pcg formulations 10 !! ! 96-11 (A. Weaver) correction to preconditioning 11 !! ! 98-02 (M. Guyon) FETI method 12 !! 8.5 ! 02-08 (G. Madec) F90: Free form 13 !! ! 02-11 (C. Talandier, A-M. Treguier) Free surface & Open boundaries 14 !! 9.0 ! 05-09 (R. Benshila) add sol_exd for extra outer halo 15 !! 9.0 ! 05-11 (V. Garnier) Surface pressure gradient organization 16 !! 9.0 ! 06-07 (S. Masson) distributed restart using iom 17 !!---------------------------------------------------------------------- 6 18 7 19 !!---------------------------------------------------------------------- … … 29 41 !! OPA 9.0 , LOCEAN-IPSL (2005) 30 42 !! $Header$ 31 !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt43 !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 32 44 !!---------------------------------------------------------------------- 33 45 … … 55 67 !! - gcdmat : preconditioning matrix (diagonal elements) 56 68 !! - gcdprc : inverse of the preconditioning matrix 57 !!58 !! History :59 !! 1.0 ! 88-04 (G. Madec) Original code60 !! ! 91-11 (G. Madec)61 !! ! 93-03 (M. Guyon) symetrical conditions62 !! ! 93-06 (M. Guyon) suppress pointers63 !! ! 96-05 (G. Madec) merge sor and pcg formulations64 !! ! 96-11 (A. Weaver) correction to preconditioning65 !! ! 98-02 (M. Guyon) FETI method66 !! 8.5 ! 02-08 (G. Madec) F90: Free form67 !! ! 02-11 (C. Talandier, A-M. Treguier) Free surface & Open boundaries68 !! 9.0 ! 05-11 (V. Garnier) Surface pressure gradient organization69 69 !!---------------------------------------------------------------------- 70 70 !! * Arguments -
trunk/NEMO/OPA_SRC/ZDF/zdftke.F90
r474 r508 5 5 !! turbulent closure parameterization 6 6 !!===================================================================== 7 !! History : 6.0 ! 91-03 (b. blanke) Original code 8 !! 7.0 ! 91-11 (G. Madec) bug fix 9 !! 7.1 ! 92-10 (G. Madec) new mixing length and eav 10 !! 7.2 ! 93-03 (M. Guyon) symetrical conditions 11 !! 7.3 ! 94-08 (G. Madec, M. Imbard) npdl flag 12 !! 7.5 ! 96-01 (G. Madec) s-coordinates 13 !! 8.0 ! 97-07 (G. Madec) lbc 14 !! 8.1 ! 99-01 (E. Stretta) new option for the mixing length 15 !! 8.5 ! 02-06 (G. Madec) add zdf_tke_init routine 16 !! 8.5 ! 02-08 (G. Madec) ri_c and Free form, F90 17 !! 9.0 ! 04-10 (C. Ethe ) 1D configuration 18 !! 9.0 ! 02-08 (G. Madec) autotasking optimization 19 !! 9.0 ! 06-07 (S. Masson) distributed restart using iom 20 !!---------------------------------------------------------------------- 7 21 #if defined key_zdftke || defined key_esopa 8 22 !!---------------------------------------------------------------------- 9 !! 'key_zdftke' TKE scheme 23 !! 'key_zdftke' TKE vertical physics 24 !!---------------------------------------------------------------------- 10 25 !!---------------------------------------------------------------------- 11 26 !! zdf_tke : update momentum and tracer Kz from a tke scheme 12 27 !! zdf_tke_init : initialization, namelist read, and parameters control 28 !! tke_rst : read/write tke restart in ocean restart file 13 29 !!---------------------------------------------------------------------- 14 !! * Modules used15 30 USE oce ! ocean dynamics and active tracers 16 31 USE dom_oce ! ocean space and time domain 17 32 USE zdf_oce ! ocean vertical physics 18 USE in_out_manager ! I/O manager19 USE lbclnk ! ocean lateral boundary conditions (or mpp link)20 33 USE phycst ! physical constants 21 34 USE taumod ! surface stress 35 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 22 36 USE prtctl ! Print control 37 USE in_out_manager ! I/O manager 38 USE iom 39 USE restart ! only for lrst_oce 23 40 24 41 IMPLICIT NONE 25 42 PRIVATE 26 43 27 !! * Routine accessibility 28 PUBLIC zdf_tke ! routine called in step module 29 PUBLIC zdf_tke_init ! routine called in zdftke_jki module 30 31 !! * Share Module variables 32 LOGICAL, PUBLIC, PARAMETER :: lk_zdftke = .TRUE. !: TKE vertical mixing flag 33 LOGICAL, PUBLIC :: & !!: ** tke namelist (namtke) ** 34 ln_rstke = .FALSE. !: =T restart with tke from a run without tke with 35 ! ! a none zero initial value for en 36 REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) :: & !: 37 en !: now turbulent kinetic energy 38 39 INTEGER, PUBLIC :: & !!! ** tke namelist (namtke) ** 40 nitke = 50 , & ! number of restart iterative loops 41 nmxl = 2 , & ! = 0/1/2/3 flag for the type of mixing length used 42 npdl = 1 , & ! = 0/1/2 flag on prandtl number on vert. eddy coeff. 43 nave = 1 , & ! = 0/1 flag for horizontal average on avt, avmu, avmv 44 navb = 0 ! = 0/1 flag for constant or profile background avt 45 REAL(wp), PUBLIC :: & !!! ** tke namlist (namtke) ** 46 ediff = 0.1_wp , & ! coeff. for vertical eddy coef.; avt=ediff*mxl*sqrt(e) 47 ediss = 0.7_wp , & ! coef. of the Kolmogoroff dissipation 48 ebb = 3.75_wp , & ! coef. of the surface input of tke 49 efave = 1._wp , & ! coef. for the tke vert. diff. coeff.; avtke=efave*avm 50 emin = 0.7071e-6_wp , & ! minimum value of tke (m2/s2) 51 emin0 = 1.e-4_wp , & ! surface minimum value of tke (m2/s2) 52 ri_c = 2._wp / 9._wp ! critic Richardson number 53 REAL(wp), PUBLIC :: & 54 eboost ! multiplicative coeff of the shear product. 55 56 !! caution vectopt_memory change the solution (last digit of the solver stat) 44 PUBLIC zdf_tke ! routine called in step module 45 PUBLIC zdf_tke_init ! routine also called in zdftke_jki module 46 PUBLIC tke_rst ! routine also called in zdftke_jki module 47 48 LOGICAL , PUBLIC, PARAMETER :: lk_zdftke = .TRUE. !: TKE vertical mixing flag 49 REAL(wp), PUBLIC :: eboost !: multiplicative coeff of the shear product. 50 REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) :: en !: now turbulent kinetic energy 57 51 # if defined key_vectopt_memory 58 REAL(wp), DIMENSION(jpi,jpj,jpk), PUBLIC :: & 59 etmean, & ! coefficient used for horizontal smoothing 60 eumean, & ! at t-, u- and v-points 61 evmean ! 52 REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) :: etmean !: coefficient used for horizontal smoothing 53 REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) :: eumean, evmean !: at t-, u- and v-points 62 54 # endif 63 55 56 !! * Namelist (namtke) 57 LOGICAL , PUBLIC :: ln_rstke = .FALSE. !: =T restart with tke from a run without tke with 58 ! ! a none zero initial value for en 59 INTEGER , PUBLIC :: nitke = 50 , & !: number of restart iterative loops 60 & nmxl = 2 , & !: = 0/1/2/3 flag for the type of mixing length used 61 & npdl = 1 , & !: = 0/1/2 flag on prandtl number on vert. eddy coeff. 62 & nave = 1 , & !: = 0/1 flag for horizontal average on avt, avmu, avmv 63 & navb = 0 !: = 0/1 flag for constant or profile background avt 64 REAL(wp), PUBLIC :: ediff = 0.1_wp , & !: coeff. for vertical eddy coef.; avt=ediff*mxl*sqrt(e) 65 & ediss = 0.7_wp , & !: coef. of the Kolmogoroff dissipation 66 & ebb = 3.75_wp , & !: coef. of the surface input of tke 67 & efave = 1._wp , & !: coef. for the tke vert. diff. coeff.; avtke=efave*avm 68 & emin = 0.7071e-6_wp , & !: minimum value of tke (m2/s2) 69 & emin0 = 1.e-4_wp , & !: surface minimum value of tke (m2/s2) 70 & ri_c = 2._wp / 9._wp !: critic Richardson number 71 NAMELIST/namtke/ ln_rstke, ediff, ediss, ebb, efave, emin, emin0, & 72 & ri_c, nitke, nmxl, npdl, nave, navb 73 64 74 # if defined key_cfg_1d 65 REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) :: & 66 e_dis, & ! dissipation turbulent lengh scale 67 e_mix, & ! mixing turbulent lengh scale 68 e_pdl, & ! prandl number 69 e_ric ! local Richardson number 75 ! ! 1D cfg only 76 REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) :: e_dis, e_mix, & ! dissipation and mixing turbulent lengh scales 77 & e_pdl, e_ric ! prandl and local Richardson numbers 70 78 #endif 71 79 … … 74 82 # include "vectopt_loop_substitute.h90" 75 83 !!---------------------------------------------------------------------- 76 !! OPA 9.0 , LOCEAN-IPSL (2005) 84 !! OPA 9.0 , LOCEAN-IPSL (2006) 85 !! $Header$ 86 !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 77 87 !!---------------------------------------------------------------------- 78 88 79 89 CONTAINS 80 90 81 SUBROUTINE zdf_tke 91 SUBROUTINE zdf_tke( kt ) 82 92 !!---------------------------------------------------------------------- 83 93 !! *** ROUTINE zdf_tke *** … … 136 146 !! update avt, avmu, avmv (before vertical eddy coef.) 137 147 !! 138 !! References : 139 !! Gaspar et al., jgr, 95, 1990, 140 !! Blanke and Delecluse, jpo, 1991 141 !! History : 142 !! 6.0 ! 91-03 (b. blanke) Original code 143 !! 7.0 ! 91-11 (G. Madec) bug fix 144 !! 7.1 ! 92-10 (G. Madec) new mixing length and eav 145 !! 7.2 ! 93-03 (M. Guyon) symetrical conditions 146 !! 7.3 ! 94-08 (G. Madec, M. Imbard) npdl flag 147 !! 7.5 ! 96-01 (G. Madec) s-coordinates 148 !! 8.0 ! 97-07 (G. Madec) lbc 149 !! 8.1 ! 99-01 (E. Stretta) new option for the mixing length 150 !! 8.5 ! 02-08 (G. Madec) ri_c and Free form, F90 151 !! 9.0 ! 04-10 (C. Ethe ) 1D configuration 148 !! References : Gaspar et al., jgr, 95, 1990, 149 !! Blanke and Delecluse, jpo, 1991 152 150 !!---------------------------------------------------------------------- 153 !! * Modules used154 151 USE oce , zwd => ua, & ! use ua as workspace 155 152 & zmxlm => ta, & ! use ta as workspace 156 153 & zmxld => sa ! use sa as workspace 157 158 !! * arguments 159 INTEGER, INTENT( in ) :: kt ! ocean time step 160 161 !! * local declarations 162 INTEGER :: ji, jj, jk ! dummy loop arguments 163 REAL(wp) :: & 164 zmlmin, zbbrau, & ! temporary scalars 165 zfact1, zfact2, zfact3, & ! 166 zrn2, zesurf, & ! 167 ztx2, zty2, zav, & ! 168 zcoef, zcof, zsh2, & ! 169 zdku, zdkv, zpdl, zri, & ! 170 zsqen, zesh2, & ! 171 zemxl, zemlm, zemlp 154 ! 155 INTEGER, INTENT(in) :: kt ! ocean time step 156 ! 157 INTEGER :: ji, jj, jk ! dummy loop arguments 158 REAL(wp) :: zmlmin, zbbrau, & ! temporary scalars 159 & zfact1, zfact2, zfact3, & ! 160 & zrn2, zesurf, & ! 161 & ztx2, zty2, zav, & ! 162 & zcoef, zcof, zsh2, & ! 163 & zdku, zdkv, zpdl, zri, & ! 164 & zsqen, zesh2, & ! 165 & zemxl, zemlm, zemlp 172 166 !!-------------------------------------------------------------------- 173 167 174 ! Initialization (first time-step only) 175 ! -------------- 176 IF( kt == nit000 ) CALL zdf_tke_init 177 178 ! Local constant initialization 168 IF( kt == nit000 ) CALL zdf_tke_init ! Initialization (first time-step only) 169 170 ! ! Local constant initialization 179 171 zmlmin = 1.e-8 180 172 zbbrau = .5 * ebb / rau0 … … 183 175 zfact3 = 0.5 * rdt * ediss 184 176 185 186 177 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 187 178 ! I. Mixing length 188 179 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 189 190 180 191 181 ! Buoyancy length scale: l=sqrt(2*e/n**2) … … 204 194 END DO 205 195 END DO 206 207 196 208 197 ! Physical limits for the mixing length … … 291 280 # endif 292 281 293 294 282 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 295 283 ! II Tubulent kinetic energy time stepping 296 284 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 297 298 285 299 286 ! 1. Vertical eddy viscosity on tke (put in zmxlm) and first estimate of avt … … 475 462 CALL lbc_lnk( en , 'W', 1. ) ; CALL lbc_lnk( avt, 'W', 1. ) 476 463 477 478 464 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 479 465 ! III. Before vertical eddy vicosity and diffusivity coefficients … … 601 587 ! ------------------------------===== 602 588 CALL lbc_lnk( avt, 'W', 1. ) 589 590 ! write en in restart file 591 ! ------------------------ 592 IF( lrst_oce ) CALL tke_rst( kt, 'WRITE' ) 603 593 604 594 IF(ln_ctl) THEN … … 624 614 !! ** Action : Increase by 1 the nstop flag is setting problem encounter 625 615 !! 626 !! history :627 !! 8.5 ! 02-06 (G. Madec) original code628 616 !!---------------------------------------------------------------------- 629 !! * Module used630 617 USE dynzdf_exp 631 618 USE trazdf_exp 632 633 !! * local declarations 634 !! caution vectopt_memory change the solution (last digit of the solver stat) 619 ! 635 620 # if defined key_vectopt_memory 636 INTEGER :: ji, jj, jk, jit ! dummy loop indices 621 ! caution vectopt_memory change the solution (last digit of the solver stat) 622 INTEGER :: ji, jj, jk ! dummy loop indices 637 623 # else 638 INTEGER :: jk , jit! dummy loop indices624 INTEGER :: jk ! dummy loop indices 639 625 # endif 640 641 NAMELIST/namtke/ ln_rstke, ediff, ediss, ebb, efave, emin, emin0, &642 ri_c, nitke, nmxl, npdl, nave, navb643 626 !!---------------------------------------------------------------------- 644 627 … … 681 664 ! Check nmxl and npdl values 682 665 IF( nmxl < 0 .OR. nmxl > 3 ) CALL ctl_stop( ' bad flag: nmxl is < 0 or > 3 ' ) 683 IF 666 IF( npdl < 0 .OR. npdl > 1 ) CALL ctl_stop( ' bad flag: npdl is < 0 or > 1 ' ) 684 667 685 668 ! Horizontal average : initialization of weighting arrays … … 691 674 IF(lwp) WRITE(numout,*) ' no horizontal average on avt, avmu, avmv' 692 675 IF(lwp) WRITE(numout,*) ' only in very high horizontal resolution !' 693 !! caution vectopt_memory change the solution (last digit of the solver stat)694 676 # if defined key_vectopt_memory 677 ! caution vectopt_memory change the solution (last digit of the solver stat) 695 678 ! weighting mean arrays etmean, eumean and evmean 696 679 ! ( 1 1 ) ( 1 ) … … 720 703 CASE ( 1 ) ! horizontal average 721 704 IF(lwp) WRITE(numout,*) ' horizontal average on avt, avmu, avmv' 722 !! caution vectopt_memory change the solution (last digit of the solver stat)723 705 # if defined key_vectopt_memory 706 ! caution vectopt_memory change the solution (last digit of the solver stat) 724 707 ! weighting mean arrays etmean, eumean and evmean 725 708 ! ( 1 1 ) ( 1/2 1/2 ) ( 1/2 1 1/2 ) … … 790 773 791 774 792 ! Initialization ofturbulent kinetic energy ( en )775 ! read or initialize turbulent kinetic energy ( en ) 793 776 ! ------------------------------------------------- 794 IF( ln_rstart ) THEN 795 ! no en field in the restart file, en set by iterative loop 796 IF( ln_rstke ) THEN 797 en (:,:,:) = emin * tmask(:,:,:) 798 DO jit = 2, nitke+1 799 CALL zdf_tke( jit ) 800 END DO 801 ENDIF 802 ! otherwise en is already read in dtrlec called by inidtr 803 ELSE 804 ! no restart: en set to emin 805 en(:,:,:) = emin * tmask(:,:,:) 806 ENDIF 807 777 CALL tke_rst( nit000, 'READ' ) 778 ! 808 779 END SUBROUTINE zdf_tke_init 780 781 782 SUBROUTINE tke_rst( kt, cdrw ) 783 !!--------------------------------------------------------------------- 784 !! *** ROUTINE ts_rst *** 785 !! 786 !! ** Purpose : Read or write filtered free surface arrays in restart file 787 !! 788 !! ** Method : 789 !! 790 !!---------------------------------------------------------------------- 791 INTEGER , INTENT(in) :: kt ! ocean time-step 792 CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag 793 ! 794 INTEGER :: jit ! dummy loop indices 795 !!---------------------------------------------------------------------- 796 ! 797 IF( TRIM(cdrw) == 'READ' ) THEN 798 IF( ln_rstart ) THEN 799 IF( iom_varid( numror, 'en' ) > 0 .AND. .NOT.(ln_rstke) ) THEN 800 CALL iom_get( numror, jpdom_local, 'en', en ) 801 ELSE 802 IF(lwp .AND. iom_varid(numror,'en') > 0 ) WRITE(numout,*) ' ===>>>> : previous run without tke scheme' 803 IF(lwp .AND. ln_rstke ) WRITE(numout,*) ' ===>>>> : We do not use en from the restart file' 804 IF(lwp) WRITE(numout,*) ' ===>>>> : en set by iterative loop' 805 IF(lwp) WRITE(numout,*) ' ======= =========' 806 en (:,:,:) = emin * tmask(:,:,:) 807 DO jit = 2, nitke+1 808 CALL zdf_tke( jit ) 809 END DO 810 ENDIF 811 ELSE 812 en(:,:,:) = emin * tmask(:,:,:) ! no restart: en set to emin 813 ENDIF 814 ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN 815 CALL iom_rstput( kt, nitrst, numrow, 'en', en ) 816 ENDIF 817 ! 818 END SUBROUTINE tke_rst 809 819 810 820 #else … … 812 822 !! Dummy module : NO TKE scheme 813 823 !!---------------------------------------------------------------------- 814 LOGICAL, PUBLIC, PARAMETER :: lk_zdftke = .FALSE. !: TKE flag824 PUBLIC, LOGICAL, PARAMETER :: lk_zdftke = .FALSE. !: TKE flag 815 825 CONTAINS 816 826 SUBROUTINE zdf_tke( kt ) ! Empty routine -
trunk/NEMO/OPA_SRC/ZDF/zdftke_jki.F90
r463 r508 5 5 !! turbulent closure parameterization 6 6 !!===================================================================== 7 !! History : 8 !! 9.0 ! 02-08 (G. Madec) autotasking optimization 9 !!---------------------------------------------------------------------- 7 10 #if defined key_zdftke || defined key_esopa 8 11 !!---------------------------------------------------------------------- … … 23 26 USE taumod ! surface stress 24 27 USE prtctl ! Print control 28 USE restart ! only for lrst_oce 25 29 26 30 IMPLICIT NONE … … 99 103 !! Gaspar et al., jgr, 95, 1990, 100 104 !! Blanke and Delecluse, jpo, 1991 101 !! History :102 !! 9.0 ! 02-08 (G. Madec) autotasking optimization103 105 !!---------------------------------------------------------------------- 104 106 !! * Modules used … … 518 520 CALL lbc_lnk( avt, 'W', 1. ) 519 521 522 ! write en in restart file 523 ! ------------------------ 524 IF( lrst_oce ) CALL tke_rst( kt, 'WRITE' ) 525 520 526 IF(ln_ctl) THEN 521 527 CALL prt_ctl(tab3d_1=en , clinfo1=' tke - e: ', tab3d_2=avt , clinfo2=' t: ', ovlap=1, kdim=jpk) -
trunk/NEMO/OPA_SRC/in_out_manager.F90
r472 r508 1 1 MODULE in_out_manager 2 !!====================================================================== 3 !! *** MODULE in_out_manager *** 4 !! Ocean physics: vertical mixing coefficient compute from the tke 5 !! turbulent closure parameterization 6 !!===================================================================== 7 !! History : 8.5 ! 02-06 (G. Madec) original code 8 !! 9.0 ! 06-07 (S. Masson) iom, add ctl_stop, ctl_warn 9 !!---------------------------------------------------------------------- 2 10 3 USE lib_print ! formated print library 11 !!---------------------------------------------------------------------- 12 !! ctl_stop : update momentum and tracer Kz from a tke scheme 13 !! ctl_warn : initialization, namelist read, and parameters control 14 !!---------------------------------------------------------------------- 4 15 USE par_kind 5 16 USE par_oce 17 USE lib_print ! formated print library 6 18 7 19 PUBLIC 8 20 9 21 !!---------------------------------------------------------------------- 10 !! namelist parameters 11 !! ------------------------------------- 12 ! namrun: parameters of the run 13 CHARACTER (len=16) :: & !: 14 cexper = "exp0" !: experiment name used for output filename 15 16 LOGICAL :: & !!: * namelist namrun * 17 ln_rstart = .FALSE. , & !: start from (F) rest or (T) a restart file 18 ln_ctl = .FALSE. !: run control for debugging 19 20 INTEGER :: & !!: * namelist namrun * 21 no = 0 , & !: job number 22 nrstdt = 0 , & !: control of the time step (0, 1 or 2) 23 nit000 = 1 , & !: index of the first time step 24 nitend = 10 , & !: index of the last time step 25 ndate0 = 961115 , & !: initial calendar date aammjj 26 nleapy = 0 , & !: Leap year calendar flag (0/1 or 30) 27 ninist = 0 , & !: initial state output flag (0/1) 28 nbench = 0 !: benchmark parameter (0/1) 22 !! namrun namelist parameters 23 !!---------------------------------------------------------------------- 24 CHARACTER (len=16) :: cexper = "exp0" !: experiment name used for output filename 25 LOGICAL :: ln_rstart = .FALSE. , & !: start from (F) rest or (T) a restart file 26 & ln_ctl = .FALSE. !: run control for debugging 27 INTEGER :: no = 0 , & !: job number 28 & nrstdt = 0 , & !: control of the time step (0, 1 or 2) 29 & nit000 = 1 , & !: index of the first time step 30 & nitend = 10 , & !: index of the last time step 31 & ndate0 = 961115 , & !: initial calendar date aammjj 32 & nleapy = 0 , & !: Leap year calendar flag (0/1 or 30) 33 & ninist = 0 , & !: initial state output flag (0/1) 34 & nbench = 0 !: benchmark parameter (0/1) 29 35 30 36 !!---------------------------------------------------------------------- 31 !! output monitoring 32 !! ----------------------------------- 33 34 INTEGER :: & !: 35 nstock = 10 , & !: restart file frequency 36 nprint = 0 , & !: level of print (0 no print) 37 nwrite = 10 , & !: restart file frequency 38 nictls = 0 , & !: Start i indice for the SUM control 39 nictle = 0 , & !: End i indice for the SUM control 40 njctls = 0 , & !: Start j indice for the SUM control 41 njctle = 0 , & !: End j indice for the SUM control 42 isplt = 1 , & !: number of processors following i 43 jsplt = 1 , & !: number of processors following j 44 ijsplt = 1 !: nb of local domain = nb of processors 37 !! output monitoring 38 !!---------------------------------------------------------------------- 39 INTEGER :: nstock = 10 , & !: restart file frequency 40 & nprint = 0 , & !: level of print (0 no print) 41 & nwrite = 10 , & !: restart file frequency 42 & nictls = 0 , & !: Start i indice for the SUM control 43 & nictle = 0 , & !: End i indice for the SUM control 44 & njctls = 0 , & !: Start j indice for the SUM control 45 & njctle = 0 , & !: End j indice for the SUM control 46 & isplt = 1 , & !: number of processors following i 47 & jsplt = 1 , & !: number of processors following j 48 & ijsplt = 1 !: nb of local domain = nb of processors 45 49 46 50 !!---------------------------------------------------------------------- 47 !! logical units 48 !! ------------------------------ 49 INTEGER :: & !: 50 numstp = 1 , & !: logical unit for time step 51 numout = 2 , & !: logical unit for output print 52 numnam = 3 , & !: logical unit for namelist 53 numnam_ice = 4 , & !: logical unit for ice namelist 54 numevo_ice = 17 , & !: logical unit for ice variables (temp. evolution) 55 numice_dmp = 18 , & !: logical unit for ice variables (damping) 56 numsol = 25 , & !: logical unit for solver statistics 57 numwri = 40 , & !: logical unit for output write 58 numisp = 41 , & !: logical unit for island statistics 59 numgap = 45 , & !: logical unit for differences diagnostic 60 numbol = 67 , & !: logical unit for "bol" diagnostics 61 numptr = 68 , & !: logical unit for Poleward TRansports 62 numflo = 69 !: logical unit for drifting floats 63 ! !: * coupled units 51 !! logical units 52 !!---------------------------------------------------------------------- 53 INTEGER :: numstp = 1 , & !: logical unit for time step 54 & numout = 2 , & !: logical unit for output print 55 & numnam = 3 , & !: logical unit for namelist 56 & numnam_ice = 4 , & !: logical unit for ice namelist 57 & numevo_ice = 17 , & !: logical unit for ice variables (temp. evolution) 58 & numsol = 25 , & !: logical unit for solver statistics 59 & numwri = 40 , & !: logical unit for output write 60 & numisp = 41 , & !: logical unit for island statistics 61 & numgap = 45 , & !: logical unit for differences diagnostic 62 & numbol = 67 , & !: logical unit for "bol" diagnostics 63 & numptr = 68 , & !: logical unit for Poleward TRansports 64 & numflo = 69 !: logical unit for drifting floats 64 65 65 66 !!---------------------------------------------------------------------- … … 67 68 !!---------------------------------------------------------------------- 68 69 69 INTEGER :: & !: 70 nstop = 0 , & !: e r r o r flag (=number of reason for a 71 ! ! prematurely stop the run) 72 nwarn = 0 !: w a r n i n g flag (=number of warning 73 ! ! found during the run) 74 75 76 CHARACTER(LEN=100) :: ctmp1, ctmp2, ctmp3 ! temporary character 77 CHARACTER (len=64) :: & !: 78 cform_err="(/,' ===>>> : E R R O R', /,' ===========',/)" , & !: 79 cform_war="(/,' ===>>> : W A R N I N G', /,' ===============',/)" !: 80 LOGICAL :: & !: 81 lwp , & !: boolean : true on the 1st processor only 82 lsp_area = .TRUE. !: to make a control print over a specific area 70 INTEGER :: nstop = 0 , & !: error flag (=number of reason for a premature stop run) 71 & nwarn = 0 !: warning flag (=number of warning found during the run) 72 CHARACTER(LEN=100) :: ctmp1, ctmp2, ctmp3 !: temporary character 73 CHARACTER (len=64) :: cform_err="(/,' ===>>> : E R R O R', /,' ===========',/)" , & !: 74 & cform_war="(/,' ===>>> : W A R N I N G', /,' ===============',/)" !: 75 LOGICAL :: lwp , & !: boolean : true on the 1st processor only 76 & lsp_area = .TRUE. !: to make a control print over a specific area 83 77 !!---------------------------------------------------------------------- 84 78 !! OPA 9.0 , LOCEAN-IPSL (2005) 85 79 !! $Header$ 86 !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt80 !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 87 81 !!---------------------------------------------------------------------- 88 82 89 90 83 CONTAINS 91 92 84 93 85 SUBROUTINE ctl_stop( cd1, cd2, cd3, cd4, cd5, & … … 96 88 !! *** ROUTINE stop_opa *** 97 89 !! 98 !! ** Purpose : ??? 99 !! 90 !! ** Purpose : ??? blah blah.... 100 91 !!----------------------------------------------------------------------- 101 CHARACTER(len=*),INTENT(in),OPTIONAL :: cd1, cd2, cd3, cd4, cd5, cd6, cd7, cd8, cd9, cd10 92 CHARACTER(len=*), INTENT(in), OPTIONAL :: cd1, cd2, cd3, cd4, cd5, & 93 & cd6, cd7, cd8, cd9, cd10 102 94 !!----------------------------------------------------------------------- 103 95 ! 104 96 nstop = nstop + 1 105 97 IF(lwp) THEN … … 117 109 ENDIF 118 110 CALL FLUSH(numout) 119 111 ! 120 112 END SUBROUTINE ctl_stop 121 113 … … 124 116 & cd6, cd7, cd8, cd9, cd10 ) 125 117 !!----------------------------------------------------------------------- 126 !! *** ROUTINE stop_ opa***118 !! *** ROUTINE stop_warn *** 127 119 !! 128 !! ** Purpose : ??? 129 !! 120 !! ** Purpose : ??? blah blah.... 130 121 !!----------------------------------------------------------------------- 131 CHARACTER(len=*),INTENT(in),OPTIONAL :: cd1, cd2, cd3, cd4, cd5, cd6, cd7, cd8, cd9, cd10 122 CHARACTER(len=*), INTENT(in), OPTIONAL :: cd1, cd2, cd3, cd4, cd5, & 123 & cd6, cd7, cd8, cd9, cd10 132 124 !!----------------------------------------------------------------------- 133 125 ! 134 126 nwarn = nwarn + 1 135 127 IF(lwp) THEN … … 147 139 ENDIF 148 140 CALL FLUSH(numout) 149 141 ! 150 142 END SUBROUTINE ctl_warn 151 143 144 !!===================================================================== 152 145 END MODULE in_out_manager -
trunk/NEMO/OPA_SRC/iom.F90
r485 r508 2 2 !!===================================================================== 3 3 !! *** MODULE iom *** 4 !!5 4 !! Input/Output manager : Library to read input files 6 !!7 !! Ongoing work : This code is here to help discussions about I/O8 !! library in the NEMO system9 5 !!==================================================================== 6 !! History : 9.0 ! 05 12 (J. Belier) Original code 7 !! 9.0 ! 06 02 (S. Masson) Adaptation to NEMO 8 !!-------------------------------------------------------------------- 9 !!gm caution add !DIR nec: improved performance to be checked as well as no result changes 10 10 11 !!-------------------------------------------------------------------- 11 12 !! iom_open : open a file read only 12 13 !! iom_close : close a file or all files opened by iom 13 !! iom_get : read a field : interface to several routines 14 !! iom_get : read a field (interfaced to several routines) 15 !! iom_gettime : read the time axis cdvar in the file !!gm : never call ?????? 14 16 !! iom_varid : get the id of a variable in a file 15 !! iom_ get_gblatt : ???17 !! iom_rstput : write a field in a restart file (interfaced to several routines) 16 18 !!-------------------------------------------------------------------- 17 !! History : 9.0 ! 05 12 (J. Belier) Original code18 !! 9.0 ! 06 02 (S. Masson) Adaptation to NEMO19 !!--------------------------------------------------------------------20 !! * Modules used21 19 USE in_out_manager ! I/O manager 22 20 USE dom_oce ! ocean space and time domain 23 USE lbclnk ! ???24 USE ioipsl ! ???21 USE lbclnk ! lateal boundary condition / mpp exchanges 22 USE ioipsl ! IOIPSL library 25 23 26 24 IMPLICIT NONE 27 25 PRIVATE 28 26 29 PUBLIC iom_open, iom_close, iom_get, iom_varid, iom_get_gblatt 30 31 !! * Interfaces 27 PUBLIC iom_open, iom_close, iom_get, iom_varid, iom_rstput, iom_gettime 28 32 29 INTERFACE iom_get 33 MODULE PROCEDURE iom_get_r_ 1d, iom_get_r_2d, iom_get_r_3d30 MODULE PROCEDURE iom_get_r_0d, iom_get_r_1d, iom_get_r_2d, iom_get_r_3d 34 31 END INTERFACE 35 36 !! * Share module variables 37 INTEGER, PARAMETER, PUBLIC :: & !: 38 jpdom_data = 1, & !: ( 1 :jpidta, 1 :jpjdta) 39 jpdom_global = 2, & !: ( 1 :jpiglo, 1 :jpjglo) 40 jpdom_local = 3, & !: One of the 3 following cases 41 jpdom_local_full = 4, & !: ( 1 :jpi , 1 :jpi ) 42 jpdom_local_noextra = 5, & !: ( 1 :nlci , 1 :nlcj ) 43 jpdom_local_noovlap = 6, & !: (nldi:nlei ,nldj:nlej ) 44 jpdom_unknown = 7 !: No dimension checking 45 46 !! * Module variables 47 INTEGER, PARAMETER :: & 48 jpmax_vars = 50, & ! maximum number of variables in one file 49 jpmax_dims = 5, & ! maximum number of dimensions for one variable 50 jpmax_digits = 5 ! maximum number of digits in the file name to reference the cpu number 51 32 INTERFACE iom_rstput 33 MODULE PROCEDURE iom_rstput_0d, iom_rstput_1d, iom_rstput_2d, iom_rstput_3d 34 END INTERFACE 35 36 INTEGER, PARAMETER, PUBLIC :: jpdom_data = 1 !: ( 1 :jpidta, 1 :jpjdta) 37 INTEGER, PARAMETER, PUBLIC :: jpdom_global = 2 !: ( 1 :jpiglo, 1 :jpjglo) 38 INTEGER, PARAMETER, PUBLIC :: jpdom_local = 3 !: One of the 3 following cases 39 INTEGER, PARAMETER, PUBLIC :: jpdom_local_full = 4 !: ( 1 :jpi , 1 :jpi ) 40 INTEGER, PARAMETER, PUBLIC :: jpdom_local_noextra = 5 !: ( 1 :nlci , 1 :nlcj ) 41 INTEGER, PARAMETER, PUBLIC :: jpdom_local_noovlap = 6 !: (nldi:nlei ,nldj:nlej ) 42 INTEGER, PARAMETER, PUBLIC :: jpdom_unknown = 7 !: No dimension checking 43 44 INTEGER, PARAMETER :: jpmax_vars = 60, & ! maximum number of variables in one file 45 & jpmax_dims = 4, & ! maximum number of dimensions for one variable 46 & jpmax_digits = 5 ! maximum number of digits for the cpu number in the file name 52 47 !$AGRIF_DO_NOT_TREAT 53 INTEGER :: iom_init = 0 54 55 TYPE :: flio_file 48 INTEGER :: iom_init = 0 49 TYPE :: flio_file 56 50 CHARACTER(LEN=240) :: name ! name of the file 57 INTEGER :: iopen ! 1 /0 is the file is open/not open51 INTEGER :: iopen ! 1(0) if the file is open(not open) 58 52 INTEGER :: nvars ! number of identified varibles in the file 59 53 INTEGER :: iduld ! id of the unlimited dimension 60 54 CHARACTER(LEN=16), DIMENSION(jpmax_vars) :: cn_var ! names of the variables 61 55 INTEGER, DIMENSION(jpmax_vars) :: ndims ! number of dimensions of the variables 62 LOGICAL, DIMENSION(jpmax_vars) :: luld ! variable includingunlimited dimension63 INTEGER, DIMENSION(jpmax_dims,jpmax_vars) :: dimsz ! size of the dimensions of the variables56 LOGICAL, DIMENSION(jpmax_vars) :: luld ! variable using the unlimited dimension 57 INTEGER, DIMENSION(jpmax_dims,jpmax_vars) :: dimsz ! size of variables dimensions 64 58 REAL(kind=wp), DIMENSION(jpmax_vars) :: scf ! scale_factor of the variables 65 59 REAL(kind=wp), DIMENSION(jpmax_vars) :: ofs ! add_offset of the variables 66 60 END TYPE flio_file 67 TYPE(flio_file), DIMENSION(flio_max_files) :: iom_file ! array containing the info for all opened files61 TYPE(flio_file), DIMENSION(flio_max_files) :: iom_file ! array containing the info for all opened files 68 62 !$AGRIF_END_DO_NOT_TREAT 69 63 … … 76 70 CONTAINS 77 71 78 SUBROUTINE iom_open( cdname, knumfl, ld img )72 SUBROUTINE iom_open( cdname, knumfl, ldwrt, kdom, ldimg ) 79 73 !!--------------------------------------------------------------------- 80 74 !! *** SUBROUTINE iom_open *** 81 75 !! 82 76 !! ** Purpose : open an input file read only (return 0 if not found) 83 !!84 !! ** Method :85 !!86 77 !!--------------------------------------------------------------------- 87 CHARACTER(len=*), INTENT(in ) :: cdname ! File name 88 INTEGER, INTENT(out) :: knumfl ! Identifier of the opened file 89 LOGICAL, INTENT(in ), OPTIONAL :: ldimg ! flg to specify that we use dimg format 90 91 CHARACTER(LEN=100) :: clname ! the name of the file based on cdname [[+clcpu]+clcpu] 92 CHARACTER(LEN=10) :: clsuffix ! ".nc" or ".dimg" 93 CHARACTER(LEN=10) :: clcpu ! the cpu number (max jpmax_digits digits) 94 LOGICAL :: llok ! check the existence 95 INTEGER :: icnt ! counter for digits in clcpu (max = jpmax_digits) 78 CHARACTER(len=*), INTENT(in ) :: cdname ! File name 79 INTEGER , INTENT( out) :: knumfl ! Identifier of the opened file 80 LOGICAL , INTENT(in ), OPTIONAL :: ldwrt ! read or write the file? 81 INTEGER , INTENT(in ), OPTIONAL :: kdom ! Type of domain to be written 82 LOGICAL , INTENT(in ), OPTIONAL :: ldimg ! use dimg format? 83 84 CHARACTER(LEN=100) :: clname ! the name of the file based on cdname [[+clcpu]+clcpu] 85 CHARACTER(LEN=100) :: cltmpn ! tempory name to store clname (in writting mode) 86 CHARACTER(LEN=10) :: clsuffix ! ".nc" or ".dimg" 87 CHARACTER(LEN=10) :: clcpu ! the cpu number (max jpmax_digits digits) 88 LOGICAL :: llok ! check the existence 89 LOGICAL :: llwrt ! 90 INTEGER :: icnt ! counter for digits in clcpu (max = jpmax_digits) 91 INTEGER :: iln, ils ! lengths of character 92 INTEGER :: idom ! type of domain 93 INTEGER :: ifliodom ! model domain identifier (see flio_dom_set) 94 INTEGER, DIMENSION(2) :: iszl ! local number of points for x,y dimensions 95 INTEGER, DIMENSION(2) :: ifst ! position of first local point for x,y dimensions 96 INTEGER, DIMENSION(2) :: ilst ! position of last local point for x,y dimensions 97 INTEGER, DIMENSION(2) :: ihst ! start halo size for x,y dimensions 98 INTEGER, DIMENSION(2) :: ihnd ! end halo size for x,y dimensions 96 99 !--------------------------------------------------------------------- 97 98 ! find the file 100 ! if iom_open is called for the first time: initialize iom_file(:)%iopen to 0 101 ! (could be done when defining iom_file in f95 but not in f90) 102 IF( iom_init == 0 ) THEN 103 iom_file(:)%iopen = 0 104 iom_init = 1 105 ENDIF 106 ! do we read or write the file? 107 ! ============= 108 IF( PRESENT(ldwrt) ) THEN ; llwrt = ldwrt 109 ELSE ; llwrt = .FALSE. 110 ENDIF 111 ! create the file name by added, if needed, TRIM(Agrif_CFixed()) and TRIM(clsuffix) 99 112 ! ============= 100 113 clname = trim(cdname) 101 114 #if defined key_agrif 102 115 if ( .NOT. Agrif_Root() ) clname = TRIM(Agrif_CFixed())//'_'//TRIM(clname) 103 #endif 116 #endif 117 ! which suffix should we use? 104 118 clsuffix = ".nc" 105 IF( PRESENT(ldimg) ) THEN 106 IF ( ldimg ) clsuffix = ".dimg" 107 ENDIF 108 ! 119 IF( PRESENT(ldimg) ) THEN ; IF( ldimg ) clsuffix = ".dimg" ; ENDIF 120 ! Add the suffix if needed 121 iln = LEN_TRIM(clname) 122 ils = LEN_TRIM(clsuffix) 123 IF( iln <= ils) clname = clname(1:iln)//TRIM(clsuffix) 124 IF( clname(iln-ils+1:iln) /= TRIM(clsuffix) ) clname = clname(1:iln)//TRIM(clsuffix) 125 cltmpn = clname ! store this name 126 ! try to find if the file to be opened already exist 109 127 INQUIRE( FILE = clname, EXIST = llok ) 110 IF( .NOT.llok ) THEN ! try to complete the name with the suffix only 111 clname = TRIM(cdname)//TRIM(clsuffix) 112 #if defined key_agrif 113 if ( .NOT. Agrif_Root() ) clname = TRIM(Agrif_CFixed())//'_'//TRIM(clname) 114 #endif 115 INQUIRE( FILE = clname, EXIST = llok ) 116 IF( .NOT.llok ) THEN ! try to complete the name with both cpu number and suffix 117 WRITE(clcpu,*) narea-1 118 clcpu = trim(adjustl(clcpu)) 119 clname = trim(cdname)//"_" 120 #if defined key_agrif 121 if ( .NOT. Agrif_Root() ) clname = TRIM(Agrif_CFixed())//'_'//TRIM(clname) 122 #endif 123 icnt = 0 124 INQUIRE( FILE = trim(clname)//trim(clcpu)//trim(clsuffix), EXIST = llok ) 125 DO WHILE( .NOT.llok .AND. icnt < jpmax_digits ) ! we try fifferent formats for the cpu number by adding 0 126 clname = trim(clname)//"0" 127 INQUIRE( FILE = trim(clname)//trim(clcpu)//trim(clsuffix), EXIST = llok ) 128 icnt = icnt + 1 129 END DO 130 IF( .NOT.llok ) THEN ! no way to find the files... 131 CALL ctl_stop( 'iom_open: file '//trim(clname)//'... not found' ) 128 IF( .NOT.llok ) THEN 129 ! we try to add the cpu number to the name 130 WRITE(clcpu,*) narea-1 131 clcpu = TRIM(ADJUSTL(clcpu)) 132 iln = INDEX(clname,TRIM(clsuffix)) 133 clname = clname(1:iln-1)//'_'//TRIM(clcpu)//TRIM(clsuffix) 134 icnt = 0 135 INQUIRE( FILE = clname, EXIST = llok ) 136 ! we try different formats for the cpu number by adding 0 137 DO WHILE( .NOT.llok .AND. icnt < jpmax_digits ) 138 clcpu = "0"//trim(clcpu) 139 clname = clname(1:iln-1)//'_'//TRIM(clcpu)//TRIM(clsuffix) 140 INQUIRE( FILE = clname, EXIST = llok ) 141 icnt = icnt + 1 142 END DO 143 ENDIF 144 ! 145 IF( llok ) THEN ! Open the file 146 ! ! ============= 147 IF( llwrt ) THEN 148 IF(lwp) WRITE(numout,*) ' iom_open ~~~ open existing file: '//TRIM(clname)//' in WRITE mode' 149 CALL flioopfd( TRIM(clname), knumfl, "WRITE" ) 150 ELSE 151 IF(lwp) WRITE(numout,*) ' iom_open ~~~ open existing file: '//TRIM(clname)//' in READ mode' 152 CALL flioopfd( TRIM(clname), knumfl ) 153 ENDIF 154 ELSE ! no way to find the file... 155 ! ! ======================= 156 IF( llwrt ) THEN 157 ! file opened in write mode 158 ! the file does not exist, we must create it... 159 ! ============= 160 llok = .TRUE. 161 ! on which domain must the file be written?? 162 ! check the domain definition 163 idom = jpdom_local_noovlap ! default definition 164 IF( PRESENT(kdom) ) idom = kdom 165 ! create the domain informations 166 ! ============= 167 SELECT CASE (idom) 168 CASE (jpdom_local_full) 169 iszl = (/ jpi , jpj /) 170 ifst = (/ nimpp , njmpp /) 171 ilst = (/ nimpp + jpi - 1 , njmpp + jpj - 1 /) 172 ihst = (/ nldi - 1 , nldj - 1 /) 173 ihnd = (/ jpi - nlei , jpj - nlej /) 174 CASE (jpdom_local_noextra) 175 iszl = (/ nlci , nlcj /) 176 ifst = (/ nimpp , njmpp /) 177 ilst = (/ nimpp + nlci - 1, njmpp + nlcj - 1 /) 178 ihst = (/ nldi - 1 , nldj - 1 /) 179 ihnd = (/ nlci - nlei , nlcj - nlej /) 180 CASE (jpdom_local_noovlap) 181 iszl = (/ nlei - nldi + 1 , nlej - nldj + 1 /) 182 ifst = (/ nimpp + nldi - 1, njmpp + nldj - 1 /) 183 ilst = (/ nimpp + nlei - 1, njmpp + nlej - 1 /) 184 ihst = (/ 0 , 0 /) 185 ihnd = (/ 0 , 0 /) 186 CASE DEFAULT 187 llok = .FALSE. 188 CALL ctl_stop( 'iom_open: wrong value of kdom, only jpdom_local* cases are accepted' ) 189 END SELECT 190 IF( llok ) THEN 191 CALL flio_dom_set( jpnij, narea-1, (/1, 2/), (/jpiglo, jpjglo/) & 192 & , iszl, ifst, ilst, ihst, ihnd, 'BOX', ifliodom ) 193 ! create the file 194 ! ============= 195 ! Note that fliocrfd may change the value of clname (add the cpu number...) 196 clname = cltmpn ! get back the file name without the cpu number in it 197 IF(lwp) WRITE(numout,*) ' iom_open ~~~ create new file: '//trim(clname)//' in WRITE mode' 198 CALL fliocrfd( clname, (/'x' , 'y' , 'z', 't'/) & 199 & , (/iszl(1), iszl(2), jpk, -1 /) & 200 & , knumfl, ifliodom ) 132 201 ENDIF 133 clname = trim(clname)//trim(clcpu)//trim(clsuffix) 202 ELSE 203 ! the file is open for read-only, it must exist... 204 iln = INDEX( cltmpn,TRIM(clsuffix) ) 205 CALL ctl_stop( 'iom_open: file '//cltmpn(1:iln-1)//'* not found' ) 134 206 ENDIF 135 207 ENDIF 136 137 ! Open the file 208 ! start to fill the information of opened files 138 209 ! ============= 139 210 IF( llok ) THEN 140 IF (lwp) WRITE(numout,*) 'iom_open ~~~ open file: '//trim(clname)141 CALL flioopfd( trim(clname), knumfl )142 IF( iom_init == 0 ) THEN143 iom_file(:)%iopen = 0144 iom_init = 1145 ENDIF146 211 iom_file(knumfl)%iopen = 1 147 212 iom_file(knumfl)%name = TRIM(clname) … … 152 217 ! does the file contain time axis (that must be unlimitted) ? 153 218 CALL flioinqf( knumfl, id_uld = iom_file(knumfl)%iduld ) 219 IF(lwp) WRITE(numout,*) ' ---> OK' 154 220 ELSE 155 knumfl = 0 156 ENDIF 157 221 knumfl = 0 ! return error flag 222 ENDIF 223 ! 158 224 END SUBROUTINE iom_open 159 225 … … 164 230 !! 165 231 !! ** Purpose : close an input file, or all files opened by iom 166 !!167 !! ** Method :168 !!169 232 !!-------------------------------------------------------------------- 170 INTEGER,INTENT(in), OPTIONAL :: knumfl ! Identifier of the file to be closed 171 ! ! If this argument is not present, 172 ! ! all the files opened by iom are closed. 173 174 INTEGER :: jf ! dummy loop indices 175 INTEGER :: i_s, i_e ! temporary integer 233 INTEGER, INTENT(in), OPTIONAL :: knumfl ! Identifier of the file to be closed 234 ! ! No argument : all the files opened by iom are closed 235 236 INTEGER :: jf ! dummy loop indices 237 INTEGER :: i_s, i_e ! temporary integer 176 238 !--------------------------------------------------------------------- 177 239 ! 178 240 IF( PRESENT(knumfl) ) THEN 179 241 i_s = knumfl … … 183 245 i_e = flio_max_files 184 246 ENDIF 185 IF ( i_s > 0 ) THEN 247 248 IF( i_s > 0 ) THEN 186 249 DO jf = i_s, i_e 187 250 IF( iom_file(jf)%iopen > 0 ) THEN 188 251 CALL flioclo( jf ) 252 IF(lwp) WRITE(numout,*) ' iom_close, close file: '//TRIM(iom_file(knumfl)%name)//' ok' 189 253 iom_file(jf)%iopen = 0 190 254 iom_file(jf)%name = 'NONE' … … 200 264 END DO 201 265 ENDIF 202 266 ! 203 267 END SUBROUTINE iom_close 204 268 205 269 206 270 !!---------------------------------------------------------------------- 207 !! INTERFACE iom_ u_getv271 !! INTERFACE iom_get_123d 208 272 !!---------------------------------------------------------------------- 209 SUBROUTINE iom_get_r_1d( knumfl, kdom , cdvar , pvar , & 210 & ktime, kstart, kcount ) 211 INTEGER , INTENT(in ) :: knumfl ! Identifier of the file 212 INTEGER , INTENT(in ) :: kdom ! Type of domain to be read 213 CHARACTER(len=*) , INTENT(in ) :: cdvar ! Name of the variable 214 REAL(wp), DIMENSION(:), INTENT(out) :: pvar ! read field 215 INTEGER , INTENT(in ) ,OPTIONAL :: ktime ! record number 216 INTEGER , DIMENSION(:), INTENT(in ), OPTIONAL :: kstart ! start position of the reading in each axis 217 INTEGER , DIMENSION(:), INTENT(in ), OPTIONAL :: kcount ! number of points to be read in each axis 218 ! 219 CHARACTER(LEN=100) :: clinfo ! info character 220 ! 221 clinfo = 'iom_get_r_1d, file: '//trim(iom_file(knumfl)%name)//', var: '//trim(cdvar) 222 IF( PRESENT(kstart) ) THEN 223 IF ( SIZE(kstart) /= 1 ) CALL ctl_stop( trim(clinfo), 'kstart must be a 1 element vector' ) 224 ENDIF 225 IF( PRESENT(kcount) ) THEN 226 IF ( SIZE(kcount) /= 1 ) CALL ctl_stop( trim(clinfo), 'kcount must be a 1 element vector' ) 227 ENDIF 228 IF ( knumfl > 0 ) CALL iom_u_getv( knumfl, kdom , cdvar , pv_r1d=pvar, & 229 & ktime=ktime, kstart=kstart, kcount=kcount ) 273 SUBROUTINE iom_get_r_0d( knumfl, cdvar, pvar ) 274 INTEGER , INTENT(in ) :: knumfl ! Identifier of the file 275 CHARACTER(len=*), INTENT(in ) :: cdvar ! Name of the variable 276 REAL(wp) , INTENT( out) :: pvar ! read field 277 ! 278 IF( knumfl > 0 .AND. iom_varid( knumfl, cdvar ) > 0 ) CALL fliogetv( knumfl, cdvar, pvar ) 279 END SUBROUTINE iom_get_r_0d 280 281 SUBROUTINE iom_get_r_1d( knumfl, kdom, cdvar, pvar, ktime, kstart, kcount ) 282 INTEGER , INTENT(in ) :: knumfl ! Identifier of the file 283 INTEGER , INTENT(in ) :: kdom ! Type of domain to be read 284 CHARACTER(len=*), INTENT(in ) :: cdvar ! Name of the variable 285 REAL(wp) , INTENT( out), DIMENSION(:) :: pvar ! read field 286 INTEGER , INTENT(in ) , OPTIONAL :: ktime ! record number 287 INTEGER , INTENT(in ), DIMENSION(1), OPTIONAL :: kstart ! start axis position of the reading 288 INTEGER , INTENT(in ), DIMENSION(1), OPTIONAL :: kcount ! number of points in each axis 289 ! 290 IF( knumfl > 0 ) CALL iom_get_123d( knumfl, kdom , cdvar , pv_r1d=pvar, & 291 & ktime=ktime, kstart=kstart, kcount=kcount ) 230 292 END SUBROUTINE iom_get_r_1d 231 SUBROUTINE iom_get_r_2d( knumfl, kdom , cdvar , pvar , & 232 & ktime, kstart, kcount ) 233 INTEGER,INTENT(in) :: knumfl 234 INTEGER,INTENT(in) :: kdom 235 CHARACTER(len=*),INTENT(in) :: cdvar 236 REAL(wp),INTENT(out),DIMENSION(:,:) :: pvar 237 INTEGER,INTENT(in),OPTIONAL :: ktime 238 INTEGER,DIMENSION(:),INTENT(in),OPTIONAL :: kstart 239 INTEGER,DIMENSION(:),INTENT(in),OPTIONAL :: kcount 240 ! 241 CHARACTER(LEN=100) :: clinfo ! info character 242 ! 243 clinfo = 'iom_get_r_2d, file: '//trim(iom_file(knumfl)%name)//', var: '//trim(cdvar) 244 IF( PRESENT(kstart) ) THEN 245 IF ( size(kstart) /= 2 ) CALL ctl_stop(trim(clinfo), 'kstart must be a 2 element vector') 246 ENDIF 247 IF( PRESENT(kcount) ) THEN 248 IF ( size(kcount) /= 2 ) CALL ctl_stop(trim(clinfo), 'kcount must be a 2 element vector') 249 ENDIF 250 IF ( knumfl > 0 ) CALL iom_u_getv( knumfl, kdom , cdvar , pv_r2d=pvar, & 251 & ktime=ktime, kstart=kstart, kcount=kcount ) 293 294 SUBROUTINE iom_get_r_2d( knumfl, kdom, cdvar, pvar, ktime, kstart, kcount ) 295 INTEGER , INTENT(in ) :: knumfl ! Identifier of the file 296 INTEGER , INTENT(in ) :: kdom ! Type of domain to be read 297 CHARACTER(len=*), INTENT(in ) :: cdvar ! Name of the variable 298 REAL(wp) , INTENT( out), DIMENSION(:,:) :: pvar ! read field 299 INTEGER , INTENT(in ) , OPTIONAL :: ktime ! record number 300 INTEGER , INTENT(in ), DIMENSION(2) , OPTIONAL :: kstart ! start axis position of the reading 301 INTEGER , INTENT(in ), DIMENSION(2) , OPTIONAL :: kcount ! number of points in each axis 302 ! 303 IF( knumfl > 0 ) CALL iom_get_123d( knumfl, kdom , cdvar , pv_r2d=pvar, & 304 & ktime=ktime, kstart=kstart, kcount=kcount ) 252 305 END SUBROUTINE iom_get_r_2d 253 SUBROUTINE iom_get_r_3d( knumfl, kdom , cdvar , pvar , & 254 & ktime, kstart, kcount ) 255 INTEGER,INTENT(in) :: knumfl 256 INTEGER,INTENT(in) :: kdom 257 CHARACTER(len=*),INTENT(in) :: cdvar 258 REAL(wp),INTENT(out),DIMENSION(:,:,:) :: pvar 259 INTEGER,INTENT(in),OPTIONAL :: ktime 260 INTEGER,DIMENSION(:),INTENT(in),OPTIONAL :: kstart 261 INTEGER,DIMENSION(:),INTENT(in),OPTIONAL :: kcount 262 ! 263 CHARACTER(LEN=100) :: clinfo ! info character 264 ! 265 clinfo = 'iom_get_r_3d, file: '//trim(iom_file(knumfl)%name)//', var: '//trim(cdvar) 266 IF ( PRESENT(kstart) ) THEN 267 IF ( size(kstart) /= 3 ) CALL ctl_stop(trim(clinfo), 'kstart must be a 3 element vector') 268 ENDIF 269 IF ( PRESENT(kcount) ) THEN 270 IF ( size(kcount) /= 3 ) CALL ctl_stop(trim(clinfo), 'kcount must be a 3 element vector') 271 ENDIF 272 IF ( knumfl > 0 ) CALL iom_u_getv( knumfl, kdom , cdvar , pv_r3d=pvar, & 273 & ktime=ktime, kstart=kstart, kcount=kcount ) 306 307 SUBROUTINE iom_get_r_3d( knumfl, kdom, cdvar, pvar, ktime, kstart, kcount ) 308 INTEGER , INTENT(in ) :: knumfl ! Identifier of the file 309 INTEGER , INTENT(in ) :: kdom ! Type of domain to be read 310 CHARACTER(len=*), INTENT(in ) :: cdvar ! Name of the variable 311 REAL(wp) , INTENT( out), DIMENSION(:,:,:) :: pvar ! read field 312 INTEGER , INTENT(in ) , OPTIONAL :: ktime ! record number 313 INTEGER , INTENT(in ), DIMENSION(3) , OPTIONAL :: kstart ! start axis position of the reading 314 INTEGER , INTENT(in ), DIMENSION(3) , OPTIONAL :: kcount ! number of points in each axis 315 ! 316 IF( knumfl > 0 ) CALL iom_get_123d( knumfl, kdom , cdvar , pv_r3d=pvar, & 317 & ktime=ktime, kstart=kstart, kcount=kcount ) 274 318 END SUBROUTINE iom_get_r_3d 275 319 !!---------------------------------------------------------------------- 276 320 277 278 SUBROUTINE iom_u_getv( knumfl, kdom , cdvar , & 279 & pv_r1d, pv_r2d, pv_r3d, & 280 & ktime , kstart, kcount ) 321 SUBROUTINE iom_get_123d( knumfl, kdom , cdvar , & 322 & pv_r1d, pv_r2d, pv_r3d, & 323 & ktime , kstart, kcount ) 281 324 !!----------------------------------------------------------------------- 282 !! *** ROUTINE iom_ u_getv***325 !! *** ROUTINE iom_get_123d *** 283 326 !! 284 327 !! ** Purpose : read a 1D/2D/3D variable 285 328 !! 286 !! ** Method : read ONE time step at each CALL 287 !! 329 !! ** Method : read ONE record at each CALL 288 330 !!----------------------------------------------------------------------- 289 INTEGER , INTENT(in) :: knumfl ! Identifier of the file290 INTEGER , INTENT(in) :: kdom ! Type of domain to be read291 CHARACTER(len=*) , INTENT(in) :: cdvar ! Name of the variable292 REAL(wp), DIMENSION(:) , INTENT(out), OPTIONAL :: pv_r1d ! read field (1D case)293 REAL(wp), DIMENSION(:,:) , INTENT(out), OPTIONAL :: pv_r2d ! read field (2D case)294 REAL(wp), DIMENSION(:,:,:) , INTENT(out), OPTIONAL :: pv_r3d ! read field (3D case)295 INTEGER , INTENT(in), OPTIONAL :: ktime ! record number296 INTEGER , DIMENSION(:) , INTENT(in), OPTIONAL :: kstart ! start position of the reading in each axis297 INTEGER , DIMENSION(:) , INTENT(in), OPTIONAL :: kcount ! number of points to be read in each axis298 299 INTEGER :: jl! loop on number of dimension300 INTEGER :: idom, &! type of domain301 & idvar, &! id of the variable302 & inbdim, &! number of dimensions of the variable303 & idmspc, &! number of spatial dimensions304 & itime, &! record number305 & istop! temporary value of nstop306 INTEGER, DIMENSION(jpmax_dims) :: istart, &! starting point to read for each axis307 & icnt, &! number of value to read along each axis308 & idimsz! size of the dimensions of the variable309 REAL(wp) :: zscf, zofs! sacle_factor and add_offset310 INTEGER :: itmp! temporary integer311 CHARACTER(LEN=100) :: clinfo! info character331 INTEGER , INTENT(in ) :: knumfl ! Identifier of the file 332 INTEGER , INTENT(in ) :: kdom ! Type of domain to be read 333 CHARACTER(len=*) , INTENT(in ) :: cdvar ! Name of the variable 334 REAL(wp), DIMENSION(:) , INTENT( out), OPTIONAL :: pv_r1d ! read field (1D case) 335 REAL(wp), DIMENSION(:,:) , INTENT( out), OPTIONAL :: pv_r2d ! read field (2D case) 336 REAL(wp), DIMENSION(:,:,:) , INTENT( out), OPTIONAL :: pv_r3d ! read field (3D case) 337 INTEGER , INTENT(in ), OPTIONAL :: ktime ! record number 338 INTEGER , DIMENSION(:) , INTENT(in ), OPTIONAL :: kstart ! start position of the reading in each axis 339 INTEGER , DIMENSION(:) , INTENT(in ), OPTIONAL :: kcount ! number of points to be read in each axis 340 ! 341 INTEGER :: jl ! loop on number of dimension 342 INTEGER :: idom, & ! type of domain 343 & idvar, & ! id of the variable 344 & inbdim, & ! number of dimensions of the variable 345 & idmspc, & ! number of spatial dimensions 346 & itime, & ! record number 347 & istop ! temporary value of nstop 348 INTEGER, DIMENSION(jpmax_dims) :: istart, & ! starting point to read for each axis 349 & icnt, & ! number of value to read along each axis 350 & idimsz ! size of the dimensions of the variable 351 REAL(wp) :: zscf, zofs ! sacle_factor and add_offset 352 INTEGER :: itmp ! temporary integer 353 CHARACTER(LEN=100) :: clinfo ! info character 312 354 !--------------------------------------------------------------------- 313 clinfo = 'iom_u_getv, file: '//trim(iom_file(knumfl)%name)//', var: '//trim(cdvar) 355 ! 356 clinfo = ' iom_get_123d, file: '//trim(iom_file(knumfl)%name)//', var: '//trim(cdvar) 314 357 ! local definition of the domain ? 315 358 idom = kdom 316 359 ! check kcount and kstart optionals parameters... 317 IF( PRESENT(kcount) .AND. (.NOT. PRESENT(kstart)) ) &360 IF( PRESENT(kcount) .AND. (.NOT. PRESENT(kstart)) ) & 318 361 CALL ctl_stop( trim(clinfo), 'kcount present needs kstart present' ) 319 IF( PRESENT(kstart) .AND. (.NOT. PRESENT(kcount)) ) &362 IF( PRESENT(kstart) .AND. (.NOT. PRESENT(kcount)) ) & 320 363 CALL ctl_stop( trim(clinfo), 'kstart present needs kcount present' ) 321 IF( PRESENT(kstart) .AND. idom /= jpdom_unknown )&364 IF( PRESENT(kstart) .AND. idom /= jpdom_unknown ) & 322 365 CALL ctl_stop( trim(clinfo), 'kstart present needs kdom = jpdom_unknown' ) 323 366 324 367 ! Search for the variable in the data base (eventually actualize data) 325 !-326 368 istop = nstop 327 369 idvar = iom_varid( knumfl, cdvar ) 328 370 ! 329 IF 371 IF( idvar > 0 ) THEN 330 372 ! to write iom_file(knumfl)%dimsz in a shorter way ! 331 373 idimsz(:) = iom_file(knumfl)%dimsz(:, idvar) 332 inbdim = iom_file(knumfl)%ndims(idvar)! number of dimensions in the file 333 idmspc = inbdim ! number of spatial dimensions in the file 334 IF( iom_file(knumfl)%luld(idvar) ) idmspc = inbdim - 1 335 IF( idmspc > 3 ) CALL ctl_stop(trim(clinfo), 'the file has more than 3 spatial dimensions', & 336 & 'this case is not coded...') 337 ! Identify the domain in case of jpdom_local definition 338 !- 339 IF( idom == jpdom_local ) THEN 374 inbdim = iom_file(knumfl)%ndims(idvar) ! number of dimensions in the file 375 idmspc = inbdim ! number of spatial dimensions in the file 376 IF( iom_file(knumfl)%luld(idvar) ) idmspc = inbdim - 1 377 IF( idmspc > 3 ) CALL ctl_stop(trim(clinfo), & 378 & 'the file has more than 3 spatial dimensions this case is not coded...' ) 379 IF( idom == jpdom_local ) THEN ! Identify the domain in case of jpdom_local definition 340 380 IF( idimsz(1) == jpi .AND. idimsz(2) == jpj ) THEN 341 381 idom = jpdom_local_full … … 348 388 ENDIF 349 389 ENDIF 350 390 ! 351 391 ! definition of istart and icnt 352 ! -392 ! 353 393 ! initializations 354 394 istart(:) = 1 355 icnt (:) = 1395 icnt (:) = 1 356 396 itime = 1 357 397 IF( PRESENT(ktime) ) itime = ktime … … 383 423 CASE (2) 384 424 ! data is 2d array (+ maybe a temporal dimension) 385 IF 425 IF( PRESENT(kstart) ) THEN 386 426 istart(1:3) = (/ kstart(1:2), itime /) 387 427 icnt(1:2) = kcount(1:2) … … 404 444 ENDIF 405 445 CASE DEFAULT 406 IF 407 CALL ctl_warn( trim(clinfo), '2D array but 3 spatial dimensions for the data...', &408 & 'As the size of the z dimension is 1 and as we try to read the first re acord, ',&409 & 'we accept this case even if there is a possible mix-up between z and time dimension ...')410 IF 446 IF( itime == 1 .AND. idimsz(3) == 1 .AND. idmspc == 3 ) THEN 447 CALL ctl_warn( trim(clinfo), '2D array but 3 spatial dimensions for the data...', & 448 & 'As the size of the z dimension is 1 and as we try to read the first record, ', & 449 & 'we accept this case even if there is a possible mix-up between z and time dimension' ) 450 IF( PRESENT(kstart) ) THEN 411 451 istart(1:2) = kstart(1:2) 412 452 icnt(1:2) = kcount(1:2) … … 428 468 ENDIF 429 469 ELSE 430 CALL ctl_stop( trim(clinfo), 'to keep iom lisibility, when reading a 2D array,', &431 & 'we do not accept data with more than 2 spatial dimension',&432 & 'Use ncwa -a to suppress the unnecessary dimensions')470 CALL ctl_stop( trim(clinfo), 'to keep iom lisibility, when reading a 2D array,', & 471 & 'we do not accept data with more than 2 spatial dimension', & 472 & 'Use ncwa -a to suppress the unnecessary dimensions' ) 433 473 ENDIF 434 474 END SELECT 435 475 ELSEIF( PRESENT(pv_r3d) ) THEN 436 476 SELECT CASE (idmspc) 437 CASE (1)438 CALL ctl_stop( trim(clinfo), 'the file has only 1 spatial dimension',&439 & 'it is impossible to read a 3d array from this file...')440 CASE (2)441 CALL ctl_stop( trim(clinfo), 'the file has only 2 spatial dimension',&442 & 'it is impossible to read a 3d array from this file...')443 CASE (3)477 CASE( 1 ) 478 CALL ctl_stop( trim(clinfo), 'the file has only 1 spatial dimension', & 479 & 'it is impossible to read a 3d array from this file...' ) 480 CASE( 2 ) 481 CALL ctl_stop( trim(clinfo), 'the file has only 2 spatial dimension', & 482 & 'it is impossible to read a 3d array from this file...' ) 483 CASE( 3 ) 444 484 ! data is 3d array (+ maybe a temporal dimension) 445 485 IF( PRESENT(kstart) ) THEN … … 469 509 ENDIF 470 510 CASE DEFAULT 471 CALL ctl_stop( trim(clinfo), 'to keep iom lisibility, when reading a 3D array,', &472 & 'we do not accept data with more than 3 spatial dimension', &473 & 'Use ncwa -a to suppress the unnecessary dimensions' )511 CALL ctl_stop( trim(clinfo), 'to keep iom lisibility, when reading a 3D array,', & 512 & 'we do not accept data with more than 3 spatial dimension', & 513 & 'Use ncwa -a to suppress the unnecessary dimensions' ) 474 514 END SELECT 475 515 ENDIF … … 491 531 itmp = size(pv_r1d) 492 532 WRITE(ctmp1,*) 'size(pv_r1d): ', itmp, ' /= icnt(1): ', icnt(1) 493 IF( itmp /= icnt(1) ) CALL ctl_stop( trim(clinfo), ctmp1 )533 IF( itmp /= icnt(1) ) CALL ctl_stop( trim(clinfo), ctmp1 ) 494 534 ELSEIF( PRESENT(pv_r2d) ) THEN 495 535 DO jl = 1, 2 … … 501 541 WRITE(ctmp1,*) 'size(pv_r2d(nldi:nlei,nldj:nlej), ',jl,'): ',itmp,' /= icnt(',jl,'): ',icnt(jl) 502 542 ENDIF 503 IF( itmp /= icnt(jl) ) CALL ctl_stop( trim(clinfo), ctmp1 )543 IF( itmp /= icnt(jl) ) CALL ctl_stop( trim(clinfo), ctmp1 ) 504 544 END DO 505 545 ELSEIF( PRESENT(pv_r3d) ) THEN … … 512 552 WRITE(ctmp1,*) 'size(pv_r3d(nldi:nlei,nldj:nlej), ',jl,'): ',itmp,' /= icnt(',jl,'): ',icnt(jl) 513 553 ENDIF 514 IF( itmp /= icnt(jl) ) CALL ctl_stop( trim(clinfo), ctmp1 )554 IF( itmp /= icnt(jl) ) CALL ctl_stop( trim(clinfo), ctmp1 ) 515 555 END DO 516 556 ENDIF … … 520 560 !- 521 561 IF( istop == nstop) THEN ! no additional errors until this point... 522 !523 istop = nstop524 562 ! 525 563 zscf = iom_file(knumfl)%scf(idvar) ! scale factor … … 529 567 CALL fliogetv( knumfl, cdvar, pv_r1d(:), start=istart(1:inbdim), count=icnt(1:inbdim) ) 530 568 !--- Apply scale_factor and offset 531 IF( zscf /= 1. ) pv_r1d(:) = pv_r1d(:) * zscf532 IF( zofs /= 0. ) pv_r1d(:) = pv_r1d(:) + zofs569 IF( zscf /= 1. ) pv_r1d(:) = pv_r1d(:) * zscf 570 IF( zofs /= 0. ) pv_r1d(:) = pv_r1d(:) + zofs 533 571 ELSEIF( PRESENT(pv_r2d) ) THEN 534 572 IF( idom /= jpdom_unknown ) THEN 535 573 CALL fliogetv( knumfl, cdvar, pv_r2d(nldi:nlei,nldj:nlej), start=istart(1:inbdim), count=icnt(1:inbdim) ) 536 574 !--- Apply scale_factor and offset 537 IF (zscf /= 1.) pv_r2d(nldi:nlei, nldj:nlej) = pv_r2d(nldi:nlei,nldj:nlej) * zscf 538 IF (zofs /= 0.) pv_r2d(nldi:nlei, nldj:nlej) = pv_r2d(nldi:nlei,nldj:nlej) + zofs 575 !CDIR NOUNROLL 576 IF( zscf /= 1.) pv_r2d(nldi:nlei, nldj:nlej) = pv_r2d(nldi:nlei,nldj:nlej) * zscf 577 !CDIR NOUNROLL 578 IF( zofs /= 0.) pv_r2d(nldi:nlei, nldj:nlej) = pv_r2d(nldi:nlei,nldj:nlej) + zofs 539 579 !--- Fill the overlap areas and extra hallows (mpp) 540 580 CALL lbc_lnk (pv_r2d,'Z',-999.,'no0') … … 542 582 CALL fliogetv( knumfl, cdvar, pv_r2d(:,:), start=istart(1:inbdim), count=icnt(1:inbdim) ) 543 583 !--- Apply scale_factor and offset 544 IF (zscf /= 1.) pv_r2d(:,:) = pv_r2d(:,:) * zscf 545 IF (zofs /= 0.) pv_r2d(:,:) = pv_r2d(:,:) + zofs 584 !CDIR COLLAPSE 585 IF( zscf /= 1.) pv_r2d(:,:) = pv_r2d(:,:) * zscf 586 !CDIR COLLAPSE 587 IF( zofs /= 0.) pv_r2d(:,:) = pv_r2d(:,:) + zofs 546 588 ENDIF 547 589 ELSEIF( PRESENT(pv_r3d) ) THEN 548 590 IF( idom /= jpdom_unknown ) THEN 549 CALL fliogetv( knumfl, cdvar, pv_r3d(nldi:nlei,nldj:nlej,:), start=istart(1:inbdim), count=icnt(1:inbdim) ) 591 CALL fliogetv( knumfl, cdvar, pv_r3d(nldi:nlei,nldj:nlej,:), start=istart(1:inbdim), & 592 & count=icnt (1:inbdim) ) 550 593 !--- Apply scale_factor and offset 551 IF( zscf /= 1. ) pv_r3d(nldi:nlei,nldj:nlej,:) = pv_r3d(nldi:nlei,nldj:nlej,:) * zscf 552 IF( zofs /= 0. ) pv_r3d(nldi:nlei,nldj:nlej,:) = pv_r3d(nldi:nlei,nldj:nlej,:) + zofs 594 !CDIR NOUNROLL 595 IF( zscf /= 1. ) pv_r3d(nldi:nlei,nldj:nlej,:) = pv_r3d(nldi:nlei,nldj:nlej,:) * zscf 596 !CDIR NOUNROLL 597 IF( zofs /= 0. ) pv_r3d(nldi:nlei,nldj:nlej,:) = pv_r3d(nldi:nlei,nldj:nlej,:) + zofs 553 598 !--- Fill the overlap areas and extra hallows (mpp) 554 599 IF( icnt(3) == jpk ) CALL lbc_lnk( pv_r3d,'Z',-999.,'no0' ) ! this if could be removed with the new lbc_lnk ... … … 556 601 CALL fliogetv( knumfl, cdvar, pv_r3d(:,:,:), start=istart(1:inbdim), count=icnt(1:inbdim) ) 557 602 !--- Apply scale_factor and offset 558 IF (zscf /= 1.) pv_r3d(:,:,:) = pv_r3d(:,:,:) * zscf 559 IF (zofs /= 0.) pv_r3d(:,:,:) = pv_r3d(:,:,:) + zofs 603 !CDIR COLLAPSE 604 IF( zscf /= 1.) pv_r3d(:,:,:) = pv_r3d(:,:,:) * zscf 605 !CDIR COLLAPSE 606 IF( zofs /= 0.) pv_r3d(:,:,:) = pv_r3d(:,:,:) + zofs 560 607 ENDIF 561 608 ENDIF 562 609 ! 563 IF( istop == nstop .AND. lwp ) &564 & WRITE(numout,*) 'read '//trim(cdvar)//' in '//trim(iom_file(knumfl)%name)//' ok'610 IF( istop == nstop .AND. lwp ) & 611 WRITE(numout,*) ' read '//trim(cdvar)//' in '//trim(iom_file(knumfl)%name)//' ok' 565 612 ENDIF 566 613 ! 567 END SUBROUTINE iom_ u_getv614 END SUBROUTINE iom_get_123d 568 615 569 616 570 617 SUBROUTINE iom_gettime( knumfl, cdvar, ptime ) 571 !!-------------------------------------------------------------------- 572 !! *** SUBROUTINE iom_close *** 573 !! 574 !! ** Purpose : read the time axis cdvar in the file 575 !! 576 !! ** Method : 577 !! 578 !!-------------------------------------------------------------------- 579 INTEGER , INTENT(in) :: knumfl ! Identifier of the file to be closed 580 CHARACTER(len=*) , INTENT(in) :: cdvar ! time axis name 581 REAL(wp), DIMENSION(:), INTENT(out) :: ptime ! the time axis 582 583 INTEGER :: idvar ! id of the variable 584 CHARACTER(LEN=100) :: clinfo ! info character 585 !--------------------------------------------------------------------- 586 clinfo = 'iom_gettime, file: '//trim(iom_file(knumfl)%name)//', var: '//trim(cdvar) 587 idvar = iom_varid( knumfl, cdvar ) 588 ! 589 ptime(:) = 0. ! default definition 590 IF ( idvar > 0 ) THEN 591 IF ( iom_file(knumfl)%ndims(idvar) == 1 ) THEN 592 IF ( iom_file(knumfl)%luld(idvar) ) THEN 593 IF ( iom_file(knumfl)%dimsz(1,idvar) == size(ptime) ) THEN 594 CALL fliogetv( knumfl, cdvar, ptime(:), start=(/ 1 /), & 595 & count=(/ iom_file(knumfl)%dimsz(1,idvar) /) ) 596 ELSE 597 WRITE(ctmp1,*) 'error with the size of ptime ',size(ptime),iom_file(knumfl)%dimsz(1,idvar) 598 CALL ctl_stop( trim(clinfo), trim(ctmp1) ) 599 ENDIF 600 ELSE 601 CALL ctl_stop( trim(clinfo), 'variable dimension is not unlimited... use iom_get' ) 602 ENDIF 603 ELSE 604 CALL ctl_stop( trim(clinfo), 'the variable has more than 1 dimension' ) 605 ENDIF 606 ELSE 607 CALL ctl_stop( trim(clinfo), 'variable not found in '//iom_file(knumfl)%name ) 608 ENDIF 609 618 !!-------------------------------------------------------------------- 619 !! *** SUBROUTINE iom_gettime *** 620 !! 621 !! ** Purpose : read the time axis cdvar in the file 622 !!-------------------------------------------------------------------- 623 INTEGER , INTENT(in ) :: knumfl ! file Identifier 624 CHARACTER(len=*) , INTENT(in ) :: cdvar ! time axis name 625 REAL(wp), DIMENSION(:), INTENT( out) :: ptime ! the time axis 626 ! 627 INTEGER :: idvar ! id of the variable 628 CHARACTER(LEN=100) :: clinfo ! info character 629 !--------------------------------------------------------------------- 630 ! 631 clinfo = 'iom_gettime, file: '//trim(iom_file(knumfl)%name)//', var: '//trim(cdvar) 632 idvar = iom_varid( knumfl, cdvar ) 633 ! 634 ptime(:) = 0. ! default definition 635 IF( idvar > 0 ) THEN 636 IF( iom_file(knumfl)%ndims(idvar) == 1 ) THEN 637 IF( iom_file(knumfl)%luld(idvar) ) THEN 638 IF( iom_file(knumfl)%dimsz(1,idvar) == size(ptime) ) THEN 639 CALL fliogetv( knumfl, cdvar, ptime(:), start=(/ 1 /), & 640 & count=(/ iom_file(knumfl)%dimsz(1,idvar) /) ) 641 ELSE 642 WRITE(ctmp1,*) 'error with the size of ptime ',size(ptime),iom_file(knumfl)%dimsz(1,idvar) 643 CALL ctl_stop( trim(clinfo), trim(ctmp1) ) 644 ENDIF 645 ELSE 646 CALL ctl_stop( trim(clinfo), 'variable dimension is not unlimited... use iom_get' ) 647 ENDIF 648 ELSE 649 CALL ctl_stop( trim(clinfo), 'the variable has more than 1 dimension' ) 650 ENDIF 651 ELSE 652 CALL ctl_stop( trim(clinfo), 'variable not found in '//iom_file(knumfl)%name ) 653 ENDIF 654 ! 610 655 END SUBROUTINE iom_gettime 611 656 612 657 613 614 658 FUNCTION iom_varid ( knumfl, cdvar, kdimsz ) 615 !!----------------------------------------------------------------------- 616 !! *** FUNCTION iom_varid *** 617 !! 618 !! ** Purpose : get the id of a variable in a file (return 0 if not found) 619 !! 620 !! ** Method : ??? 621 !! 622 !!----------------------------------------------------------------------- 623 INTEGER , INTENT(in) :: knumfl ! file Identifier 624 CHARACTER(len=*) , INTENT(in) :: cdvar ! name of the variable 625 INTEGER, DIMENSION(:), INTENT(out), OPTIONAL :: kdimsz ! size of the dimensions 626 ! 627 INTEGER :: iom_varid, idvar, i_nvd, ji 628 INTEGER, DIMENSION(jpmax_dims) :: idimid 629 LOGICAL :: ll_fnd 630 CHARACTER(LEN=100) :: clinfo ! info character 631 !!----------------------------------------------------------------------- 632 clinfo = 'iom_varid, file: '//trim(iom_file(knumfl)%name)//', var: '//trim(cdvar) 633 iom_varid = 0 ! default definition 634 IF ( PRESENT(kdimsz) ) kdimsz(:) = 0 ! default definition 635 ! 636 IF ( knumfl > 0 ) THEN 637 IF( iom_file(knumfl)%iopen == 0 ) THEN 638 CALL ctl_stop( trim(clinfo), 'the file is not open' ) 639 ELSE 640 ! 641 ll_fnd = .FALSE. 642 idvar = 0 643 ! 644 DO WHILE ( .NOT.ll_fnd .AND. idvar < iom_file(knumfl)%nvars ) 645 idvar = idvar + 1 646 ll_fnd = ( TRIM(cdvar) == TRIM(iom_file(knumfl)%cn_var(idvar)) ) 647 END DO 648 ! 649 IF( .NOT.ll_fnd ) THEN 650 idvar = idvar + 1 651 IF( idvar <= jpmax_vars ) THEN 652 CALL flioinqv( knumfl, cdvar, ll_fnd, nb_dims = i_nvd ) 653 IF( ll_fnd ) THEN 654 IF( i_nvd <= jpmax_dims ) THEN 655 iom_file(knumfl)%nvars = idvar 656 iom_file(knumfl)%cn_var(idvar) = trim(cdvar) 657 iom_file(knumfl)%ndims(idvar) = i_nvd 658 CALL flioinqv( knumfl, cdvar, ll_fnd, & 659 & len_dims = iom_file(knumfl)%dimsz(1:i_nvd,idvar), & 660 & id_dims = idimid(1:i_nvd) ) 661 DO ji = 1, i_nvd 662 IF ( idimid(ji) == iom_file(knumfl)%iduld ) iom_file(knumfl)%luld(idvar) = .TRUE. 663 END DO 664 !---------- 665 !---------- Deal with scale_factor and offset 666 CALL flioinqa( knumfl, cdvar, 'scale_factor', ll_fnd ) 667 IF (ll_fnd) THEN 668 CALL fliogeta( knumfl, cdvar, 'scale_factor', iom_file(knumfl)%scf(idvar) ) 669 ELSE 670 iom_file(knumfl)%scf(idvar) = 1. 671 END IF 672 CALL flioinqa( knumfl, cdvar, 'offset', ll_fnd ) 673 IF( ll_fnd ) THEN 674 CALL fliogeta( knumfl, cdvar, 'offset', iom_file(knumfl)%ofs(idvar) ) 675 ELSE 676 iom_file(knumfl)%ofs(idvar) = 0. 677 END IF 678 ! 679 iom_varid = idvar 680 IF ( PRESENT(kdimsz) ) THEN 681 IF ( i_nvd == size(kdimsz) ) THEN 682 kdimsz(:) = iom_file(knumfl)%dimsz(1:i_nvd,idvar) 683 ELSE 684 WRITE(ctmp1,*) i_nvd, size(kdimsz) 685 CALL ctl_stop( trim(clinfo), 'error in kdimsz size'//trim(ctmp1) ) 686 ENDIF 687 ENDIF 688 ELSE 689 CALL ctl_stop( trim(clinfo), 'Too many dimensions in the file '//iom_file(knumfl)%name, & 690 & 'increase the parameter jpmax_vars') 691 ENDIF 692 ELSE 693 CALL ctl_warn( trim(clinfo), 'Variable '//trim(cdvar)// & 694 & ' is not found in the file '//trim(iom_file(knumfl)%name) ) 695 ENDIF 696 ELSE 697 CALL ctl_stop( trim(clinfo), 'Too many variables in the file '//iom_file(knumfl)%name, & 698 & 'increase the parameter jpmax_vars') 699 ENDIF 700 ELSE 701 iom_varid = idvar 702 IF ( PRESENT(kdimsz) ) THEN 703 i_nvd = iom_file(knumfl)%ndims(idvar) 704 IF ( i_nvd == size(kdimsz) ) THEN 705 kdimsz(:) = iom_file(knumfl)%dimsz(1:i_nvd,idvar) 706 ELSE 707 WRITE(ctmp1,*) i_nvd, size(kdimsz) 708 CALL ctl_stop( trim(clinfo), 'error in kdimsz size'//trim(ctmp1) ) 709 ENDIF 710 ENDIF 711 ENDIF 712 ENDIF 713 ENDIF 714 659 !!----------------------------------------------------------------------- 660 !! *** FUNCTION iom_varid *** 661 !! 662 !! ** Purpose : get the id of a variable in a file (return 0 if not found) 663 !!----------------------------------------------------------------------- 664 INTEGER , INTENT(in ) :: knumfl ! file Identifier 665 CHARACTER(len=*) , INTENT(in ) :: cdvar ! name of the variable 666 INTEGER, DIMENSION(:), INTENT( out), OPTIONAL :: kdimsz ! size of the dimensions 667 ! 668 INTEGER :: ji ! dummy loop index 669 INTEGER :: iom_varid, idvar, i_nvd 670 INTEGER, DIMENSION(jpmax_dims) :: idimid 671 LOGICAL :: ll_fnd 672 CHARACTER(LEN=100) :: clinfo ! info character 673 !!----------------------------------------------------------------------- 674 clinfo = 'iom_varid, file: '//trim(iom_file(knumfl)%name)//', var: '//trim(cdvar) 675 iom_varid = 0 ! default definition 676 IF( PRESENT(kdimsz) ) kdimsz(:) = 0 ! default definition 677 ! 678 IF( knumfl > 0 ) THEN 679 IF( iom_file(knumfl)%iopen == 0 ) THEN 680 CALL ctl_stop( trim(clinfo), 'the file is not open' ) 681 ELSE 682 ll_fnd = .FALSE. 683 idvar = 0 684 ! 685 DO WHILE ( .NOT.ll_fnd .AND. idvar < iom_file(knumfl)%nvars ) 686 idvar = idvar + 1 687 ll_fnd = ( TRIM(cdvar) == TRIM(iom_file(knumfl)%cn_var(idvar)) ) 688 END DO 689 ! 690 IF( .NOT.ll_fnd ) THEN 691 idvar = idvar + 1 692 IF( idvar <= jpmax_vars ) THEN 693 CALL flioinqv( knumfl, cdvar, ll_fnd, nb_dims = i_nvd ) 694 IF( ll_fnd ) THEN 695 IF( i_nvd <= jpmax_dims ) THEN 696 iom_file(knumfl)%nvars = idvar 697 iom_file(knumfl)%cn_var(idvar) = trim(cdvar) 698 iom_file(knumfl)%ndims(idvar) = i_nvd 699 CALL flioinqv( knumfl, cdvar, ll_fnd, & 700 & len_dims = iom_file(knumfl)%dimsz(1:i_nvd,idvar), & 701 & id_dims = idimid(1:i_nvd) ) 702 DO ji = 1, i_nvd 703 IF( idimid(ji) == iom_file(knumfl)%iduld ) iom_file(knumfl)%luld(idvar) = .TRUE. 704 END DO 705 !---------- 706 !---------- Deal with scale_factor and offset 707 CALL flioinqa( knumfl, cdvar, 'scale_factor', ll_fnd ) 708 IF( ll_fnd) THEN 709 CALL fliogeta( knumfl, cdvar, 'scale_factor', iom_file(knumfl)%scf(idvar) ) 710 ELSE 711 iom_file(knumfl)%scf(idvar) = 1. 712 END IF 713 CALL flioinqa( knumfl, cdvar, 'offset', ll_fnd ) 714 IF( ll_fnd ) THEN 715 CALL fliogeta( knumfl, cdvar, 'offset', iom_file(knumfl)%ofs(idvar) ) 716 ELSE 717 iom_file(knumfl)%ofs(idvar) = 0. 718 END IF 719 ! 720 iom_varid = idvar 721 IF( PRESENT(kdimsz) ) THEN 722 IF( i_nvd == size(kdimsz) ) THEN 723 kdimsz(:) = iom_file(knumfl)%dimsz(1:i_nvd,idvar) 724 ELSE 725 WRITE(ctmp1,*) i_nvd, size(kdimsz) 726 CALL ctl_stop( trim(clinfo), 'error in kdimsz size'//trim(ctmp1) ) 727 ENDIF 728 ENDIF 729 ELSE 730 CALL ctl_stop( trim(clinfo), 'Too many dimensions in the file '//iom_file(knumfl)%name, & 731 & 'increase the parameter jpmax_vars') 732 ENDIF 733 !!$ ELSE 734 !!$ CALL ctl_warn( trim(clinfo), 'Variable '//trim(cdvar)// & 735 !!$ & ' is not found in the file '//trim(iom_file(knumfl)%name) ) 736 ENDIF 737 ELSE 738 CALL ctl_stop( trim(clinfo), 'Too many variables in the file '//iom_file(knumfl)%name, & 739 & 'increase the parameter jpmax_vars') 740 ENDIF 741 ELSE 742 iom_varid = idvar 743 IF( PRESENT(kdimsz) ) THEN 744 i_nvd = iom_file(knumfl)%ndims(idvar) 745 IF( i_nvd == size(kdimsz) ) THEN 746 kdimsz(:) = iom_file(knumfl)%dimsz(1:i_nvd,idvar) 747 ELSE 748 WRITE(ctmp1,*) i_nvd, size(kdimsz) 749 CALL ctl_stop( trim(clinfo), 'error in kdimsz size'//trim(ctmp1) ) 750 ENDIF 751 ENDIF 752 ENDIF 753 ENDIF 754 ENDIF 755 ! 715 756 END FUNCTION iom_varid 716 757 717 718 FUNCTION iom_get_gblatt( knumfl, kinfonum ) 719 !!----------------------------------------------------------------------- 720 !! *** FUNCTION iom_get_gblatt *** 758 !!---------------------------------------------------------------------- 759 !! INTERFACE iom_rstput 760 !!---------------------------------------------------------------------- 761 SUBROUTINE iom_rstput_0d( kt, kwrite, knumfl, cdvar, pvar ) 762 INTEGER , INTENT(in) :: kt ! ocean time-step 763 INTEGER , INTENT(in) :: kwrite ! writing time-step 764 INTEGER , INTENT(in) :: knumfl ! Identifier of the file 765 CHARACTER(len=*), INTENT(in) :: cdvar ! time axis name 766 REAL(wp) , INTENT(in) :: pvar ! read field 767 IF( knumfl > 0 ) CALL iom_rstput_0123d( kt, kwrite, knumfl, cdvar, pv_r0d = pvar ) 768 END SUBROUTINE iom_rstput_0d 769 770 SUBROUTINE iom_rstput_1d( kt, kwrite, knumfl, cdvar, pvar ) 771 INTEGER , INTENT(in) :: kt ! ocean time-step 772 INTEGER , INTENT(in) :: kwrite ! writing time-step 773 INTEGER , INTENT(in) :: knumfl ! Identifier of the file 774 CHARACTER(len=*), INTENT(in) :: cdvar ! time axis name 775 REAL(wp) , INTENT(in), DIMENSION( jpk) :: pvar ! read field 776 IF( knumfl > 0 ) CALL iom_rstput_0123d( kt, kwrite, knumfl, cdvar, pv_r1d = pvar ) 777 END SUBROUTINE iom_rstput_1d 778 779 SUBROUTINE iom_rstput_2d( kt, kwrite, knumfl, cdvar, pvar ) 780 INTEGER , INTENT(in) :: kt ! ocean time-step 781 INTEGER , INTENT(in) :: kwrite ! writing time-step 782 INTEGER , INTENT(in) :: knumfl ! Identifier of the file 783 CHARACTER(len=*), INTENT(in) :: cdvar ! time axis name 784 REAL(wp) , INTENT(in), DIMENSION(jpi,jpj ) :: pvar ! read field 785 IF( knumfl > 0 ) CALL iom_rstput_0123d( kt, kwrite, knumfl, cdvar, pv_r2d = pvar ) 786 END SUBROUTINE iom_rstput_2d 787 788 SUBROUTINE iom_rstput_3d( kt, kwrite, knumfl, cdvar, pvar ) 789 INTEGER , INTENT(in) :: kt ! ocean time-step 790 INTEGER , INTENT(in) :: kwrite ! writing time-step 791 INTEGER , INTENT(in) :: knumfl ! Identifier of the file 792 CHARACTER(len=*), INTENT(in) :: cdvar ! time axis name 793 REAL(wp) , INTENT(in), DIMENSION(jpi,jpj,jpk) :: pvar ! read field 794 IF( knumfl > 0 ) CALL iom_rstput_0123d( kt, kwrite, knumfl, cdvar, pv_r3d = pvar ) 795 END SUBROUTINE iom_rstput_3d 796 !!---------------------------------------------------------------------- 797 798 SUBROUTINE iom_rstput_0123d( kt, kwrite, knumfl, cdvar , & 799 & pv_r0d, pv_r1d, pv_r2d, pv_r3d ) 800 !!-------------------------------------------------------------------- 801 !! *** SUBROUTINE iom_rstput *** 721 802 !! 722 !! ** Purpose : ??? 723 !! 724 !! ** Method : ??? 725 !! 726 !!----------------------------------------------------------------------- 727 INTEGER,INTENT(in) :: knumfl 728 INTEGER, intent(in) :: kinfonum 729 ! 730 CHARACTER(LEN=10) :: clinfo 731 REAL(wp) :: iom_get_gblatt 732 !!----------------------------------------------------------------------- 733 734 WRITE(clinfo,*) kinfonum 735 clinfo = 'info'//trim(adjustl(clinfo)) 736 CALL fliogeta (knumfl, "?", clinfo, iom_get_gblatt) 737 738 END FUNCTION iom_get_gblatt 739 803 !! ** Purpose : read the time axis cdvar in the file 804 !!-------------------------------------------------------------------- 805 INTEGER , INTENT(in) :: kt ! ocean time-step 806 INTEGER , INTENT(in) :: kwrite ! writing time-step 807 INTEGER , INTENT(in) :: knumfl ! Identifier of the file 808 CHARACTER(len=*) , INTENT(in) :: cdvar ! time axis name 809 REAL(wp) , INTENT(in), OPTIONAL :: pv_r0d ! read field 810 REAL(wp), DIMENSION(:) , INTENT(in), OPTIONAL :: pv_r1d ! read field 811 REAL(wp), DIMENSION(:,:) , INTENT(in), OPTIONAL :: pv_r2d ! read field 812 REAL(wp), DIMENSION(:,:,:), INTENT(in), OPTIONAL :: pv_r3d ! read field 813 ! 814 INTEGER :: idims, idvar 815 INTEGER :: ix1, ix2, iy1, iy2 816 INTEGER, DIMENSION(4) :: idimsz, idimid 817 CHARACTER(LEN=100) :: clinfo ! info character 818 !--------------------------------------------------------------------- 819 ! 820 clinfo = ' iom_rstput_0123d, file: '//TRIM(iom_file(knumfl)%name)//', var: '//TRIM(cdvar) 821 822 ! define dimension variables if it is not already done 823 ! ========================== 824 IF( iom_file(knumfl)%nvars == 0 ) THEN 825 ! define the dimension variables if it is not already done 826 CALL fliodefv( knumfl,'nav_lon', (/1,2/), v_t=flio_r4 , axis='X', & 827 & long_name="Longitude", units="degrees_east" ) 828 CALL fliodefv( knumfl,'nav_lat', (/1,2/), v_t=flio_r4 , axis='Y', & 829 & long_name="Latitude", units="degrees_north" ) 830 CALL fliodefv( knumfl,'nav_lev', (/3/) , v_t=flio_i4 , axis='Z', & 831 & long_name="Model levels",units="model_levels") 832 CALL fliodefv( knumfl,'time_counter', (/4/), v_t=flio_r4, axis='T', & 833 & long_name="Time axis", units='seconds since 0001-01-01 00:00:00' ) 834 ! update informations structure related the dimension variable we just added... 835 iom_file(knumfl)%nvars = 4 836 iom_file(knumfl)%luld(1:4) = (/ .FALSE., .FALSE., .FALSE., .TRUE. /) 837 iom_file(knumfl)%cn_var(1:3) = (/ 'nav_lon', 'nav_lat', 'nav_lev' /) 838 iom_file(knumfl)%cn_var(4) = 'time_counter' 839 iom_file(knumfl)%ndims(1:4) = (/ 2, 2, 1, 1 /) 840 IF(lwp) WRITE(numout,*) TRIM(clinfo)//' define dimension variables done' 841 ENDIF 842 843 ! define the data if it is not already done 844 ! =============== 845 idvar = iom_varid( knumfl, cdvar ) 846 IF( idvar <= 0 ) THEN 847 ! variable definition 848 IF( PRESENT(pv_r0d) ) THEN ; idims = 0 849 ELSEIF( PRESENT(pv_r1d) ) THEN ; idims = 2 ; idimid(1:idims) = (/ 3,4/) 850 ELSEIF( PRESENT(pv_r2d) ) THEN ; idims = 3 ; idimid(1:idims) = (/1,2 ,4/) 851 ELSEIF( PRESENT(pv_r3d) ) THEN ; idims = 4 ; idimid(1:idims) = (/1,2,3,4/) 852 ENDIF 853 IF( PRESENT(pv_r0d) ) THEN ; CALL fliodefv (knumfl, cdvar , v_t = flio_r8) 854 ELSE ; CALL fliodefv (knumfl, cdvar, idimid(1:idims), v_t = flio_r8) 855 ENDIF 856 ! update informations structure related the new variable we want to add... 857 idvar = iom_file(knumfl)%nvars + 1 858 iom_file(knumfl)%nvars = idvar 859 iom_file(knumfl)%cn_var(idvar) = TRIM(cdvar) 860 iom_file(knumfl)%scf(idvar) = 1. 861 iom_file(knumfl)%ofs(idvar) = 0. 862 iom_file(knumfl)%ndims(idvar) = idims 863 IF( .NOT. PRESENT(pv_r0d) ) THEN 864 iom_file(knumfl)%luld(idvar) = .TRUE. 865 CALL flioinqf( knumfl, ln_dim = idimsz ) 866 iom_file(knumfl)%dimsz(1:idims-1,idvar) = idimsz(idimid(1:idims-1)) 867 ENDIF 868 IF(lwp) WRITE(numout,*) TRIM(clinfo)//' defined ok' 869 ENDIF 870 871 ! time step kwrite : write the variable 872 IF( kt == kwrite ) THEN 873 ! on what kind of domain must the data be written? 874 IF( PRESENT(pv_r2d) .OR. PRESENT(pv_r3d) ) THEN 875 idimsz(1:2) = iom_file(knumfl)%dimsz(1:2,idvar) 876 IF( idimsz(1) == (nlei - nldi + 1) .AND. idimsz(2) == (nlej - nldj + 1) ) THEN 877 ix1 = nldi ; ix2 = nlei ; iy1 = nldj ; iy2 = nlej 878 ELSEIF( idimsz(1) == nlci .AND. idimsz(2) == nlcj ) THEN 879 ix1 = 1 ; ix2 = nlci ; iy1 = 1 ; iy2 = nlcj 880 ELSEIF( idimsz(1) == jpi .AND. idimsz(2) == jpj ) THEN 881 ix1 = 1 ; ix2 = jpi ; iy1 = 1 ; iy2 = jpj 882 ELSE 883 CALL ctl_stop( 'iom_rstput_0123d: should have been an impossible case...' ) 884 ENDIF 885 886 ! write dimension variables if it is not already done 887 ! ============= 888 IF( iom_file(knumfl)%dimsz(1, 1) == 0 ) THEN 889 CALL flioputv( knumfl, 'nav_lon' , glamt(ix1:ix2, iy1:iy2) ) 890 CALL flioputv( knumfl, 'nav_lat' , gphit(ix1:ix2, iy1:iy2) ) 891 CALL flioputv( knumfl, 'nav_lev' , gdept_0 ) 892 CALL flioputv( knumfl, 'time_counter', kt ) ! +++ WRONG VALUE: to be improved but not really useful... 893 ! update informations structure related the dimension variable 894 iom_file(knumfl)%dimsz(1:2, 1) = idimsz(1:2) 895 iom_file(knumfl)%dimsz(1:2, 2) = idimsz(1:2) 896 iom_file(knumfl)%dimsz(1, 3:4) = (/idimsz(3), 1/) 897 IF(lwp) WRITE(numout,*) TRIM(clinfo)//' write dimension variables done' 898 ENDIF 899 ENDIF 900 901 ! write the data 902 ! ============= 903 IF( PRESENT(pv_r0d) ) THEN ; CALL flioputv( knumfl, cdvar, pv_r0d ) 904 ELSEIF( PRESENT(pv_r1d) ) THEN ; CALL flioputv( knumfl, cdvar, pv_r1d( :) ) 905 ELSEIF( PRESENT(pv_r2d) ) THEN ; CALL flioputv( knumfl, cdvar, pv_r2d(ix1:ix2, iy1:iy2 ) ) 906 ELSEIF( PRESENT(pv_r3d) ) THEN ; CALL flioputv( knumfl, cdvar, pv_r3d(ix1:ix2, iy1:iy2, :) ) 907 ENDIF 908 ! add 1 to the size of the temporal dimension (not really useful...) 909 IF( iom_file(knumfl)%luld(idvar) ) iom_file(knumfl)%dimsz(iom_file(knumfl)%ndims(idvar), idvar) & 910 & = iom_file(knumfl)%dimsz(iom_file(knumfl)%ndims(idvar), idvar) + 1 911 IF(lwp) WRITE(numout,*) TRIM(clinfo)//' written ok' 912 ENDIF 913 ! 914 END SUBROUTINE iom_rstput_0123d 915 740 916 !!====================================================================== 741 917 END MODULE iom -
trunk/NEMO/OPA_SRC/istate.F90
r479 r508 4 4 !! Ocean state : initial state setting 5 5 !!===================================================================== 6 !! History : 4.0 ! 89-12 (P. Andrich) Original code 7 !! 5.0 ! 91-11 (G. Madec) rewritting 8 !! 6.0 ! 96-01 (G. Madec) terrain following coordinates 9 !! 8.0 ! 01-09 (M. Levy, M. Ben Jelloul) istate_eel 10 !! 8.0 ! 01-09 (M. Levy, M. Ben Jelloul) istate_uvg 11 !! 9.0 ! 03-08 (G. Madec) F90: Free form, modules 12 !! 9.0 ! 03-09 (G. Madec, C. Talandier) add EEL R5 13 !! 9.0 ! 04-05 (A. Koch-Larrouy) istate_gyre 14 !! 9.0 ! 06-07 (S. Masson) distributed restart using iom 15 !!---------------------------------------------------------------------- 6 16 7 17 !!---------------------------------------------------------------------- … … 13 23 !! istate_uvg : initial velocity in geostropic balance 14 24 !!---------------------------------------------------------------------- 15 !! * Modules used16 25 USE oce ! ocean dynamics and active tracers 17 26 USE dom_oce ! ocean space and time domain … … 19 28 USE ldftra_oce ! ocean active tracers: lateral physics 20 29 USE zdf_oce ! ocean vertical physics 21 USE in_out_manager ! I/O manager22 30 USE phycst ! physical constants 23 31 USE wzvmod ! verctical velocity (wzv routine) … … 26 34 USE restart ! ocean restart (rst_read routine) 27 35 USE solisl ! ??? 28 36 USE in_out_manager ! I/O manager 37 USE iom 38 29 39 IMPLICIT NONE 30 40 PRIVATE 31 41 32 !! * Routine accessibility 33 PUBLIC istate_init ! routine called by step.F90 42 PUBLIC istate_init ! routine called by step.F90 34 43 35 44 !! * Substitutions … … 37 46 # include "vectopt_loop_substitute.h90" 38 47 !!---------------------------------------------------------------------- 39 !! OPA 9.0 , LOCEAN-IPSL (200 5)48 !! OPA 9.0 , LOCEAN-IPSL (2006) 40 49 !! $Header$ 41 !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt50 !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 42 51 !!---------------------------------------------------------------------- 43 52 … … 48 57 !! *** ROUTINE istate_init *** 49 58 !! 50 !! ** Purpose : Initialization of the dynamics and tracers. 51 !! 52 !! ** Method : 53 !! 54 !! History : 55 !! 4.0 ! 91-03 () Original code 56 !! ! 91-11 (G. Madec) 57 !! 9.0 ! 03-09 (G. Madec) F90: Free form, modules, orthogonality 58 !!---------------------------------------------------------------------- 59 USE iom 60 !! * Local declarations 61 !CT INTEGER :: inum 62 !!---------------------------------------------------------------------- 63 64 65 ! Initialization to zero 66 ! ---------------------- 67 68 ! before fields ! now fields ! after fields ! 69 ; ub (:,:,:) = 0.e0 ; un (:,:,:) = 0.e0 ; ua (:,:,:) = 0.e0 70 ; vb (:,:,:) = 0.e0 ; vn (:,:,:) = 0.e0 ; va (:,:,:) = 0.e0 71 ; ; wn (:,:,:) = 0.e0 ; 72 ; rotb (:,:,:) = 0.e0 ; rotn (:,:,:) = 0.e0 ; 73 ; hdivb(:,:,:) = 0.e0 ; hdivn(:,:,:) = 0.e0 ; 74 75 ; tb (:,:,:) = 0.e0 ; tn (:,:,:) = 0.e0 ; ta (:,:,:) = 0.e0 76 ; sb (:,:,:) = 0.e0 ; sn (:,:,:) = 0.e0 ; sa (:,:,:) = 0.e0 59 !! ** Purpose : Initialization of the dynamics and tracer fields. 60 !!---------------------------------------------------------------------- 61 62 IF(lwp) WRITE(numout,*) 63 IF(lwp) WRITE(numout,*) 'istate_ini : Initialization of the dynamics and tracers' 64 IF(lwp) WRITE(numout,*) '~~~~~~~~~~' 77 65 78 66 rhd (:,:,:) = 0.e0 79 67 rhop (:,:,:) = 0.e0 80 68 rn2 (:,:,:) = 0.e0 81 82 #if defined key_dynspg_rl83 ! rigid-lid formulation84 bsfb(:,:) = 0.e0 ! before barotropic stream-function85 bsfn(:,:) = 0.e0 ! now barotropic stream-function86 bsfd(:,:) = 0.e0 ! barotropic stream-function trend87 #endif88 ! free surface formulation89 sshb(:,:) = 0.e0 ! before sea-surface height90 sshn(:,:) = 0.e0 ! now sea-surface height91 92 69 93 70 IF( ln_rstart ) THEN ! Restart from a file … … 100 77 neuler = 0 ! Set time-step indicator at nit000 (euler forward) 101 78 adatrj = 0._wp 79 ! ! Initialization of ocean to zero 80 ! before fields ! now fields 81 ; ub (:,:,:) = 0.e0 ; un (:,:,:) = 0.e0 82 ; vb (:,:,:) = 0.e0 ; vn (:,:,:) = 0.e0 83 ; rotb (:,:,:) = 0.e0 ; rotn (:,:,:) = 0.e0 84 ; hdivb(:,:,:) = 0.e0 ; hdivn(:,:,:) = 0.e0 85 ! 102 86 IF( cp_cfg == 'eel' ) THEN 103 87 CALL istate_eel ! EEL configuration : start from pre-defined … … 107 91 ! ! and salinity fields 108 92 ELSE 109 ! ! Other configurations: Initial temperature and salinity fields 110 111 !CT CALL iom_open ('ssh', inum) 112 !CT CALL iom_get( inum, jpdom_local, 'sshb', sshb ) ! free surface formulation (ssh) 113 !CT sshn(:,:) = sshb(:,:) 114 !CT CALL iom_close (inum) 115 93 ! ! Other configurations: Initial temperature and salinity fields 116 94 #if defined key_dtatem 117 95 CALL dta_tem( nit000 ) ! read 3D temperature data … … 120 98 #else 121 99 IF(lwp) WRITE(numout,*) ! analytical temperature profile 122 IF(lwp) WRITE(numout,*)' Temperature initialization using an analytic profile'100 IF(lwp) WRITE(numout,*)' Temperature initialization using an analytic profile' 123 101 CALL istate_tem 124 102 #endif … … 130 108 ! No salinity data 131 109 IF(lwp)WRITE(numout,*) ! analytical salinity profile 132 IF(lwp)WRITE(numout,*)' Salinity initialisation using a constant value'110 IF(lwp)WRITE(numout,*)' Salinity initialisation using a constant value' 133 111 CALL istate_sal 134 112 #endif … … 139 117 ! ! ----------------- 140 118 CALL wzv( nit000 ) ! from horizontal divergence 141 119 ! 142 120 END SUBROUTINE istate_init 143 121 … … 153 131 !! 154 132 !! References : Philander ??? 155 !! 156 !! History : 157 !! 4.0 ! 89-12 (P. Andrich) Original code 158 !! 6.0 ! 96-01 (G. Madec) terrain following coordinates 159 !! 9.0 ! 02-09 (G. Madec) F90: Free form 160 !!---------------------------------------------------------------------- 161 !! * Local declarations 133 !!---------------------------------------------------------------------- 162 134 INTEGER :: ji, jj, jk 163 135 !!---------------------------------------------------------------------- 164 136 ! 165 137 IF(lwp) WRITE(numout,*) 166 138 IF(lwp) WRITE(numout,*) 'istate_tem : initial temperature profile' … … 181 153 & 1 , jpi , 5 , 1 , jpk , & 182 154 & 1 , 1. , numout ) 183 155 ! 184 156 END SUBROUTINE istate_tem 185 157 … … 194 166 !! 195 167 !! ** Action : Initialize sn and sb 196 !! 197 !! History : 198 !! 4.0 ! 89-12 (P. Andrich) Original code 199 !! 8.5 ! 02-09 (G. Madec) F90: Free form 200 !!---------------------------------------------------------------------- 201 !! * Local declarations 168 !!---------------------------------------------------------------------- 202 169 REAL(wp) :: zsal = 35.50_wp 203 170 !!---------------------------------------------------------------------- … … 224 191 !! - set velocity field including horizontal divergence 225 192 !! and relative vorticity fields 226 !! 227 !! History : 228 !! 8.0 ! 01-09 (M. Levy, M. Ben Jelloul) read file for EEL 2 & 6 229 !! 9.0 ! 03-09 (G. Madec, C. Talandier) EEL 5 230 !! 9.0 ! 05-11 (V. Garnier) Surface pressure gradient organization 231 !!---------------------------------------------------------------------- 232 !! * Modules used 193 !!---------------------------------------------------------------------- 233 194 USE eosbn2 ! eq. of state, Brunt Vaisala frequency (eos routine) 234 195 USE divcur ! hor. divergence & rel. vorticity (div_cur routine) 235 196 USE iom 236 197 237 !! * Local declarations238 198 INTEGER :: inum ! temporary logical unit 239 199 INTEGER :: ji, jj, jk ! dummy loop indices 240 200 INTEGER :: ijloc 241 REAL(wp) :: & 242 zh1, zh2, zslope, zcst, & ! temporary scalars 243 zfcor 244 REAL(wp) :: & 245 zt1 = 12._wp, & ! surface temperature value (EEL R5) 246 zt2 = 2._wp, & ! bottom temperature value (EEL R5) 247 zsal = 35.5_wp, & ! constant salinity (EEL R2, R5 and R6) 248 zueel = 0.1_wp ! constant uniform zonal velocity (EEL R5) 201 REAL(wp) :: zh1, zh2, zslope, zcst, zfcor ! temporary scalars 202 REAL(wp) :: zt1 = 12._wp, & ! surface temperature value (EEL R5) 203 & zt2 = 2._wp, & ! bottom temperature value (EEL R5) 204 & zsal = 35.5_wp, & ! constant salinity (EEL R2, R5 and R6) 205 & zueel = 0.1_wp ! constant uniform zonal velocity (EEL R5) 249 206 # if ! defined key_dynspg_rl 250 REAL(wp), DIMENSION(jpiglo,jpjglo) :: & 251 zssh ! initial ssh over the global domain 207 REAL(wp), DIMENSION(jpiglo,jpjglo) :: zssh ! initial ssh over the global domain 252 208 # endif 253 209 !!---------------------------------------------------------------------- … … 389 345 !! ** Method : - set temprature field 390 346 !! - set salinity field 391 !! 392 !! ** History : 393 !! 9.0 ! 04-05 (A. Koch-Larrouy) Original code 394 !!---------------------------------------------------------------------- 395 !! * Modules used 396 USE iom 397 398 !! * Local variables 399 INTEGER :: inum ! temporary logical unit 400 INTEGER, PARAMETER :: & 401 ntsinit = 0 ! (0/1) (analytical/input data files) T&S initialization 402 347 !!---------------------------------------------------------------------- 403 348 INTEGER :: ji, jj, jk ! dummy loop indices 349 INTEGER :: inum ! temporary logical unit 350 INTEGER, PARAMETER :: ntsinit = 0 ! (0/1) (analytical/input data files) T&S initialization 404 351 !!---------------------------------------------------------------------- 405 352 … … 469 416 ENDIF 470 417 471 472 418 END SUBROUTINE istate_gyre 473 474 419 475 420 … … 484 429 !! surface to the bottom. 485 430 !! p=integral [ rau*g dz ] 486 !! 487 !! History : 488 !! 8.1 ! 01-09 (M. Levy, M. Ben Jelloul) Original code 489 !! 8.5 ! 02-09 (G. Madec) F90: Free form 490 !! 9.0 ! 05-11 (V. Garnier) Surface pressure gradient organization 491 !!---------------------------------------------------------------------- 492 !! * Modules used 431 !!---------------------------------------------------------------------- 493 432 USE eosbn2 ! eq. of state, Brunt Vaisala frequency (eos routine) 494 433 USE dynspg ! surface pressure gradient (dyn_spg routine) … … 496 435 USE lbclnk ! ocean lateral boundary condition (or mpp link) 497 436 498 !! * Local declarations499 437 INTEGER :: ji, jj, jk ! dummy loop indices 500 438 INTEGER :: indic ! ??? 501 REAL(wp) :: & 502 zmsv, zphv, zmsu, zphu, & ! temporary scalars 503 zalfg 504 REAL(wp), DIMENSION (jpi,jpj,jpk) :: & 505 zprn ! workspace 439 REAL(wp) :: zmsv, zphv, zmsu, zphu, zalfg ! temporary scalars 440 REAL(wp), DIMENSION (jpi,jpj,jpk) :: zprn ! workspace 506 441 !!---------------------------------------------------------------------- 507 442 … … 514 449 515 450 zalfg = 0.5 * grav * rau0 516 ! Surface value 517 zprn(:,:,1) = zalfg * fse3w(:,:,1) * ( 1 + rhd(:,:,1) ) 518 519 ! Vertical integration from the surface 520 DO jk = 2, jpkm1 451 452 zprn(:,:,1) = zalfg * fse3w(:,:,1) * ( 1 + rhd(:,:,1) ) ! Surface value 453 454 DO jk = 2, jpkm1 ! Vertical integration from the surface 521 455 zprn(:,:,jk) = zprn(:,:,jk-1) & 522 456 & + zalfg * fse3w(:,:,jk) * ( 2. + rhd(:,:,jk) + rhd(:,:,jk-1) ) … … 525 459 ! Compute geostrophic balance 526 460 ! --------------------------- 527 528 461 DO jk = 1, jpkm1 529 462 DO jj = 2, jpjm1 … … 571 504 ! after initializing u and v, we need to calculate the initial streamfunction bsf. 572 505 ! Otherwise, only the trend will be computed and the model will blow up (inconsistency). 573 574 506 ! to do that, we call dyn_spg with a special trick: 575 ! we fill ua and va with the velocities divided by dt, 576 ! and the streamfunction will be brought to the right 577 ! value assuming the velocities have been set up in 578 ! one time step. 579 ! we then set bsfd to zero (first guess for next step 580 ! is d(psi)/dt = 0.) 581 582 ! sets up s false trend to calculate the barotropic 583 ! streamfunction. 507 ! we fill ua and va with the velocities divided by dt, and the streamfunction will be brought to the 508 ! right value assuming the velocities have been set up in one time step. 509 ! we then set bsfd to zero (first guess for next step is d(psi)/dt = 0.) 510 ! sets up s false trend to calculate the barotropic streamfunction. 584 511 585 512 ua(:,:,:) = ub(:,:,:) / rdt … … 612 539 hdivb(:,:,:) = hdivn(:,:,:) ! set the before to the now value 613 540 rotb (:,:,:) = rotn (:,:,:) ! set the before to the now value 614 541 ! 615 542 END SUBROUTINE istate_uvg 616 543 -
trunk/NEMO/OPA_SRC/opa.F90
r503 r508 48 48 USE obc_par ! open boundary cond. parameters 49 49 USE obcini ! open boundary cond. initialization (obc_ini routine) 50 USE solver ! solver initialization (solver_init routine)51 50 USE istate ! initial state setting (istate_init routine) 52 51 USE eosbn2 ! equation of state (eos bn2 routine) … … 69 68 70 69 USE step ! OPA time-stepping (stp routine) 71 USE dynspg_oce ! Control choice of surface pressure gradient schemes72 70 USE prtctl ! Print control (prt_ctl_init routine) 73 71 USE ini1d ! re-initialization of u-v mask for the 1D configuration … … 245 243 CALL istate_init ! ocean initial state (Dynamics and tracers) 246 244 247 IF( lk_dynspg_flt .OR. lk_dynspg_rl ) THEN248 CALL solver_init( nit000 ) ! Elliptic solver249 ENDIF250 251 245 !!add 252 246 CALL eos( tb, sb, rhd, rhop ) ! before potential and in situ densities -
trunk/NEMO/OPA_SRC/par_oce.F90
r467 r508 26 26 27 27 INTEGER, PUBLIC, PARAMETER :: & !: 28 jpni = 1, & !: number of processors following i29 jpnj = 1, & !: number of processors following j30 jpnij = 1, & !: nb of local domain = nb of processors28 jpni = 2, & !: number of processors following i 29 jpnj = 2, & !: number of processors following j 30 jpnij = 4, & !: nb of local domain = nb of processors 31 31 ! ! ( <= jpni x jpnj ) 32 32 jpr2di = 0, & !: number of columns for extra outer halo -
trunk/NEMO/OPA_SRC/restart.F90
r473 r508 3 3 !! *** MODULE restart *** 4 4 !! Ocean restart : write the ocean restart file 5 !!===================================================================== 6 7 !!---------------------------------------------------------------------- 8 !! rst_write : write of the restart file 9 !! rst_read : read the restart file 10 !!---------------------------------------------------------------------- 11 !! * Modules used 5 !!====================================================================== 6 !! History : ! 99-11 (M. Imbard) Original code 7 !! 8.5 ! 02-08 (G. Madec) F90: Free form 8 !! 9.0 ! 05-11 (V. Garnier) Surface pressure gradient organization 9 !! 9.0 ! 06-07 (S. Masson) use IOM for restart 10 !!---------------------------------------------------------------------- 11 12 !!---------------------------------------------------------------------- 13 !! rst_opn : open the ocean restart file 14 !! rst_write : write the ocean restart file 15 !! rst_read : read the ocean restart file 16 !!---------------------------------------------------------------------- 12 17 USE dom_oce ! ocean space and time domain 13 18 USE oce ! ocean dynamics and tracers 14 19 USE phycst ! physical constants 15 USE in_out_manager ! I/O manager16 20 USE daymod ! calendar 17 USE sol_oce ! ocean elliptic solver18 USE zdf_oce ! ???19 USE zdftke ! turbulent kinetic energy scheme20 21 USE ice_oce ! ice variables 21 22 USE blk_oce ! bulk variables 22 USE flx_oce ! sea-ice/ocean forcings variables23 USE dynspg_oce ! free surface time splitting scheme variables24 23 USE cpl_oce, ONLY : lk_cpl ! 24 USE in_out_manager ! I/O manager 25 USE iom ! I/O module 25 26 26 27 IMPLICIT NONE 27 28 PRIVATE 28 29 29 !! * Routine accessibility 30 PUBLIC rst_write ! routine called by step.F90 31 PUBLIC rst_read ! routine called by inidtr.F90 32 33 !! * Module variables 34 CHARACTER (len=48) :: & 35 crestart = 'initial.nc' ! restart file name 36 !!---------------------------------------------------------------------- 37 !! OPA 9.0 , LOCEAN-IPSL (2005) 30 PUBLIC rst_opn ! routine called by step module 31 PUBLIC rst_write ! routine called by step module 32 PUBLIC rst_read ! routine called by opa module 33 34 LOGICAL, PUBLIC :: lrst_oce !: logical to control the oce restart write 35 INTEGER, PUBLIC :: nitrst !: time step at which restart file should be written 36 INTEGER, PUBLIC :: numror, numrow !: logical unit for cean restart (read and write) 37 38 !! * Substitutions 39 # include "vectopt_loop_substitute.h90" 40 !!---------------------------------------------------------------------- 41 !! OPA 9.0 , LOCEAN-IPSL (2006) 38 42 !! $Header$ 39 !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt 40 !!---------------------------------------------------------------------- 41 43 !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 44 !!---------------------------------------------------------------------- 42 45 43 46 CONTAINS 47 48 SUBROUTINE rst_opn( kt ) 49 !!--------------------------------------------------------------------- 50 !! *** ROUTINE rst_opn *** 51 !! 52 !! ** Purpose : + initialization (should be read in the namelist) of nitrst 53 !! + open the restart when we are one time step before nitrst 54 !! - restart header is defined when kt = nitrst-1 55 !! - restart data are written when kt = nitrst 56 !! + define lrst_oce to .TRUE. when we need to define or write the restart 57 !!---------------------------------------------------------------------- 58 INTEGER, INTENT(in) :: kt ! ocean time-step 59 !! 60 CHARACTER(LEN=20) :: clkt ! ocean time-step deine as a character 61 CHARACTER(LEN=50) :: clname ! ice output restart file name 62 !!---------------------------------------------------------------------- 63 ! 64 IF( kt == nit000 ) THEN ! default initialization, to do: should be read in the namelist... 65 nitrst = nitend ! to do: should be read in the namelist in a cleaver way... 66 lrst_oce = .FALSE. 67 ENDIF 68 69 IF ( kt == nitrst-1 .AND. lrst_oce ) THEN 70 CALL ctl_stop( 'rst_opn: we cannot create an ocean restart at every time step' ) 71 numrow = 0 72 ELSEIF( kt == nitrst-1 .OR. nitend == nit000 ) THEN ! beware if model runs only one time step 73 ! beware of the format used to write kt (default is i8.8, that should be large enough) 74 IF( nitrst > 1.0e9 ) THEN 75 WRITE(clkt,*) nitrst 76 ELSE 77 WRITE(clkt,'(i8.8)') nitrst 78 ENDIF 79 ! create the file 80 IF(lwp) WRITE(numout,*) 81 clname = TRIM(cexper)//"_"//TRIM(ADJUSTL(clkt))//"_restart" 82 IF(lwp) WRITE(numout,*) ' open ocean restart.output NetCDF file: '//clname 83 CALL iom_open( clname, numrow, ldwrt = .TRUE. ) 84 lrst_oce = .TRUE. 85 ENDIF 86 ! 87 END SUBROUTINE rst_opn 88 44 89 45 90 #if ( defined key_mpp_mpi || defined key_mpp_shmem ) && defined key_dimgout … … 66 111 !! ** Purpose : Write restart fields in NetCDF format 67 112 !! 68 !! ** Method : Write in num wrs file each nstock time stepin NetCDF113 !! ** Method : Write in numrow when kt == nitrst in NetCDF 69 114 !! file, save fields which are necessary for restart 70 !! 71 !! History : 72 !! ! 99-11 (M. Imbard) Original code 73 !! 8.5 ! 02-08 (G. Madec) F90: Free form 74 !! 9.0 ! 05-11 (V. Garnier) Surface pressure gradient organization 75 !!---------------------------------------------------------------------- 76 !! * Modules used 77 USE ioipsl 78 79 !! * Arguments 80 INTEGER, INTENT( in ) :: kt ! ocean time-step 81 82 !! * Local declarations 83 LOGICAL :: llbon 84 CHARACTER (len=50) :: clname, cln 85 INTEGER :: ic, jc, itime 86 INTEGER :: inumwrs 87 REAL(wp) :: zdate0 88 REAL(wp), DIMENSION( 1) :: zfice, zfblk ! used only in case of ice & bulk 89 REAL(wp), DIMENSION(10) :: zinfo(10) 90 REAL(wp), DIMENSION(jpi,jpj) :: ztab 91 #if defined key_agrif 92 Integer :: knum 93 #endif 94 !!---------------------------------------------------------------------- 95 96 IF( kt == nit000 ) THEN 97 IF(lwp) WRITE(numout,*) 98 IF(lwp) WRITE(numout,*) 'rst_wri : write restart.output NetCDF file' 99 IF(lwp) WRITE(numout,*) '~~~~~~~' 100 zfice(1) = 1.e0 ; zfblk(1) = 1.e0 101 ENDIF 102 103 104 IF( MOD( kt, nstock ) == 0 .OR. kt == nitend ) THEN 105 106 ! 0. Initializations 107 ! ------------------ 108 109 IF(lwp) WRITE(numout,*) ' ' 110 IF(lwp) WRITE(numout,*) 'rst_write : write the restart file in NetCDF format ', & 111 'at it= ',kt,' date= ',ndastp 112 IF(lwp) WRITE(numout,*) '~~~~~~~~~' 113 114 ! Job informations 115 zinfo(:) = 0.e0 116 zinfo(1) = FLOAT( no ) ! job number 117 zinfo(2) = FLOAT( kt ) ! time-step 118 zinfo(3) = FLOAT( 2 - nsolv ) ! pcg solver 119 zinfo(4) = FLOAT( nsolv - 1 ) ! sor solver 120 IF( lk_zdftke ) THEN 121 zinfo(5) = 1.e0 ! TKE 122 ELSE 123 zinfo(5) = 0.e0 ! no TKE 124 ENDIF 125 zinfo(6) = FLOAT( ndastp ) ! date 126 zinfo(7) = adatrj ! ??? 127 128 ! delete the restart file if it exists 129 INQUIRE( FILE=crestart, EXIST=llbon ) 130 IF(llbon) THEN 131 #if defined key_agrif 132 knum =Agrif_Get_Unit() 133 OPEN( UNIT=knum, FILE=crestart, STATUS='old' ) 134 CLOSE( knum, STATUS='delete' ) 135 #else 136 OPEN( UNIT=inumwrs, FILE=crestart, STATUS='old' ) 137 CLOSE( inumwrs, STATUS='delete' ) 138 #endif 139 ENDIF 140 141 ! Name of the new restart file 142 ic = 1 143 DO jc = 1, 16 144 IF( cexper(jc:jc) /= ' ' ) ic = jc 145 END DO 146 WRITE(cln,'("_",i4.4,i2.2,i2.2,"_restart")') nyear, nmonth, nday 147 clname = cexper(1:ic)//cln 148 ic = 1 149 DO jc = 1, 48 150 IF( clname(jc:jc) /= ' ' ) ic = jc 151 END DO 152 crestart = clname(1:ic)//".nc" 153 itime = 0 154 CALL ymds2ju( nyear, nmonth, nday, 0.e0, zdate0 ) 155 CALL restini( 'NONE', jpi, jpj, glamt, gphit, jpk, gdept_0, clname, & 156 itime, zdate0, rdt*nstock ,inumwrs, domain_id=nidom ) 157 158 CALL restput( inumwrs, 'info' , 1 , 1 , 10 , 0, zinfo ) ! restart informations 159 160 CALL restput( inumwrs, 'ub' , jpi, jpj, jpk, 0, ub ) ! prognostic variables 161 CALL restput( inumwrs, 'vb' , jpi, jpj, jpk, 0, vb ) 162 CALL restput( inumwrs, 'tb' , jpi, jpj, jpk, 0, tb ) 163 CALL restput( inumwrs, 'sb' , jpi, jpj, jpk, 0, sb ) 164 CALL restput( inumwrs, 'rotb' , jpi, jpj, jpk, 0, rotb ) 165 CALL restput( inumwrs, 'hdivb' , jpi, jpj, jpk, 0, hdivb ) 166 CALL restput( inumwrs, 'un' , jpi, jpj, jpk, 0, un ) 167 CALL restput( inumwrs, 'vn' , jpi, jpj, jpk, 0, vn ) 168 CALL restput( inumwrs, 'tn' , jpi, jpj, jpk, 0, tn ) 169 CALL restput( inumwrs, 'sn' , jpi, jpj, jpk, 0, sn ) 170 CALL restput( inumwrs, 'rotn' , jpi, jpj, jpk, 0, rotn ) 171 CALL restput( inumwrs, 'hdivn' , jpi, jpj, jpk, 0, hdivn ) 172 173 ztab(:,:) = gcx(1:jpi,1:jpj) 174 CALL restput( inumwrs, 'gcx' , jpi, jpj, 1 , 0, ztab ) ! Read elliptic solver arrays 175 ztab(:,:) = gcxb(1:jpi,1:jpj) 176 CALL restput( inumwrs, 'gcxb' , jpi, jpj, 1 , 0, ztab ) 177 # if defined key_dynspg_rl 178 CALL restput( inumwrs, 'bsfb' , jpi, jpj, 1 , 0, bsfb ) ! Rigid-lid formulation (bsf) 179 CALL restput( inumwrs, 'bsfn' , jpi, jpj, 1 , 0, bsfn ) 180 CALL restput( inumwrs, 'bsfd' , jpi, jpj, 1 , 0, bsfd ) 181 # else 182 CALL restput( inumwrs, 'sshb' , jpi, jpj, 1 , 0, sshb ) ! free surface formulation (ssh) 183 CALL restput( inumwrs, 'sshn' , jpi, jpj, 1 , 0, sshn ) 184 # if defined key_dynspg_ts 185 CALL restput( inumwrs, 'sshb_b' , jpi, jpj, 1 , 0, sshb_b ) ! free surface formulation (ssh) 186 CALL restput( inumwrs, 'sshn_b' , jpi, jpj, 1 , 0, sshn_b ) ! issued from barotropic loop 187 CALL restput( inumwrs, 'un_b' , jpi, jpj, 1 , 0, un_b ) ! horizontal transports 188 CALL restput( inumwrs, 'vn_b' , jpi, jpj, 1 , 0, vn_b ) ! issued from barotropic loop 115 !!---------------------------------------------------------------------- 116 INTEGER, INTENT(in) :: kt ! ocean time-step 117 !!---------------------------------------------------------------------- 118 119 IF(lwp) THEN 120 WRITE(numout,*) 121 WRITE(numout,*) 'rst_write : write ocean NetCDF restart file kt =', kt,' date= ', ndastp 122 WRITE(numout,*) '~~~~~~~~~' 123 ENDIF 124 125 ! calendar control 126 CALL iom_rstput( kt, nitrst, numrow, 'kt' , REAL( kt , wp) ) ! time-step 127 CALL iom_rstput( kt, nitrst, numrow, 'ndastp' , REAL( ndastp, wp) ) ! date 128 CALL iom_rstput( kt, nitrst, numrow, 'adatrj' , adatrj ) ! number of elapsed days since 129 ! ! the begining of the run [s] 130 131 ! prognostic variables 132 CALL iom_rstput( kt, nitrst, numrow, 'ub' , ub ) 133 CALL iom_rstput( kt, nitrst, numrow, 'vb' , vb ) 134 CALL iom_rstput( kt, nitrst, numrow, 'tb' , tb ) 135 CALL iom_rstput( kt, nitrst, numrow, 'sb' , sb ) 136 CALL iom_rstput( kt, nitrst, numrow, 'rotb' , rotb ) 137 CALL iom_rstput( kt, nitrst, numrow, 'hdivb' , hdivb ) 138 CALL iom_rstput( kt, nitrst, numrow, 'un' , un ) 139 CALL iom_rstput( kt, nitrst, numrow, 'vn' , vn ) 140 CALL iom_rstput( kt, nitrst, numrow, 'tn' , tn ) 141 CALL iom_rstput( kt, nitrst, numrow, 'sn' , sn ) 142 CALL iom_rstput( kt, nitrst, numrow, 'rotn' , rotn ) 143 CALL iom_rstput( kt, nitrst, numrow, 'hdivn' , hdivn ) 144 145 # if defined key_ice_lim 146 CALL iom_rstput( kt, nitrst, numrow, 'nfice' , REAL( nfice, wp) ) ! ice computation frequency 147 CALL iom_rstput( kt, nitrst, numrow, 'sst_io' , sst_io ) 148 CALL iom_rstput( kt, nitrst, numrow, 'sss_io' , sss_io ) 149 CALL iom_rstput( kt, nitrst, numrow, 'u_io' , u_io ) 150 CALL iom_rstput( kt, nitrst, numrow, 'v_io' , v_io ) 151 # if defined key_coupled 152 CALL iom_rstput( kt, nitrst, numrow, 'alb_ice', alb_ice ) 189 153 # endif 190 154 # endif 191 # if defined key_zdftke || defined key_esopa192 IF( lk_zdftke ) THEN193 CALL restput( inumwrs, 'en' , jpi, jpj, jpk, 0, en ) ! TKE arrays194 ENDIF195 # endif196 # if defined key_ice_lim197 zfice(1) = FLOAT( nfice ) ! Louvain La Neuve Sea Ice Model198 CALL restput( inumwrs, 'nfice' , 1, 1, 1 , 0, zfice )199 CALL restput( inumwrs, 'sst_io' , jpi, jpj, 1 , 0, sst_io )200 CALL restput( inumwrs, 'sss_io' , jpi, jpj, 1 , 0, sss_io )201 CALL restput( inumwrs, 'u_io' , jpi, jpj, 1 , 0, u_io )202 CALL restput( inumwrs, 'v_io' , jpi, jpj, 1 , 0, v_io )203 # if defined key_coupled204 CALL restput( inumwrs, 'alb_ice', jpi, jpj, 1 , 0, alb_ice )205 # endif206 # endif207 155 # if defined key_flx_bulk_monthly || defined key_flx_bulk_daily 208 zfblk(1) = FLOAT( nfbulk ) ! Bulk209 CALL restput( inumwrs, 'nfbulk' , 1, 1, 1 , 0, zfblk)210 CALL restput( inumwrs, 'gsst' , jpi, jpj, 1 , 0, gsst ) 211 # endif 212 213 CALL restclo( inumwrs ) ! close the restart file214 215 ENDIF 216 156 CALL iom_rstput( kt, nitrst, numrow, 'nfbulk' , REAL( nfbulk, wp) ) ! bulk computation frequency 157 CALL iom_rstput( kt, nitrst, numrow, 'gsst' , gsst ) 158 # endif 159 160 IF( kt == nitrst ) THEN 161 CALL iom_close( numrow ) ! close the restart file (only at last time step) 162 lrst_oce = .FALSE. 163 ENDIF 164 ! 217 165 END SUBROUTINE rst_write 218 166 … … 246 194 !! nrstdt = 2 the duration of the experiment in days (adatrj) 247 195 !! has been stored in the restart file. 248 !! 249 !! History : 250 !! ! 99-05 (M. Imbard) Original code 251 !! 8.5 ! 02-09 (G. Madec) F90: Free form 252 !! 9.0 ! 05-11 (V. Garnier) Surface pressure gradient organization 253 !!---------------------------------------------------------------------- 254 !! * Modules used 255 USE iom 256 257 !! * Local declarations 258 INTEGER :: & 259 inum ! temporary logical unit 260 REAL(wp), DIMENSION(1, 1, 10) :: zinfo 261 REAL(wp), DIMENSION(1, 1, 1) :: zzz 262 INTEGER :: ios 263 # if defined key_ice_lim 196 !!---------------------------------------------------------------------- 197 REAL(wp) :: zcoef, zkt, zndastp, znfice, znfbulk 198 # if defined key_ice_lim 264 199 INTEGER :: ji, jj 265 # endif 266 !!---------------------------------------------------------------------- 267 268 IF(lwp) WRITE(numout,*) 269 IF(lwp) WRITE(numout,*) 'rst_read : read the NetCDF restart file' 270 IF(lwp) WRITE(numout,*) '~~~~~~~~' 271 272 IF(lwp) WRITE(numout,*) ' Info on the present job : ' 273 IF(lwp) WRITE(numout,*) ' job number : ', no 274 IF(lwp) WRITE(numout,*) ' time-step : ', nit000 275 IF(lwp) WRITE(numout,*) ' solver type : ', nsolv 276 IF( lk_zdftke ) THEN 277 IF(lwp) WRITE(numout,*) ' tke option : 1 ' 278 ELSE 279 IF(lwp) WRITE(numout,*) ' tke option : 0 ' 280 ENDIF 281 IF(lwp) WRITE(numout,*) ' date ndastp : ', ndastp 282 IF(lwp) WRITE(numout,*) 283 284 ! Time domain : restart 285 ! ------------------------- 286 287 IF(lwp) WRITE(numout,*) 288 IF(lwp) WRITE(numout,*) 289 IF(lwp) WRITE(numout,*) ' *** restart option' 290 SELECT CASE ( nrstdt ) 291 CASE ( 0 ) 292 IF(lwp) WRITE(numout,*) ' nrstdt = 0 no control of nit000' 293 CASE ( 1 ) 294 IF(lwp) WRITE(numout,*) ' nrstdt = 1 we control the date of nit000' 295 CASE ( 2 ) 296 IF(lwp) WRITE(numout,*) ' nrstdt = 2 the date adatrj is read in restart file' 297 CASE DEFAULT 298 IF(lwp) WRITE(numout,*) ' ===>>>> nrstdt not equal 0, 1 or 2 : no control of the date' 299 IF(lwp) WRITE(numout,*) ' ======= =========' 300 END SELECT 301 302 CALL iom_open ( 'restart', inum ) 303 304 CALL iom_get ( inum, jpdom_unknown, 'info', zinfo ) 305 306 IF(lwp) WRITE(numout,*) 307 IF(lwp) WRITE(numout,*) ' Info on the restart file read : ' 308 IF(lwp) WRITE(numout,*) ' job number : ', NINT( zinfo(1, 1, 1) ) 309 IF(lwp) WRITE(numout,*) ' time-step : ', NINT( zinfo(1, 1, 2) ) 310 IF(lwp) WRITE(numout,*) ' solver type : ', NINT( zinfo(1, 1, 4) ) + 1 311 IF(lwp) WRITE(numout,*) ' tke option : ', NINT( zinfo(1, 1, 5) ) 312 IF(lwp) WRITE(numout,*) ' date ndastp : ', NINT( zinfo(1, 1, 6) ) 313 IF(lwp) WRITE(numout,*) 314 200 # endif 201 !!---------------------------------------------------------------------- 202 203 IF(lwp) THEN ! Contol prints 204 WRITE(numout,*) 205 WRITE(numout,*) 'rst_read : read oce NetCDF restart file' 206 WRITE(numout,*) '~~~~~~~~' 207 208 WRITE(numout,*) ' *** Info on the present job : ' 209 WRITE(numout,*) ' time-step : ', nit000 210 !!$ WRITE(numout,*) ' solver type : ', nsolv 211 !!$ IF( lk_zdftke ) THEN 212 !!$ WRITE(numout,*) ' tke option : 1 ' 213 !!$ ELSE 214 !!$ WRITE(numout,*) ' tke option : 0 ' 215 !!$ ENDIF 216 WRITE(numout,*) ' date ndastp : ', ndastp 217 WRITE(numout,*) 218 WRITE(numout,*) ' *** restart option' 219 SELECT CASE ( nrstdt ) 220 CASE ( 0 ) 221 WRITE(numout,*) ' nrstdt = 0 no control of nit000' 222 CASE ( 1 ) 223 WRITE(numout,*) ' nrstdt = 1 we control the date of nit000' 224 CASE ( 2 ) 225 WRITE(numout,*) ' nrstdt = 2 the date adatrj is read in restart file' 226 CASE DEFAULT 227 WRITE(numout,*) ' ===>>>> nrstdt not equal 0, 1 or 2 : no control of the date' 228 WRITE(numout,*) ' ======= =========' 229 END SELECT 230 WRITE(numout,*) 231 ENDIF 232 233 CALL iom_open( 'restart', numror ) ! Open 234 235 ! Calendar informations 236 CALL iom_get( numror, 'kt' , zkt ) ! time-step 237 CALL iom_get( numror, 'ndastp', zndastp ) ! date 238 ! Additional contol prints 239 IF(lwp) THEN 240 WRITE(numout,*) 241 WRITE(numout,*) ' *** Info on the restart file read : ' 242 WRITE(numout,*) ' time-step : ', NINT( zkt ) 243 !!$ WRITE(numout,*) ' solver type : ', +++ 244 !!$ WRITE(numout,*) ' tke option : ', +++ 245 WRITE(numout,*) ' date ndastp : ', NINT( zndastp ) 246 WRITE(numout,*) 247 ENDIF 315 248 ! Control of date 316 IF( nit000 - NINT( z info(1, 1, 2)) /= 1 .AND. nrstdt /= 0 ) &249 IF( nit000 - NINT( zkt ) /= 1 .AND. nrstdt /= 0 ) & 317 250 & CALL ctl_stop( ' ===>>>> : problem with nit000 for the restart', & 318 251 & ' verify the restart file or rerun with nrstdt = 0 (namelist)' ) 319 320 252 ! re-initialisation of adatrj0 321 adatrj0 = ( FLOAT( nit000-1 ) * rdttra(1) ) / rday 322 253 adatrj0 = ( REAL( nit000-1, wp ) * rdttra(1) ) / rday 323 254 IF ( nrstdt == 2 ) THEN 324 255 ! by default ndatsp has been set to ndate0 in dom_nam 325 256 ! ndate0 has been read in the namelist (standard OPA 8) 326 257 ! here when nrstdt=2 we keep the final date of previous run 327 ndastp = NINT( zinfo(1, 1, 6) ) 328 adatrj0 = zinfo(1, 1, 7) 329 ENDIF 330 331 CALL iom_get( inum, jpdom_local, 'ub' , ub ) ! Read prognostic variables 332 CALL iom_get( inum, jpdom_local, 'vb' , vb ) 333 CALL iom_get( inum, jpdom_local, 'tb' , tb ) 334 CALL iom_get( inum, jpdom_local, 'sb' , sb ) 335 CALL iom_get( inum, jpdom_local, 'rotb' , rotb ) 336 CALL iom_get( inum, jpdom_local, 'hdivb', hdivb ) 337 CALL iom_get( inum, jpdom_local, 'un' , un ) 338 CALL iom_get( inum, jpdom_local, 'vn' , vn ) 339 CALL iom_get( inum, jpdom_local, 'tn' , tn ) 340 CALL iom_get( inum, jpdom_local, 'sn' , sn ) 341 CALL iom_get( inum, jpdom_local, 'rotn' , rotn ) 342 CALL iom_get( inum, jpdom_local, 'hdivn', hdivn ) 343 ! Caution : extrahallow 344 ! gcx and gcxb are defined as: DIMENSION(1-jpr2di:jpi+jpr2di,1-jpr2dj:jpj+jpr2dj) 345 CALL iom_get( inum, jpdom_local, 'gcx' , gcx (1:jpi,1:jpj) ) 346 CALL iom_get( inum, jpdom_local, 'gcxb', gcxb(1:jpi,1:jpj) ) ! Read elliptic solver arrays 347 # if defined key_dynspg_rl 348 CALL iom_get( inum, jpdom_local, 'bsfb', bsfb ) ! Rigid-lid formulation (bsf) 349 CALL iom_get( inum, jpdom_local, 'bsfn', bsfn ) 350 CALL iom_get( inum, jpdom_local, 'bsfd', bsfd ) 351 # else 352 CALL iom_get( inum, jpdom_local, 'sshb', sshb ) ! free surface formulation (ssh) 353 CALL iom_get( inum, jpdom_local, 'sshn', sshn ) 354 # if defined key_dynspg_ts 355 CALL iom_get( inum, jpdom_local, 'sshb_b', sshb_b ) ! free surface formulation (ssh) 356 CALL iom_get( inum, jpdom_local, 'sshn_b', sshn_b ) ! issued from barotropic loop 357 CALL iom_get( inum, jpdom_local, 'un_b' , un_b ) ! horizontal transports 358 CALL iom_get( inum, jpdom_local, 'vn_b' , vn_b ) ! issued from barotropic loop 359 # endif 360 # endif 361 # if defined key_zdftke || defined key_esopa 362 IF( lk_zdftke ) THEN 363 IF( NINT( zinfo(1, 1, 5) ) == 1 ) THEN ! Read tke arrays 364 CALL iom_get( inum, jpdom_local, 'en', en ) 365 ln_rstke = .FALSE. 366 ELSE 367 IF(lwp) WRITE(numout,*) ' ===>>>> : the previous restart file did not used tke scheme' 368 IF(lwp) WRITE(numout,*) ' ======= =======' 369 nrstdt = 2 370 ln_rstke = .TRUE. 371 ENDIF 372 ENDIF 373 # endif 258 ndastp = NINT( zndastp ) 259 CALL iom_get( numror, 'adatrj', adatrj ) ! number of elapsed days since the begining of last run 260 ENDIF 261 262 ! ! Read prognostic variables 263 CALL iom_get( numror, jpdom_local, 'ub' , ub ) ! before i-component velocity 264 CALL iom_get( numror, jpdom_local, 'vb' , vb ) ! before j-component velocity 265 CALL iom_get( numror, jpdom_local, 'tb' , tb ) ! before temperature 266 CALL iom_get( numror, jpdom_local, 'sb' , sb ) ! before salinity 267 CALL iom_get( numror, jpdom_local, 'rotb' , rotb ) ! before curl 268 CALL iom_get( numror, jpdom_local, 'hdivb', hdivb ) ! before horizontal divergence 269 CALL iom_get( numror, jpdom_local, 'un' , un ) ! now i-component velocity 270 CALL iom_get( numror, jpdom_local, 'vn' , vn ) ! now j-component velocity 271 CALL iom_get( numror, jpdom_local, 'tn' , tn ) ! now temperature 272 CALL iom_get( numror, jpdom_local, 'sn' , sn ) ! now salinity 273 CALL iom_get( numror, jpdom_local, 'rotn' , rotn ) ! now curl 274 CALL iom_get( numror, jpdom_local, 'hdivn', hdivn ) ! now horizontal divergence 275 276 277 IF( neuler == 0 ) THEN ! Euler restart (neuler=0) 278 tb (:,:,:) = tn (:,:,:) ! all before fields set to now field values 279 sb (:,:,:) = sn (:,:,:) 280 ub (:,:,:) = un (:,:,:) 281 vb (:,:,:) = vn (:,:,:) 282 rotb (:,:,:) = rotn (:,:,:) 283 hdivb(:,:,:) = hdivn(:,:,:) 284 ENDIF 285 286 !!sm: TO BE MOVED IN NEW SURFACE MODULE... 287 374 288 # if defined key_ice_lim 375 289 ! Louvain La Neuve Sea Ice Model 376 ios = iom_varid( inum, 'nfice' ) 377 IF( ios > 0 ) then 378 CALL iom_get( inum, jpdom_unknown, 'nfice' , zzz ) 379 zinfo(1, 1, 8) = zzz(1, 1, 1) 380 CALL iom_get( inum, jpdom_local, 'sst_io', sst_io ) 381 CALL iom_get( inum, jpdom_local, 'sss_io', sss_io ) 382 CALL iom_get( inum, jpdom_local, 'u_io' , u_io ) 383 CALL iom_get( inum, jpdom_local, 'v_io' , v_io ) 290 IF( iom_varid( numror, 'nfice' ) > 0 ) then 291 CALL iom_get( numror , 'nfice' , znfice ) ! ice computation frequency 292 CALL iom_get( numror, jpdom_local, 'sst_io' , sst_io ) 293 CALL iom_get( numror, jpdom_local, 'sss_io' , sss_io ) 294 CALL iom_get( numror, jpdom_local, 'u_io' , u_io ) 295 CALL iom_get( numror, jpdom_local, 'v_io' , v_io ) 384 296 #if defined key_coupled 385 CALL iom_get( inum, jpdom_local, 'alb_ice', alb_ice )297 CALL iom_get( numror, jpdom_local, 'alb_ice', alb_ice ) 386 298 #endif 387 ENDIF 388 IF( zinfo(1, 1, 8) /= FLOAT(nfice) .OR. ios == 0 ) THEN 299 IF( znfice /= REAL( nfice, wp ) ) THEN ! if nfice changed between 2 runs 300 zcoef = REAL( nfice-1, wp ) / znfice 301 sst_io(:,:) = zcoef * sst_io(:,:) 302 sss_io(:,:) = zcoef * sss_io(:,:) 303 u_io (:,:) = zcoef * u_io (:,:) 304 v_io (:,:) = zcoef * v_io (:,:) 305 ENDIF 306 ELSE 389 307 IF(lwp) WRITE(numout,*) 390 308 IF(lwp) WRITE(numout,*) 'rst_read : LLN sea Ice Model => Ice initialization' 391 309 IF(lwp) WRITE(numout,*) 392 sst_io(:,:) = ( nfice-1 )*( tn(:,:,1) + rt0 ) !!bug a explanation is needed here! 393 sss_io(:,:) = ( nfice-1 )* sn(:,:,1) 310 zcoef = REAL( nfice-1, wp ) 311 sst_io(:,:) = zcoef *( tn(:,:,1) + rt0 ) !!bug a explanation is needed here! 312 sss_io(:,:) = zcoef * sn(:,:,1) 313 zcoef = 0.5 * REAL( nfice-1, wp ) 394 314 DO jj = 2, jpj 395 DO ji = 2, jpi396 u_io(ji,jj) = ( nfice-1 ) * 0.5* ( un(ji-1,jj ,1) + un(ji-1,jj-1,1) )397 v_io(ji,jj) = ( nfice-1 ) * 0.5* ( vn(ji ,jj-1,1) + vn(ji-1,jj-1,1) )315 DO ji = fs_2, jpi ! vector opt. 316 u_io(ji,jj) = zcoef * ( un(ji-1,jj ,1) + un(ji-1,jj-1,1) ) 317 v_io(ji,jj) = zcoef * ( vn(ji ,jj-1,1) + vn(ji-1,jj-1,1) ) 398 318 END DO 399 319 END DO … … 405 325 # if defined key_flx_bulk_monthly || defined key_flx_bulk_daily 406 326 ! Louvain La Neuve Sea Ice Model 407 ios = iom_varid( inum, 'nfbulk' ) 408 IF( ios > 0 ) then 409 CALL iom_get( inum, jpdom_unknown, 'nfbulk' , zzz ) 410 CALL iom_get( inum, jpdom_local, 'gsst' , gsst ) 411 zinfo(1, 1, 9) = zzz(1, 1, 1) 412 ENDIF 413 IF( zinfo(1, 1, 9) /= FLOAT(nfbulk) .OR. ios == 0 ) THEN 327 IF( iom_varid( numror, 'nfbulk' ) > 0 ) THEN 328 CALL iom_get( numror , 'nfbulk', znfbulk ) ! bulk computation frequency 329 CALL iom_get( numror, jpdom_local, 'gsst' , gsst ) 330 IF( znfbulk /= REAL(nfbulk, wp) ) THEN ! if you change nfbulk between 2 runs 331 zcoef = REAL( nfbulk-1, wp ) / znfbulk 332 gsst(:,:) = zcoef * gsst(:,:) 333 ENDIF 334 ELSE 414 335 IF(lwp) WRITE(numout,*) 415 336 IF(lwp) WRITE(numout,*) 'rst_read : LLN sea Ice Model => Ice initialization' 416 337 IF(lwp) WRITE(numout,*) 417 gsst(:,:) = 0. 418 gsst(:,:) = gsst(:,:) + ( nfbulk-1 )*( tn(:,:,1) + rt0 ) 338 gsst(:,:) = REAL( nfbulk - 1, wp )*( tn(:,:,1) + rt0 ) 419 339 ENDIF 420 340 # endif 421 341 422 CALL iom_close( inum ) 423 424 ! In case of restart with neuler = 0 then put all before fields = to now fields 425 IF ( neuler == 0 ) THEN 426 tb(:,:,:)=tn(:,:,:) 427 sb(:,:,:)=sn(:,:,:) 428 ub(:,:,:)=un(:,:,:) 429 vb(:,:,:)=vn(:,:,:) 430 rotb(:,:,:)=rotn(:,:,:) 431 hdivb(:,:,:)=hdivn(:,:,:) 432 #if defined key_dynspg_rl 433 ! rigid lid 434 bsfb(:,:)=bsfn(:,:) 435 #else 436 ! free surface formulation (eta) 437 sshb(:,:)=sshn(:,:) 342 !!sm: end of TO BE MOVED IN NEW SURFACE MODULE... 343 ! 344 END SUBROUTINE rst_read 345 438 346 #endif 439 ENDIF 440 441 END SUBROUTINE rst_read 442 443 #endif 347 444 348 !!===================================================================== 445 349 END MODULE restart -
trunk/NEMO/OPA_SRC/step.F90
r503 r508 4 4 !! Time-stepping : manager of the ocean, tracer and ice time stepping 5 5 !!====================================================================== 6 !! History : 7 !! ! 91-03 () Original code 8 !! ! 91-11 (G. Madec) 9 !! ! 92-06 (M. Imbard) add a first output record 10 !! ! 96-04 (G. Madec) introduction of dynspg 11 !! ! 96-04 (M.A. Foujols) introduction of passive tracer 12 !! 8.0 ! 97-06 (G. Madec) new architecture of call 13 !! 8.2 ! 97-06 (G. Madec, M. Imbard, G. Roullet) free surface 14 !! 8.2 ! 99-02 (G. Madec, N. Grima) hpg implicit 15 !! 8.2 ! 00-07 (J-M Molines, M. Imbard) Open Bondary Conditions 16 !! 9.0 ! 02-06 (G. Madec) free form, suppress macro-tasking 17 !! " ! 04-08 (C. Talandier) New trends organization 18 !! " ! 05-01 (C. Ethe) Add the KPP closure scheme 19 !! " ! 05-11 (V. Garnier) Surface pressure gradient organization 20 !! " ! 05-11 (G. Madec) Reorganisation of tra and dyn calls 6 !! History : ! 91-03 () Original code 7 !! ! 91-11 (G. Madec) 8 !! ! 92-06 (M. Imbard) add a first output record 9 !! ! 96-04 (G. Madec) introduction of dynspg 10 !! ! 96-04 (M.A. Foujols) introduction of passive tracer 11 !! 8.0 ! 97-06 (G. Madec) new architecture of call 12 !! 8.2 ! 97-06 (G. Madec, M. Imbard, G. Roullet) free surface 13 !! 8.2 ! 99-02 (G. Madec, N. Grima) hpg implicit 14 !! 8.2 ! 00-07 (J-M Molines, M. Imbard) Open Bondary Conditions 15 !! 9.0 ! 02-06 (G. Madec) free form, suppress macro-tasking 16 !! " " ! 04-08 (C. Talandier) New trends organization 17 !! " " ! 05-01 (C. Ethe) Add the KPP closure scheme 18 !! " " ! 05-11 (V. Garnier) Surface pressure gradient organization 19 !! " " ! 05-11 (G. Madec) Reorganisation of tra and dyn calls 20 !! " " ! 06-01 (L. Debreu, C. Mazauric) Agrif implementation 21 !! " " ! 06-07 (S. Masson) restart using iom 21 22 !!---------------------------------------------------------------------- 22 23 !! stp : OPA system time-stepping 23 24 !!---------------------------------------------------------------------- 24 !! * Modules used25 25 USE oce ! ocean dynamics and tracers variables 26 26 USE dom_oce ! ocean space and time domain variables … … 30 30 USE cpl_oce ! coupled ocean-atmosphere variables 31 31 USE in_out_manager ! I/O manager 32 USE iom 32 33 USE lbclnk 33 34 … … 126 127 PRIVATE 127 128 128 !! * Routine accessibility129 129 PUBLIC stp ! called by opa.F90 130 130 … … 135 135 !! OPA 9.0 , LOCEAN-IPSL (2005) 136 136 !! $Header$ 137 !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt137 !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 138 138 !!---------------------------------------------------------------------- 139 139 140 140 CONTAINS 141 141 142 #if !defined key_agrif 142 #if defined key_agrif 143 SUBROUTINE stp( ) 144 #else 143 145 SUBROUTINE stp( kstp ) 144 #else145 SUBROUTINE stp( )146 146 #endif 147 147 !!---------------------------------------------------------------------- … … 160 160 !! -7- Compute the diagnostics variables (rd,N2, div,cur,w) 161 161 !! -8- Outputs and diagnostics 162 !!163 162 !!---------------------------------------------------------------------- 164 163 !! * Arguments 165 #if !defined key_agrif 164 #if defined key_agrif 165 INTEGER :: kstp ! ocean time-step index 166 #else 166 167 INTEGER, INTENT( in ) :: kstp ! ocean time-step index 167 #else168 INTEGER :: kstp ! ocean time-step index169 168 #endif 170 169 … … 176 175 kstp = nit000 + Agrif_Nb_Step() 177 176 ! IF ( Agrif_Root() .and. lwp) Write(*,*) '---' 178 ! IF (lwp) Write(*,*) 'Grid N °',Agrif_Fixed(),' time step ',kstp177 ! IF (lwp) Write(*,*) 'Grid Number',Agrif_Fixed(),' time step ',kstp 179 178 #endif 180 179 indic = 1 ! reset to no error condition 180 181 181 adatrj = adatrj + rdt/86400._wp 182 182 183 183 CALL day( kstp ) ! Calendar 184 185 CALL rst_opn( kstp ) ! Open the restart file 184 186 185 187 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> … … 326 328 CALL tra_ldf ( kstp ) ! lateral mixing 327 329 #if defined key_agrif 328 IF (.NOT. Agrif_Root()) CALL Agrif_Sponge_tra( kstp )! tracers sponge330 IF(.NOT. Agrif_Root()) CALL Agrif_Sponge_tra( kstp ) ! tracers sponge 329 331 #endif 330 332 CALL tra_zdf ( kstp ) ! vertical mixing … … 363 365 CALL dyn_ldf( kstp ) ! lateral mixing 364 366 #if defined key_agrif 365 IF (.NOT. Agrif_Root()) CALL Agrif_Sponge_dyn( kstp )! momemtum sponge367 IF(.NOT. Agrif_Root()) CALL Agrif_Sponge_dyn( kstp ) ! momemtum sponge 366 368 #endif 367 369 CALL dyn_hpg( kstp ) ! horizontal gradient of Hydrostatic pressure … … 400 402 401 403 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 402 ! Control, diagnostics and outputs 403 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 404 ! N.B. ua, va, ta, sa arrays are used as workspace in this section 405 !----------------------------------------------------------------------- 406 407 ! ! Time loop: control and print 408 CALL stp_ctl( kstp, indic ) 409 IF ( indic < 0 ) CALL ctl_stop( 'step: indic < 0' ) 410 411 IF ( nstop == 0 ) THEN 412 ! ! Diagnostics: 404 ! Control, and restarts 405 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 406 ! N.B. ua, va, ta, sa arrays are used as workspace in this section 407 !----------------------------------------------------------------------- 408 ! ! Time loop: control and print 409 CALL stp_ctl( kstp, indic ) 410 IF( indic < 0 ) CALL ctl_stop( 'step: indic < 0' ) 411 412 IF( kstp == nit000 ) CALL iom_close( numror ) ! close input ocean restart file 413 IF( lrst_oce ) CALL rst_write ( kstp ) ! write output ocean restart file 414 IF( lk_obc ) CALL obc_rst_wri( kstp ) ! write open boundary restart file 415 IF( lk_trdmld ) CALL trd_mld_rst_write( kstp ) ! ocean model: restart file output for trends diagnostics 416 417 418 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 419 ! diagnostics and outputs 420 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 421 ! N.B. ua, va, ta, sa arrays are used as workspace in this section 422 !----------------------------------------------------------------------- 423 424 IF ( nstop == 0 ) THEN ! Diagnostics 413 425 IF( lk_floats ) CALL flo_stp( kstp ) ! drifting Floats 414 426 IF( lk_trddyn ) CALL trd_dwr( kstp ) ! trends: dynamics … … 423 435 IF( ln_diaptr ) CALL dia_ptr( kstp ) ! Poleward TRansports diagnostics 424 436 425 ! ! save and outputs 426 CALL rst_write ( kstp ) ! ocean model: restart file output 427 IF( lk_obc ) CALL obc_rst_wri( kstp ) ! ocean model: open boundary restart file output 428 IF( lk_trdmld ) CALL trd_mld_rst_write( kstp ) ! ocean model: restart file output for trends diagnostics 437 ! ! Outputs 429 438 CALL dia_wri ( kstp, indic ) ! ocean model: outputs 430 431 439 ENDIF 432 440 … … 436 444 437 445 IF( lk_cpl ) CALL cpl_stp( kstp ) ! coupled mode : field exchanges 438 446 ! 439 447 END SUBROUTINE stp 440 448
Note: See TracChangeset
for help on using the changeset viewer.