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.
Changeset 15459 for NEMO/trunk/src/TOP/PISCES/P4Z/p4zsms.F90 – NEMO

Ignore:
Timestamp:
2021-10-29T10:19:18+02:00 (3 years ago)
Author:
cetlod
Message:

Various bug fixes and more comments in PISCES routines ; sette test OK in debug mode, nn_hls=1/2, with tiling ; run.stat unchanged ; of course tracer.stat different

File:
1 edited

Legend:

Unmodified
Added
Removed

  • NEMO/trunk/src/TOP/PISCES/P4Z/p4zsms.F90

    r15182 r15459  
    3030   PRIVATE 
    3131 
    32    PUBLIC   p4z_sms_init   ! called in p4zsms.F90 
    33    PUBLIC   p4z_sms        ! called in p4zsms.F90 
     32   PUBLIC   p4z_sms_init   ! called in trcini_pisces.F90 
     33   PUBLIC   p4z_sms        ! called in trcsms_pisces.F90 
    3434 
    3535   INTEGER ::    numco2, numnut, numnit      ! logical unit for co2 budget 
    36    REAL(wp) ::   alkbudget, no3budget, silbudget, ferbudget, po4budget 
     36   REAL(wp) ::   alkbudget, no3budget, silbudget, ferbudget, po4budget ! total budget of the different conservative elements 
    3737   REAL(wp) ::   xfact, xfact1, xfact2, xfact3 
    3838 
    … …  
    5656      !!              routines of PISCES bio-model 
    5757      !! 
    58       !! ** Method  : - at each new day ... 
    59       !!              - several calls of bio and sed ??? 
    60       !!              - ... 
     58      !! ** Method  : - calls the various SMS subroutines 
     59      !!              - calls the sediment module (if ln_sediment)   
     60      !!              - several calls of bio and sed (possible time-splitting) 
     61      !!              - handles the potential negative concentrations (xnegtr) 
    6162      !!--------------------------------------------------------------------- 
    6263      ! 
    … …  
    9192      IF( ln_pisdmp .AND. MOD( kt - 1, nn_pisdmp ) == 0 )   CALL p4z_dmp( kt, Kbb, Kmm )      ! Relaxation of some tracers 
    9293      ! 
    93       rfact = rDt_trc 
     94      rfact = rDt_trc  ! time step of PISCES 
    9495      ! 
    9596      IF( ( ln_top_euler .AND. kt == nittrc000 )  .OR. ( .NOT.ln_top_euler .AND. kt <= nittrc000 + 1 ) ) THEN 
    96          rfactr  = 1. / rfact 
    97          rfact2  = rfact / REAL( nrdttrc, wp ) 
    98          rfact2r = 1. / rfact2 
    99          xstep = rfact2 / rday         ! Time step duration for biology 
     97         rfactr  = 1. / rfact  ! inverse of the time step 
     98         rfact2  = rfact / REAL( nrdttrc, wp )  ! time step of the biological SMS 
     99         rfact2r = 1. / rfact2  ! Inverse of the biological time step 
     100         xstep = rfact2 / rday         ! Time step duration for biology relative to a day 
    100101         xfact = 1.e+3 * rfact2r 
    101102         IF(lwp) WRITE(numout,*)  
    … …  
    131132         CALL p4z_flx( kt, jnt, Kbb, Kmm, Krhs )   ! Compute surface fluxes 
    132133         ! 
     134         ! Handling of the negative concentrations 
     135         ! The biological SMS may generate negative concentrations 
     136         ! Trends are tested at each grid cell. If a negative concentrations  
     137         ! is created at a grid cell, all the sources and sinks at that grid  
     138         ! cell are scale to avoid that negative concentration. This approach  
     139         ! is quite simplistic but it conserves mass. 
     140         ! ------------------------------------------------------------------ 
    133141         xnegtr(:,:,:) = 1.e0 
    134142         DO jn = jp_pcs0, jp_pcs1 
    … …  
    142150         !                                ! where at least 1 tracer concentration becomes negative 
    143151         !                                !  
     152         ! Concentrations are updated 
    144153         DO jn = jp_pcs0, jp_pcs1 
    145154           tr(:,:,:,jn,Kbb) = tr(:,:,:,jn,Kbb) + xnegtr(:,:,:) * tr(:,:,:,jn,Krhs) 
    … …  
    194203        ENDIF 
    195204        ! 
     205        ! Trends are are reset to 0 
    196206         DO jn = jp_pcs0, jp_pcs1 
    197207            tr(:,:,:,jn,Krhs) = 0._wp 
    … …  
    202212#endif 
    203213      ! 
     214      ! If ln_sediment is set to .true. then the sediment module is called 
    204215      IF( ln_sediment ) THEN  
    205216         ! 
    … …  
    213224         ztrbbio(:,:,:,jn) = 0._wp 
    214225      END DO 
     226      ! 
    215227      ! 
    216228      IF( l_trdtrc ) THEN 
    … …  
    236248      !!                     ***  p4z_sms_init  ***   
    237249      !! 
    238       !! ** Purpose :   read PISCES namelist 
    239       !! 
    240       !! ** input   :   file 'namelist.trc.s' containing the following 
    241       !!             namelist: natext, natbio, natsms 
     250      !! ** Purpose :   read the general PISCES namelist 
     251      !! 
     252      !! ** input   :   file 'namelist_pisces' containing the following 
     253      !!                namelist: nampisbio, nampisdmp, nampismass  
    242254      !!---------------------------------------------------------------------- 
    243255      INTEGER :: ios                 ! Local integer output status for namelist read 
    244256      !! 
    245       NAMELIST/nampisbio/ nrdttrc, wsbio, xkmort, ferat3, wsbio2, wsbio2max, wsbio2scale,    & 
    246          &                  ldocp, ldocz, lthet, no3rat3, po4rat3 
     257      NAMELIST/nampisbio/ nrdttrc, wsbio, xkmort, feratz, feratm, wsbio2, wsbio2max,    & 
     258         &                wsbio2scale, ldocp, ldocz, lthet, no3rat3, po4rat3 
    247259         ! 
    248260      NAMELIST/nampisdmp/ ln_pisdmp, nn_pisdmp 
    … …  
    268280         WRITE(numout,*) '      half saturation constant for mortality    xkmort      =', xkmort  
    269281         IF( ln_p5z ) THEN 
    270             WRITE(numout,*) '      N/C in zooplankton                        no3rat3     =', no3rat3 
    271             WRITE(numout,*) '      P/C in zooplankton                        po4rat3     =', po4rat3 
    272          ENDIF 
    273          WRITE(numout,*) '      Fe/C in zooplankton                       ferat3      =', ferat3 
     282            WRITE(numout,*) '      N/C in zooplankton                     no3rat3     =', no3rat3 
     283            WRITE(numout,*) '      P/C in zooplankton                     po4rat3     =', po4rat3 
     284         ENDIF 
     285         WRITE(numout,*) '      Fe/C in microzooplankton                  feratz      =', feratz 
     286         WRITE(numout,*) '      Fe/C in microzooplankton                  feratz      =', feratm 
    274287         WRITE(numout,*) '      Big particles sinking speed               wsbio2      =', wsbio2 
    275288         WRITE(numout,*) '      Big particles maximum sinking speed       wsbio2max   =', wsbio2max 
    … …  
    317330      !!                   ***  ROUTINE p4z_rst  *** 
    318331      !! 
    319       !!  ** Purpose : Read or write variables in restart file: 
     332      !!  ** Purpose : Read or write specific PISCES variables in restart file: 
    320333      !! 
    321334      !!  WRITE(READ) mode: 
    … …  
    330343      IF( TRIM(cdrw) == 'READ' ) THEN 
    331344         ! 
     345         ! Read the specific variable of PISCES 
    332346         IF(lwp) WRITE(numout,*) 
    333347         IF(lwp) WRITE(numout,*) ' p4z_rst : Read specific variables from pisces model ' 
    334348         IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~~~~' 
    335349         !  
     350         ! Read the pH. If not in the restart file, then it is initialized from 
     351         ! the initial conditions 
    336352         IF( iom_varid( numrtr, 'PH', ldstop = .FALSE. ) > 0 ) THEN 
    337353            CALL iom_get( numrtr, jpdom_auto, 'PH' , hi(:,:,:)  ) 
    … …  
    341357         ENDIF 
    342358         CALL iom_get( numrtr, jpdom_auto, 'Silicalim', xksi(:,:) ) 
     359 
     360         ! Read the Si half saturation constant and the maximum Silica concentration 
    343361         IF( iom_varid( numrtr, 'Silicamax', ldstop = .FALSE. ) > 0 ) THEN 
    344362            CALL iom_get( numrtr, jpdom_auto, 'Silicamax' , xksimax(:,:)  ) 
    … …  
    346364            xksimax(:,:) = xksi(:,:) 
    347365         ENDIF 
    348          ! 
     366 
     367         ! Read the Fe3 consumption term by phytoplankton 
     368         IF( iom_varid( numrtr, 'Consfe3', ldstop = .FALSE. ) > 0 ) THEN 
     369            CALL iom_get( numrtr, jpdom_auto, 'Consfe3' , consfe3(:,:,:)  ) 
     370         ELSE 
     371            consfe3(:,:,:) = 0._wp 
     372         ENDIF 
     373 
     374 
     375         ! Read the cumulative total flux. If not in the restart file, it is set to 0           
    349376         IF( iom_varid( numrtr, 'tcflxcum', ldstop = .FALSE. ) > 0 ) THEN  ! cumulative total flux of carbon 
    350377            CALL iom_get( numrtr, 'tcflxcum' , t_oce_co2_flx_cum  ) 
    … …  
    353380         ENDIF 
    354381         ! 
     382         ! PISCES size proxy 
     383         IF( iom_varid( numrtr, 'sized', ldstop = .FALSE. ) > 0 ) THEN 
     384            CALL iom_get( numrtr, jpdom_auto, 'sized' , sized(:,:,:)  ) 
     385            sized(:,:,:) = MAX( 1.0, sized(:,:,:) ) 
     386         ELSE 
     387            sized(:,:,:) = 1. 
     388         ENDIF 
     389         ! 
     390         IF( iom_varid( numrtr, 'sizen', ldstop = .FALSE. ) > 0 ) THEN 
     391            CALL iom_get( numrtr, jpdom_auto, 'sizen' , sizen(:,:,:)  ) 
     392            sizen(:,:,:) = MAX( 1.0, sizen(:,:,:) ) 
     393         ELSE 
     394            sizen(:,:,:) = 1. 
     395         ENDIF 
     396 
     397         ! PISCES-QUOTA specific part 
    355398         IF( ln_p5z ) THEN 
    356             IF( iom_varid( numrtr, 'sized', ldstop = .FALSE. ) > 0 ) THEN 
     399            ! Read the size of the different phytoplankton groups 
     400            ! If not in the restart file, they are set to 1 
     401            IF( iom_varid( numrtr, 'sizep', ldstop = .FALSE. ) > 0 ) THEN 
    357402               CALL iom_get( numrtr, jpdom_auto, 'sizep' , sizep(:,:,:)  ) 
    358                CALL iom_get( numrtr, jpdom_auto, 'sizen' , sizen(:,:,:)  ) 
    359                CALL iom_get( numrtr, jpdom_auto, 'sized' , sized(:,:,:)  ) 
     403               sizep(:,:,:) = MAX( 1.0, sizep(:,:,:) ) 
    360404            ELSE 
    361405               sizep(:,:,:) = 1. 
    362                sizen(:,:,:) = 1. 
    363                sized(:,:,:) = 1. 
    364406            ENDIF 
    365407        ENDIF 
    366408        ! 
    367409      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN 
     410         ! write the specific variables of PISCES 
    368411         IF( kt == nitrst ) THEN 
    369412            IF(lwp) WRITE(numout,*) 
    … …  
    375418         CALL iom_rstput( kt, nitrst, numrtw, 'Silicamax', xksimax(:,:) ) 
    376419         CALL iom_rstput( kt, nitrst, numrtw, 'tcflxcum', t_oce_co2_flx_cum ) 
    377          IF( ln_p5z ) THEN 
    378             CALL iom_rstput( kt, nitrst, numrtw, 'sizep', sizep(:,:,:) ) 
    379             CALL iom_rstput( kt, nitrst, numrtw, 'sizen', sizen(:,:,:) ) 
    380             CALL iom_rstput( kt, nitrst, numrtw, 'sized', sized(:,:,:) ) 
    381          ENDIF 
     420         CALL iom_rstput( kt, nitrst, numrtw, 'Consfe3', consfe3(:,:,:) ) ! Si max concentration 
     421         CALL iom_rstput( kt, nitrst, numrtw, 'tcflxcum', t_oce_co2_flx_cum ) ! Cumulative CO2 flux 
     422         CALL iom_rstput( kt, nitrst, numrtw, 'sizen', sizen(:,:,:) )  ! Size of nanophytoplankton 
     423         CALL iom_rstput( kt, nitrst, numrtw, 'sized', sized(:,:,:) )  ! Size of diatoms 
     424         IF( ln_p5z ) CALL iom_rstput( kt, nitrst, numrtw, 'sizep', sizep(:,:,:) )  ! Size of picophytoplankton 
    382425      ENDIF 
    383426      ! 
    … …  
    389432      !!                    ***  p4z_dmp  *** 
    390433      !! 
    391       !! ** purpose  : Relaxation of some tracers 
     434      !! ** purpose  : Relaxation of the total budget of some elements 
     435      !!               This routine avoids the model to drift far from the  
     436      !!               observed content in various elements 
     437      !!               Elements that may be relaxed : Alk, P, N, Si 
    392438      !!---------------------------------------------------------------------- 
    393439      ! 
    … …  
    396442      ! 
    397443      REAL(wp) ::  alkmean = 2426.     ! mean value of alkalinity ( Glodap ; for Goyet 2391. ) 
    398       REAL(wp) ::  po4mean = 2.165     ! mean value of phosphates 
    399       REAL(wp) ::  no3mean = 30.90     ! mean value of nitrate 
    400       REAL(wp) ::  silmean = 91.51     ! mean value of silicate 
     444      REAL(wp) ::  po4mean = 2.174     ! mean value of phosphate 
     445      REAL(wp) ::  no3mean = 31.00     ! mean value of nitrate 
     446      REAL(wp) ::  silmean = 90.33     ! mean value of silicate 
    401447      ! 
    402448      REAL(wp) :: zarea, zalksumn, zpo4sumn, zno3sumn, zsilsumn 
    … …  
    419465            zsilsumn = glob_sum( 'p4zsms', tr(:,:,:,jpsil,Kmm) * cvol(:,:,:)  ) * zarea 
    420466  
     467            ! Correct the trn mean content of alkalinity 
    421468            IF(lwp) WRITE(numout,*) '       TALKN mean : ', zalksumn 
    422469            tr(:,:,:,jptal,Kmm) = tr(:,:,:,jptal,Kmm) * alkmean / zalksumn 
    423470 
     471            ! Correct the trn mean content of PO4 
    424472            IF(lwp) WRITE(numout,*) '       PO4N  mean : ', zpo4sumn 
    425473            tr(:,:,:,jppo4,Kmm) = tr(:,:,:,jppo4,Kmm) * po4mean / zpo4sumn 
    426474 
     475            ! Correct the trn mean content of NO3 
    427476            IF(lwp) WRITE(numout,*) '       NO3N  mean : ', zno3sumn 
    428477            tr(:,:,:,jpno3,Kmm) = tr(:,:,:,jpno3,Kmm) * no3mean / zno3sumn 
    429478 
     479            ! Correct the trn mean content of SiO3 
    430480            IF(lwp) WRITE(numout,*) '       SiO3N mean : ', zsilsumn 
    431481            tr(:,:,:,jpsil,Kmm) = MIN( 400.e-6,tr(:,:,:,jpsil,Kmm) * silmean / zsilsumn ) 
    … …  
    439489  
    440490               IF(lwp) WRITE(numout,*) ' ' 
     491               ! Correct the trb mean content of alkalinity 
    441492               IF(lwp) WRITE(numout,*) '       TALKB mean : ', zalksumb 
    442493               tr(:,:,:,jptal,Kbb) = tr(:,:,:,jptal,Kbb) * alkmean / zalksumb 
    443494 
     495               ! Correct the trb mean content of PO4 
    444496               IF(lwp) WRITE(numout,*) '       PO4B  mean : ', zpo4sumb 
    445497               tr(:,:,:,jppo4,Kbb) = tr(:,:,:,jppo4,Kbb) * po4mean / zpo4sumb 
    446498 
     499               ! Correct the trb mean content of NO3 
    447500               IF(lwp) WRITE(numout,*) '       NO3B  mean : ', zno3sumb 
    448501               tr(:,:,:,jpno3,Kbb) = tr(:,:,:,jpno3,Kbb) * no3mean / zno3sumb 
    449502 
     503               ! Correct the trb mean content of SiO3 
    450504               IF(lwp) WRITE(numout,*) '       SiO3B mean : ', zsilsumb 
    451505               tr(:,:,:,jpsil,Kbb) = MIN( 400.e-6,tr(:,:,:,jpsil,Kbb) * silmean / zsilsumb ) 
    … …  
    487541      ENDIF 
    488542 
     543      ! Compute the budget of NO3 
    489544      IF( iom_use( "pno3tot" ) .OR. ( ln_check_mass .AND. kt == nitend )  ) THEN 
    490          !   Compute the budget of NO3, ALK, Si, Fer 
    491545         IF( ln_p4z ) THEN 
    492546            zwork(:,:,:) =    tr(:,:,:,jpno3,Kmm) + tr(:,:,:,jpnh4,Kmm)                      & 
    … …  
    506560      ENDIF 
    507561      ! 
     562      ! Compute the budget of PO4 
    508563      IF( iom_use( "ppo4tot" ) .OR. ( ln_check_mass .AND. kt == nitend )  ) THEN 
    509564         IF( ln_p4z ) THEN 
    … …  
    524579      ENDIF 
    525580      ! 
     581      ! Compute the budget of SiO3 
    526582      IF( iom_use( "psiltot" ) .OR. ( ln_check_mass .AND. kt == nitend )  ) THEN 
    527583         zwork(:,:,:) =  tr(:,:,:,jpsil,Kmm) + tr(:,:,:,jpgsi,Kmm) + tr(:,:,:,jpdsi,Kmm)  
    … …  
    540596      ENDIF 
    541597      ! 
     598      ! Compute the budget of Iron 
    542599      IF( iom_use( "pfertot" ) .OR. ( ln_check_mass .AND. kt == nitend )  ) THEN 
    543600         zwork(:,:,:) =   tr(:,:,:,jpfer,Kmm) + tr(:,:,:,jpnfe,Kmm) + tr(:,:,:,jpdfe,Kmm)   & 
    544601            &         +   tr(:,:,:,jpbfe,Kmm) + tr(:,:,:,jpsfe,Kmm)                      & 
    545             &         + ( tr(:,:,:,jpzoo,Kmm) + tr(:,:,:,jpmes,Kmm) )  * ferat3     
     602            &         + ( tr(:,:,:,jpzoo,Kmm) * feratz + tr(:,:,:,jpmes,Kmm) ) * feratm    
    546603         ! 
    547604         ferbudget = glob_sum( 'p4zsms', zwork(:,:,:) * cvol(:,:,:)  )   
Note: See TracChangeset for help on using the changeset viewer.