source: LMDZ6/branches/Amaury_dev/libf/phylmd/phyaqua_mod.F90 @ 5157

Last change on this file since 5157 was 5144, checked in by abarral, 7 weeks ago

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