source: LMDZ5/branches/testing/libf/phylmd/phyaqua_mod.F90 @ 5403

Last change on this file since 5403 was 2839, checked in by Laurent Fairhead, 8 years ago

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