New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
nemogcm.F90 in branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC – NEMO

source: branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90 @ 4409

Last change on this file since 4409 was 4409, checked in by trackstand2, 10 years ago

Changes to allow jpk to be modified to deepest level within a subdomain. jpkorig holds original value.

  • Property svn:keywords set to Id
File size: 41.1 KB
RevLine 
[2496]1MODULE nemogcm
[2442]2   !!======================================================================
[2496]3   !!                       ***  MODULE nemogcm   ***
4   !! Ocean system   : NEMO GCM (ocean dynamics, on-line tracers, biochemistry and sea-ice)
[2442]5   !!======================================================================
[1593]6   !! History :  OPA  ! 1990-10  (C. Levy, G. Madec)  Original code
7   !!            7.0  ! 1991-11  (M. Imbard, C. Levy, G. Madec)
8   !!            7.1  ! 1993-03  (M. Imbard, C. Levy, G. Madec, O. Marti, M. Guyon, A. Lazar,
[2104]9   !!                             P. Delecluse, C. Perigaud, G. Caniaux, B. Colot, C. Maes) release 7.1
[1593]10   !!             -   ! 1992-06  (L.Terray)  coupling implementation
11   !!             -   ! 1993-11  (M.A. Filiberti) IGLOO sea-ice
12   !!            8.0  ! 1996-03  (M. Imbard, C. Levy, G. Madec, O. Marti, M. Guyon, A. Lazar,
[2104]13   !!                             P. Delecluse, L.Terray, M.A. Filiberti, J. Vialar, A.M. Treguier, M. Levy) release 8.0
[1593]14   !!            8.1  ! 1997-06  (M. Imbard, G. Madec)
15   !!            8.2  ! 1999-11  (M. Imbard, H. Goosse)  LIM sea-ice model
16   !!                 ! 1999-12  (V. Thierry, A-M. Treguier, M. Imbard, M-A. Foujols)  OPEN-MP
17   !!                 ! 2000-07  (J-M Molines, M. Imbard)  Open Boundary Conditions  (CLIPPER)
18   !!   NEMO     1.0  ! 2002-08  (G. Madec)  F90: Free form and modules
19   !!             -   ! 2004-06  (R. Redler, NEC CCRLE, Germany) add OASIS[3/4] coupled interfaces
20   !!             -   ! 2004-08  (C. Talandier) New trends organization
21   !!             -   ! 2005-06  (C. Ethe) Add the 1D configuration possibility
22   !!             -   ! 2005-11  (V. Garnier) Surface pressure gradient organization
23   !!             -   ! 2006-03  (L. Debreu, C. Mazauric)  Agrif implementation
24   !!             -   ! 2006-04  (G. Madec, R. Benshila)  Step reorganization
25   !!             -   ! 2007-07  (J. Chanut, A. Sellar) Unstructured open boundaries (BDY)
26   !!            3.2  ! 2009-08  (S. Masson)  open/write in the listing file in mpp
[2236]27   !!            3.3  ! 2010-05  (K. Mogensen, A. Weaver, M. Martin, D. Lea) Assimilation interface
28   !!             -   ! 2010-10  (C. Ethe, G. Madec) reorganisation of initialisation phase
[2715]29   !!            4.0  ! 2011-01  (A. R. Porter, STFC Daresbury) dynamical allocation
[1593]30   !!----------------------------------------------------------------------
[3]31
32   !!----------------------------------------------------------------------
[2496]33   !!   nemo_gcm       : solve ocean dynamics, tracer, biogeochemistry and/or sea-ice
34   !!   nemo_init      : initialization of the NEMO system
35   !!   nemo_ctl       : initialisation of the contol print
36   !!   nemo_closefile : close remaining open files
[2715]37   !!   nemo_alloc     : dynamical allocation
38   !!   nemo_partition : calculate MPP domain decomposition
[3187]39   !!   sqfact         : calculate factors of the no. of MPI processes
[3]40   !!----------------------------------------------------------------------
[2382]41   USE step_oce        ! module used in the ocean time stepping module
[888]42   USE sbc_oce         ! surface boundary condition: ocean
[2392]43   USE cla             ! cross land advection               (tra_cla routine)
[3]44   USE domcfg          ! domain configuration               (dom_cfg routine)
45   USE mppini          ! shared/distributed memory setting (mpp_init routine)
46   USE domain          ! domain initialization             (dom_init routine)
47   USE obcini          ! open boundary cond. initialization (obc_ini routine)
[911]48   USE bdyini          ! unstructured open boundary cond. initialization (bdy_init routine)
[3]49   USE istate          ! initial state setting          (istate_init routine)
50   USE ldfdyn          ! lateral viscosity setting      (ldfdyn_init routine)
51   USE ldftra          ! lateral diffusivity setting    (ldftra_init routine)
[2392]52   USE zdfini          ! vertical physics setting          (zdf_init routine)
[3]53   USE phycst          ! physical constant                  (par_cst routine)
[2236]54   USE trdmod          ! momentum/tracers trends       (trd_mod_init routine)
55   USE asminc          ! assimilation increments       (asm_inc_init routine)
56   USE asmtrj          ! writing out state trajectory
57   USE sshwzv          ! vertical velocity used in asm
58   USE diaptr          ! poleward transports           (dia_ptr_init routine)
59   USE diaobs          ! Observation diagnostics       (dia_obs_init routine)
[2496]60   USE step            ! NEMO time-stepping                 (stp     routine)
[532]61#if defined key_oasis3
[1359]62   USE cpl_oasis3      ! OASIS3 coupling
[599]63#elif defined key_oasis4
[1359]64   USE cpl_oasis4      ! OASIS4 coupling (not working)
[532]65#endif
[900]66   USE c1d             ! 1D configuration
67   USE step_c1d        ! Time stepping loop for the 1D configuration
[1594]68#if defined key_top
[1593]69   USE trcini          ! passive tracer initialisation
[1594]70#endif
[1593]71   USE lib_mpp         ! distributed memory computing
[1412]72#if defined key_iomput
73   USE mod_ioclient
[1359]74#endif
[3432]75   USE partition_mod   ! irregular domain partitioning
76   USE timing, ONLY: timing_init, timing_finalize, timing_disable, timing_enable
[268]77
[3837]78#define ARPDEBUG
[3432]79
[2715]80   IMPLICIT NONE
[3]81   PRIVATE
82
[2496]83   PUBLIC   nemo_gcm    ! called by model.F90
84   PUBLIC   nemo_init   ! needed by AGRIF
[467]85
[2498]86   CHARACTER(lc) ::   cform_aaa="( /, 'AAAAAAAA', / ) "     ! flag for output listing
[1593]87
[3]88   !!----------------------------------------------------------------------
[2715]89   !! NEMO/OPA 4.0 , NEMO Consortium (2011)
[2392]90   !! $Id$
[2329]91   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
[3]92   !!----------------------------------------------------------------------
93CONTAINS
94
[2496]95   SUBROUTINE nemo_gcm
[3]96      !!----------------------------------------------------------------------
[2496]97      !!                     ***  ROUTINE nemo_gcm  ***
[3]98      !!
[2496]99      !! ** Purpose :   NEMO solves the primitive equations on an orthogonal
[1593]100      !!              curvilinear mesh on the sphere.
[3]101      !!
102      !! ** Method  : - model general initialization
103      !!              - launch the time-stepping (stp routine)
[1593]104      !!              - finalize the run by closing files and communications
[3]105      !!
[2715]106      !! References : Madec, Delecluse, Imbard, and Levy, 1997:  internal report, IPSL.
[1593]107      !!              Madec, 2008, internal report, IPSL.
[3]108      !!----------------------------------------------------------------------
109      INTEGER ::   istp       ! time step index
[389]110      !!----------------------------------------------------------------------
[2382]111      !
[392]112#if defined key_agrif
[1593]113      CALL Agrif_Init_Grids()      ! AGRIF: set the meshes
[389]114#endif
115
[1593]116      !                            !-----------------------!
[2496]117      CALL nemo_init               !==  Initialisations  ==!
[1593]118      !                            !-----------------------!
[2715]119#if defined key_agrif
120      CALL Agrif_Declare_Var       ! AGRIF: set the meshes
121# if defined key_top
122      CALL Agrif_Declare_Var_Top   ! AGRIF: set the meshes
123# endif
124#endif
[682]125      ! check that all process are still there... If some process have an error,
126      ! they will never enter in step and other processes will wait until the end of the cpu time!
[900]127      IF( lk_mpp )   CALL mpp_max( nstop )
[682]128
[1593]129      IF(lwp) WRITE(numout,cform_aaa)   ! Flag AAAAAAA
130
[3432]131      CALL timing_enable()
[1593]132      !                            !-----------------------!
133      !                            !==   time stepping   ==!
134      !                            !-----------------------!
[900]135      istp = nit000
[2236]136#if defined key_c1d
[389]137         DO WHILE ( istp <= nitend .AND. nstop == 0 )
[900]138            CALL stp_c1d( istp )
[389]139            istp = istp + 1
140         END DO
[2236]141#else
142          IF( lk_asminc ) THEN
143             IF( ln_bkgwri ) CALL asm_bkg_wri( nit000 - 1 )    ! Output background fields
144             IF( ln_trjwri ) CALL asm_trj_wri( nit000 - 1 )    ! Output trajectory fields
145             IF( ln_asmdin ) THEN                        ! Direct initialization
146                IF( ln_trainc ) CALL tra_asm_inc( nit000 - 1 )    ! Tracers
147                IF( ln_dyninc ) THEN
148                   CALL dyn_asm_inc( nit000 - 1 )    ! Dynamics
149                   IF ( ln_asmdin ) CALL ssh_wzv ( nit000 - 1 )      ! update vertical velocity
150                ENDIF
151                IF( ln_sshinc ) CALL ssh_asm_inc( nit000 - 1 )    ! SSH
152             ENDIF
153          ENDIF
154       
[389]155         DO WHILE ( istp <= nitend .AND. nstop == 0 )
[392]156#if defined key_agrif
[1593]157            CALL Agrif_Step( stp )           ! AGRIF: time stepping
[389]158#else
[1593]159            CALL stp( istp )                 ! standard time stepping
[389]160#endif
161            istp = istp + 1
[900]162            IF( lk_mpp )   CALL mpp_max( nstop )
[389]163         END DO
[2236]164#endif
165
166      IF( lk_diaobs ) CALL dia_obs_wri
[1593]167       
168      !                            !------------------------!
169      !                            !==  finalize the run  ==!
170      !                            !------------------------!
171      IF(lwp) WRITE(numout,cform_aaa)   ! Flag AAAAAAA
172      !
173      IF( nstop /= 0 .AND. lwp ) THEN   ! error print
[682]174         WRITE(numout,cform_err)
175         WRITE(numout,*) nstop, ' error have been found' 
[389]176      ENDIF
[1593]177      !
[3432]178      CALL timing_finalize              ! Timing report
179
[2496]180      CALL nemo_closefile
[532]181#if defined key_oasis3 || defined key_oasis4
[1976]182      CALL cpl_prism_finalize           ! end coupling and mpp communications with OASIS
[532]183#else
[1593]184      IF( lk_mpp )   CALL mppstop       ! end mpp communications
[532]185#endif
[900]186      !
[2496]187   END SUBROUTINE nemo_gcm
[389]188
189
[2496]190   SUBROUTINE nemo_init
[389]191      !!----------------------------------------------------------------------
[2496]192      !!                     ***  ROUTINE nemo_init  ***
[389]193      !!
[2496]194      !! ** Purpose :   initialization of the NEMO GCM
[389]195      !!----------------------------------------------------------------------
[2715]196      INTEGER ::   ji            ! dummy loop indices
197      INTEGER ::   ilocal_comm   ! local integer
[3432]198      CHARACTER(len=80), DIMENSION(24) ::   cltxt
[1593]199      !!
[1601]200      NAMELIST/namctl/ ln_ctl  , nn_print, nn_ictls, nn_ictle,   &
[1976]201         &             nn_isplt, nn_jsplt, nn_jctls, nn_jctle, nn_bench
[3]202      !!----------------------------------------------------------------------
[1593]203      !
[3432]204      cltxt(:) = ''
[2496]205      !
[1593]206      !                             ! open Namelist file
[1581]207      CALL ctl_opn( numnam, 'namelist', 'OLD', 'FORMATTED', 'SEQUENTIAL', -1, 6, .FALSE. )
[1593]208      !
209      READ( numnam, namctl )        ! Namelist namctl : Control prints & Benchmark
210      !
211      !                             !--------------------------------------------!
212      !                             !  set communicator & select the local node  !
213      !                             !--------------------------------------------!
[1412]214#if defined key_iomput
[2200]215      IF( Agrif_Root() ) THEN
[1412]216# if defined key_oasis3 || defined key_oasis4
[2715]217         CALL cpl_prism_init( ilocal_comm )                 ! nemo local communicator given by oasis
[1412]218# endif
[2715]219         CALL  init_ioclient( ilocal_comm )                 ! exchange io_server nemo local communicator with the io_server
[2200]220      ENDIF
[2715]221      narea = mynode( cltxt, numnam, nstop, ilocal_comm )   ! Nodes selection
[532]222#else
[1412]223# if defined key_oasis3 || defined key_oasis4
[2236]224      IF( Agrif_Root() ) THEN
[2715]225         CALL cpl_prism_init( ilocal_comm )                 ! nemo local communicator given by oasis
[2236]226      ENDIF
[2715]227      narea = mynode( cltxt, numnam, nstop, ilocal_comm )   ! Nodes selection (control print return in cltxt)
[1412]228# else
[2082]229      ilocal_comm = 0
[2715]230      narea = mynode( cltxt, numnam, nstop )                 ! Nodes selection (control print return in cltxt)
[1412]231# endif
[532]232#endif
[2715]233      narea = narea + 1                                     ! mynode return the rank of proc (0 --> jpnij -1 )
[3]234
[2715]235      lwp = (narea == 1) .OR. ln_ctl                        ! control of all listing output print
[1579]236
[3837]237      ! Calculate domain z dimensions as needed when partitioning.
238      ! This used to be done in par_oce.F90 when they were parameters rather
239      ! than variables
240      IF( Agrif_Root() ) THEN
241         jpk = jpkdta                                             ! third dim
242         jpkm1 = jpk-1                                            ! inner domain indices
[4409]243         jpkorig = jpk                        ! Copy of jpk that is NOT modified
244         jpkf    = jpk                        ! Max depth of this sub-domain. Initially set to jpk here
245                                              ! but altered later in domzgr
[3837]246      ENDIF
247
[3432]248      CALL timing_init                                      ! Init timing module
249      CALL timing_disable                                   ! but disable during startup
250
[2715]251      ! If dimensions of processor grid weren't specified in the namelist file
252      ! then we calculate them here now that we have our communicator size
253      IF( (jpni < 1) .OR. (jpnj < 1) )THEN
254#if   defined key_mpp_mpi
[3432]255#if   defined key_mpp_rkpart
256         IF( Agrif_Root() ) CALL nemo_recursive_partition(mppsize)
257#else
[2715]258         IF( Agrif_Root() ) CALL nemo_partition(mppsize)
[3432]259#endif
[2715]260#else
261         jpni  = 1
262         jpnj  = 1
263         jpnij = jpni*jpnj
264#endif
[3837]265
266#if   defined key_mpp_rkpart
267      ELSE
[4409]268         CALL ctl_stop( 'STOP', &
269                        'nemo_init : invalid inputs in namelist - cannot specify jpn{i,j}>0 &
270                        & when using recursive k-section paritioning!' )
[3837]271#endif
[2715]272      END IF
273
274      ! Calculate domain dimensions given calculated jpni and jpnj
275      ! This used to be done in par_oce.F90 when they were parameters rather
276      ! than variables
277      IF( Agrif_Root() ) THEN
[3432]278#if ! defined key_mpp_rkpart
[2715]279         jpi = ( jpiglo-2*jpreci + (jpni-1) ) / jpni + 2*jpreci   ! first  dim.
280         jpj = ( jpjglo-2*jprecj + (jpnj-1) ) / jpnj + 2*jprecj   ! second dim.
281         jpim1 = jpi-1                                            ! inner domain indices
282         jpjm1 = jpj-1                                            !   "           "
283         jpij  = jpi*jpj                                          !  jpi x j
[3432]284#endif
[2715]285      ENDIF
286
[1593]287      IF(lwp) THEN                            ! open listing units
288         !
[1581]289         CALL ctl_opn( numout, 'ocean.output', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, 6, .FALSE., narea )
[1593]290         !
[1579]291         WRITE(numout,*)
[1593]292         WRITE(numout,*) '         CNRS - NERC - Met OFFICE - MERCATOR-ocean'
293         WRITE(numout,*) '                       NEMO team'
[1579]294         WRITE(numout,*) '            Ocean General Circulation Model'
[2104]295         WRITE(numout,*) '                  version 3.3  (2010) '
[1579]296         WRITE(numout,*)
297         WRITE(numout,*)
[3432]298         DO ji = 1, SIZE(cltxt,1) 
[1593]299            IF( TRIM(cltxt(ji)) /= '' )   WRITE(numout,*) cltxt(ji)      ! control print of mynode
[1579]300         END DO
[1593]301         WRITE(numout,cform_aaa)                                         ! Flag AAAAAAA
302         !
[473]303      ENDIF
[2715]304
305      ! Now we know the dimensions of the grid and numout has been set we can
306      ! allocate arrays
307      CALL nemo_alloc()
308
[2496]309      !                             !-------------------------------!
310      !                             !  NEMO general initialization  !
311      !                             !-------------------------------!
[473]312
[2496]313      CALL nemo_ctl                          ! Control prints & Benchmark
[531]314
[2082]315      !                                      ! Domain decomposition
[3432]316#if defined key_mpp_rkpart
317                                          CALL mpp_init3     ! Remainder of set-up for
318                                                             ! recursive partitioning
319#else
[1593]320      IF( jpni*jpnj == jpnij ) THEN   ;   CALL mpp_init      ! standard cutting out
321      ELSE                            ;   CALL mpp_init2     ! eliminate land processors
[3]322      ENDIF
[3432]323#endif
[2382]324      !
[2082]325      !                                      ! General initialization
[3432]326!                            CALL     timing_init! Timing module
[2027]327                            CALL     phy_cst    ! Physical constants
328                            CALL     eos_init   ! Equation of state
329                            CALL     dom_cfg    ! Domain configuration
330                            CALL     dom_init   ! Domain
[413]331
[2027]332      IF( ln_ctl        )   CALL prt_ctl_init   ! Print control
333
334      IF( lk_obc        )   CALL     obc_init   ! Open boundaries
335      IF( lk_bdy        )   CALL     bdy_init   ! Unstructured open boundaries
336
337                            CALL  istate_init   ! ocean initial state (Dynamics and tracers)
338
[3]339      !                                     ! Ocean physics
[2104]340                            CALL     sbc_init   ! Forcings : surface module
[2082]341      !                                         ! Vertical physics
342                            CALL     zdf_init      ! namelist read
343                            CALL zdf_bfr_init      ! bottom friction
344      IF( lk_zdfric     )   CALL zdf_ric_init      ! Richardson number dependent Kz
[2329]345      IF( lk_zdftke     )   CALL zdf_tke_init      ! TKE closure scheme
346      IF( lk_zdfgls     )   CALL zdf_gls_init      ! GLS closure scheme
347      IF( lk_zdfkpp     )   CALL zdf_kpp_init      ! KPP closure scheme
[2082]348      IF( lk_zdftmx     )   CALL zdf_tmx_init      ! tidal vertical mixing
[2027]349      IF( lk_zdfddm .AND. .NOT. lk_zdfkpp )   & 
[2082]350         &                  CALL zdf_ddm_init      ! double diffusive mixing
351      !                                         ! Lateral physics
352                            CALL ldf_tra_init      ! Lateral ocean tracer physics
353                            CALL ldf_dyn_init      ! Lateral ocean momentum physics
[2392]354      IF( lk_ldfslp     )   CALL ldf_slp_init      ! slope of lateral mixing
[2082]355
[2027]356      !                                     ! Active tracers
357                            CALL tra_qsr_init   ! penetrative solar radiation qsr
[2325]358                            CALL tra_bbc_init   ! bottom heat flux
[2027]359      IF( lk_trabbl     )   CALL tra_bbl_init   ! advective (and/or diffusive) bottom boundary layer scheme
360      IF( lk_tradmp     )   CALL tra_dmp_init   ! internal damping trends
361                            CALL tra_adv_init   ! horizontal & vertical advection
362                            CALL tra_ldf_init   ! lateral mixing
363                            CALL tra_zdf_init   ! vertical mixing and after tracer fields
364
365      !                                     ! Dynamics
366                            CALL dyn_adv_init   ! advection (vector or flux form)
[2104]367                            CALL dyn_vor_init   ! vorticity term including Coriolis
[2027]368                            CALL dyn_ldf_init   ! lateral mixing
[2104]369                            CALL dyn_hpg_init   ! horizontal gradient of Hydrostatic pressure
[2027]370                            CALL dyn_zdf_init   ! vertical diffusion
371                            CALL dyn_spg_init   ! surface pressure gradient
[2392]372                           
373      !                                     ! Misc. options
374      IF( nn_cla == 1   )   CALL cla_init       ! Cross Land Advection
375     
[1594]376#if defined key_top
[2027]377      !                                     ! Passive tracers
[2082]378                            CALL     trc_init
[1594]379#endif
[2382]380      !                                     ! Diagnostics
[2027]381                            CALL     iom_init   ! iom_put initialization
[2392]382      IF( lk_floats     )   CALL     flo_init   ! drifting Floats
383      IF( lk_diaar5     )   CALL dia_ar5_init   ! ar5 diag
[2027]384                            CALL dia_ptr_init   ! Poleward TRansports initialization
[2148]385                            CALL dia_hsb_init   ! heat content, salt content and volume budgets
[2027]386                            CALL trd_mod_init   ! Mixed-layer/Vorticity/Integral constraints trends
[2392]387      IF( lk_diaobs     ) THEN                  ! Observation & model comparison
[2382]388                            CALL dia_obs_init            ! Initialize observational data
389                            CALL dia_obs( nit000 - 1 )   ! Observation operator for restart
[2236]390      ENDIF     
[2382]391      !                                     ! Assimilation increments
[2392]392      IF( lk_asminc     )   CALL asm_inc_init   ! Initialize assimilation increments
[2382]393      IF(lwp) WRITE(numout,*) 'Euler time step switch is ', neuler
[1593]394      !
[2496]395   END SUBROUTINE nemo_init
[467]396
397
[2496]398   SUBROUTINE nemo_ctl
[467]399      !!----------------------------------------------------------------------
[2496]400      !!                     ***  ROUTINE nemo_ctl  ***
[467]401      !!
[2442]402      !! ** Purpose :   control print setting
[467]403      !!
[2442]404      !! ** Method  : - print namctl information and check some consistencies
[467]405      !!----------------------------------------------------------------------
[2442]406      !
[2496]407      IF(lwp) THEN                  ! control print
[531]408         WRITE(numout,*)
[2496]409         WRITE(numout,*) 'nemo_ctl: Control prints & Benchmark'
[531]410         WRITE(numout,*) '~~~~~~~ '
[1593]411         WRITE(numout,*) '   Namelist namctl'
[1601]412         WRITE(numout,*) '      run control (for debugging)     ln_ctl     = ', ln_ctl
413         WRITE(numout,*) '      level of print                  nn_print   = ', nn_print
414         WRITE(numout,*) '      Start i indice for SUM control  nn_ictls   = ', nn_ictls
415         WRITE(numout,*) '      End i indice for SUM control    nn_ictle   = ', nn_ictle
416         WRITE(numout,*) '      Start j indice for SUM control  nn_jctls   = ', nn_jctls
417         WRITE(numout,*) '      End j indice for SUM control    nn_jctle   = ', nn_jctle
418         WRITE(numout,*) '      number of proc. following i     nn_isplt   = ', nn_isplt
419         WRITE(numout,*) '      number of proc. following j     nn_jsplt   = ', nn_jsplt
420         WRITE(numout,*) '      benchmark parameter (0/1)       nn_bench   = ', nn_bench
[531]421      ENDIF
[2442]422      !
[1601]423      nprint    = nn_print          ! convert DOCTOR namelist names into OLD names
424      nictls    = nn_ictls
425      nictle    = nn_ictle
426      njctls    = nn_jctls
427      njctle    = nn_jctle
428      isplt     = nn_isplt
429      jsplt     = nn_jsplt
430      nbench    = nn_bench
[2442]431      !                             ! Parameter control
[1593]432      !
433      IF( ln_ctl ) THEN                 ! sub-domain area indices for the control prints
434         IF( lk_mpp ) THEN
[2496]435            isplt = jpni   ;   jsplt = jpnj   ;   ijsplt = jpni*jpnj   ! the domain is forced to the real split domain
[531]436         ELSE
437            IF( isplt == 1 .AND. jsplt == 1  ) THEN
[1593]438               CALL ctl_warn( ' - isplt & jsplt are equal to 1',   &
439                  &           ' - the print control will be done over the whole domain' )
[531]440            ENDIF
[1593]441            ijsplt = isplt * jsplt            ! total number of processors ijsplt
[531]442         ENDIF
443         IF(lwp) WRITE(numout,*)'          - The total number of processors over which the'
444         IF(lwp) WRITE(numout,*)'            print control will be done is ijsplt : ', ijsplt
[1593]445         !
446         !                              ! indices used for the SUM control
447         IF( nictls+nictle+njctls+njctle == 0 )   THEN    ! print control done over the default area
448            lsp_area = .FALSE.                       
449         ELSE                                             ! print control done over a specific  area
[531]450            lsp_area = .TRUE.
451            IF( nictls < 1 .OR. nictls > jpiglo )   THEN
452               CALL ctl_warn( '          - nictls must be 1<=nictls>=jpiglo, it is forced to 1' )
453               nictls = 1
454            ENDIF
455            IF( nictle < 1 .OR. nictle > jpiglo )   THEN
456               CALL ctl_warn( '          - nictle must be 1<=nictle>=jpiglo, it is forced to jpiglo' )
457               nictle = jpiglo
458            ENDIF
459            IF( njctls < 1 .OR. njctls > jpjglo )   THEN
460               CALL ctl_warn( '          - njctls must be 1<=njctls>=jpjglo, it is forced to 1' )
461               njctls = 1
462            ENDIF
463            IF( njctle < 1 .OR. njctle > jpjglo )   THEN
464               CALL ctl_warn( '          - njctle must be 1<=njctle>=jpjglo, it is forced to jpjglo' )
465               njctle = jpjglo
466            ENDIF
[1593]467         ENDIF
468      ENDIF
[2442]469      !
470      IF( nbench == 1 ) THEN              ! Benchmark
[531]471         SELECT CASE ( cp_cfg )
[1593]472         CASE ( 'gyre' )   ;   CALL ctl_warn( ' The Benchmark is activated ' )
473         CASE DEFAULT      ;   CALL ctl_stop( ' The Benchmark is based on the GYRE configuration:',   &
474            &                                 ' key_gyre must be used or set nbench = 0' )
[531]475         END SELECT
476      ENDIF
[1593]477      !
[2496]478      IF( lk_c1d .AND. .NOT.lk_iomput )   CALL ctl_stop( 'nemo_ctl: The 1D configuration must be used ',   &
479         &                                               'with the IOM Input/Output manager. '         ,   &
[2442]480         &                                               'Compile with key_iomput enabled' )
[2409]481      !
[2496]482   END SUBROUTINE nemo_ctl
[467]483
484
[2496]485   SUBROUTINE nemo_closefile
[467]486      !!----------------------------------------------------------------------
[2496]487      !!                     ***  ROUTINE nemo_closefile  ***
[467]488      !!
489      !! ** Purpose :   Close the files
490      !!----------------------------------------------------------------------
[1593]491      !
492      IF( lk_mpp )   CALL mppsync
493      !
[1685]494      CALL iom_close                                 ! close all input/output files managed by iom_*
[1593]495      !
[1685]496      IF( numstp     /= -1 )   CLOSE( numstp     )   ! time-step file
497      IF( numsol     /= -1 )   CLOSE( numsol     )   ! solver file
498      IF( numnam     /= -1 )   CLOSE( numnam     )   ! oce namelist
499      IF( numnam_ice /= -1 )   CLOSE( numnam_ice )   ! ice namelist
500      IF( numevo_ice /= -1 )   CLOSE( numevo_ice )   ! ice variables (temp. evolution)
501      IF( numout     /=  6 )   CLOSE( numout     )   ! standard model output file
[1593]502      !
[2442]503      numout = 6                                     ! redefine numout in case it is used after this point...
504      !
[2496]505   END SUBROUTINE nemo_closefile
[467]506
[2715]507
508   SUBROUTINE nemo_alloc
509      !!----------------------------------------------------------------------
510      !!                     ***  ROUTINE nemo_alloc  ***
511      !!
512      !! ** Purpose :   Allocate all the dynamic arrays of the OPA modules
513      !!
514      !! ** Method  :
515      !!----------------------------------------------------------------------
516      USE diawri    , ONLY: dia_wri_alloc
517      USE dom_oce   , ONLY: dom_oce_alloc
518      USE ldfdyn_oce, ONLY: ldfdyn_oce_alloc
519      USE ldftra_oce, ONLY: ldftra_oce_alloc
520      USE trc_oce   , ONLY: trc_oce_alloc
521      USE wrk_nemo  , ONLY: wrk_alloc
[3432]522      USE exchmod   , ONLY: exchmod_alloc
[2715]523      !
524      INTEGER :: ierr
525      !!----------------------------------------------------------------------
526      !
527      ierr =        oce_alloc       ()          ! ocean
528      ierr = ierr + dia_wri_alloc   ()
529      ierr = ierr + dom_oce_alloc   ()          ! ocean domain
530      ierr = ierr + ldfdyn_oce_alloc()          ! ocean lateral  physics : dynamics
531      ierr = ierr + ldftra_oce_alloc()          ! ocean lateral  physics : tracers
532      ierr = ierr + zdf_oce_alloc   ()          ! ocean vertical physics
533      !
534      ierr = ierr + lib_mpp_alloc   (numout)    ! mpp exchanges
535      ierr = ierr + trc_oce_alloc   ()          ! shared TRC / TRA arrays
536      !
537      ierr = ierr + wrk_alloc(numout, lwp)      ! workspace
538      !
[3432]539      ierr = ierr + exchmod_alloc()             ! New mpp msg framework
540      !
[2715]541      IF( lk_mpp    )   CALL mpp_sum( ierr )
542      IF( ierr /= 0 )   CALL ctl_stop( 'STOP', 'nemo_alloc : unable to allocate standard ocean arrays' )
543      !
544   END SUBROUTINE nemo_alloc
545
546
547   SUBROUTINE nemo_partition( num_pes )
[3432]548      USE mapcomm_mod, ONLY: trimmed
[2715]549      !!----------------------------------------------------------------------
550      !!                 ***  ROUTINE nemo_partition  ***
551      !!
[3187]552      !! ** Purpose : Work out a sensible factorisation of the number of
553      !!              processors for the x and y dimensions.
[2715]554      !! ** Method  :
555      !!----------------------------------------------------------------------
556      INTEGER, INTENT(in) :: num_pes ! The number of MPI processes we have
557      !
[3187]558      INTEGER             :: ifact1, ifact2 ! factors of num_pes, ifact1 <= ifact2
[2715]559      !!----------------------------------------------------------------------
560
[3187]561      ! Factorise the number of processors into ifact1*ifact2, such that
562      ! ifact1 and ifact2 are as nearly equal as possible.
[2715]563
[3187]564      CALL sqfact( num_pes, ifact1, ifact2 )
[2715]565
[3187]566      ! Make sure that the smaller dimension of the processor grid
567      ! is given the smaller dimension of the global domain
568      IF( jpiglo <= jpjglo) THEN
569         jpni = ifact1
570         jpnj = ifact2
[2715]571      ELSE
[3187]572         jpni = ifact2
573         jpnj = ifact1
[2715]574      ENDIF
[3187]575
576      ! This should never happen
577      IF( (jpni*jpnj) /= num_pes) THEN
578         WRITE (numout, *) 'WARNING: internal error - factorisation of number of PEs failed'
579      ENDIF
580
581      ! This should only happen if num_pes is prime
582      IF( ifact1 <= 1 ) THEN
583         WRITE (numout, *) 'WARNING: factorisation of number of PEs failed'
584         WRITE (numout, *) '       : using grid of ',jpni,' x ',jpnj
585      ENDIF
[2715]586      !
587      jpnij = jpni*jpnj
588      !
[3432]589
590      ! Array that stores whether domain boundaries have been trimmed. Not used in
591      ! this case (regular domain decomp.) so set all to false.
592      ALLOCATE(trimmed(4,jpnij))
593      trimmed(:,:) = .FALSE.
594
[2715]595   END SUBROUTINE nemo_partition
596
[3432]597
598   SUBROUTINE nemo_recursive_partition( num_pes )
[3837]599      USE in_out_manager, ONLY: numnam
[4409]600      USE dom_oce,        ONLY: ln_zco
[3837]601      USE dom_oce,        ONLY: gdepw_0, gdept_0, e3w_0, e3t_0, &
602                                mig, mjg, mi0, mi1, mj0, mj1,  mbathy, bathy
603      USE domzgr,         ONLY: zgr_z, zgr_bat, namzgr, zgr_zco, zgr_zps
604      USE closea,         ONLY: dom_clo
605      USE domain,         ONLY: dom_nam
606      USE iom,            ONLY: jpiglo, jpjglo, wp, jpdom_unknown, &
607                                iom_open, iom_get, iom_close
[3849]608      USE mapcomm_mod, ONLY: ielb, ieub, pielb, pjelb, pieub, pjeub,          &
609                             iesub, jesub, jeub, ilbext, iubext, jubext,      &
[3432]610                             jlbext, pnactive, piesub, pjesub, jelb, pilbext, &
[3849]611                             piubext, pjlbext, pjubext, LAND, trimmed,        &
612                             msgtrim_z, set_num_subdomains
613      USE partition_mod, ONLY: partition_rk, partition_mca_rk, read_partition, &
614                               imask, ibotlevel, partition_mask_alloc,         &
[3837]615                               smooth_global_bathy, global_bot_level
[3432]616      USE par_oce,       ONLY: do_exchanges
617#if defined key_mpp_mpi
618      USE mpi
619#endif
620      !!----------------------------------------------------------------------
621      !!                 ***  ROUTINE nemo_recursive_partition  ***
622      !!
623      !! ** Purpose : Work out a sensible factorisation of the number of
624      !!              processors for the x and y dimensions.
625      !! ** Method  :
626      !!----------------------------------------------------------------------
627      IMPLICIT none
628      INTEGER, INTENT(in) :: num_pes ! The number of MPI processes we have
629      ! Local vars
630      INTEGER :: ierr                          ! Error flag
631      INTEGER :: inum                          ! temporary logical unit
632      INTEGER :: ii,jj,iproc                   ! Loop index
[3837]633      CHARACTER(LEN=8) :: lstr                 ! Local string for reading env. var.
634      INTEGER          :: lztrim               ! Local int for      "      "    "
[3432]635      REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdta  ! temporary data workspace
636      !!----------------------------------------------------------------------
637
[3837]638      ! Allocate masking arrays used in partitioning
639      CALL partition_mask_alloc(jpiglo,jpjglo,ierr)
[3432]640      IF(ierr /= 0)THEN
[3837]641         CALL ctl_stop('nemo_recursive_partition: failed to allocate masking arrays')
642         RETURN
643      END IF
644
645      ! Allocate local workspace array for this routine
646      ALLOCATE(zdta(jpiglo,jpjglo), Stat=ierr)
647      IF(ierr /= 0)THEN
[3432]648         CALL ctl_stop('nemo_recursive_partition: failed to allocate workspace arrays')
649         RETURN
650      END IF
651
[3849]652      ! Check whether user has specified halo trimming in z via environment
653      ! variable.
[3837]654      ! Halo trimming in z is on by default
655      msgtrim_z = .TRUE.
[3849]656      CALL GET_ENVIRONMENT_VARIABLE(NAME='NEMO_MSGTRIM_Z', VALUE=lstr, &
657                                    STATUS=ierr)
[3837]658      IF( ierr == 0)THEN
[4409]659         READ(lstr,FMT="(I10)",IOSTAT=ierr) lztrim
[3837]660         IF(ierr == 0)THEN
661            IF (lztrim == 0) msgtrim_z = .FALSE.
662         ELSE
663            CALL ctl_warn('nemo_recursive_partition: failed to parse value of NEMO_MSGTRIM_Z environment variable: '//TRIM(lstr))
664         END IF
665      END IF
666
[3849]667      IF(lwp) WRITE(*,*) 'ARPDBG: msgtrim_z = ',msgtrim_z
[3837]668
669      ! ============================
670      ! Generate a global mask from the model bathymetry
671      ! ============================
[3432]672
[3837]673      ! Read the z-coordinate options from the namelist file
674      REWIND(numnam)
675      READ  (numnam, namzgr)
[3432]676
[3837]677      ! Read domain options from namelist file
678      CALL dom_nam()
[3432]679
[3849]680      ! Allocate these arrays so we can use domzgr::zgr_z routine; free them
[3837]681      ! when we're done so as not to upset the 'official' allocation once
682      ! the domain decomposition is done.
683      ALLOCATE(gdepw_0(jpk), gdept_0(jpk), e3w_0(jpk), e3t_0(jpk), &
684               mig(jpiglo), mjg(jpjglo), &
685               mbathy(jpiglo,jpjglo), bathy(jpiglo,jpjglo), Stat=ierr)
686      IF(ierr /= 0)THEN
[3849]687         CALL ctl_stop('STOP', &
688                       'nemo_recursive_partition: failed to allocate zgr_z() arrays')
[3837]689         RETURN
690      END IF
691
692      ! Set-up reference depth coordinates
693      CALL zgr_z()
694
695      ! Set-up sub-domain limits as global domain for zgr_bat()
696      nldi = 2 ; nlci = jpiglo - 1
697      nldj = 2 ; nlcj = jpjglo - 1
698      jpi = jpiglo
699      jpj = jpjglo
700
701      ! Set-up fake m{i,j}g arrays for zgr_bat() call
702      DO ii = 1, jpiglo, 1
703         mig(ii) = ii
704         mi0(ii) = ii
705         mi1(ii) = ii
706      END DO
707      DO jj = 1, jpjglo, 1
708         mjg(jj) = jj
709         mj0(jj) = jj
710         mj1(jj) = jj
711      END DO
712
713      ! Initialise closed seas so loop over closed seas in zgr_bat works
714      CALL dom_clo()
715
716      ! Read-in bathy (if required) of global domain
717      CALL zgr_bat(.TRUE.)
718
[3432]719      ! land/sea mask (zero on land, 1 otherwise) over the global/zoom domain
720      imask(:,:)=1
721
[3837]722      ! Copy bathymetry in case we need to smooth it
723      zdta(:,:) = bathy(:,:)
724
725      IF(ln_sco)THEN
726         ! If ln_sco defined then the bathymetry gets smoothed before the
727         ! simulation begins and that process can alter the coastlines (bug!)
728         ! therefore we do it here too before calculating our mask.
729         CALL smooth_global_bathy(zdta, mbathy)
730      ELSE IF(ln_zps)THEN
731         CALL zgr_zps(.TRUE.)
732      ELSE IF(ln_zco)THEN
733         ! Not certain this is required since mbathy computed in zgr_bat()
734         ! in this case.
735         !CALL zgr_zco()
[3432]736      END IF
737
[3837]738      ! Compute the deepest/last ocean level for every point on the grid
739      ibotlevel(:,:) = mbathy(:,:)
740      CALL global_bot_level(ibotlevel)
741
742      ! Comment-out line below to achieve a regular partition
743      WHERE ( zdta(:,:) <= 1.0E-20 ) imask = LAND
744
[3432]745      ! Allocate partitioning arrays.
746
[3849]747      IF ( .NOT. ALLOCATED(pielb) ) THEN
[3432]748         ALLOCATE (pielb(num_pes),   pieub(num_pes), piesub(num_pes),     &
749                   pilbext(num_pes), piubext(num_pes),                    &
750                   pjelb(num_pes),   pjeub(num_pes), pjesub(num_pes),     &
751                   pjlbext(num_pes), pjubext(num_pes), pnactive(num_pes), &
[3849]752                   trimmed(4,num_pes), Stat = ierr)
[3432]753         IF(ierr /= 0)THEN
754            CALL ctl_stop('STOP', &
755                          'nemo_recursive_partition: failed to allocate partitioning arrays')
756            RETURN
757         END IF
758      ENDIF
759
[3849]760      ! Set error flag so that we calculate domain decomp if not reading
761      ! existing decomposition or if read fails.
762      ierr = 1
[3432]763
[3849]764      IF( nn_readpart )THEN
765         ! Read the partitioning to use from disk
766         CALL read_partition(ierr)
767         IF ( ierr /= 0 ) THEN
768            CALL ctl_warn('Read of pre-calculated domain decomposition failed - will calculate one instead.')
769         END IF
770      END IF
[3432]771
[3849]772      ! Set the number of sub-domains for which we are to partition
773      ! (module var in mapcomm_mod)
774      CALL set_num_subdomains(num_pes)
775
776      IF(ierr /= 0)THEN
777         ! Multi-core aware version of recursive k-section partitioning.
778         ! Currently only accounts for whether a grid point is wet or dry.
779         ! It has no knowledge of the number of wet levels at a point.
780         CALL partition_mca_rk ( imask, 1, jpiglo, 1, jpjglo, ierr )
781
782         ! Now we can do recursive k-section partitioning
783         ! ARPDBG - BUG if limits on array below are set to anything other than
784         ! 1 and jp{i,j}glo then check for external boundaries in a few lines
785         ! time WILL FAIL!
786         ! CALL partition_rk ( imask, 1, jpiglo, 1, jpjglo, ierr )
787     END IF
788
[3432]789      ! Check the error code from partitioning.
790      IF ( ierr /= 0 ) THEN
791         CALL ctl_stop('STOP','nemo_recursive_partition: Partitioning failed')
792         RETURN
793      ENDIF
794
[3837]795      ! If we used generate_fake_land() above then we must set
796      ! the mask correctly now we've partitioned. This is only
797      ! necessary when testing.
[3432]798      !WHERE ( zdta(:,:) <= 0. ) imask = 0
799
[3837]800      ! ARPDBG Quick and dirty dump to stdout in gnuplot form
801      IF(narea == 1)THEN
802         OPEN(UNIT=998, FILE="imask.dat", &
803              STATUS='REPLACE', ACTION='WRITE', IOSTAT=jj)
804         IF( jj == 0 )THEN
805            WRITE (998,*) '# Depth map'
806            WRITE (998,*) '# i   j  bathy  imask   ibotlevel   mbathy'
807            DO jj = 1, jpjglo, 1
808               DO ii = 1, jpiglo, 1
809                  WRITE (998,"(I4,1x,I4,1x,E16.6,1x,I4,1x,I4,1x,I4)") &
810                  ii, jj, zdta(ii,jj), imask(ii,jj), ibotlevel(ii,jj), mbathy(ii,jj)
811               END DO
812               WRITE (998,*)
813            END DO
814            CLOSE(998)
815         END IF
816      END IF
[3432]817
818      jpkm1 = jpk - 1
819
820      ! This chunk taken directly from original mpp_ini - not sure why nbondi
821      ! is reset? However, if it isn't reset then bad things happen in dommsk
822      ! so I'm doing what the original code does...
823      nperio = 0
824      nbondi = 0
825      IF( jperio == 1 .OR. jperio == 4 .OR. jperio == 6 ) THEN
826         IF( jpni == 1 )THEN
827            nbondi = 2
828            nperio = 1
829         END IF
830      END IF
831
832#if defined ARPDEBUG
[3837]833      ! This output is REQUIRED by the check_nemo_comms.pl test script
[3432]834      WRITE (*,FMT="(I4,' : ARPDBG: ielb, ieub, iesub = ',3I5)") narea-1,&
835            ielb, ieub, iesub
836      WRITE (*,FMT="(I4,' : ARPDBG: jelb, jeub, jesub = ',3I5)") narea-1,&
837            jelb, jeub, jesub
838      WRITE (*,FMT="(I4,' : ARPDBG: nldi, nlei, nlci = ',3I5)") narea-1, &
839            nldi, nlei, nlci
840      WRITE (*,FMT="(I4,' : ARPDBG: nldj, nlej, nlcj = ',3I5)") narea-1, &
841            nldj, nlej, nlcj
842      WRITE (*,FMT="(I4,' : ARPDBG: jpi, jpj = ',2I5)") narea-1, jpi, jpj
843      WRITE (*,FMT="(I4,' : ARPDBG: nimpp, njmpp = ',2I5)") narea-1, &
844            nimpp, njmpp
845#endif
846
847      ! Debugging option - can turn off all halo exchanges by setting this to
848      ! false.
849      do_exchanges = .TRUE.
850
[3837]851      ! Free the domzgr/_oce member arrays that we used earlier in zgr_z() and
852      ! zgr_bat().
853      DEALLOCATE(gdepw_0, gdept_0, e3w_0, e3t_0, mig, mjg,  &
854                 mbathy, bathy)
855
[3432]856   END SUBROUTINE nemo_recursive_partition
857
858
[3187]859   SUBROUTINE sqfact ( kn, kna, knb )
[2715]860      !!----------------------------------------------------------------------
[3187]861      !!                     ***  ROUTINE sqfact  ***
[2715]862      !!
[3187]863      !! ** Purpose :   return factors (kna, knb) of kn, such that
864      !!                (1) kna*knb=kn
865      !!                (2) kna and knb are as near equal as possible
866      !!                (3) kna < knb
867      !! ** Method  :   Search backwards from the square root of kn,
868      !!                until we find an integer that cleanly divides kn
869      !! ** Preconditions : kn must be positive
[2715]870      !!----------------------------------------------------------------------
[3187]871      INTEGER, INTENT(in   ) ::   kn
872      INTEGER, INTENT(  out) ::   kna, knb
873       
874      ! Search backwards from the square root of n.
[2715]875
[3849]876      fact_loop: DO kna=INT(SQRT(REAL(kn))),1,-1
[3432]877         IF ( kn/kna*kna == kn ) THEN
878            EXIT fact_loop
879         ENDIF
[3187]880      END DO fact_loop
[2715]881
[3187]882      IF( kna < 1 ) kna = 1 
[2715]883
[3187]884      ! kna divides kn cleanly. Work out the other factor.
885      knb = kn/kna
[2715]886
[3187]887   END SUBROUTINE sqfact
[2715]888
[3432]889
890   SUBROUTINE generate_fake_land(imask)
891      !!----------------------------------------------------------------------
892      !! Generate a fake land mass to test the decomposition code
893      !!----------------------------------------------------------------------
894      USE par_oce, ONLY: jpiglo, jpjglo
895      USE partition_mod, ONLY: write_partition_map
896      IMPLICIT none
897      INTEGER, DIMENSION(jpiglo,jpjglo), INTENT(inout) :: imask
898      ! Locals
899      INTEGER :: ii, jj
900      INTEGER :: icentre, jcentre
901      INTEGER :: iwidth, iheight
902      INTEGER :: istart, istop
903
904      ! imask is zero on land points , unity on ocean points
905      iwidth = jpiglo/8
906      iheight = jpjglo/8
907
908      icentre = jpiglo/2
909      jcentre = jpjglo/2
910
911      istart = icentre - iwidth
912      istop = icentre + iwidth
913      DO jj = jcentre, jcentre - iheight, -1
914         imask(istart:istop,jj) = 0
915         istart = istart + 1
916         istop = istop - 1
917      END DO
918      istart = icentre - iwidth
919      istop = icentre + iwidth
920      DO jj = jcentre+1, jcentre + iheight, 1
921         imask(istart:istop,jj) = 0
922         istart = istart + 1
923         istop = istop - 1
924      END DO
925
926! Quick and dirty dump to stdout in gnuplot form
927!!$      WRITE (*,*) 'GNUPLOT MAP'
928!!$      DO jj = 1, jpjglo, 1
929!!$         DO ii = 1, jpiglo, 1
930!!$            WRITE (*,*) ii, jj, imask(ii,jj)
931!!$         END DO
932!!$         WRITE (*,*)
933!!$      END DO
934!!$      WRITE (*,*) 'END GNUPLOT MAP'
935
936   END SUBROUTINE generate_fake_land
937
[3]938   !!======================================================================
[2496]939END MODULE nemogcm
Note: See TracBrowser for help on using the repository browser.