source: LMDZ6/trunk/libf/phylmd/phyaqua_mod.f90 @ 5627

Last change on this file since 5627 was 5627, checked in by amaison, 2 months ago

Representation of heterogeneous continental subsurfaces with parameter or flux aggregation in the simplified surface model (bucket) for 1D case studies.

  • Property copyright set to
    Name of program: LMDZ
    Creation date: 1984
    Version: LMDZ5
    License: CeCILL version 2
    Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
    See the license file in the root directory
  • Property svn:keywords set to Id
File size: 30.3 KB
RevLine 
[3540]1!
2! $Id: phyaqua_mod.f90 5627 2025-04-25 12:21:02Z amaison $
3!
[1992]4MODULE phyaqua_mod
5  ! Routines complementaires pour la physique planetaire.
[5282]6  USE clesphys_mod_h
7    IMPLICIT NONE
[1529]8
[1992]9CONTAINS
[1529]10
[3579]11  SUBROUTINE iniaqua(nlon,year_len,iflag_phys)
[1529]12
[1992]13    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
14    ! Creation d'un etat initial et de conditions aux limites
15    ! (resp startphy.nc et limit.nc) pour des configurations idealisees
16    ! du modele LMDZ dans sa version terrestre.
17    ! iflag_phys est un parametre qui controle
18    ! iflag_phys = N
19    ! de 100 a 199 : aqua planetes avec SST forcees
20    ! N-100 determine le type de SSTs
21    ! de 200 a 299 : terra planetes avec Ts calcule
[1529]22
[1992]23    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
[1529]24
[1992]25    USE dimphy, ONLY: klon
[2351]26    USE geometry_mod, ONLY : latitude
[1992]27    USE surface_data, ONLY: type_ocean, ok_veget
28    USE pbl_surface_mod, ONLY: pbl_surface_init
29    USE fonte_neige_mod, ONLY: fonte_neige_init
30    USE phys_state_var_mod
[2344]31    USE time_phylmdz_mod, ONLY: day_ref, ndays, pdtphys, &
32                                day_ini,day_end
[1992]33    USE indice_sol_mod
[2344]34    USE nrtype, ONLY: pi
[3435]35!    USE ioipsl
36    USE mod_phys_lmdz_para, ONLY: is_master
37    USE mod_phys_lmdz_transfert_para, ONLY: bcast
38    USE mod_grid_phy_lmdz
39    USE ioipsl_getin_p_mod, ONLY : getin_p
[5273]40    USE phys_cal_mod , ONLY: calend, year_len_phy => year_len
41    USE dimsoil_mod_h, ONLY: nsoilmx
[5285]42    USE yomcst_mod_h
[5274]43IMPLICIT NONE
[1529]44
[5274]45
[1671]46
[3579]47    INTEGER, INTENT (IN) :: nlon, year_len, iflag_phys
[1992]48    ! IM ajout latfi, lonfi
[2351]49!    REAL, INTENT (IN) :: lonfi(nlon), latfi(nlon)
[1529]50
[1992]51    INTEGER type_profil, type_aqua
[1529]52
[1992]53    ! Ajouts initialisation des surfaces
54    REAL :: run_off_lic_0(nlon)
55    REAL :: qsolsrf(nlon, nbsrf), snsrf(nlon, nbsrf)
56    REAL :: tsoil(nlon, nsoilmx, nbsrf)
57    REAL :: tslab(nlon), seaice(nlon)
[2243]58    REAL fder(nlon)
[1529]59
60
61
[1992]62    ! Arguments :
63    ! -----------
[1529]64
[1992]65    ! integer radpas
66    INTEGER it, unit, i, k, itap
[1529]67
[1992]68    REAL rugos, albedo
69    REAL tsurf
70    REAL time, timestep, day, day0
[2243]71    REAL qsol_f
[1992]72    REAL rugsrel(nlon)
73    LOGICAL alb_ocean
[1529]74
[1992]75    CHARACTER *80 ans, file_forctl, file_fordat, file_start
76    CHARACTER *100 file, var
77    CHARACTER *2 cnbl
[1529]78
[3540]79    REAL phy_nat(nlon, year_len)
80    REAL phy_alb(nlon, year_len)
81    REAL phy_sst(nlon, year_len)
82    REAL phy_bil(nlon, year_len)
83    REAL phy_rug(nlon, year_len)
84    REAL phy_ice(nlon, year_len)
85    REAL phy_fter(nlon, year_len)
86    REAL phy_foce(nlon, year_len)
87    REAL phy_fsic(nlon, year_len)
88    REAL phy_flic(nlon, year_len)
[1529]89
[1992]90    INTEGER, SAVE :: read_climoz = 0 ! read ozone climatology
[3435]91!$OMP THREADPRIVATE(read_climoz)
[1529]92
[1992]93    ! -------------------------------------------------------------------------
94    ! declaration pour l'appel a phyredem
95    ! -------------------------------------------------------------------------
[1529]96
[1992]97    ! real pctsrf(nlon,nbsrf),ftsol(nlon,nbsrf)
98    REAL falbe(nlon, nbsrf), falblw(nlon, nbsrf)
99    ! real pbl_tke(nlon,llm,nbsrf)
100    ! real rain_fall(nlon),snow_fall(nlon)
101    ! real solsw(nlon), sollw(nlon),radsol(nlon)
102    ! real t_ancien(nlon,llm),q_ancien(nlon,llm),rnebcon(nlon,llm)
103    ! real ratqs(nlon,llm)
104    ! real clwcon(nlon,llm)
[1529]105
[1992]106    INTEGER longcles
107    PARAMETER (longcles=20)
108    REAL clesphy0(longcles)
[1529]109
110
[1992]111    ! -----------------------------------------------------------------------
112    ! dynamial tendencies :
113    ! ---------------------
[1529]114
[1992]115    INTEGER l, ierr, aslun
[1529]116
[1992]117    REAL paire
[1529]118
[3531]119    ! Local
120    CHARACTER (LEN=20) :: modname='phyaqua'
121    CHARACTER (LEN=80) :: abort_message
[1529]122
[3531]123
[1992]124    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
125    ! INITIALISATIONS
126    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
[1529]127
[1992]128    ! -----------------------------------------------------------------------
129    ! Initialisations  des constantes
130    ! -------------------------------
[1529]131
[3579]132    !IF (calend .EQ. "earth_360d") Then
133      year_len_phy = year_len
134    !END IF
135   
[5084]136    if (year_len.ne.360) then
[3579]137      write (*,*) year_len
[4463]138      call abort_physic("iniaqua", 'iniaqua: 360 day calendar is required !', 1)
[3540]139    endif
140
[1992]141    type_aqua = iflag_phys/100
142    type_profil = iflag_phys - type_aqua*100
143    PRINT *, 'iniaqua:type_aqua, type_profil', type_aqua, type_profil
[1529]144
[1992]145    IF (klon/=nlon) THEN
146      WRITE (*, *) 'iniaqua: klon=', klon, ' nlon=', nlon
[3531]147      abort_message= 'probleme de dimensions dans iniaqua'
148      CALL abort_physic(modname,abort_message,1)
[1992]149    END IF
150    CALL phys_state_var_init(read_climoz)
[1529]151
152
[1992]153    read_climoz = 0
154    day0 = 217.
155    day = day0
156    it = 0
157    time = 0.
[1529]158
[1992]159    ! -----------------------------------------------------------------------
160    ! initialisations de la physique
161    ! -----------------------------------------------------------------------
[1529]162
[2344]163    day_ini = day_ref
164    day_end = day_ini + ndays
[1759]165
[3435]166    nbapp_rad = 24
167    CALL getin_p('nbapp_rad', nbapp_rad)
168
[1992]169    ! ---------------------------------------------------------------------
170    ! Creation des conditions aux limites:
171    ! ------------------------------------
172    ! Initialisations des constantes
173    ! Ajouter les manquants dans planete.def... (albedo etc)
[3435]174    co2_ppm = 348.
175    CALL getin_p('co2_ppm', co2_ppm)
176
177    solaire = 1365.
178    CALL getin_p('solaire', solaire)
179 
[1992]180    ! CALL getin('albedo',albedo) ! albedo is set below, depending on
181    ! type_aqua
[3435]182    alb_ocean = .TRUE.
183    CALL getin_p('alb_ocean', alb_ocean)
184
[1992]185    WRITE (*, *) 'iniaqua: co2_ppm=', co2_ppm
186    WRITE (*, *) 'iniaqua: solaire=', solaire
187    WRITE (*, *) 'iniaqua: alb_ocean=', alb_ocean
[1529]188
[1992]189    radsol = 0.
190    qsol_f = 10.
[1529]191
[1992]192    ! Conditions aux limites:
193    ! -----------------------
[1529]194
[1992]195    qsol(:) = qsol_f
196    rugsrel = 0.0 ! (rugsrel = rugoro)
197    rugoro = 0.0
198    u_ancien = 0.0
199    v_ancien = 0.0
200    agesno = 50.0
201    ! Relief plat
202    zmea = 0.
203    zstd = 0.
204    zsig = 0.
205    zgam = 0.
206    zthe = 0.
207    zpic = 0.
208    zval = 0.
[1529]209
[1992]210    ! Une seule surface
211    pctsrf = 0.
212    IF (type_aqua==1) THEN
213      rugos = 1.E-4
214      albedo = 0.19
215      pctsrf(:, is_oce) = 1.
216    ELSE IF (type_aqua==2) THEN
217      rugos = 0.03
218      albedo = 0.1
219      pctsrf(:, is_ter) = 1.
220    END IF
[1529]221
[3435]222    CALL getin_p('rugos', rugos)
223
[1992]224    WRITE (*, *) 'iniaqua: rugos=', rugos
[2209]225    zmasq(:) = pctsrf(:, is_ter)
[1529]226
[1992]227    ! pctsrf_pot(:,is_oce) = 1. - zmasq(:)
228    ! pctsrf_pot(:,is_sic) = 1. - zmasq(:)
[1529]229
[1992]230    ! Si alb_ocean on calcule un albedo oceanique moyen
231    ! if (alb_ocean) then
232    ! Voir pourquoi on avait ca.
233    ! CALL ini_alb_oce(phy_alb)
234    ! else
235    phy_alb(:, :) = albedo ! albedo land only (old value condsurf_jyg=0.3)
236    ! endif !alb_ocean
[1529]237
[3540]238    DO i = 1, year_len
[1992]239      ! IM Terraplanete   phy_sst(:,i) = 260.+50.*cos(rlatd(:))**2
240      ! IM ajout calcul profil sst selon le cas considere (cf. FBr)
[1529]241
[1992]242      phy_nat(:, i) = 1.0 ! 0=ocean libre, 1=land, 2=glacier, 3=banquise
243      phy_bil(:, i) = 1.0 ! ne sert que pour les slab_ocean
244      phy_rug(:, i) = rugos ! longueur rugosite utilisee sur land only
245      phy_ice(:, i) = 0.0 ! fraction de glace (?)
246      phy_fter(:, i) = pctsrf(:, is_ter) ! fraction de glace (?)
247      phy_foce(:, i) = pctsrf(:, is_oce) ! fraction de glace (?)
248      phy_fsic(:, i) = pctsrf(:, is_sic) ! fraction de glace (?)
249      phy_flic(:, i) = pctsrf(:, is_lic) ! fraction de glace (?)
250    END DO
251    ! IM calcul profil sst
[2351]252    CALL profil_sst(nlon, latitude, type_profil, phy_sst)
[1529]253
[3435]254    IF (grid_type==unstructured) THEN
255      CALL writelim_unstruct(klon, phy_nat, phy_alb, phy_sst, phy_bil, phy_rug, phy_ice, &
256                             phy_fter, phy_foce, phy_flic, phy_fsic)
257    ELSE
258     
259       CALL writelim(klon, phy_nat, phy_alb, phy_sst, phy_bil, phy_rug, phy_ice, &
260                     phy_fter, phy_foce, phy_flic, phy_fsic)
261    ENDIF
[1529]262
[1992]263    ! ---------------------------------------------------------------------
264    ! Ecriture de l'etat initial:
265    ! ---------------------------
[1529]266
267
[1992]268    ! Ecriture etat initial physique
[1529]269
[2344]270    timestep = pdtphys
271    radpas = nint(rday/timestep/float(nbapp_rad))
[1529]272
[1992]273    DO i = 1, longcles
274      clesphy0(i) = 0.
275    END DO
276    clesphy0(1) = float(iflag_con)
277    clesphy0(2) = float(nbapp_rad)
278    ! IF( cycle_diurne  ) clesphy0(3) =  1.
279    clesphy0(3) = 1. ! cycle_diurne
280    clesphy0(4) = 1. ! soil_model
281    clesphy0(5) = 1. ! new_oliq
282    clesphy0(6) = 0. ! ok_orodr
283    clesphy0(7) = 0. ! ok_orolf
284    clesphy0(8) = 0. ! ok_limitvrai
[1529]285
286
[1992]287    ! =======================================================================
288    ! Profils initiaux
289    ! =======================================================================
[1529]290
[1992]291    ! On initialise les temperatures de surfaces comme les sst
292    DO i = 1, nlon
293      ftsol(i, :) = phy_sst(i, 1)
294      tsoil(i, :, :) = phy_sst(i, 1)
295      tslab(i) = phy_sst(i, 1)
296    END DO
[1529]297
[1992]298    falbe(:, :) = albedo
299    falblw(:, :) = albedo
300    rain_fall(:) = 0.
301    snow_fall(:) = 0.
302    solsw(:) = 0.
303    sollw(:) = 0.
304    radsol(:) = 0.
[1529]305
[1992]306    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
307    ! intialisation bidon mais pas grave
308    t_ancien(:, :) = 0.
309    q_ancien(:, :) = 0.
310    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
311    rnebcon = 0.
312    ratqs = 0.
313    clwcon = 0.
314    pbl_tke = 1.E-8
[1529]315
[1992]316    ! variables supplementaires pour appel a plb_surface_init
317    fder(:) = 0.
318    seaice(:) = 0.
319    run_off_lic_0 = 0.
[2243]320    fevap = 0.
[1529]321
322
[1992]323    ! Initialisations necessaires avant phyredem
324    type_ocean = 'force'
325    CALL fonte_neige_init(run_off_lic_0)
326    qsolsrf(:, :) = qsol(1) ! humidite du sol des sous surface
327    snsrf(:, :) = 0. ! couverture de neige des sous surface
[2243]328    z0m(:, :) = rugos ! couverture de neige des sous surface
329    z0h=z0m
[1530]330
331
[2243]332    CALL pbl_surface_init(fder, snsrf, qsolsrf, tsoil)
[1529]333
[1992]334    PRINT *, 'iniaqua: before phyredem'
[1529]335
[3435]336    pbl_tke(:,:,:) = 1.e-8
[1992]337    falb1 = albedo
338    falb2 = albedo
339    zmax0 = 0.
340    f0 = 0.
341    sig1 = 0.
342    w01 = 0.
343    wake_deltat = 0.
344    wake_deltaq = 0.
345    wake_s = 0.
[3435]346    wake_dens = 0.
[1992]347    wake_cstar = 0.
348    wake_pe = 0.
349    wake_fip = 0.
350    fm_therm = 0.
351    entr_therm = 0.
352    detr_therm = 0.
[2827]353    ale_bl = 0.
354    ale_bl_trig =0.
355    alp_bl =0.
[2979]356    treedrg(:,:,:)=0.
[5627]357    tsurf_tersrf(:,:) = 0.
358    qsurf_tersrf(:,:) = 0.
359    cdragm_tersrf(:,:) = 0.
360    cdragh_tersrf(:,:) = 0.
361    swnet_tersrf(:,:) = 0.
362    lwnet_tersrf(:,:) = 0.
363    fluxsens_tersrf(:,:) = 0.
364    fluxlat_tersrf(:,:) = 0.
[1529]365
[3579]366    u10m = 0.
367    v10m = 0.
368
369    ql_ancien   = 0.
370    qs_ancien   = 0.
[4523]371    qbs_ancien  = 0.
[3579]372    u_ancien    = 0.
373    v_ancien    = 0.
374    prw_ancien  = 0.
375    prlw_ancien = 0.
376    prsw_ancien = 0. 
[4523]377    prbsw_ancien= 0.
[3579]378
379    ale_wake    = 0.
380    ale_bl_stat = 0. 
381
382
[3435]383!ym error : the sub surface dimension is the third not second : forgotten for iniaqua
384!    falb_dir(:,is_ter,:)=0.08; falb_dir(:,is_lic,:)=0.6
385!    falb_dir(:,is_oce,:)=0.5;  falb_dir(:,is_sic,:)=0.6
386    falb_dir(:,:,is_ter)=0.08; falb_dir(:,:,is_lic)=0.6
387    falb_dir(:,:,is_oce)=0.5;  falb_dir(:,:,is_sic)=0.6
[1529]388
[3435]389!ym falb_dif has been forgotten, initialize with defaukt value found in phyetat0 or 0 ?
390!ym probably the uninitialized value was 0 for standard (regular grid) case
391    falb_dif(:,:,:)=0
392
393
[1992]394    CALL phyredem('startphy.nc')
[1529]395
[1992]396    PRINT *, 'iniaqua: after phyredem'
397    CALL phys_state_var_end
[1529]398
[1992]399    RETURN
400  END SUBROUTINE iniaqua
[1529]401
402
[1992]403  ! ====================================================================
404  ! ====================================================================
405  SUBROUTINE zenang_an(cycle_diurne, gmtime, rlat, rlon, rmu0, fract)
[5285]406    USE yomcst_mod_h
[1992]407    USE dimphy
408    IMPLICIT NONE
409    ! ====================================================================
410    ! =============================================================
411    ! CALL zenang(cycle_diurne,gmtime,rlat,rlon,rmu0,fract)
412    ! Auteur : A. Campoy et F. Hourdin
413    ! Objet  : calculer les valeurs moyennes du cos de l'angle zenithal
414    ! et l'ensoleillement moyen entre gmtime1 et gmtime2
415    ! connaissant la declinaison, la latitude et la longitude.
[1529]416
[1992]417    ! Dans cette version particuliere, on calcule le rayonnement
418    ! moyen sur l'année à chaque latitude.
419    ! angle zenithal calculé pour obtenir un
420    ! Fit polynomial de  l'ensoleillement moyen au sommet de l'atmosphere
421    ! en moyenne annuelle.
422    ! Spécifique de la terre. Utilisé pour les aqua planetes.
[1529]423
[1992]424    ! Rque   : Different de la routine angle en ce sens que zenang
425    ! fournit des moyennes de pmu0 et non des valeurs
426    ! instantanees, du coup frac prend toutes les valeurs
427    ! entre 0 et 1.
428    ! Date   : premiere version le 13 decembre 1994
429    ! revu pour  GCM  le 30 septembre 1996
430    ! ===============================================================
431    ! longi----INPUT : la longitude vraie de la terre dans son plan
432    ! solaire a partir de l'equinoxe de printemps (degre)
433    ! gmtime---INPUT : temps universel en fraction de jour
434    ! pdtrad---INPUT : pas de temps du rayonnement (secondes)
435    ! lat------INPUT : latitude en degres
436    ! long-----INPUT : longitude en degres
437    ! pmu0-----OUTPUT: angle zenithal moyen entre gmtime et gmtime+pdtrad
438    ! frac-----OUTPUT: ensoleillement moyen entre gmtime et gmtime+pdtrad
439    ! ================================================================
440    ! ================================================================
441    LOGICAL cycle_diurne
442    REAL gmtime
443    REAL rlat(klon), rlon(klon), rmu0(klon), fract(klon)
444    ! ================================================================
445    INTEGER i
446    REAL gmtime1, gmtime2
447    REAL pi_local
[1529]448
449
[1992]450    REAL rmu0m(klon), rmu0a(klon)
[1529]451
452
[1992]453    pi_local = 4.0*atan(1.0)
[1529]454
[1992]455    ! ================================================================
456    ! Calcul de l'angle zenithal moyen sur la journee
457    ! ================================================================
458
459    DO i = 1, klon
460      fract(i) = 1.
461      ! Calcule du flux moyen
462      IF (abs(rlat(i))<=28.75) THEN
463        rmu0m(i) = (210.1924+206.6059*cos(0.0174533*rlat(i))**2)/1365.
464      ELSE IF (abs(rlat(i))<=43.75) THEN
465        rmu0m(i) = (187.4562+236.1853*cos(0.0174533*rlat(i))**2)/1365.
466      ELSE IF (abs(rlat(i))<=71.25) THEN
467        rmu0m(i) = (162.4439+284.1192*cos(0.0174533*rlat(i))**2)/1365.
[1529]468      ELSE
[1992]469        rmu0m(i) = (172.8125+183.7673*cos(0.0174533*rlat(i))**2)/1365.
470      END IF
471    END DO
[1529]472
[1992]473    ! ================================================================
474    ! Avec ou sans cycle diurne
475    ! ================================================================
[1529]476
[1992]477    IF (cycle_diurne) THEN
[1529]478
[1992]479      ! On redecompose flux  au sommet suivant un cycle diurne idealise
480      ! identique a toutes les latitudes.
[1671]481
[1992]482      DO i = 1, klon
483        rmu0a(i) = 2.*rmu0m(i)*sqrt(2.)*pi_local/(4.-pi_local)
484        rmu0(i) = rmu0a(i)*abs(sin(pi_local*gmtime+pi_local*rlon(i)/360.)) - &
485          rmu0a(i)/sqrt(2.)
486      END DO
[1671]487
[1992]488      DO i = 1, klon
489        IF (rmu0(i)<=0.) THEN
490          rmu0(i) = 0.
491          fract(i) = 0.
492        ELSE
493          fract(i) = 1.
494        END IF
495      END DO
[1671]496
[1992]497      ! Affichage de l'angel zenitale
498      ! print*,'************************************'
499      ! print*,'************************************'
500      ! print*,'************************************'
501      ! print*,'latitude=',rlat(i),'longitude=',rlon(i)
502      ! print*,'rmu0m=',rmu0m(i)
503      ! print*,'rmu0a=',rmu0a(i)
504      ! print*,'rmu0=',rmu0(i)
[1529]505
[1992]506    ELSE
[1671]507
[1992]508      DO i = 1, klon
509        fract(i) = 0.5
510        rmu0(i) = rmu0m(i)*2.
511      END DO
512
513    END IF
514
515    RETURN
516  END SUBROUTINE zenang_an
517
518  ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
519
[3435]520  SUBROUTINE writelim_unstruct(klon, phy_nat, phy_alb, phy_sst, phy_bil, phy_rug, &
521      phy_ice, phy_fter, phy_foce, phy_flic, phy_fsic)
522
523    USE mod_phys_lmdz_para, ONLY: is_omp_master, klon_mpi
524    USE mod_phys_lmdz_transfert_para, ONLY: gather_omp
[4619]525    USE lmdz_xios
[3435]526    IMPLICIT NONE
527
528    INTEGER, INTENT (IN) :: klon
529    REAL, INTENT (IN) :: phy_nat(klon, 360)
530    REAL, INTENT (IN) :: phy_alb(klon, 360)
531    REAL, INTENT (IN) :: phy_sst(klon, 360)
532    REAL, INTENT (IN) :: phy_bil(klon, 360)
533    REAL, INTENT (IN) :: phy_rug(klon, 360)
534    REAL, INTENT (IN) :: phy_ice(klon, 360)
535    REAL, INTENT (IN) :: phy_fter(klon, 360)
536    REAL, INTENT (IN) :: phy_foce(klon, 360)
537    REAL, INTENT (IN) :: phy_flic(klon, 360)
538    REAL, INTENT (IN) :: phy_fsic(klon, 360)
539
540    REAL :: phy_mpi(klon_mpi, 360) ! temporary variable, to store phy_***(:)
541      ! on the whole physics grid
542 
[4619]543    IF (using_xios) THEN
544      PRINT *, 'writelim: Ecriture du fichier limit'
[3435]545
[4619]546      CALL gather_omp(phy_foce, phy_mpi)
547      IF (is_omp_master) CALL xios_send_field('foce_limout',phy_mpi)
[3435]548
[4619]549      CALL gather_omp(phy_fsic, phy_mpi)
550      IF (is_omp_master) CALL xios_send_field('fsic_limout',phy_mpi)
[3435]551     
[4619]552      CALL gather_omp(phy_fter, phy_mpi)
553      IF (is_omp_master) CALL xios_send_field('fter_limout',phy_mpi)
[3435]554     
[4619]555      CALL gather_omp(phy_flic, phy_mpi)
556      IF (is_omp_master) CALL xios_send_field('flic_limout',phy_mpi)
[3435]557
[4619]558      CALL gather_omp(phy_sst, phy_mpi)
559      IF (is_omp_master) CALL xios_send_field('sst_limout',phy_mpi)
[3435]560
[4619]561      CALL gather_omp(phy_bil, phy_mpi)
562      IF (is_omp_master) CALL xios_send_field('bils_limout',phy_mpi)
[3435]563
[4619]564      CALL gather_omp(phy_alb, phy_mpi)
565      IF (is_omp_master) CALL xios_send_field('alb_limout',phy_mpi)
[3435]566
[4619]567      CALL gather_omp(phy_rug, phy_mpi)
568      IF (is_omp_master) CALL xios_send_field('rug_limout',phy_mpi)
569    ENDIF
[3435]570  END SUBROUTINE writelim_unstruct
571
572
573
[1992]574  SUBROUTINE writelim(klon, phy_nat, phy_alb, phy_sst, phy_bil, phy_rug, &
575      phy_ice, phy_fter, phy_foce, phy_flic, phy_fsic)
576
[3435]577    USE mod_phys_lmdz_para, ONLY: is_master
[1992]578    USE mod_grid_phy_lmdz, ONLY: klon_glo
579    USE mod_phys_lmdz_transfert_para, ONLY: gather
[3540]580    USE phys_cal_mod, ONLY: year_len
[5270]581    USE netcdf, ONLY: nf90_def_var, nf90_put_var, nf90_get_var, nf90_strerror, nf90_close, &
582            nf90_enddef, nf90_put_att, nf90_unlimited, nf90_noerr, nf90_global, nf90_clobber, &
583            nf90_64bit_offset, nf90_def_dim, nf90_create
[5249]584    USE lmdz_cppkeys_wrapper, ONLY: nf90_format
[1992]585    IMPLICIT NONE
586
587    INTEGER, INTENT (IN) :: klon
[3540]588    REAL, INTENT (IN) :: phy_nat(klon, year_len)
589    REAL, INTENT (IN) :: phy_alb(klon, year_len)
590    REAL, INTENT (IN) :: phy_sst(klon, year_len)
591    REAL, INTENT (IN) :: phy_bil(klon, year_len)
592    REAL, INTENT (IN) :: phy_rug(klon, year_len)
593    REAL, INTENT (IN) :: phy_ice(klon, year_len)
594    REAL, INTENT (IN) :: phy_fter(klon, year_len)
595    REAL, INTENT (IN) :: phy_foce(klon, year_len)
596    REAL, INTENT (IN) :: phy_flic(klon, year_len)
597    REAL, INTENT (IN) :: phy_fsic(klon, year_len)
[1992]598
[3540]599    REAL :: phy_glo(klon_glo, year_len) ! temporary variable, to store phy_***(:)
[1992]600      ! on the whole physics grid
601    INTEGER :: k
602    INTEGER ierr
603    INTEGER dimfirst(3)
604    INTEGER dimlast(3)
605
606    INTEGER nid, ndim, ntim
607    INTEGER dims(2), debut(2), epais(2)
608    INTEGER id_tim
609    INTEGER id_nat, id_sst, id_bils, id_rug, id_alb
610    INTEGER id_fter, id_foce, id_fsic, id_flic
611
[3435]612    IF (is_master) THEN
[1992]613
614      PRINT *, 'writelim: Ecriture du fichier limit'
615
[5270]616      ierr = nf90_create('limit.nc', IOR(nf90_clobber,nf90_64bit_offset), nid)
[1992]617
[5270]618      ierr = nf90_put_att(nid, nf90_global, 'title', 'Fichier conditions aux limites')
619      ! !        ierr = nf90_def_dim (nid, "points_physiques", klon, ndim)
620      ierr = nf90_def_dim(nid, 'points_physiques', klon_glo, ndim)
621      ierr = nf90_def_dim(nid, 'time', nf90_unlimited, ntim)
[1992]622
623      dims(1) = ndim
624      dims(2) = ntim
625
[5249]626      ierr = nf90_def_var(nid, 'TEMPS', nf90_format, [ntim], id_tim)
[5270]627      ierr = nf90_put_att(nid, id_tim, 'title', 'Jour dans l annee')
[2198]628
[5249]629      ierr = nf90_def_var(nid, 'NAT', nf90_format, dims, id_nat)
[5270]630      ierr = nf90_put_att(nid, id_nat, 'title', 'Nature du sol (0,1,2,3)')
[2198]631
[5249]632      ierr = nf90_def_var(nid, 'SST', nf90_format, dims, id_sst)
[5270]633      ierr = nf90_put_att(nid, id_sst, 'title', 'Temperature superficielle de la mer')
[2198]634
[5249]635      ierr = nf90_def_var(nid, 'BILS', nf90_format, dims, id_bils)
[5270]636      ierr = nf90_put_att(nid, id_bils, 'title', 'Reference flux de chaleur au sol')
[2198]637
[5249]638      ierr = nf90_def_var(nid, 'ALB', nf90_format, dims, id_alb)
[5270]639      ierr = nf90_put_att(nid, id_alb, 'title', 'Albedo a la surface')
[2198]640
[5249]641      ierr = nf90_def_var(nid, 'RUG', nf90_format, dims, id_rug)
[5270]642      ierr = nf90_put_att(nid, id_rug, 'title', 'Rugosite')
[1992]643
[5249]644      ierr = nf90_def_var(nid, 'FTER', nf90_format, dims, id_fter)
[5270]645      ierr = nf90_put_att(nid, id_fter, 'title','Frac. Land')
[5249]646      ierr = nf90_def_var(nid, 'FOCE', nf90_format, dims, id_foce)
[5270]647      ierr = nf90_put_att(nid, id_foce, 'title','Frac. Ocean')
[5249]648      ierr = nf90_def_var(nid, 'FSIC', nf90_format, dims, id_fsic)
[5270]649      ierr = nf90_put_att(nid, id_fsic, 'title','Frac. Sea Ice')
[5249]650      ierr = nf90_def_var(nid, 'FLIC', nf90_format, dims, id_flic)
[5270]651      ierr = nf90_put_att(nid, id_flic, 'title','Frac. Land Ice')
[1992]652
[5270]653      ierr = nf90_enddef(nid)
654      IF (ierr/=nf90_noerr) THEN
[1992]655        WRITE (*, *) 'writelim error: failed to end define mode'
[5270]656        WRITE (*, *) nf90_strerror(ierr)
[1992]657      END IF
658
659
660      ! write the 'times'
[3540]661      DO k = 1, year_len
[5249]662        ierr = nf90_put_var(nid, id_tim, k, [k])
[5270]663        IF (ierr/=nf90_noerr) THEN
[1992]664          WRITE (*, *) 'writelim error with temps(k),k=', k
[5270]665          WRITE (*, *) nf90_strerror(ierr)
[1992]666        END IF
667      END DO
[1529]668
[3435]669    END IF ! of if (is_master)
[1671]670
[1992]671    ! write the fields, after having collected them on master
[1671]672
[1992]673    CALL gather(phy_nat, phy_glo)
[3435]674    IF (is_master) THEN
[5249]675      ierr = nf90_put_var(nid, id_nat, phy_glo)
[5270]676      IF (ierr/=nf90_noerr) THEN
[1992]677        WRITE (*, *) 'writelim error with phy_nat'
[5270]678        WRITE (*, *) nf90_strerror(ierr)
[1992]679      END IF
680    END IF
[1671]681
[1992]682    CALL gather(phy_sst, phy_glo)
[3435]683    IF (is_master) THEN
[5249]684      ierr = nf90_put_var(nid, id_sst, phy_glo)
[5270]685      IF (ierr/=nf90_noerr) THEN
[1992]686        WRITE (*, *) 'writelim error with phy_sst'
[5270]687        WRITE (*, *) nf90_strerror(ierr)
[1992]688      END IF
689    END IF
[1671]690
[1992]691    CALL gather(phy_bil, phy_glo)
[3435]692    IF (is_master) THEN
[5249]693      ierr = nf90_put_var(nid, id_bils, phy_glo)
[5270]694      IF (ierr/=nf90_noerr) THEN
[1992]695        WRITE (*, *) 'writelim error with phy_bil'
[5270]696        WRITE (*, *) nf90_strerror(ierr)
[1992]697      END IF
698    END IF
[1671]699
[1992]700    CALL gather(phy_alb, phy_glo)
[3435]701    IF (is_master) THEN
[5249]702      ierr = nf90_put_var(nid, id_alb, phy_glo)
[5270]703      IF (ierr/=nf90_noerr) THEN
[1992]704        WRITE (*, *) 'writelim error with phy_alb'
[5270]705        WRITE (*, *) nf90_strerror(ierr)
[1992]706      END IF
707    END IF
[1671]708
[1992]709    CALL gather(phy_rug, phy_glo)
[3435]710    IF (is_master) THEN
[5249]711      ierr = nf90_put_var(nid, id_rug, phy_glo)
[5270]712      IF (ierr/=nf90_noerr) THEN
[1992]713        WRITE (*, *) 'writelim error with phy_rug'
[5270]714        WRITE (*, *) nf90_strerror(ierr)
[1992]715      END IF
716    END IF
[1671]717
[1992]718    CALL gather(phy_fter, phy_glo)
[3435]719    IF (is_master) THEN
[5249]720      ierr = nf90_put_var(nid, id_fter, phy_glo)
[5270]721      IF (ierr/=nf90_noerr) THEN
[1992]722        WRITE (*, *) 'writelim error with phy_fter'
[5270]723        WRITE (*, *) nf90_strerror(ierr)
[1992]724      END IF
725    END IF
[1671]726
[1992]727    CALL gather(phy_foce, phy_glo)
[3435]728    IF (is_master) THEN
[5249]729      ierr = nf90_put_var(nid, id_foce, phy_glo)
[5270]730      IF (ierr/=nf90_noerr) THEN
[1992]731        WRITE (*, *) 'writelim error with phy_foce'
[5270]732        WRITE (*, *) nf90_strerror(ierr)
[1992]733      END IF
734    END IF
[1671]735
[1992]736    CALL gather(phy_fsic, phy_glo)
[3435]737    IF (is_master) THEN
[5249]738      ierr = nf90_put_var(nid, id_fsic, phy_glo)
[5270]739      IF (ierr/=nf90_noerr) THEN
[1992]740        WRITE (*, *) 'writelim error with phy_fsic'
[5270]741        WRITE (*, *) nf90_strerror(ierr)
[1992]742      END IF
743    END IF
[1671]744
[1992]745    CALL gather(phy_flic, phy_glo)
[3435]746    IF (is_master) THEN
[5249]747      ierr = nf90_put_var(nid, id_flic, phy_glo)
[5270]748      IF (ierr/=nf90_noerr) THEN
[1992]749        WRITE (*, *) 'writelim error with phy_flic'
[5270]750        WRITE (*, *) nf90_strerror(ierr)
[1992]751      END IF
752    END IF
[1671]753
[1992]754    ! close file:
[3435]755    IF (is_master) THEN
[5270]756      ierr = nf90_close(nid)
[1992]757    END IF
[1671]758
[1992]759  END SUBROUTINE writelim
[1529]760
[1992]761  ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
[1671]762
[1992]763  SUBROUTINE profil_sst(nlon, rlatd, type_profil, phy_sst)
764    USE dimphy
[3540]765    USE phys_cal_mod , ONLY: year_len
[1992]766    IMPLICIT NONE
[1529]767
[1992]768    INTEGER nlon, type_profil, i, k, j
[3540]769    REAL :: rlatd(nlon), phy_sst(nlon, year_len)
[1992]770    INTEGER imn, imx, amn, amx, kmn, kmx
771    INTEGER p, pplus, nlat_max
772    PARAMETER (nlat_max=72)
773    REAL x_anom_sst(nlat_max)
[3531]774    CHARACTER (LEN=20) :: modname='profil_sst'
775    CHARACTER (LEN=80) :: abort_message
[1529]776
[3531]777    IF (klon/=nlon) THEN
778       abort_message='probleme de dimensions dans profil_sst'
779       CALL abort_physic(modname,abort_message,1)
780    ENDIF
[1992]781    WRITE (*, *) ' profil_sst: type_profil=', type_profil
[3540]782    DO i = 1, year_len
[1992]783      ! phy_sst(:,i) = 260.+50.*cos(rlatd(:))**2
[1529]784
[1992]785      ! Rajout fbrlmd
[1529]786
[1992]787      IF (type_profil==1) THEN
788        ! Méthode 1 "Control" faible plateau à l'Equateur
789        DO j = 1, klon
790          phy_sst(j, i) = 273. + 27.*(1-sin(1.5*rlatd(j))**2)
791          ! PI/3=1.047197551
792          IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
793            phy_sst(j, i) = 273.
794          END IF
795        END DO
796      END IF
797      IF (type_profil==2) THEN
798        ! Méthode 2 "Flat" fort plateau à l'Equateur
799        DO j = 1, klon
800          phy_sst(j, i) = 273. + 27.*(1-sin(1.5*rlatd(j))**4)
801          IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
802            phy_sst(j, i) = 273.
803          END IF
804        END DO
805      END IF
[1529]806
807
[1992]808      IF (type_profil==3) THEN
809        ! Méthode 3 "Qobs" plateau réel à l'Equateur
810        DO j = 1, klon
811          phy_sst(j, i) = 273. + 0.5*27.*(2-sin(1.5*rlatd(j))**2-sin(1.5* &
812            rlatd(j))**4)
813          IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
814            phy_sst(j, i) = 273.
815          END IF
816        END DO
817      END IF
[1529]818
[1992]819      IF (type_profil==4) THEN
820        ! Méthode 4 : Méthode 3 + SST+2 "Qobs" plateau réel à l'Equateur
821        DO j = 1, klon
822          phy_sst(j, i) = 273. + 0.5*29.*(2-sin(1.5*rlatd(j))**2-sin(1.5* &
823            rlatd(j))**4)
824          IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
825            phy_sst(j, i) = 273.
826          END IF
827        END DO
828      END IF
[1529]829
[1992]830      IF (type_profil==5) THEN
831        ! Méthode 5 : Méthode 3 + +2K "Qobs" plateau réel à l'Equateur
832        DO j = 1, klon
833          phy_sst(j, i) = 273. + 2. + 0.5*27.*(2-sin(1.5*rlatd(j))**2-sin(1.5 &
834            *rlatd(j))**4)
835          IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
836            phy_sst(j, i) = 273. + 2.
837          END IF
[1529]838
[1992]839        END DO
840      END IF
[1529]841
[1992]842      IF (type_profil==6) THEN
843        ! Méthode 6 "cst" valeur constante de SST
844        DO j = 1, klon
845          phy_sst(j, i) = 288.
846        END DO
847      END IF
[1529]848
849
[1992]850      IF (type_profil==7) THEN
851        ! Méthode 7 "cst" valeur constante de SST +2
852        DO j = 1, klon
853          phy_sst(j, i) = 288. + 2.
854        END DO
855      END IF
[1529]856
[1992]857      p = 0
858      IF (type_profil==8) THEN
859        ! Méthode 8 profil anomalies SST du modèle couplé AR4
860        DO j = 1, klon
861          IF (rlatd(j)==rlatd(j-1)) THEN
862            phy_sst(j, i) = 273. + x_anom_sst(pplus) + &
863              0.5*27.*(2-sin(1.5*rlatd(j))**2-sin(1.5*rlatd(j))**4)
864            IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
865              phy_sst(j, i) = 273. + x_anom_sst(pplus)
866            END IF
867          ELSE
868            p = p + 1
869            pplus = 73 - p
870            phy_sst(j, i) = 273. + x_anom_sst(pplus) + &
871              0.5*27.*(2-sin(1.5*rlatd(j))**2-sin(1.5*rlatd(j))**4)
872            IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
873              phy_sst(j, i) = 273. + x_anom_sst(pplus)
874            END IF
875            WRITE (*, *) rlatd(j), x_anom_sst(pplus), phy_sst(j, i)
876          END IF
877        END DO
878      END IF
[1529]879
[1992]880      IF (type_profil==9) THEN
881        ! Méthode 5 : Méthode 3 + -2K "Qobs" plateau réel à l'Equateur
882        DO j = 1, klon
883          phy_sst(j, i) = 273. - 2. + 0.5*27.*(2-sin(1.5*rlatd(j))**2-sin(1.5 &
884            *rlatd(j))**4)
885          IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
886            phy_sst(j, i) = 273. - 2.
887          END IF
888        END DO
889      END IF
[1529]890
891
[1992]892      IF (type_profil==10) THEN
893        ! Méthode 10 : Méthode 3 + +4K "Qobs" plateau réel à l'Equateur
894        DO j = 1, klon
895          phy_sst(j, i) = 273. + 4. + 0.5*27.*(2-sin(1.5*rlatd(j))**2-sin(1.5 &
896            *rlatd(j))**4)
897          IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
898            phy_sst(j, i) = 273. + 4.
899          END IF
900        END DO
901      END IF
[1529]902
[1992]903      IF (type_profil==11) THEN
904        ! Méthode 11 : Méthode 3 + 4CO2 "Qobs" plateau réel à l'Equateur
905        DO j = 1, klon
906          phy_sst(j, i) = 273. + 0.5*27.*(2-sin(1.5*rlatd(j))**2-sin(1.5* &
907            rlatd(j))**4)
908          IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
909            phy_sst(j, i) = 273.
910          END IF
911        END DO
912      END IF
[1529]913
[1992]914      IF (type_profil==12) THEN
915        ! Méthode 12 : Méthode 10 + 4CO2 "Qobs" plateau réel à l'Equateur
916        DO j = 1, klon
917          phy_sst(j, i) = 273. + 4. + 0.5*27.*(2-sin(1.5*rlatd(j))**2-sin(1.5 &
918            *rlatd(j))**4)
919          IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
920            phy_sst(j, i) = 273. + 4.
921          END IF
922        END DO
923      END IF
[1529]924
[1992]925      IF (type_profil==13) THEN
926        ! Méthode 13 "Qmax" plateau réel à l'Equateur augmenté !
927        DO j = 1, klon
928          phy_sst(j, i) = 273. + 0.5*29.*(2-sin(1.5*rlatd(j))**2-sin(1.5* &
929            rlatd(j))**4)
930          IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
931            phy_sst(j, i) = 273.
932          END IF
933        END DO
934      END IF
[1529]935
[1992]936      IF (type_profil==14) THEN
937        ! Méthode 13 "Qmax2K" plateau réel à l'Equateur augmenté +2K !
938        DO j = 1, klon
939          phy_sst(j, i) = 273. + 2. + 0.5*29.*(2-sin(1.5*rlatd(j))**2-sin(1.5 &
940            *rlatd(j))**4)
941          IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
942            phy_sst(j, i) = 273.
943          END IF
944        END DO
945      END IF
[1529]946
[5084]947      if (type_profil.EQ.20) then
[2107]948      print*,'Profile SST 20'
949!     Méthode 13 "Qmax2K" plateau réel �|  l'Equateur augmenté +2K
950
951      do j=1,klon
952        phy_sst(j,i)=248.+55.*(1-sin(rlatd(j))**2)
953      enddo
954      endif
955
[5084]956      if (type_profil.EQ.21) then
[2107]957      print*,'Profile SST 21'
958!     Méthode 13 "Qmax2K" plateau réel �|  l'Equateur augmenté +2K
959      do j=1,klon
960        phy_sst(j,i)=252.+55.*(1-sin(rlatd(j))**2)
961      enddo
962      endif
963
964
965
[1992]966    END DO
967
968    ! IM beg : verif profil SST: phy_sst
969    amn = min(phy_sst(1,1), 1000.)
970    amx = max(phy_sst(1,1), -1000.)
971    imn = 1
972    kmn = 1
973    imx = 1
974    kmx = 1
[3540]975    DO k = 1, year_len
[1992]976      DO i = 2, nlon
977        IF (phy_sst(i,k)<amn) THEN
978          amn = phy_sst(i, k)
979          imn = i
980          kmn = k
981        END IF
982        IF (phy_sst(i,k)>amx) THEN
983          amx = phy_sst(i, k)
984          imx = i
985          kmx = k
986        END IF
987      END DO
988    END DO
989
990    PRINT *, 'profil_sst: imn, kmn, phy_sst(imn,kmn) ', imn, kmn, amn
991    PRINT *, 'profil_sst: imx, kmx, phy_sst(imx,kmx) ', imx, kmx, amx
992    ! IM end : verif profil SST: phy_sst
993
994    RETURN
995  END SUBROUTINE profil_sst
996
997END MODULE phyaqua_mod
Note: See TracBrowser for help on using the repository browser.